text:
/* gsl_histogram_stat.c
 * Copyright (C) 2000  Simone Piccardi
 *
 * This library 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 library; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
/***************************************************************
 *
 * File gsl_histogram_stat.c:
 * Routines for statisticalcomputations on histograms.
 * Need GSL library and header.
 * Contains the routines:
 * gsl_histogram_mean    compute histogram mean
 * gsl_histogram_sigma   compute histogram sigma
 *
 * Author: S. Piccardi
 * Jan. 2000
 *
 ***************************************************************/
#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

/* FIXME: We skip negative values in the histogram h->bin[i] < 0,
   since those correspond to negative weights (BJG) */

double
gsl_histogram_mean (const gsl_histogram * h)
{
  const size_t n = h->n;
  size_t i;

  /* Compute the bin-weighted arithmetic mean M of a histogram using the
     recurrence relation

     M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n)))
     W(n) = W(n-1) + w(n)

   */

  long double wmean = 0;
  long double W = 0;

  for (i = 0; i < n; i++)
    {
      double xi = (h->range[i + 1] + h->range[i]) / 2;
      double wi = h->bin[i];

      if (wi > 0)
        {
          W += wi;
          wmean += (xi - wmean) * (wi / W);
        }
    }

  return wmean;
}

double
gsl_histogram_sigma (const gsl_histogram * h)
{
  const size_t n = h->n;
  size_t i;

  long double wvariance = 0 ;
  long double wmean = 0;
  long double W = 0;

  /* Use a two-pass algorithm for stability.  Could also use a single
     pass formula, as given in N.J.Higham 'Accuracy and Stability of
     Numerical Methods', p.12 */

  /* Compute the mean */

  for (i = 0; i < n; i++)
    {
      double xi = (h->range[i + 1] + h->range[i]) / 2;
      double wi = h->bin[i];

      if (wi > 0)
        {
          W += wi;
          wmean += (xi - wmean) * (wi / W);
        }
    }

  /* Compute the variance */

  W = 0.0;

  for (i = 0; i < n; i++)
    {
      double xi = ((h->range[i + 1]) + (h->range[i])) / 2;
      double wi = h->bin[i];

      if (wi > 0) {
        const long double delta = (xi - wmean);
        W += wi ;
        wvariance += (delta * delta - wvariance) * (wi / W);
      }
    }

  {
    double sigma = sqrt (wvariance) ;
    return sigma;
  }
}

/*
  sum up all bins of histogram
 */

double
gsl_histogram_sum(const gsl_histogram * h)
{
  double sum=0;
  size_t i=0;
  size_t n;
  n=h->n;

  while(i < n)
    sum += h->bin[i++];

  return sum;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* spmatrix/util.c
 *
 * Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018 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 <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>

/*
gsl_spmatrix_cumsum()

Compute the cumulative sum:

p[j] = Sum_{k=0...j-1} c[k]

0 <= j < n + 1

Alternatively,
p[0] = 0
p[j] = p[j - 1] + c[j - 1]

Inputs: n - length of input array
        c - (input/output) array of size n + 1
            on input, contains the n values c[k]
            on output, contains the n + 1 values p[j]

Return: success or error
*/

void
gsl_spmatrix_cumsum(const size_t n, int * c)
{
  int sum = 0;
  size_t k;

  for (k = 0; k < n; ++k)
    {
      int ck = c[k];
      c[k] = sum;
      sum += ck;
    }

  c[n] = sum;
}

ext: c

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <cblas.h>

#include "tasks.h"

void trsm_task_seq(void *ptr)
{
    struct trsm_task_arg *arg = (struct trsm_task_arg*) ptr;

    int n       = arg->n;
    int m       = arg->m;
    double *A21 = arg->A21;
    double *A11 = arg->A11;
    int ldA     = arg->ldA;

    // Compute A21 := A21 * inv(A11').
    cblas_dtrsm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit,
                m, n,
                1.0, A11, ldA,
                     A21, ldA);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
#ifndef S3D_CV_IMAGE_OPERATION_INPUT_OUTPUT_ADAPTER_H
#define S3D_CV_IMAGE_OPERATION_INPUT_OUTPUT_ADAPTER_H

#include <s3d/multiview/stan_results.h>

#include <opencv2/core/mat.hpp>

#include <gsl/gsl>

namespace s3d {
namespace image_operation {

class ImageOperations;

class InputOutputAdapter {
public:
  InputOutputAdapter(gsl::not_null<ImageOperations*> imageOperations);

  void setInputImages(const cv::Mat& leftImage, const cv::Mat& rightImage);
  std::tuple<cv::Mat&, cv::Mat&> getOutputImages();
  bool applyAllOperations();
  void copyInputToOutputs();
  bool canApplyOperations();

  StanResults results;

private:
  cv::Mat inputLeftImage{};
  cv::Mat inputRightImage{};
  cv::Mat outputLeftImage{};
  cv::Mat outputRightImage{};

  ImageOperations* operations_{};
};

} // namespace s3d
} // namespace image_operation

#endif // S3D_CV_IMAGE_OPERATION_INPUT_OUTPUT_ADAPTER_H

ext: h


text:
#include <stdio.h>
#include <math.h>
#include <nlopt.h>

double myfunc(unsigned n, const double *x, double *grad, void *my_func_data)
{
    //printf("%f\t%f\t", grad[0], grad[1]);
    //printf("%f\t%f\n", x[0], x[1]);
    if (grad) {
        grad[0] = 0.0;
        grad[1] = 0.5 / sqrt(x[1]);
    }
    //printf("%f\n", sqrt(x[1]));
    //sleep(1);
    return sqrt(x[1]);
}

typedef struct {
    double a, b;
} my_constraint_data;

double myconstraint(unsigned n, const double *x, double *grad, void *data)
{
    my_constraint_data *d = (my_constraint_data *) data;
    double a = d->a, b = d->b;
    if (grad) {
        grad[0] = 3 * a * (a*x[0] + b) * (a*x[0] + b);
        grad[1] = -1.0;
    }
    return ((a*x[0] + b) * (a*x[0] + b) * (a*x[0] + b) - x[1]);
 }

int main(){
    double lb[2] = { -HUGE_VAL, 0 }; /* lower bounds */
    nlopt_opt opt;

    //opt = nlopt_create(NLOPT_LD_MMA, 2); /* algorithm and dimensionality */
    opt = nlopt_create(NLOPT_LN_COBYLA, 2); /* algorithm and dimensionality */
    nlopt_set_lower_bounds(opt, lb);
    nlopt_set_min_objective(opt, myfunc, NULL);

    my_constraint_data data[2] = { {2,0}, {-1,1} };

    nlopt_add_inequality_constraint(opt, myconstraint, &data[0], 1e-8);
    nlopt_add_inequality_constraint(opt, myconstraint, &data[1], 1e-8);

    //nlopt_set_xtol_rel(opt, 1e-4);
    nlopt_set_stopval(opt, (sqrt(8./27.) + 1e-3));

    double x[2] = { 1.234, 5.678 };  /* `*`some` `initial` `guess`*` */
    double minf; /* `*`the` `minimum` `objective` `value,` `upon` `return`*` */
    if (nlopt_optimize(opt, x, &minf) < 0) {
	printf("nlopt failed!\n");
    }
    else {
	printf("found minimum at f(%g,%g) = %0.10g\n", x[0], x[1], minf);
    }

    nlopt_destroy(opt);
}

max_stars_repo_name: Vindaar/nimnlopt

max_issues_repo_name: Vindaar/nimnlopt

max_forks_repo_name: Vindaar/nimnlopt


text:
#ifndef CollectorMap_h
#define CollectorMap_h
/** This is free and unencumbered software released into the public domain.
The authors of ISIS do not claim copyright on the contents of this file.
For more details about the LICENSE terms and the AUTHORS, you will
find files of those names at the top level of this repository. **/

/* SPDX-License-Identifier: CC0-1.0 */
#include <cctype>
#include <QString>
#include <map>
#include <algorithm>
#include <functional>
#include <sstream>
#include <iomanip>
#include <exception>
#include "IException.h"
#include "IString.h"
#include <gsl/gsl_math.h>

namespace Isis {

  /**
   * @brief Provides a simple comparison between two values
   *
   * This simple comparison function object is provided with no special frills
   * that does pretty much exactly what std::less does.
   */
  template <typename K> struct SimpleCompare {

    /**
     * Returns true if v1 is less than v2
     *
     * @param v1 Input constant
     * @param v2 Input constant
     *
     * @return bool Returns true if v1 is less than v2
     */
    bool operator()(const K &v1, const K &v2) const {
      return (v1 < v2);
    }

  };

  /**
  * @brief Provides a case insensitive string comparison
  *
  * This string comparison functor object compares two strings ignoring case. Use
  * this policy when your key into the collector map is a string and you want to
  * ignore case when finding elements in the collection.
  */
  template <typename K> struct NoCaseStringCompare {

    /**
     * Compares v1 and v2 as case insensitive strings, and returns true of v1 is
     * less than v2 (as those strings).
     *
     * @param v1 Input constant
     * @param v2 Input constant
     *
     * @return bool Returns true if v1 is less than v2 in string format
     */
    bool operator()(const K &v1, const K &v2) const {
      return (IString::DownCase(v1) < IString::DownCase(v2));
    }

  };

  /**
   * @brief Provides a robust comparison of double/float values
   *
   * This functor compares floating point values using a default epsilon of
   * 1.0e-6. It can be used for doubles or floats, however floats will be promoted
   * to double for the comparison.
   */
  template <typename K> struct RobustFloatCompare {

    /**
     * Compares v1 and v2 as floating point values.
     *
     * @param v1
     * @param v2
     *
     * @return bool
     */
    bool operator()(const K &v1, const K &v2) const {
      return (gsl_fcmp(v1, v2, -1.0E-6) < 0);
    }

  };

  /**
   * @brief Supplies a NOOP default for removal of a CollectorMap entry
   *
   * This simple declaration is basically a NOOP that implements removal
   * of a CollectionMap entry.  It is most useful (and the default behavior)
   * when the storage element of the CollectorMap is anything but a pointer.
   * Pointers that require deletion should use the PointerRemoval policy unless
   * the pointers are owned by another "entity".
   */
  template <typename T> struct NoopRemoval {

    // defines used for cross platform suppression of unused parameter warning. This will prevent
    // distroy() from causing unused parameter warnings in clang. This is done for the purpose
    // of porting ISIS to OSX 10.11
    // #define MON_Internal_UnusedStringify(macro_arg_string_literal) #macro_arg_string_literal
    // #define MONUnusedParameter(macro_arg_parameter) _Pragma(MON_Internal_UnusedStringify(unused(macro_arg_parameter)))

    protected:

      /**
       * Destroys the CollectorMap entry
       *
       * @param element The CollectorMap to be destroyed
       */
      void destroy(T *element) {
        // arbitrary cast type to suppress unused parameter warnings generated
        // by clang when building on mac
        (void)element;

        return;
      }

  };

  /**
   * @brief Supplies a policy for deleting pointers that CollectorMap owns
   *
   * Defines a method to delete pointers when removed from a CollectorMap.
   * This is necessary to prevent memory leaks and defer the deletion to
   * removal from CollectorMap class.
   */
  template <typename T> struct PointerRemoval {
    protected:

      /**
       * Destroys the CollectorMap pointer's CollectorMap
       *
       * @param element The pointer pointing to the CollectorMap to be destroyed
       */
      void destroy(T *element) {
        delete(*element);
        return;
      }

  };

  /**
   * @brief Policy for deleting arrays that CollectorMap owns
   *
   * Defines a method to delete arrays when removed from a CollectorMap.
   * This is necessary to prevent memory leaks and defer the deletion to
   * removal from CollectorMap class.
   */
  template <typename T> struct ArrayRemoval {
    protected:

      /**
       * Destroys the array of CollectorMaps
       *
       * @param element The array of CollectorMaps to be destroyed
       */
      void destroy(T *element) {
        delete [](*element);
        return;
      }

  };

  /**
   * @brief (Default) Policy for copying map elements
   *
   * Defines a method to copy simple elements from an existing map to a
   * destination map.  This policy just makes a direct copy of the element to the
   * destination.
   *
   * This policy assumes the assignment operator handles the proper copying of
   * each element T in the collection.
   */
  template <typename T> struct DefaultCopy {
    protected:

      /**
       * Returns a copy of the input
       *
       * @param src The map element to be copied
       *
       * @return const T& The copy of the input
       */
      const T &copy(const T &src) const {
        return (src);
      }

  };

  /**
   * @brief Pointer to object policy for copying map elements
   *
   * Defines a copy method to properly handle pointers to objects (assumed) when
   * copying the complete CollectorMap.  This implementation assumes the copy
   * constructor properly handles the creation of a new element from a different
   * one.
   *
   * This policy assumes the assignment operator handles the proper copying of
   * each element T* in the collection.
   *
   * This employs an intersting technique of redirection.  Because the type T is
   * actually T*, the templated allocate() method exists to get down to the T
   * class base level.  Looks strange but it works.
   */
  template <typename T> struct PointerCopy {
    protected:
      /**
       * @brief Allocate new object using copy construtor and new pointer
       *
       * This copy method takes a pointer to a pointer (T is actually a T*) and
       * allocates a new object using the copy constructor.
       *
       * @param src  Pointer to pointer of new class to allocate
       *
       * @return T  Pointer to new object type T
       */
      T copy(const T &src)  const {
        return (allocate(*(src)));
      }

    private:
      /**
       * @brief Allocate new object using copy constructor
       *
       * @param obj Source object to create new one from
       *
       * @return P* Pointer to newly allocated object
       */
      template <typename P>
      P *allocate(const P &obj) const {
        return (new P(obj));
      }
  };

  /**
   *  @brief Collector/container for arbitrary items
   *
   *  Used to contain types with iterators of const and non-const conditions.
   *  This is a multimap that contains arbitrary keys with arbitrary elements.  It
   *  is intended to be used for pointers and copyable objects.  They should be
   *  rather efficient in the copy out operation so large objects may not be
   *  suitable or classes that do not have a good copy operator.  During testing
   *  it was noted that an object is copied up to four times and destroyed three
   *  times upon an add() operation.
   *
   *  This class is implemented using policies.  The ComparePolicy is used to test
   *  key elements such as strings and double values.  The NoCaseStringCompare
   *  policy is provided that expedites case insensitive string key comparisons.
   *  The RobustFloatCompare implements the comparison of double or float key
   *  types.  Direct comparisons of floats can be problematic due to round off and
   *  storage manifestations of these values in conputers.   The default policy,
   *  SimpleCompare, does a simple parameter to key equality test.
   *
   *  The RemovalPolicy is provided when a map value is removed from the list.
   *  This allows pointers and arrays to be stored in the map as well.  To store
   *  pointers, use PointerRemoval and for arrays there is the ArrayRemoval
   *  policy.  The default is the NoopRemoval policy which simply lets the
   *  destructor handle removals.
   *
   *  The CopyPolicy is necessary to properly handle the copying of elements.
   *  This is especially important for pointers and arrays.  In order to minimize
   *  difficult passing strategies, map elements are passed by address and the
   *  return type is the element type.  DefaultCopy simply copies the elements as
   *  is relying on the element T assigment operator to do the right thing. For
   *  pointers to objects, the PointerCopy allocates the object using the copy
   *  constructor.  One could provide a similar operator assuming a clone()
   *  method existed for the type T element. The ArrayCopy policy is left to the
   *  user to provide their own as it cannot support arrays of varying length.
   *  (One should use std::vector instead!) Users can supply their own CopyPolicy
   *  that need only expose a copy(cont T *src) method.
   *
   *  Here are some examples that demonstrate how this policy-based template class
   *  can be used:
   *
   *  @code
   * // Create a unique string key list that stores double floating point values.
   * // Use the default removal and copy policies but allow testing for character
   * // keys without regard to case via the NoCaseStringCompare.
   *  #include "CollectorMap.h"
   *
   *  CollectorMap<QString, double, NoCaseStringCompare > dmap;
   *  cout << "\nSize of double map = " << dmap.size() << endl;
   *  dmap.add("one", 1.0);
   *  dmap.add("two", 2.0);
   *  cout << "Size of double map = " << dmap.size() << endl;
   *
   *  cout << "One = " << dmap.get("one") << endl;
   *  cout << "Two = " << dmap.get("Two") << endl;
   *
   *  const double &one = dmap.get("one");
   *  cout << "\nTest Const one = " << one << endl;
   *
   *  dmap.remove("one");
   *  @endcode
   *
   *  Using this class internal to classes is perhaps where it may be applied more
   *  frequently.  The example below shows how to declare an integer key using
   *  pointers to classes:
   *
   *  @code
   * #include "CollectorMap.h"
   *
   *  class ClassTest {
   *   public:
   *    ClassTest(int n = 0) : _n(n) {  }
   *    ~ClassTest() {  }
   *    int Ident() const { return (_n); }
   *   private:
   *    int _n;
   *  };
   *
   *
   * //  Typedefs are sometimes convenient in these cases
   *  typedef CollectorMap<int, ClassTest *, SimpleCompare,
   *                       PointerRemoval, PointerCopy> PointerMap;
   *
   *  PointerMap ctest2;
   *  ctest2.add(4,new ClassTest(4));
   *  ctest2.add(5,new ClassTest(5));
   *  ctest2.add(6,new ClassTest(6));
   *  ctest2.add(7,new ClassTest(7));
   *
   *  cout << "Remove ClassTest 6\n";
   *  ctest2.remove(6);
   *
   * //  Creates a copy of ctest2 using the PointerCopy policy
   *  PointerMap map2(ctest2);
   *
   *  cout << "Find element 7: "  << map2.find(7)->Ident() << endl;
   *
   *  @endcode
   *
   *  And, finally, an example of how to use duplicate keys:
   *
   *  @code
   *  #include "CollectorMap.h"
   *
   *  typedef CollectorMap<int,QString> IntStr;
   *  IntStr dupstr(IntStr::DuplicateKeys);
   *  dupstr.add(1,"One");
   *  dupstr.add(1, "One #2");
   *  dupstr.add(1,"One #3");
   *  dupstr.add(2,"Two");
   *  dupstr.add(2,"Two #2");
   *  dupstr.add(3,"Three");
   *
   *  cout << "Size of Dup object: " << dupstr.size() << endl;
   *  cout << "Number Ones:   " << dupstr.count(1) << endl;
   *  cout << "Number Twos:   " << dupstr.count(2) << endl;
   *  cout << "Number Threes: " << dupstr.count(3) << endl;
   *  cout << "Number Fours:  " << dupstr.count(4) << endl;
   *
   *  IntStr::CollectorConstIter isIter;
   *  int j = 0;
   *  for (isIter = dupstr.begin() ; isIter != dupstr.end() ; ++isIter, j++) {
   *       cout << "IntStr[" << j << "] = {" << isIter->first << ", "
   *            << isIter->second << "}, Index: " << dupstr.index(isIter->first)
   *            << endl;
   *       cout << "Nth Test Ident       = " << dupstr.getNth(j) << endl;
   *  }
   *  @endcode
   *
   *  The output of the above example is:
   *  @code
   * Size of Dup object: 6
   * Number Ones:   3
   * Number Twos:   2
   * Number Threes: 1
   * Number Fours:  0
   * IntStr[0] = {1, One}, Index: 0
   * Nth Test Ident       = One
   * IntStr[1] = {1, One #2}, Index: 0
   * Nth Test Ident       = One #2
   * IntStr[2] = * {1, One #3}, Index: 0
   * Nth Test Ident       = One #3
   * IntStr[3] * = {2, Two}, Index: 3
   * Nth Test Ident       = Two
   * IntStr[4] = * {2, Two #2}, Index: 3
   * Nth Test Ident       = Two #2
   * IntStr[5] = * {3, Three}, Index: 5
   * Nth Test Ident = Three
   *  @endcode
   *
   *  @ingroup Utility
   *
   *  @author 2006-06-21 Kris Becker
   *
   *  @internal
   *    @history 2006-07-03 Kris Becker Added the ability to stored duplicate keys
   *             if needed (using a multimap instead of a map).  Initial default
   *             behavior of unique keys is retained.  See KeyPolicy.
   *    @history 2006-07-28 Kris Becker Fixed a bug in the NoCaseStringCompare
   *             implementation.  Prior to this fix, it would not function
   *             properly at all for case-insenstive keys.
   *    @history 2006-08-30 Kris Becker Fixed bug in copy constructors that
   *             attempted to use a virtual method in the object being created
   *             when it *must* use the method from the one it is being created
   *             from. (g++ 4.1 on Suse 10.1 didn't like this bug at all!)
   *    @history 2008-06-18 Christopher Austin Fixed Documentation
   *    @history 2017-08-30 Summer Stapleton - Updated documentation. References #4807.
   *
   */
  template < typename K, typename T,
           template <class> class ComparePolicy = SimpleCompare,
           template <class> class RemovalPolicy = NoopRemoval,
           template <class> class CopyPolicy = DefaultCopy
           >
  class CollectorMap : public RemovalPolicy<T>, public CopyPolicy<T> {
    public:
      typedef T                                      CollectorType; //!< Data type
      /** A multimap attacking a key to a CollectorType and a ComparePolicy<CollectorType>*/
      typedef std::multimap<K, CollectorType, ComparePolicy<K> > CollectorList;
      //! CollectorList iterator type declaration
      typedef typename CollectorList::iterator       CollectorIter;
      //! CollectorList constant iterator type declaration
      typedef typename CollectorList::const_iterator CollectorConstIter;

      /**
       * @brief Enumerated selection of key behaviour
       *
       * Using this enumeration during construction allows the user of this class
       * to specify if the keys used to identify elements are unique or can be
       * duplicated.
       */
      enum KeyPolicy { UniqueKeys,     //!<  Constrain keys to be unique
                       DuplicateKeys   //!<  Allow duplication of keys
                     };

      /** Constructor */
      CollectorMap() : _keyPolicy(UniqueKeys) { }

      /**
       * @brief Allows the user to choose if keys can be duplicated
       *
       * This constructor is provided to the user that wants to explicity define how
       * the keys, namely insertions are managed.  The default is unique keys in the
       * noop constructor...this one allows instantiation of either policy.
       *
       * @param keyPolicy   Can be UniqueKeys or DuplicateKeys
       */
      CollectorMap(const KeyPolicy &keyPolicy) : _keyPolicy(keyPolicy) { }

      /**  Destructor handles removal of the elements within the collection
       *
       * This must take into account the removal strategy and apply to any
       * remaining elements.
       */
      virtual ~CollectorMap() {
        selfDestruct();
      }

      /**
       * @brief Copy constructor invokes the copy policy as provided by the users
       *
       * This copy constructor will transfer the map of an incoming CollectorMap to
       * a newly created one.  This process employs the user selectable CopyPolicy.
       * It invokes the copy() method exposed in the copy policy.
       *
       * @param cmap The CollectorMap to be copied
       */
      CollectorMap(const CollectorMap &cmap) {
        _keyPolicy = cmap._keyPolicy;
        CollectorConstIter cItr;
        for(cItr = cmap._list.begin() ; cItr != cmap._list.end() ; cItr++) {
          _list.insert(std::make_pair(cItr->first, cmap.copy(cItr->second)));
        }
      }

      /**
       * @brief Assignment operator for the CollectorMap class object
       *
       * This object assignment operator is provided to properly handle the copying
       * of CollectorMap elements to a new instantiation.  This implements the
       * CopyPolicy for each element in the @b cmap object to the current one.  This
       * is a two step operation:  first destroy any elements that exist in the
       * destination object (using the RemovalPolicy) and then copy all elements
       * from the @b cmap object to the current one using the copy() method exposed
       * in the CopyPolicy.
       *
       * @param cmap The CollectorMap to be copied
       *
       * @returns Pointer to this CollectorMap
       */
      CollectorMap &operator=(const CollectorMap &cmap) {
        if(&cmap != this) {
          selfDestruct();
          _keyPolicy = cmap._keyPolicy;
          CollectorConstIter cItr;
          for(cItr = cmap._list.begin() ; cItr != cmap._list.end() ; cItr++) {
            _list.insert(std::make_pair(cItr->first, cmap.copy(cItr->second)));
          }
        }
        return (*this);
      }

      /**
       * Returns the size of the collection
       *
       * @return int Number of elements in collection
       */
      int size() const {
        return (_list.size());
      }

      /**
       * @brief Returns the number of keys found in the list
       *
       * For unique keys, this will always be 1.  If duplicate keys are allowed,
       * this will return the number of keys in the container.
       *
       * @param key  Key to return count for
       *
       * @return int Number keys in container
       */
      int count(const K &key) const {
        return (_list.count(key));
      }

      /**
       * Adds the element to the list.
       *
       * If the element exists and the key policy is restricted to uniqueness, it is
       * replaced after the removal strategy is applied.  If it doesn't exist, it is
       * inserted into the list.  For duplicate keys, it is simply inserted.
       *
       * @param key  Key in the associative map for the value
       * @param value Value to be associated with the key
       */
      void add(const K &key, const T &value) {
        if(_keyPolicy == UniqueKeys) remove(key);
        _list.insert(std::make_pair(key, value));
        return;
      }

      /**
       * Checks the existance of a particular key in the list
       * @param key Key to search for in the list
       * @return bool True if the key exists, false otherwise
       */
      bool exists(const K &key) const {
        CollectorConstIter cItr = _list.find(key);
        return (cItr != _list.end());
      }

      /**
       * @brief Returns the value associated with the name provided
       *
       * If the specifed name and value does not exist in the list, an out_of_range
       * exception is thrown.  Use @b exists to predetermine of the value is in the
       * list.
       *
       * @param key Key to fetch the value for
       * @return T Value associated with name
       * @throws IException if the value is not found
       */
      T &get(const K &key) {
        CollectorIter cItr = _list.find(key);
        if(cItr == _list.end()) {
          QString mess = "Requested value does not exist!";
          throw IException(IException::Programmer, mess, _FILEINFO_);
        }
        return (cItr->second);
      }

      /**
       * @brief Const version returning the value associated with the given name
       *
       * @param key Key to fetch the value for
       *
       * @returns The value associated with the given name
       */
      const T &get(const K &key) const {
        CollectorConstIter cItr = _list.find(key);
        if(cItr == _list.end()) {
          QString mess = "Requested value does not exist!";
          throw IException(IException::Programmer, mess, _FILEINFO_);
        }
        return (cItr->second);
      }

      /**
       * @brief Returns the index of the first occuring element in the list
       *
       * This returns the index such that the getNth() methods would retrieve the
       * element with key.  For duplicate keys, it is garaunteed to return the first
       * element.  It will return -1 if the element is not in the list.
       *
       * @param key Key to fetch the value for
       *
       * @return int Zero-based index of (first) element with key.  If it doesn't
       *         exist, -1 is returned.
       */
      int index(const K &key) const {
        CollectorConstIter cItr = _list.lower_bound(key);
        if(cItr == _list.end()) {
          return (-1);
        }
        else {
          return (std::distance(_list.begin(), cItr));
        }
      }

      /**
       * @brief Returns the nth value in the collection
       *
       * If the specifed value does not exist in the list, an out_of_range exception
       * is thrown.  Use @b size() to predetermine if the range is valid.
       *
       * @param nth Return the Nth value in the list
       *
       * @return T Value associated with name
       */
      T &getNth(int nth) {
        CollectorIter cItr;
        int i;
        for(cItr = _list.begin(), i = 0 ; cItr != _list.end() ; ++cItr, i++) {
          if(i == nth) break;
        }

        if(cItr == _list.end()) {
          std::ostringstream mess;
          mess << "Requested index (" << nth << ") out of range" << std::endl;
          throw IException(IException::Programmer, mess.str(), _FILEINFO_);
        }
        return (cItr->second);
      }

      /**
        * @brief Returns the nth value in the collection
        *
        * If the specifed value does not exist in the list, an out_of_range exception
        * is thrown.  Use @b size() to predetermine if the range is valid.
        *
        * @param nth Return the Nth value in the list
        *
        * @return T Value associated with name
        */
      const T &getNth(int nth) const {
        CollectorConstIter cItr;
        int i;
        for(cItr = _list.begin(), i = 0 ; cItr != _list.end() ; ++cItr, i++) {
          if(i == nth) break;
        }
        if(cItr == _list.end()) {
          std::ostringstream mess;
          mess << "Requested index (" << nth << ") out of range" << std::endl;
          throw IException(IException::Programmer, mess.str(), _FILEINFO_);
        }
        return (cItr->second);
      }

      /**
        * @brief Returns the nth key in the collection
        *
        * If the specifed key does not exist in the list, an out_of_range exception
        * is thrown.  Use @b size() to predetermine if the range is valid.
        *
        * @param nth Return the Nth key in the list
        *
        * @return K Key associated with name
        */
      const K &key(int nth) const {
        CollectorConstIter cItr;
        int i;
        for(cItr = _list.begin(), i = 0 ; cItr != _list.end() ; ++cItr, i++) {
          if(i == nth) break;
        }
        if(cItr == _list.end()) {
          std::ostringstream mess;
          mess << "Requested key index (" << nth << ") out of range" << std::endl;
          throw IException(IException::Programmer, mess.str(), _FILEINFO_);
        }
        return (cItr->first);
      }

      /**
       * Removes and entry from the list
       *
       * @param key Name of key/value pair to remove from the list
       *
       * @return int Number of elements erased
       */
      int remove(const K &key) {
        CollectorIter Itr1 = _list.lower_bound(key);
        if(Itr1 == _list.end()) return (0);

        CollectorIter Itr2 = _list.upper_bound(key);
        while(Itr1 != Itr2) {
          this->destroy(&Itr1->second);
          ++Itr1;
        }
        return (_list.erase(key));
      }

      /**
       * Const iterator into list
       *
       * @return CollectorConstIter Returns a const iterator to the list
       */
      CollectorConstIter begin() const {
        return _list.begin();
      }

      /**
       * Const iterator to end of list
       *
       * @return CollectorConstIter  Returns the const end of the list
       */
      CollectorConstIter end() const {
        return _list.end();
      }

      /**
       * Returns the start of the list for iterating purposes
       *
       * @return CollectorIter Returns an iterator on the collection
       */
      CollectorIter begin() {
        return _list.begin();
      }

      /**
       * Returns the end of the list
       *
       * @return CollectorIter Returns the end of the list for determining the end
       *         of the iteration loop
       */
      CollectorIter end() {
        return _list.end();
      }

    private:
      KeyPolicy      _keyPolicy;  //!<  Unique or duplicate key constraint
      CollectorList  _list;       //!< The list

      /**
       * @brief Thourough destruction of list
       *
       * This method iterates through each element in the list applying the
       * RemovalPolicy to each value in the map.  It then clears the internal list
       * for subsequent reuse if needed.
       */
      void selfDestruct() {
        CollectorIter itr;
        for(itr = _list.begin() ; itr != _list.end() ; itr++) {
          this->destroy(&itr->second);
        }
        _list.clear();
      }

  };

};
#endif


text:
#ifndef __LINALG_LIB_WRAPPER_H__
#define __LINALG_LIB_WRAPPER_H__

// Wrapper for linear algebra library (BLAS, LAPACK)

#if !defined(USE_MKL) && !defined(USE_OPENBLAS)
#define USE_OPENBLAS
#endif

#ifdef USE_MKL
#include <mkl.h>
#define BLAS_SET_NUM_THREADS mkl_set_num_threads
#endif

#ifdef USE_OPENBLAS
#include <cblas.h>
#include <lapacke.h>
#define BLAS_SET_NUM_THREADS openblas_set_num_threads
#endif

#endif


text:
/*
   Copyright [2021] [IBM Corporation]
   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.
*/

#ifndef _MCAS_COMMON_BYTE_
#define _MCAS_COMMON_BYTE_

#include <cstddef>
#include <gsl/gsl_byte>

#ifndef MCAS_BYTE_USES_STD
/* For compilation with C++14 use gsl::byte, not C++17 std::byte */
#define MCAS_BYTE_USES_STD 0
#endif

namespace common
{
#if MCAS_BYTE_USES_STD
	using byte = std::byte;
#else
	using byte = gsl::byte; /* can be std::byte in C++17 */
#endif
}

#endif

max_stars_repo_name: IBM/artemis

max_issues_repo_name: IBM/artemis

max_forks_repo_name: IBM/artemis


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void read_matrix(int** index, int** matrix, double scaling, int N_kw, char* input_fileName)
{
    FILE *fp = fopen(input_fileName, "r");
    fscanf(fp, "%*[^\n]\n");

    for (int ii = 0; ii < N_kw; ii++)
        fscanf(fp, "%d,", &((*index)[ii]));
    fscanf(fp, "%*[^\n]\n");

    int tmp;
    for (int ii = 0; ii < N_kw; ii++)
    {
        for (int jj = 0; jj < N_kw; jj++)
        {
            fscanf(fp, "%d,", &tmp);
            (*matrix)[ii*N_kw + jj] = (int) (scaling * tmp);
        }
        fscanf(fp, "%*[^\n]\n");
    }
    fclose(fp);
}

void pad_matrix(int** matrix_padded, int** matrix, int N_kw, int N_doc, int freq_max)
{
    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    // perform padding on the keywords
    int ii, jj;
    #pragma omp parallel for private(ii)
    for (int ii = 0; ii < N_kw; ii++)
        (*matrix_padded)[ii*N_kw + ii] = 2*freq_max;

    // perform padding
    #pragma omp parallel for private(ii, jj)
    for (ii = 0; ii < N_kw; ii++)
    {
        for (jj = 0; jj < N_kw; jj++)
        {
            if (ii > jj)
            {
                int n1 = 2*freq_max - (*matrix)[ii*N_kw + ii];
                int n2 = 2*freq_max - (*matrix)[jj*N_kw + jj];
                int x1 = gsl_ran_hypergeometric(r, n1, 2*N_doc-n1, (*matrix_padded)[jj*N_kw + jj]);
                int x2 = gsl_ran_hypergeometric(r, n2, 2*N_doc-n2, (*matrix_padded)[ii*N_kw + ii]);

                (*matrix_padded)[ii*N_kw + jj] = (*matrix)[ii*N_kw + jj] + x1 + x2;
                (*matrix_padded)[jj*N_kw + ii] = (*matrix_padded)[ii*N_kw + jj];
            }
        }
    }
    gsl_rng_free(r);
}

void observe_matrix(gsl_matrix* matrix_obs, int** matrix_padded, int N_kw)
{
    // perform observed count generation
    for (int ii = 0; ii < N_kw; ii++)
        for (int jj = 0; jj < N_kw; jj++)
            gsl_matrix_set(matrix_obs, ii, jj, (double) ((*matrix_padded)[ii*N_kw + jj]));
}

void permutation_generation(int* idx1, int* idx2, int** permutation_tmp, int** permutation, int** permutation_inv, int N_kw, int N_obs)
{
    *idx1 = rand() % N_obs;
    *idx2 = -1;
    int idx_old = (*permutation)[*idx1];
    int idx_new = rand() % N_kw;

    (*permutation_tmp)[*idx1] = idx_new;

    if ((*permutation_inv)[idx_new] >= 0)
    {
        *idx2 = (*permutation_inv)[idx_new];
        (*permutation_tmp)[*idx2] = idx_old;
    }
}

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
//Gets deltas (first order differences) of X,
//which is assumed to be pre-allocated of double size
//so that the appended rows/cols will hold the deltas.
//Thus, for dim==0, C must be even, with C/2 the original ncols of X.
//And,  for dim==1, R must be even, with R/2 the original nrows of X.

//I implement this just like FIR for speed and shorter code,
//except non-causal and mid-sample of B is 0,
//so I don't explicitly make B (e.g., B[n] just equals sc*n).

//Note that this may treat edge samples differently than other code.
//But this could be changed with a few lines of additional code here,
//without changing the super-efficient FIR implementation.
//That is, since out-of-range samps were assumed to be 0, nothing was
//added to Y for them, so can just add something later.

#include <stdio.h>
#include <cblas.h>

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int add_deltas_s (float *X, const int iscolmajor, const int R, const int C, const int dim, const int N);
int add_deltas_d (double *X, const int iscolmajor, const int R, const int C, const int dim, const int N);

int add_deltas_s (float *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const float z = 0.0f;
    const int No = R*C/2;
    int r, c, n;
    float sc = 1.0f;

    //Checks
    if (R<1) { fprintf(stderr,"error in add_deltas_s: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in add_deltas_s: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in add_deltas_s: N (delta winlength) must be positive\n"); return 1; }
    if (dim==0 && C%2!=0) { fprintf(stderr,"error in add_deltas_s: C (ncols X) must be even for dim==0\n"); return 1; }
    if (dim==1 && R%2!=0) { fprintf(stderr,"error in add_deltas_s: R (nrows X) must be even for dim==1\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5f/sc;

    if (dim==0)
    {
        if (iscolmajor)
        {
            cblas_scopy(No,&z,0,&X[No],1);
            for (c=0; c<C/2; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c*R],0,&X[No+c*R],1);        //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c*R],1,&X[No+c*R+n],1);    //past samps
                    cblas_saxpy(R-n,sc*n,&X[c*R+n],1,&X[No+c*R],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[c*R+R-1],0,&X[No+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C/2; c++)
            {
                cblas_scopy(R,&z,0,&X[C/2+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c],0,&X[C/2+c],C);                //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c],C,&X[C/2+c+n*C],C);          //past samps
                    cblas_saxpy(R-n,sc*n,&X[c+n*C],C,&X[C/2+c],C);           //future samps
                    cblas_saxpy(n,sc*n,&X[c+C*(R-1)],0,&X[C/2+c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R/2; r++)
            {
                cblas_scopy(C,&z,0,&X[R/2+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r],0,&X[R/2+r],R);                //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r],R,&X[R/2+r+n*R],R);          //past samps
                    cblas_saxpy(C-n,sc*n,&X[r+n*R],R,&X[R/2+r],R);           //future samps
                    cblas_saxpy(n,sc*n,&X[r+R*(C-1)],0,&X[R/2+r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            cblas_scopy(No,&z,0,&X[No],1);
            for (r=0; r<R/2; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r*C],0,&X[No+r*C],1);        //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r*C],1,&X[No+r*C+n],1);    //past samps
                    cblas_saxpy(C-n,sc*n,&X[r*C+n],1,&X[No+r*C],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[r*C+C-1],0,&X[No+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in add_deltas_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int add_deltas_d (double *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const double z = 0.0;
    const int No = R*C/2;
    int r, c, n;
    double sc = 1.0;

    //Checks
    if (R<1) { fprintf(stderr,"error in add_deltas_d: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in add_deltas_d: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in add_deltas_d: N (delta winlength) must be positive\n"); return 1; }
    if (dim==0 && C%2!=0) { fprintf(stderr,"error in add_deltas_d: C (ncols X) must be even for dim==0\n"); return 1; }
    if (dim==1 && R%2!=0) { fprintf(stderr,"error in add_deltas_d: R (nrows X) must be even for dim==1\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5/sc;

    if (dim==0)
    {
        if (iscolmajor)
        {
            cblas_dcopy(No,&z,0,&X[No],1);
            for (c=0; c<C/2; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c*R],0,&X[No+c*R],1);        //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c*R],1,&X[No+c*R+n],1);    //past samps
                    cblas_daxpy(R-n,sc*n,&X[c*R+n],1,&X[No+c*R],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[c*R+R-1],0,&X[No+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C/2; c++)
            {
                cblas_dcopy(R,&z,0,&X[C/2+c],C);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c],0,&X[C/2+c],C);                //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c],C,&X[C/2+c+n*C],C);          //past samps
                    cblas_daxpy(R-n,sc*n,&X[c+n*C],C,&X[C/2+c],C);           //future samps
                    cblas_daxpy(n,sc*n,&X[c+C*(R-1)],0,&X[C/2+c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R/2; r++)
            {
                cblas_dcopy(C,&z,0,&X[R/2+r],R);
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r],0,&X[R/2+r],R);                //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r],R,&X[R/2+r+n*R],R);          //past samps
                    cblas_daxpy(C-n,sc*n,&X[r+n*R],R,&X[R/2+r],R);           //future samps
                    cblas_daxpy(n,sc*n,&X[r+R*(C-1)],0,&X[R/2+r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            cblas_dcopy(No,&z,0,&X[No],1);
            for (r=0; r<R/2; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r*C],0,&X[No+r*C],1);        //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r*C],1,&X[No+r*C+n],1);    //past samps
                    cblas_daxpy(C-n,sc*n,&X[r*C+n],1,&X[No+r*C],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[r*C+C-1],0,&X[No+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in add_deltas_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#ifndef QDM_ISPLINE_H
#define QDM_ISPLINE_H 1

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

void
qdm_ispline_vector(
    gsl_vector *result,
    const double tau,
    const size_t spline_df,
    const gsl_vector *knots
);

void
qdm_ispline_matrix(
    gsl_matrix *result,
    const gsl_vector *taus,
    const size_t spline_df,
    const gsl_vector *knots
);

#endif /* QDM_ISPLINE_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

#pragma comment(lib,"gsl")

int main()
{
	printf("%g\n",gsl_sf_bessel_J0(5.0));
	return 0;
}


text:
#pragma once

#include <Windows.h>

#include <cstdint>
#include <string>
#include <vector>

#include <gsl/span>
#include "hid_handle.h"

namespace hid
{
	struct HidOpenFlags
	{
		enum T : uint32_t
		{
			none,
			exclusive = 1u << 0u,
			async     = 1u << 1u
		};
	};

	using HidOpenFlags_t = uint32_t;

	struct HidCaps
	{
		uint16_t usage;
		uint16_t usagePage;
		uint16_t inputReportSize;
		uint16_t outputReportSize;
		uint16_t featureReportSize;
		uint16_t linkCollectionNodes;
		uint16_t inputButtonCaps;
		uint16_t inputValueCaps;
		uint16_t inputDataIndices;
		uint16_t outputButtonCaps;
		uint16_t outputValueCaps;
		uint16_t outputDataIndices;
		uint16_t featureButtonCaps;
		uint16_t featureValueCaps;
		uint16_t featureDataIndices;
	};

	struct HidAttributes
	{
		uint16_t vendorId;
		uint16_t productId;
		uint16_t versionNumber;
	};

	class HidInstance
	{
		HidOpenFlags_t flags = 0;

		Handle handle = Handle(nullptr, true);

		HidCaps caps_ {};
		HidAttributes attributes_ {};

		OVERLAPPED overlappedIn = {};
		OVERLAPPED overlappedOut = {};

		bool pendingRead_ = false;
		bool pendingWrite_ = false;

		size_t nativeError_ = 0;

	public:
		std::wstring path;
		std::wstring instanceId;
		std::wstring serialString;

		std::vector<uint8_t> inputBuffer;
		std::vector<uint8_t> outputBuffer;

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

		HidInstance() = default;

		HidInstance(std::wstring path, std::wstring instanceId);
		explicit HidInstance(std::wstring path);
		HidInstance(HidInstance&& other) noexcept;

		~HidInstance();

		HidInstance& operator=(HidInstance&& other) noexcept;

		bool isOpen() const;
		bool isExclusive() const;
		bool isAsync() const;
		const HidCaps& caps() const;
		const HidAttributes& attributes() const;

		bool readMetadata();
		bool readCaps();
		bool readSerial();
		bool readAttributes();
		bool getFeature(const gsl::span<uint8_t>& buffer) const;
		bool setFeature(const gsl::span<uint8_t>& buffer) const;

		bool open(HidOpenFlags_t openFlags);
		void close();

		inline auto nativeError() const
		{
			return nativeError_;
		}

		bool read(void* buffer, size_t size) const;
		bool read(const gsl::span<uint8_t>& buffer) const;
		bool read();

		bool readAsync();

		bool write(const void* buffer, size_t size) const;
		bool write(const gsl::span<const uint8_t>& buffer) const;
		bool write() const;

		bool writeAsync();

		bool asyncReadPending() const;
		bool asyncReadInProgress();
		bool asyncWritePending() const;
		bool asyncWriteInProgress();

		void cancelAsyncReadAndWait();
		void cancelAsyncWriteAndWait();

		bool setOutputReport(const gsl::span<uint8_t>& buffer) const;
		bool setOutputReport();

	private:
		void cancelAsyncAndWait(OVERLAPPED* overlapped);
		bool asyncInProgress(OVERLAPPED* overlapped);

		bool readCaps(HANDLE h);
		bool readSerial(HANDLE h);
		bool readAttributes(HANDLE h);
	};
}


text:
/* interpolation/gsl_interp.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman
 */
#ifndef __GSL_INTERP_H__
#define __GSL_INTERP_H__
#include <stdlib.h>
#include <gsl/gsl_inline.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

/* evaluation accelerator */
typedef struct {
  size_t  cache;        /* cache of index   */
  size_t  miss_count;   /* keep statistics  */
  size_t  hit_count;
}
gsl_interp_accel;

/* interpolation object type */
typedef struct {
  const char * name;
  unsigned int min_size;
  void *  (*alloc) (size_t size);
  int     (*init)    (void *, const double xa[], const double ya[], size_t size);
  int     (*eval)    (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y);
  int     (*eval_deriv)  (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
  int     (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
  int     (*eval_integ)  (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
  void    (*free)         (void *);

} gsl_interp_type;

/* general interpolation object */
typedef struct {
  const gsl_interp_type * type;
  double  xmin;
  double  xmax;
  size_t  size;
  void * state;
} gsl_interp;

/* available types */
GSL_VAR const gsl_interp_type * gsl_interp_linear;
GSL_VAR const gsl_interp_type * gsl_interp_polynomial;
GSL_VAR const gsl_interp_type * gsl_interp_cspline;
GSL_VAR const gsl_interp_type * gsl_interp_cspline_periodic;
GSL_VAR const gsl_interp_type * gsl_interp_akima;
GSL_VAR const gsl_interp_type * gsl_interp_akima_periodic;

gsl_interp_accel *
gsl_interp_accel_alloc(void);

int
gsl_interp_accel_reset (gsl_interp_accel * a);

void
gsl_interp_accel_free(gsl_interp_accel * a);

gsl_interp *
gsl_interp_alloc(const gsl_interp_type * T, size_t n);

int
gsl_interp_init(gsl_interp * obj, const double xa[], const double ya[], size_t size);

const char * gsl_interp_name(const gsl_interp * interp);
unsigned int gsl_interp_min_size(const gsl_interp * interp);

int
gsl_interp_eval_e(const gsl_interp * obj,
                  const double xa[], const double ya[], double x,
                  gsl_interp_accel * a, double * y);

double
gsl_interp_eval(const gsl_interp * obj,
                const double xa[], const double ya[], double x,
                gsl_interp_accel * a);

int
gsl_interp_eval_deriv_e(const gsl_interp * obj,
                        const double xa[], const double ya[], double x,
                        gsl_interp_accel * a,
                        double * d);

double
gsl_interp_eval_deriv(const gsl_interp * obj,
                      const double xa[], const double ya[], double x,
                      gsl_interp_accel * a);

int
gsl_interp_eval_deriv2_e(const gsl_interp * obj,
                         const double xa[], const double ya[], double x,
                         gsl_interp_accel * a,
                         double * d2);

double
gsl_interp_eval_deriv2(const gsl_interp * obj,
                       const double xa[], const double ya[], double x,
                       gsl_interp_accel * a);

int
gsl_interp_eval_integ_e(const gsl_interp * obj,
                        const double xa[], const double ya[],
                        double a, double b,
                        gsl_interp_accel * acc,
                        double * result);

double
gsl_interp_eval_integ(const gsl_interp * obj,
                      const double xa[], const double ya[],
                      double a, double b,
                      gsl_interp_accel * acc);

void
gsl_interp_free(gsl_interp * interp);

INLINE_DECL size_t
gsl_interp_bsearch(const double x_array[], double x,
                   size_t index_lo, size_t index_hi);

#ifdef HAVE_INLINE

/* Perform a binary search of an array of values.
 *
 * The parameters index_lo and index_hi provide an initial bracket,
 * and it is assumed that index_lo < index_hi. The resulting index
 * is guaranteed to be strictly less than index_hi and greater than
 * or equal to index_lo, so that the implicit bracket [index, index+1]
 * always corresponds to a region within the implicit value range of
 * the value array.
 *
 * Note that this means the relationship of 'x' to x_array[index]
 * and x_array[index+1] depends on the result region, i.e. the
 * behaviour at the boundaries may not correspond to what you
 * expect. We have the following complete specification of the
 * behaviour.
 * Suppose the input is x_array[] = { x0, x1, ..., xN }
 *    if ( x == x0 )           then  index == 0
 *    if ( x > x0 && x <= x1 ) then  index == 0, and sim. for other interior pts
 *    if ( x == xN )           then  index == N-1
 *    if ( x > xN )            then  index == N-1
 *    if ( x < x0 )            then  index == 0
 */

INLINE_FUN size_t
gsl_interp_bsearch(const double x_array[], double x,
                   size_t index_lo, size_t index_hi)
{
  size_t ilo = index_lo;
  size_t ihi = index_hi;
  while(ihi > ilo + 1) {
    size_t i = (ihi + ilo)/2;
    if(x_array[i] > x)
      ihi = i;
    else
      ilo = i;
  }

  return ilo;
}
#endif

INLINE_DECL size_t
gsl_interp_accel_find(gsl_interp_accel * a, const double x_array[], size_t size, double x);

#ifdef HAVE_INLINE
INLINE_FUN size_t
gsl_interp_accel_find(gsl_interp_accel * a, const double xa[], size_t len, double x)
{
  size_t x_index = a->cache;

  if(x < xa[x_index]) {
    a->miss_count++;
    a->cache = gsl_interp_bsearch(xa, x, 0, x_index);
  }
  else if(x >= xa[x_index + 1]) {
    a->miss_count++;
    a->cache = gsl_interp_bsearch(xa, x, x_index, len-1);
  }
  else {
    a->hit_count++;
  }

  return a->cache;
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_INTERP_H__ */

max_forks_repo_name: gersteinlab/LESSeq


text:
#ifndef RANDOM_H
#define RANDOM_H

extern "C" {
#include <gsl/gsl_rng.h>
#include <gsl/gsl_const_mksa.h>
}

#include <math.h>
#include <mkl.h>
#include <iostream>
#include "matrix.h"
#include "file_IO.h"

// note that the RANDOM package should not has the dependency with lib_traj.h
class RNG_BOOST
{
 public:
  MKL_LONG N_THREADS_BD;
  MKL_LONG N_THREADS_SS;
  gsl_rng **BOOST_BD;
  gsl_rng **BOOST_SS;
  bool INITIALIZATION_BD;
  bool INITIALIZATION_SS;
  RNG_BOOST()
    {
      std::cout<< "ERR: basic constructor for RNG_BOOST class is not allowed\n";
    }
  RNG_BOOST
    (COND& given_condition)
    {
      INITIALIZATION_BD = TRUE;
      const gsl_rng_type *T_boost;
      N_THREADS_BD = atol(given_condition("N_THREADS_BD").c_str());
      BOOST_BD = new gsl_rng* [N_THREADS_BD];
      gsl_rng_env_setup();
      T_boost = gsl_rng_default;
      for(MKL_LONG i=0; i<N_THREADS_BD; i++)
        {
          BOOST_BD[i] = gsl_rng_alloc(T_boost);
          MKL_LONG seed_BD = atol(given_condition("basic_random_seed").c_str());
          gsl_rng_set(BOOST_BD[i], seed_BD + i); // it set the seed with index i
        }
      INITIALIZATION_SS = FALSE; // default
      if(given_condition("Method") == "NAPLE_ASSOCIATION")
	if(given_condition("Step") != "EQUILIBRATION")
	 /* given_condition("Step")!="EQUILIBRATION" && given_condition("Method")=="NAPLE_ASSOCIATION") */
        // identify the NAPLE_ASSOCIATION topological evolution with equilibration condtion or not
	  {
	    INITIALIZATION_SS = TRUE;
	    N_THREADS_SS = atol(given_condition("N_THREADS_SS").c_str());
	    /* BOOST_SS = (gsl_rng**)mkl_malloc(N_THREADS_SS*sizeof(gsl_rng*), BIT); */
	    BOOST_SS = new gsl_rng* [N_THREADS_SS];
	    for(MKL_LONG i=0; i<N_THREADS_SS; i++)
	      {
		BOOST_SS[i] = gsl_rng_alloc(T_boost);
		MKL_LONG seed_SS = atol(given_condition("basic_random_seed_SS").c_str());
		gsl_rng_set(BOOST_SS[i], seed_SS + i);
	      }
        }
    }

  ~RNG_BOOST()
    {
      if(INITIALIZATION_BD)
        {
          for(MKL_LONG i=0; i<N_THREADS_BD; i++)
            {
              gsl_rng_free(BOOST_BD[i]);
            }
          delete[] BOOST_BD;
          if(INITIALIZATION_SS)
            {
              for(MKL_LONG i=0; i<N_THREADS_SS; i++)
                {
                  gsl_rng_free(BOOST_SS[i]);
                }
              delete[] BOOST_SS;
            }
        }
    }
};

namespace RANDOM
{
  MKL_LONG
    random_vector_generator
    (MATRIX& R_VEC_TRANS);
  MKL_LONG
    single_random_vector_generator
    (MATRIX& given_vec);
  MKL_LONG
    single_random_vector_generator_variance
    (MATRIX& given_vec, double s_2);
  MKL_LONG
    single_random_vector_generator_boost
    (MATRIX& given_vec, gsl_rng* r_boost);
  MKL_LONG
    single_random_vector_generator_variance_boost
    (MATRIX& given_vec, double s_2, gsl_rng* r_boost);
  MKL_LONG
    single_unit_random_vector_generator
    (MATRIX& given_vec);
  MKL_LONG
    unit_random_vector_generator
    (MATRIX& R_VEC_TRANS);
  MKL_LONG
    unit_random_vector_generator_2D
    (MATRIX& R_VEC_TRANS);

  MKL_LONG
    return_LONG_INT_rand
    (MKL_LONG SUP);
  MKL_LONG
    return_LONG_INT_rand_boost
    (gsl_rng* r, MKL_LONG SUP);
  double
    return_double_rand_SUP1();
  double
    return_double_rand_SUP1_boost
    (gsl_rng* r);
  MKL_LONG
    get_LONG_ARR_rand_boost
    (gsl_rng* r, MKL_LONG SUP, MKL_LONG* given_long_arr, MKL_LONG N_arr);
  MKL_LONG
    get_DOUBLE_ARR_rand_boost
    (gsl_rng* r, double* given_double_arr, MKL_LONG N_arr);

}

#endif


text:
#ifndef RANDOM_H
#define RANDOM_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

extern gsl_rng *rng_state;

void rand_init(int);
int rand_store(void*);
void rand_load(void*, int);

#define SQRT2INV 0.70710678118654752440084436210
#define rand_uniform() gsl_rng_uniform_pos(rng_state)
#define rand_gaussian() gsl_ran_gaussian_ziggurat(rng_state, SQRT2INV)
#define rand_z2() ((rand_uniform() < 0.5) ? SQRT2INV : -SQRT2INV)

#endif

max_stars_repo_name: mrlhansen/YaHMC

max_issues_repo_name: mrlhansen/YaHMC

max_forks_repo_name: mrlhansen/YaHMC


text:
/*
 * optimize_tnc.h
 *
 *  Created on: Feb 9, 2010
 *      Author: smitty
 */

#ifndef _OPTIMIZE_STATE_RECONSTRUCTOR_PERIODS_NLOPT_H_
#define _OPTIMIZE_STATE_RECONSTRUCTOR_PERIODS_NLOPT_H_

#include <nlopt.h>

#include "state_reconstructor.h"
#include "rate_model.h"

#include <armadillo>
using namespace arma;

void optimize_sr_periods_nlopt(vector<RateModel> * _rm,StateReconstructor * _sr,
    vector<mat> * _free_mask, int _nfree);

#endif /* _OPTIMIZE_STATE_RECONSTRUCTOR_PERIODS_NLOPT_H_ */


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

/*!
 * @brief Converts v-w coordinates to lune coordinates.
 *
 * @param[in] nv      Number of v coordinates.
 * @param[in] v       v coordinates s.t. \f$ v \in [-1/3, 1/3] \f$.  This
 *                    an array of dimension [nv].
 * @param[in] nw      Number of w coordinates.
 * @param[in] w       w coordinates (similar to lune latitudes) s.t.
 *                    \f$ w \in [-3\pi/8, 3/pi/8] \f$.  This is an
 *                    array of dimension [nw].
 *
 * @param[out] gamma  Longitude (degrees) \f$ \gamma \in [-30,30] \f$.
 *                    This is an array of dimension [nv].
 * @param[out] delta  lune latitude (degrees) \f$ \delta \in [-90,90] \f$.
 *                    This is an array of dimension [nw].
 *
 * @result 0 indicates success.
 *
 * @author Carl Tape and converted to C by Ben Baker
 *
 * @date 2016
 *
 * @copyright MIT
 *
 */
int compearth_rect2lune(const int nv, const double *__restrict__ v,
                        const int nw, const double *__restrict__ w,
                        double *__restrict__ gamma,
                        double *__restrict__ delta)

{
    double u[CE_CHUNKSIZE] __attribute__((aligned(64)));
    int i, j, ierr, nuLoc;
/*
    const int maxit = 20;
    const int useHalley = 2;
    const double tol = 1.e-12;
*/
    const double pi38 = 3.0*M_PI/8.0;
    const double pi180i = 180.0/M_PI;
    ierr = 0;
    // convert v -> gamma
    compearth_v2gamma(nv, v, gamma);
    // convert latitude to colatitude
    for (i=0; i<nw; i=i+CE_CHUNKSIZE)
    {
        nuLoc = MIN(CE_CHUNKSIZE, nw - i);
        for (j=0; j<nuLoc; j++)
        {
            u[j] = pi38 - w[i+j];
        }
        // convert u -> beta
        ierr = compearth_u2beta(nuLoc, u, &delta[i]);
        //ierr = compearth_u2beta(nuLoc, maxit, useHalley, u, tol, &delta[i]);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Computing u2beta\n", __func__);
            return -1;
        }
    }
    // convert to gamma to degrees
    cblas_dscal(nv, pi180i, gamma, 1);
    // convert delta from colatitude to latitude and radians to degrees
    #pragma omp simd
    for (i=0; i<nw; i++)
    {
        delta[i] = 90.0 - pi180i*delta[i];
    }
    return ierr;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
/* vector/gsl_vector_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, 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 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 __GSL_VECTOR_DOUBLE_H__
#define __GSL_VECTOR_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_double.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

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block *block;
  int owner;
}
gsl_vector;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_view;

typedef _gsl_vector_view gsl_vector_view;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_const_view;

typedef const _gsl_vector_const_view gsl_vector_const_view;

/* Allocation */

GSL_EXPORT gsl_vector *gsl_vector_alloc (const size_t n);
GSL_EXPORT gsl_vector *gsl_vector_calloc (const size_t n);

GSL_EXPORT gsl_vector *gsl_vector_alloc_from_block (gsl_block * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_EXPORT gsl_vector *gsl_vector_alloc_from_vector (gsl_vector * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_EXPORT void gsl_vector_free (gsl_vector * v);

/* Views */

GSL_EXPORT
_gsl_vector_view
gsl_vector_view_array (double *v, size_t n);

GSL_EXPORT
_gsl_vector_view
gsl_vector_view_array_with_stride (double *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_const_view
gsl_vector_const_view_array (const double *v, size_t n);

GSL_EXPORT
_gsl_vector_const_view
gsl_vector_const_view_array_with_stride (const double *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_view
gsl_vector_subvector (gsl_vector *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_view
gsl_vector_subvector_with_stride (gsl_vector *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_const_view
gsl_vector_const_subvector (const gsl_vector *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_const_view
gsl_vector_const_subvector_with_stride (const gsl_vector *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT double gsl_vector_get (const gsl_vector * v, const size_t i);
GSL_EXPORT void gsl_vector_set (gsl_vector * v, const size_t i, double x);

GSL_EXPORT double *gsl_vector_ptr (gsl_vector * v, const size_t i);
GSL_EXPORT const double *gsl_vector_const_ptr (const gsl_vector * v, const size_t i);

GSL_EXPORT void gsl_vector_set_zero (gsl_vector * v);
GSL_EXPORT void gsl_vector_set_all (gsl_vector * v, double x);
GSL_EXPORT int gsl_vector_set_basis (gsl_vector * v, size_t i);

GSL_EXPORT int gsl_vector_fread (FILE * stream, gsl_vector * v);
GSL_EXPORT int gsl_vector_fwrite (FILE * stream, const gsl_vector * v);
GSL_EXPORT int gsl_vector_fscanf (FILE * stream, gsl_vector * v);
GSL_EXPORT int gsl_vector_fprintf (FILE * stream, const gsl_vector * v,
                                   const char *format);

GSL_EXPORT int gsl_vector_memcpy (gsl_vector * dest, const gsl_vector * src);

GSL_EXPORT int gsl_vector_reverse (gsl_vector * v);

GSL_EXPORT int gsl_vector_swap (gsl_vector * v, gsl_vector * w);
GSL_EXPORT int gsl_vector_swap_elements (gsl_vector * v, const size_t i, const size_t j);

GSL_EXPORT double gsl_vector_max (const gsl_vector * v);
GSL_EXPORT double gsl_vector_min (const gsl_vector * v);
GSL_EXPORT void gsl_vector_minmax (const gsl_vector * v, double * min_out, double * max_out);

GSL_EXPORT size_t gsl_vector_max_index (const gsl_vector * v);
GSL_EXPORT size_t gsl_vector_min_index (const gsl_vector * v);
GSL_EXPORT void gsl_vector_minmax_index (const gsl_vector * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_add (gsl_vector * a, const gsl_vector * b);
GSL_EXPORT int gsl_vector_sub (gsl_vector * a, const gsl_vector * b);
GSL_EXPORT int gsl_vector_mul (gsl_vector * a, const gsl_vector * b);
GSL_EXPORT int gsl_vector_div (gsl_vector * a, const gsl_vector * b);
GSL_EXPORT int gsl_vector_scale (gsl_vector * a, const double x);
GSL_EXPORT int gsl_vector_add_constant (gsl_vector * a, const double x);

GSL_EXPORT int gsl_vector_isnull (const gsl_vector * v);

#ifdef HAVE_INLINE

extern inline
double
gsl_vector_get (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_set (gsl_vector * v, const size_t i, double x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
double *
gsl_vector_ptr (gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (double *) (v->data + i * v->stride);
}

extern inline
const double *
gsl_vector_const_ptr (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const double *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/* multimin/steepest_descent.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 *
 * 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.
 */

/* steepest_descent.c -- the steepest descent algorithm */

#include <config.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_blas.h>

typedef struct
{
  double step;
  double max_step;
  gsl_vector *x1;
  gsl_vector *g1;
}
steepest_descent_state_t;

static int
steepest_descent_alloc (void *vstate, size_t n)
{
  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;

  state->x1 = gsl_vector_alloc (n);

  if (state->x1 == NULL)
    {
      GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
    }

  state->g1 = gsl_vector_alloc (n);

  if (state->g1 == NULL)
    {
      gsl_vector_free (state->x1);
      GSL_ERROR ("failed to allocate space for g1", GSL_ENOMEM);
    }

  return GSL_SUCCESS;
}

static int
steepest_descent_set (void *vstate, gsl_multimin_function_fdf * fdf,
                      const gsl_vector * x, double *f,
                      gsl_vector * gradient, double step_size)
{
  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;

  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);

  state->step = step_size;
  state->max_step = step_size;

  return GSL_SUCCESS;
}

static void
steepest_descent_free (void *vstate)
{
  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;

  gsl_vector_free (state->x1);
  gsl_vector_free (state->g1);
}

static int
steepest_descent_restart (void *vstate)
{
  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;

  state->step = state->max_step;

  return GSL_SUCCESS;
}

static int
steepest_descent_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
                          gsl_vector * x, double *f,
                          gsl_vector * gradient, gsl_vector * dx)
{
  steepest_descent_state_t *state = (steepest_descent_state_t *) vstate;

  gsl_vector *x1 = state->x1;
  gsl_vector *g1 = state->g1;

  double f0 = *f;
  double f1, fm, f_trial;
  double step0 = 0.0, step1, stepm, step_trial, step = state->step;
  size_t iter = 0;

  /* compute new trial point at x1= x - step * dir, where dir is the
     normalized gradient */

  double gnorm = gsl_blas_dnrm2 (gradient);

  gsl_vector_set_zero (dx);
  gsl_blas_daxpy (-step / gnorm, gradient, dx);

  gsl_vector_memcpy (x1, x);
  gsl_blas_daxpy (1.0, dx, x1);

  /* evaluate function and gradient at new point x1 */

  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x1, &f1, g1);

  if (f1 < f0)
    {
      state->step = 2.0 * step;

      gsl_vector_memcpy (x, x1);
      gsl_vector_memcpy (gradient, g1);

      *f = f1;

      return GSL_SUCCESS;
    }

  step1 = step;

trial:

  stepm = 0.5 * step1;

  gsl_vector_set_zero (dx);
  gsl_blas_daxpy (-stepm / gnorm, gradient, dx);

  gsl_vector_memcpy (x1, x);
  gsl_blas_daxpy (1.0, dx, x1);

  fm = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
  printf("trying stepm = %.18e fm=%g\n", stepm, fm);
  if (fm >= f0)
    {
      /* downhill step failed, reduce step-size and try again */
      f1 = fm;
      step1 = stepm;
      goto trial;
    }

  /* We now have a triplet (0,f0) (stepm, fm) (step1, f1) */

minimize:
  iter++;

  if ((stepm - step0) > (step1 - stepm))
    {
      step_trial = stepm - 0.38 * (stepm - step0);
    }
  else
    {
      step_trial = stepm + 0.38 * (step1 - stepm);
    }

  gsl_vector_set_zero (dx);
  gsl_blas_daxpy (-step_trial / gnorm, gradient, dx);

  gsl_vector_memcpy (x1, x);
  gsl_blas_daxpy (1.0, dx, x1);

  f_trial = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);

  if (f_trial > fm)
    {
      if (step_trial < stepm)
        {
          step0 = step_trial;
          f0 = f_trial;
        }
      else
        {
          step1 = step_trial;
          f1 = f_trial;
        }
    }
  else
    {
      if (step_trial < stepm)
        {
          step1 = stepm;
          f1 = fm;
        }
      else
        {
          step0 = stepm;
          f0 = fm;
        }

      stepm = step_trial;
      fm = f_trial;
    }

  printf("f = %.18e\n", fm);

  if (iter > 100)
    {
      gsl_vector_set_zero (dx);
      gsl_blas_daxpy (-stepm / gnorm, gradient, dx);
      gsl_blas_daxpy (1.0, dx, x);

      GSL_MULTIMIN_FN_EVAL_DF (fdf, x, gradient);
      *f = fm;
      state->step = stepm;

      return GSL_SUCCESS;
    }

  goto minimize;
}

static const gsl_multimin_fdfminimizer_type steepest_descent_type =
  { "steepest_descent",         /* name */
  sizeof (steepest_descent_state_t),
  &steepest_descent_alloc,
  &steepest_descent_set,
  &steepest_descent_iterate,
  &steepest_descent_restart,
  &steepest_descent_free
};

const gsl_multimin_fdfminimizer_type
  *gsl_multimin_fdfminimizer_steepest_descent = &steepest_descent_type;

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

/*
 * =====================================================================================
 *
 *       Filename:  Expi.c
 *
 *    Description:  Calculation of the exponential integrals (w/o the GSL special
 *    			functions).
 *
 *        Version:  1.0
 *        Created:  15/04/2013 11:54:19
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:
 *
 * =====================================================================================
 */

#include "funcs.h"
#include <gsl/gsl_integration.h>
/* #include <gsl/gsl_sf_expint.h> */

double fu ( double t, void* params )
{
	double f = -exp(-t)/t ;
	return f ;
}

double ex ( double t, void* params )
{
	double e = exp(-t) ;
	return e ;
}

int expi ( double x, double* result, double* abserr )
{
	double r, err ;

	gsl_integration_workspace* expi_ws =
	       	gsl_integration_workspace_alloc (WS_SZ) ;

	gsl_function F ;
	F.function = &fu ;

	int status ;
	status = gsl_integration_qagiu ( &F, x, 10e-9, .001 , WS_SZ, expi_ws, &r,
				&err) ;

	*result = r ; *abserr = err ;

/*	Using the GSL special functions, it is simply:
 *
 *      *result = - gsl_sf_expint_E1(x) ;
 */

	gsl_integration_workspace_free (expi_ws) ;

	return status;
}

int expi_plus ( double x, double* result, double* abserr )
{
	double r, err ;

	gsl_integration_workspace *expi_ws =
	       	gsl_integration_workspace_alloc (WS_SZ) ;

	gsl_function F ;
	F.function = &ex ;

	int status ;
	status = gsl_integration_qawc ( &F, -x, x, 0, 1e-9, .001, WS_SZ,
			expi_ws, &r, &err ) ;

	double R , EXPI , ERREXPI ;
	int s = expi ( x , &EXPI, &ERREXPI ) ;
	R = - (r - EXPI) ;                          /* The minus (-) sign because of the def.
							of ex 				*/
	double ERR ;
	ERR = err + ERREXPI ;

	*result = R ;
	*abserr = ERR ;

/* 	Using the GSL exponential integral, it is simply:
 *
 * 	R = gsl_sf_expint_Ei(x)  ;
 *
 */
	gsl_integration_workspace_free (expi_ws) ;

	return status + s ;
}


text:
#ifndef __mcas_SERVER_TASK_KEY_FIND_H__
#define __mcas_SERVER_TASK_KEY_FIND_H__

#include <common/logging.h>
#include <gsl/pointers>
#include <unistd.h>
#include <string>

#include "task.h"

namespace mcas
{
/**
 * Key search task.  We limit the number of hops we search so as to bound
 * the worst case execution time.
 *
 */
class Key_find_task : public Shard_task,
                      private common::log_source
{
  static constexpr unsigned MAX_COMPARES_PER_WORK = 5;

 public:
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"  // uninitialized _expr, _out_key, _type
  Key_find_task(const std::string& expression,
                const offset_t offset,
                Connection_handler* handler,
                gsl::not_null<component::IKVIndex*> index,
                const unsigned debug_level)
      : Shard_task(handler),
        log_source(debug_level),
        _offset(offset),
        _index(index)
  {
    using namespace component;
    _index->add_ref();

    CPLOG(1, "offset=%lu", offset);
    CPLOG(1,"expr: (%s)", expression.c_str());

    if (expression == "next:") {
      _type = IKVIndex::FIND_TYPE_NEXT;
      _expr = expression.substr(5);
    }
    else if (expression.substr(0, 6) == "regex:") {
      _type = IKVIndex::FIND_TYPE_REGEX;
      _expr = expression.substr(6);
    }
    else if (expression.substr(0, 6) == "exact:") {
      _type = IKVIndex::FIND_TYPE_EXACT;
      _expr = expression.substr(6);
    }
    else if (expression.substr(0, 7) == "prefix:") {
      _type = IKVIndex::FIND_TYPE_PREFIX;
      _expr = expression.substr(7);
    }
    else
      throw Logic_exception("unhandled expression");

  }
#pragma GCC diagnostic pop

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

  status_t do_work() override
  {
    using namespace component;

    status_t hr;
    try {
      hr = _index->find(_expr, _offset, _type, _offset, _out_key, MAX_COMPARES_PER_WORK);

      if (hr == E_MAX_REACHED) {
        _offset++;
        return component::IKVStore::S_MORE;
      }
      else if (hr == S_OK) {
        CPLOG(2, "matched: (%s)", _out_key.c_str());
        return S_OK;
      }
      else {
        _out_key.clear();
        return hr;
      }
    }
    catch (...) {
      return E_FAIL;
    }

    throw Logic_exception("unexpected code path (hr=%d)", hr);
  }

  const void* get_result() const override { return _out_key.data(); }

  size_t get_result_length() const override { return _out_key.length(); }

  offset_t matched_position() const override { return _offset; }

 private:
  std::string                             _expr;
  std::string                             _out_key;
  component::IKVIndex::find_t             _type;
  offset_t                                _offset;
  component::Itf_ref<component::IKVIndex> _index;
};

}  // namespace mcas
#endif  // __mcas_SERVER_TASK_KEY_FIND_H__


text:
#include <stdio.h>
#include<stdlib.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

void *xcalloc(int items, int size)
{
	void *ptr = calloc(items, size);
	if (ptr == NULL)
	{
		printf("Unable to allocate memory\n");
		exit(EXIT_FAILURE);
	}
	return ptr;
}

void GetLGLPoints(int N, double *x)
{
	if (N == 0)
	{
		printf("Support for constants doesn't exist\n");
		exit(1);
	}
	else if (N == 1)
	{
		x[0] = -1.0;
		x[1] = 1.0;
	}
	else if (N==2)
	{
		x[0] = -1.0;
		x[1] = 0.0;
		x[2] = 1.0;
	}
	else if (N==3)
	{
		x[0] = -1.0;
		x[1] = -0.4472;
		x[2] = 0.4472;
		x[3] = 1.0;
	}
	else if (N==4)
	{
		x[0] = -1.0;
		x[1] = -0.6547;
		x[2] = 0;
		x[3] = 0.6547;
		x[4] = 1.0;
	}
	else
	{
		printf("Currently only have up to fourth order polynomial encoded for LGL points\n");
		exit(1);
	}
}

void GetLGLWeights(int N, double *w)
{
	if (N == 0)
	{
		printf("Support for constants doesn't exist\n");
		exit(1);
	}
	else if (N == 1)
	{
		w[0] = 1.0;
		w[1] = 1.0;
	}
	else if (N==2)
	{
		w[0] = 0.3333333333333;
		w[1] = 1.3333333333333;
		w[2] = 0.3333333333333;
	}
	else if (N==3)
	{
		w[0] = 0.16666666666667;
		w[1] = 0.83333333333333;
		w[2] = 0.83333333333333;
		w[3] = 0.16666666666667;
	}
	else if (N==4)
	{
		w[0] = 0.1;
		w[1] = 0.5444444;
		w[2] = 0.7111111;
		w[3] = 0.5444444;
		w[4] = 0.1;
	}
	else
	{
		printf("Currently only have up to fourth order polynomial encoded for LGL points\n");
		exit(1);
	}

}

void calculateNodalCoordinates(int N, int NumEl, double* X, double* Xnodes)
{
	double r[N+1];
	GetLGLPoints(N, r);
	for(int i = 0; i < NumEl; i++)
	{
		for (int j = 0; j < N+1; j++)
		{
			Xnodes[i*(N+1)+j] = X[i] + 0.5*(r[j]+1)*(X[i+1]-X[i]);
		}
	}

}

void LegendrePoly(double *x, int Np, int N, double* P)
{
	// Purpose: Evaluates the normalized Legendre polynomials
	// 					at points x for order N and returns the values in the
	// 					array P that is Np long

	double aold = 0.0, anew=0.0;

	double **PL = xcalloc(N+1, sizeof(double*));
	for (int i =0; i < N+1; i++)
	{
		PL[i] =  xcalloc(Np, sizeof(double));
	}

	// Initial values P_0(x) and P_1(x)
	if (N==0)
	{
		double constVal = 1.0/sqrt(2);
		for (int j = 0; j < Np; j++)
		{
			P[j] = constVal;
		}

		free(PL[0]);
		free(PL);
		return;
	}
	else
	{
		for (int j = 0; j < Np; j++)
		{
			PL[0][j] = 1.0/sqrt(2);
		}
	}

	double *prow = xcalloc(Np, sizeof(double));
	double coeff = sqrt(1.5);
	for (int i = 0; i < Np; i++)
	{
		prow[i] = coeff*x[i];
	}

	if(N==1)
	{
		for (int i = 0; i < Np; i++)
		{
			P[i] = prow[i];
		}
		free(prow);
		free(PL[0]);
		free(PL[1]);
		free(PL);
		return;
	}
	else
	{
		for (int j =0; j < Np; j++)
		{
			PL[1][j] = coeff*x[j];
		}

	}

	// Repeat value in recurrence
	aold = sqrt(1.0/3);

	// Forward recurrence
	for (int i = 2; i <= N; i++)
	{
		anew = sqrt(i*i/((2.0*i+1)*(2*i-1)));
		for (int j = 0; j < Np; j++)
		{
			PL[i][j] = (x[j]*PL[i-1][j]-aold*PL[i-2][j])/anew;
		}
		aold = anew;
	}

	for (int j = 0; j < Np; j++)
	{
		P[j] = PL[N][j];
	}

	for (int i =0; i < N; i++)
	{
		free(PL[i]);
	}
	free(PL);

	return;
}

void Vandermonde1D(double *x, int Np, int N, gsl_matrix** V1D)
{
	// Purpose: Evaluates the Vandermonde matrix V_{ij} = phi_i(r_j);

	*V1D = gsl_matrix_alloc(Np, N+1);
	double V1DT[Np];
	for(int i = 0; i < N+1; i++)
	{
			LegendrePoly(x, Np, i, V1DT);
			for (int j = 0; j < Np; j++)
			{
				gsl_matrix_set(*V1D, i, j, V1DT[j]);
				//gsl_matrix_set(*V1D, j, i, V1DT[j]);
			}
	}
}

void GradLegendrePoly(double *x, int Np, int N, double* DP)
{
	// Purpose: Evaluates the derivatives of normalized Legendre polynomial
	// 					at points x for order N and returns the values in the
	// 					array P that is Np long
	//

	double aold, anew;

	double **DPL = xcalloc(N+1, sizeof(double*));
	for (int i =0; i < N+1; i++)
	{
		DPL[i] =  xcalloc(Np, sizeof(double));
	}

	// Initial values P_0(x) and P_1(x)
	if (N==0)
	{
		for (int j = 0; j < Np; j++)
		{
			DP[0] = 0;
		}
		free(DPL[0]);
		free(DPL);
		return;
	}
	else
	{
		for (int j = 0; j < Np; j++)
		{
			DPL[0][j] = 0;
		}
	}

	double *dprow = xcalloc(Np, sizeof(double));
	double deriv = sqrt(1.5);
	for (int i = 0; i < Np; i++)
	{
		dprow[i] = deriv;
	}

	if(N==1)
	{
		for (int i = 0; i < Np; i++)
		{
			DP[i] = deriv;
		}
		free(dprow);
		free(DPL[0]);
		free(DPL[1]);
		free(DPL);
		return;
	}
	else
	{
		for (int j =0; j < Np; j++)
		{
			DPL[1][j] = deriv;
		}

	}

	// Repeat value in recurrence
	aold = sqrt(1.0/3);

	// Forward recurrence
	for (int i = 2; i <= N; i++)
	{
		anew = sqrt(i*i/((2.0*i+1)*(2*i-1)));
		double P[Np];
		LegendrePoly(x, Np, i-1, P);
		for (int j = 0; j < Np; j++)
		{
			DPL[i][j] = 1.0/(anew)*(P[j] + x[j]*DPL[i-1][j] - aold*DPL[i-2][j]);
		}
		aold = anew;
	}

	for (int j = 0; j < Np; j++)
	{
		DP[j] = DPL[N][j];
	}

	for (int i =0; i < N; i++)
	{
		free(DPL[i]);
	}
	free(DPL);

	return;

}

void GradVandermonde1D(double *r, int Np, int N, gsl_matrix** GV)
{
	// Purpose: Evaluates the gradient of modal basis (i) at (r)
	// 					at order N and return them in GV whose size is
	// 					ixNp
	//
	*GV = gsl_matrix_alloc(Np, N+1);
	double GVT[Np];
	for (int i =0; i < N+1; i++)
	{
		GradLegendrePoly(r, Np, i, GVT);
		for (int j = 0; j < Np; j++)
			gsl_matrix_set(*GV,i,j, GVT[j]);
			//gsl_matrix_set(*GV,j,i,GVT[j]);
	}

}

void VolIntMat1D(double *x, int N, gsl_matrix* V, gsl_matrix *Vr, gsl_matrix** VolMat)
{

	gsl_permutation *p = gsl_permutation_alloc(N+1);
	gsl_matrix *Vcopy = gsl_matrix_alloc(N+1, N+1);
	gsl_matrix_memcpy(Vcopy, V);

	gsl_matrix *Vinv = gsl_matrix_alloc(N+1, N+1);

	int s;
	gsl_linalg_LU_decomp(Vcopy, p, &s);
	gsl_linalg_LU_invert(Vcopy, p, Vinv);

	gsl_matrix *tmp1 = gsl_matrix_alloc(N+1, N+1);
	gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, Vinv, Vinv, 0.0, tmp1);

	gsl_matrix *tmp2 = gsl_matrix_alloc(N+1, N+1);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Vr, tmp1, 0.0, tmp2);

	*VolMat = gsl_matrix_alloc(N+1, N+1);
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, V, tmp2, 0.0, *
VolMat);
	//gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Vr, Vinv, 0.0, *Dr);

	gsl_permutation_free(p);
	gsl_matrix_free(Vcopy);
	gsl_matrix_free(Vinv);
	gsl_matrix_free(tmp1);
	gsl_matrix_free(tmp2);

}

void Lift1D(int Np, gsl_matrix *V, gsl_matrix **LIFT)
{
	*LIFT = gsl_matrix_alloc(Np,2);

	gsl_matrix *EMAT = gsl_matrix_calloc(Np, 2);
	gsl_matrix_set(EMAT, 0, 0, 1.0);
	gsl_matrix_set(EMAT, Np-1, 1, -1.0);

	gsl_matrix *tmp = gsl_matrix_alloc(Np,2);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, V, EMAT, 0.0, tmp);
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, V, tmp, 0.0, *LIFT);

	gsl_matrix_free(EMAT);
	gsl_matrix_free(tmp);

}

void calculateMassMatrix(int Np, gsl_matrix* V, gsl_matrix** MassMatrix)
{
	*MassMatrix = gsl_matrix_alloc(Np,Np);
	gsl_matrix *MassMatrixInv = gsl_matrix_alloc(Np,Np);
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, V, V, 0.0, MassMatrixInv);

	int s = 0;
	gsl_permutation *p = gsl_permutation_alloc(Np);
	gsl_linalg_LU_decomp(MassMatrixInv, p, &s);
	gsl_linalg_LU_invert(MassMatrixInv, p, *MassMatrix);

	gsl_matrix_free(MassMatrixInv);
	gsl_permutation_free(p);

}

void calculateLIFTVolMat(int N, int Np, gsl_matrix **LIFT, gsl_matrix **VolMat,gsl_matrix **MassMatrix)
{
	double r[Np];
	GetLGLPoints(N,r);

	gsl_matrix *V;
	gsl_matrix *GV;

	Vandermonde1D(r,Np,N,&V);
	GradVandermonde1D(r, Np, N, &GV);
	VolIntMat1D(r, N, V, GV, VolMat);
	Lift1D(Np, V, LIFT);
	calculateMassMatrix(Np, V,MassMatrix);

	gsl_matrix_free(V);
	gsl_matrix_free(GV);

}

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions computing coeffcients of Not-A-Knot and Quadratic splines in matrix form.
 *
 */

#ifndef _SPLINECOEFFS_H
#define _SPLINECOEFFS_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

void BuildNotAKnotSpline(
  gsl_matrix* splinecoeffs,   /* Output: matrix containing all the spline coeffs (already allocated) */
  gsl_vector* vectx,          /* Input: vector x*/
  gsl_vector* vecty,          /* Input: vector y */
  int n);                      /* Size of x, y, and of output matrix */

void BuildQuadSpline(
  gsl_matrix* splinecoeffs,   /* Output: matrix containing all the spline coeffs (already allocated) */
  gsl_vector* vectx,          /* Input: vector x*/
  gsl_vector* vecty,          /* Input: vector y */
  int n);                      /* Size of x, y, and of output matrix */

void BuildSplineCoeffs(
  CAmpPhaseSpline** splines,                  /*  */
  CAmpPhaseFrequencySeries* freqseries);      /*  */

void BuildListmodesCAmpPhaseSpline(
  ListmodesCAmpPhaseSpline** listspline,              /* Output: list of modes of splines in matrix form */
  ListmodesCAmpPhaseFrequencySeries* listh);          /* Input: list of modes in amplitude/phase form */

/* Functions for spline evaluation */

/* Note: for the spines in matrix form, the first column contains the x values, so the coeffs start at 1 */
double EvalCubic(
  gsl_vector* coeffs,  /**/
  double eps,          /**/
  double eps2,         /**/
  double eps3);         /**/

double EvalQuad(
  gsl_vector* coeffs,  /**/
  double eps,          /**/
  double eps2);         /**/

void EvalCAmpPhaseSpline(
  CAmpPhaseSpline* splines,                     //input
  CAmpPhaseFrequencySeries* freqseries);  //in/out defines CAmpPhase from defined freqs

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _SPLINECOEFFS_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#pragma once

#include <Windows.h>

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstddef>
#include <cassert>
#include <cstdarg>
#include <climits>
#include <cinttypes>
#include <cmath>
#include <ctime>

#include <algorithm>
#include <any>
#include <array>
#include <chrono>
#include <complex>
#include <deque>
#include <exception>
#include <filesystem>
#include <functional>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <list>
#include <locale>
#include <memory>
#include <mutex>
#include <new>
#include <numeric>
#include <optional>
#include <random>
#include <ratio>
#include <sstream>
#include <stack>
#include <string>
#include <string_view>
#include <system_error>
#include <thread>
#include <tuple>
#include <type_traits>
#include <typeinfo>
#include <unordered_set>
#include <utility>
#include <valarray>
#include <vector>

#include <gsl/gsl>

max_stars_repo_name: lopezfjose/Digger

max_issues_repo_name: lopezfjose/Digger

max_forks_repo_name: lopezfjose/Digger