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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>

#include "test_sf.h"

#define WKB_TOL (1.0e+04 * TEST_SQRT_TOL0)

/*
 * Test the identities:
 *
 * Sum_{k=0}^n (n choose k) (2 y)^{n - k} H_k(x) = H_n(x + y)
 *
 * and
 *
 * Sum_{k=0}^n (n choose k) y^{n-k} He_k(x) = He_n(x + y)
 *
 * see: http://mathworld.wolfram.com/HermitePolynomial.html (Eq. 55)
 */
void
test_hermite_id1(const int n, const double x, const double y)
{
  double *a = malloc((n + 1) * sizeof(double));
  double *b = malloc((n + 1) * sizeof(double));
  double lhs, rhs;
  int k;

  a[0] = gsl_pow_int(2.0 * y, n);
  b[0] = gsl_pow_int(y, n);
  for (k = 1; k <= n; ++k)
    {
      double fac = (n - k + 1.0) / (k * y);
      a[k] = 0.5 * fac * a[k - 1];
      b[k] = fac * b[k - 1];
    }

  lhs = gsl_sf_hermite_series(n, x, a);
  rhs = gsl_sf_hermite(n, x + y);
  gsl_test_rel(lhs, rhs, TEST_TOL4, "identity1 phys n=%d x=%g y=%g", n, x, y);

  lhs = gsl_sf_hermite_prob_series(n, x, b);
  rhs = gsl_sf_hermite_prob(n, x + y);
  gsl_test_rel(lhs, rhs, TEST_TOL3, "identity1 prob n=%d x=%g y=%g", n, x, y);

  free(a);
  free(b);
}

int
test_hermite(void)
{
  gsl_sf_result r;

  int s = 0;
  int m, n, sa;
  double res[256];
  double x;
  const double aizero1 = -2.3381074104597670384891972524467; /* first zero of the Airy function Ai */

  /* test some known identities */
  test_hermite_id1(10, 0.75, 0.33);
  test_hermite_id1(8, -0.75, 1.20);
  test_hermite_id1(7, 2.88, -3.2);

  TEST_SF(s, gsl_sf_hermite_prob_e, (0, 0.75, &r),  1.,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_e, (1, 0.75, &r),  0.75,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_prob_e, (25, 0., &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_e, (28, 0., &r),  2.13458046676875e14,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_e, (30, 0., &r), -6.190283353629375e15,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_e, (25, 0.75, &r),  -1.08128685847680748265939328423e12,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_e, (28, 0.75, &r),  -1.60620252094658918105511125135e14,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF_RETURN(s, gsl_sf_hermite_prob_e, (2800, 0.75, &r), GSL_EOVRFLW);
  TEST_SF_RETURN(s, gsl_sf_hermite_prob_e, (2800, 1.75, &r), GSL_EOVRFLW);

  TEST_SF(s, gsl_sf_hermite_prob_deriv_e, (225, 128, 0.75, &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_prob_deriv_e, (5, 128, 0.75, &r),  -3.0288278964712702882066404e112,  TEST_TOL1, GSL_SUCCESS);

  x = 0.75;
  sa = 0;
  gsl_sf_hermite_prob_array(0, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array(0, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array(1, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   x,   TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array(1, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   x, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   823.810509681701660156250, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  1.03749254986255755872498e78, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_prob_array(100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_deriv(0, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   823.810509681701660156250, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  1.03749254986255755872498e78, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_prob_array_deriv(0, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_deriv(1000, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array_deriv(1000, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_deriv(99, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[99], +0.0,    9.332621544394415268169923886e155, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   6.999466158295811451127442914e157, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array_deriv(99, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_deriv(100, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   9.332621544394415268169923886e157, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array_deriv(100, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_array_deriv(23, 100, 0.75, res);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[23], +0.0,  2.5852016738884976640000000000e22, TEST_TOL0);
  TEST_SF_VAL(sa, res[24], +0.0,  4.6533630129992957952000000000e23, TEST_TOL0);
  TEST_SF_VAL(sa, res[37], +0.0,  2.3592417210568968566591219172e37, TEST_TOL0);
  TEST_SF_VAL(sa, res[60], +0.0, -9.2573208827175536243052086845e59, TEST_TOL0);
  TEST_SF_VAL(sa, res[61], +0.0,  5.6259625429686219137261735305e59, TEST_TOL1);
  TEST_SF_VAL(sa, res[100], +0.0, 2.6503570965896336273549197e100, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_array_deriv(23, 37, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_prob_deriv_array(100, 50, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   -3.88338863813139372375411561e31, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   -4.04757862431646677625108652e32, TEST_TOL1);
  TEST_SF_VAL(sa, res[10],  +0.0,    7.9614368698398116765703194e38, TEST_TOL0);
  TEST_SF_VAL(sa, res[30],  +0.0,   -9.1416928183915197188795338e54, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,    0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_prob_deriv_array(100, 50, 0.75)");
  s += sa;

  n = 128;
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_prob_series_e, (n, x, res, &r),  -4.0451066556993485405907339548e68,  TEST_TOL0, GSL_SUCCESS);

  /* phys */

  x = 0.75;

  TEST_SF(s, gsl_sf_hermite_e, (0, 0.75, &r),  1.0,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_e, (1, 0.75, &r),  1.5,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_e, (25, 0., &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_e, (28, 0., &r),  3.497296636753920000e18,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_e, (30, 0., &r), -2.028432049317273600e20,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_e, (25, 0.75, &r),  -9.7029819451106077507781088352e15,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_e, (28, 0.75, &r),   3.7538457078067672096408339776e18,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF_RETURN(s, gsl_sf_hermite_e, (2800, 0.75, &r), GSL_EOVRFLW);
  TEST_SF_RETURN(s, gsl_sf_hermite_e, (2800, 10.1, &r), GSL_EOVRFLW);

  TEST_SF(s, gsl_sf_hermite_deriv_e, (225, 128, 0.75, &r),  0.,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_deriv_e, (5, 128, 0.75, &r),  2.89461215568095657569833e132,  TEST_TOL0, GSL_SUCCESS);

  x = 0.75;
  sa = 0;
  gsl_sf_hermite_array(0, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array(0, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array(1, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   2.*x,   TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array(1, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   2.0*x, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   38740.4384765625, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  -1.4611185395125104593177790757e93, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array(100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array_deriv(0, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   1.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   38740.4384765625, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,  -1.4611185395125104593177790757e93, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array_deriv(0, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array_deriv(1000, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array_deriv(1000, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array_deriv(99, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[99], +0.0,    5.915251651227242890408570128e185, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   8.872877476840864335612855192e187, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array_deriv(99, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array_deriv(100, 100, 0.75, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   1.183050330245448578081714026e188, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array_deriv(100, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_array_deriv(23, 100, 0.75, res);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.0, TEST_TOL0);
  TEST_SF_VAL(sa, res[23],  +0.0,   2.168624344319444261221171200e29, TEST_TOL0);
  TEST_SF_VAL(sa, res[24],  +0.0,   7.807047639549999340396216320e30, TEST_TOL0);
  TEST_SF_VAL(sa, res[37],  +0.0,   1.930387357696033719818118732e46, TEST_TOL0);
  TEST_SF_VAL(sa, res[60],  +0.0,   6.775378005383186748501182409e71, TEST_TOL0);
  TEST_SF_VAL(sa, res[61],  +0.0,  -4.451215867508936256056845902e73, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   2.957966000491202678467161e118, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_array_deriv(23, 100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_deriv_array(100, 50, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,  -8.26632218305863100726861832e38, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   2.40954750392844799126151557e40, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   1.52281030265187793605875604e49, TEST_TOL0);
  TEST_SF_VAL(sa, res[30],  +0.0,   9.52199419132990437892101915e65, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   0.0, TEST_TOL0);
  gsl_test(sa, "gsl_sf_hermite_deriv_array(100, 50, 0.75)");
  s += sa;

  n = 128;
  /* arbitrary weights */
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_series_e, (n, x, res, &r),  1.07772223811696567390619566842e88,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (0,  1.3, &r),  0.322651504564963746879400858624,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (0,  1.3, &r),  0.322651504564963746879400858624,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (1,  1.3, &r),  0.593187573778613235895531272243,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (1,  1.3, &r),  0.593187573778613235895531272243,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (1, -1.3, &r), -0.593187573778613235895531272243,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (1, -1.3, &r), -0.593187573778613235895531272243,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (27, 0, &r),  0.0,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (27, 0, &r),  0.0,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (28, 0, &r),  0.290371943657199641200016132937,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (28, 0, &r),  0.290371943657199641200016132937,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (28, 0.75, &r),  0.23526280808621240649319140441,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (28, 0.75, &r),  0.23526280808621240649319140441,  TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (200, 0.75, &r), -0.13725356483699291817038427801, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (200, 0.75, &r), -0.13725356483699291817038427801, TEST_TOL3, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (100028, 0.75, &r),  -0.02903467369856961147236598086,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (100028, 0.75, &r),  -0.02903467369856961147236598086,  TEST_TOL5, GSL_SUCCESS);

  n = 10025;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  -0.05301278920004176459797680403,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  -0.05301278920004176459797680403,  TEST_TOL4, GSL_SUCCESS);

  n = 10028;
  x = ((sqrt(2*n+1.)+aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.06992968509693993526829596970,  TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  0.06992968509693993526829596970,  TEST_TOL3, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.08049000991742150521366671021,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  0.08049000991742150521366671021,  TEST_TOL4, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-(aizero1/pow(8.*n,1/6.))/2.5);
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  0.08048800667512084252723933250,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  0.08048800667512084252723933250,  TEST_TOL4, GSL_SUCCESS);

  n = 10025;
  x = (sqrt(2*n+1.)-2.5*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  7.97206830806663013555068100e-6,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  7.97206830806663013555068100e-6,  1.0e-9, GSL_SUCCESS);

  n = 10028;
  x = (sqrt(2*n+1.)-2.5*(aizero1/pow(8.*n,1/6.)));
  TEST_SF(s, gsl_sf_hermite_func_e, (n, x, &r),  7.97188517397786729928465829e-6,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (n, x, &r),  7.97188517397786729928465829e-6,  1.0e-8, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_e, (10000, 60.0, &r), 0.03162606955427450540143292572, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_fast_e, (10000, 60.0, &r), 0.03162606955427450540143292572, TEST_TOL4, GSL_SUCCESS);

  x = 0.75;

  sa = 0;
  gsl_sf_hermite_func_array(100, x, res);
  TEST_SF_VAL(sa, res[0],   +0.0,   0.566979307027693616978839335983, TEST_TOL0);
  TEST_SF_VAL(sa, res[1],   +0.0,   0.601372369187597546203014795470, TEST_TOL0);
  TEST_SF_VAL(sa, res[10],  +0.0,   0.360329854170806945032958735574, TEST_TOL0);
  TEST_SF_VAL(sa, res[100], +0.0,   -0.07616422890563462489003733382, TEST_TOL1);
  gsl_test(sa, "gsl_sf_hermite_func_array(100, 0.75)");
  s += sa;

  sa = 0;
  gsl_sf_hermite_func_array(250, 20.0, res);
  TEST_SF_VAL(sa, res[199], +0.0,   0.254621970261340422399150964119, TEST_TOL2);
  TEST_SF_VAL(sa, res[200], +0.0,   0.288364948026205642290411878357, TEST_TOL2);
  TEST_SF_VAL(sa, res[210], +0.0,   0.266625427572822255774117774166, TEST_TOL2);
  TEST_SF_VAL(sa, res[220], +0.0,  -0.296028413764592652216845612832, TEST_TOL2);
  TEST_SF_VAL(sa, res[230], +0.0,   0.229657495510699141971182819560, TEST_TOL2);
  TEST_SF_VAL(sa, res[240], +0.0,  -0.024027870622288819185159269632, TEST_TOL2);
  TEST_SF_VAL(sa, res[250], +0.0,  -0.236101289420398152417654177998, TEST_TOL2);
  gsl_test(sa, "gsl_sf_hermite_func_array(250, 20.0)");
  s += sa;

  n = 128;
  /* arbitrary weights */
  res[0] = 1.;
  for(m=1; m<=n; m++){
    res[m] = res[m-1]/2.;
  }
  TEST_SF(s, gsl_sf_hermite_func_series_e, (n, x, res, &r),  0.81717103529960997134154552556,  TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_der_e, (0, 28, 0.75, &r),  0.235262808086212406493191404,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (1, 28, 0.75, &r),  1.289485094958329643927802330,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (2, 28, 0.75, &r), -13.27764473136561269145948989,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (3, 28, 0.75, &r), -72.42242083458141066943555691,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (4, 28, 0.75, &r),  753.6960554274941800190147503,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (5, 28, 0.75, &r),  4035.32788513029308540826835,   TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_hermite_func_der_e, (0, 380, 0.75, &r), -0.0400554661321992411631174,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (1, 380, 0.75, &r), -4.0417244263030600591206553,    TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (2, 380, 0.75, &r),  30.4596785269042604519780923,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hermite_func_der_e, (3, 380, 0.75, &r),  3073.4187352276349348458186556, TEST_TOL2, GSL_SUCCESS);

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 17 */
    double He17z[8] = { 0.751842600703896170737870774614, 1.50988330779674075905491513417, 2.28101944025298889535537879396,
                        3.07379717532819355851658337833, 3.90006571719800990903311840097, 4.77853158962998382710540812497,
                        5.74446007865940618125547815768,6.88912243989533223256205432938 };

    n = 17;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He17z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 18 */
    double He18z[9] = { 0.365245755507697595916901619097, 1.09839551809150122773848360538, 1.83977992150864548966395498992,
                        2.59583368891124032910545091458, 3.37473653577809099529779309480, 4.18802023162940370448450911428,
                        5.05407268544273984538327527397, 6.00774591135959752029303858752, 7.13946484914647887560975631213 };

    n = 18;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He18z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 23 */
    double He23z[11] = { 0.648471153534495816722576841197, 1.29987646830397886997876116860, 1.95732755293342410739100839243,
                         2.62432363405918177067330340783, 3.30504002175296456723204112903, 4.00477532173330406712238633738,
                         4.73072419745147329426707133987, 5.49347398647179412855289367747, 6.31034985444839982842886078177,
                         7.21465943505186138595859194492, 8.29338602741735258945596770157 };

    n = 23;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He23z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the probabilists' Hermite polynomial of order 24 */
    double He24z[12] = { 0.317370096629452319318170455994, 0.953421922932109084904629632351, 1.59348042981642010695074168129,
                         2.24046785169175236246653790858, 2.89772864322331368932008199475, 3.56930676407356024709649151613,
                         4.26038360501990548884317727406, 4.97804137463912033462166468006, 5.73274717525120114834341822330,
                         6.54167500509863444148277523364, 7.43789066602166310850331715501, 8.50780351919525720508386233432 };

    n = 24;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_prob_zero_e, (n, m, &r),  He24z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 17 */
    double H17z[8] = { 0.531633001342654731349086553718, 1.06764872574345055363045773799, 1.61292431422123133311288254454,
                       2.17350282666662081927537907149, 2.75776291570388873092640349574, 3.37893209114149408338327069289,
                       4.06194667587547430689245559698, 4.87134519367440308834927655662 };

    n = 17;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_zero_e, (n, m, &r),  H17z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 18 */
    double H18z[9] = { 0.258267750519096759258116098711, 0.776682919267411661316659462284, 1.30092085838961736566626555439,
                       1.83553160426162889225383944409, 2.38629908916668600026459301424, 2.96137750553160684477863254906,
                       3.57376906848626607950067599377, 4.24811787356812646302342016090, 5.04836400887446676837203757885 };

    n = 18;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_zero_e, (n, m, &r),  H18z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 23 */
    double H23z[11] = { 0.458538350068104797757887329284, 0.919151465442563765431719239593, 1.38403958568249523732634717118,
                        1.85567703767137106251504753718, 2.33701621147445578644623502174, 2.83180378712615690144806140734,
                        3.34512715994122457247439814585, 3.88447270810610186607248760288, 4.46209117374000667673186157071,
                        5.10153461047667712968749766165, 5.86430949898457256538748413474 };

    n = 23;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_zero_e, (n, m, &r),  H23z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  {
    /* positive zeros of the physicists' Hermite polynomial of order 24 */
    double H24z[12] = { 0.224414547472515585151136715527, 0.674171107037212236000245923730, 1.12676081761124507213306126773,
                        1.58425001096169414850563336202, 2.04900357366169891178708399532, 2.52388101701142697419907602333,
                        3.01254613756556482565453858421, 3.52000681303452471128987227609, 4.05366440244814950394766297923,
                        4.62566275642378726504864923776, 5.25938292766804436743072304398, 6.01592556142573971734857350899 };

    n = 24;
    for (m=1; m<=n/2; m++) {
      TEST_SF(s, gsl_sf_hermite_zero_e, (n, m, &r),  H24z[m-1],  TEST_TOL0, GSL_SUCCESS);
    }
  }

  n = 2121;
  x = -1.*n;
  res[0] = (double) n;
  sa = 0;
  for (m=1; m<=n/2; m++) {
    gsl_sf_hermite_prob_zero_e(n, m, &r);
    if (x>=r.val) {
      sa += TEST_SF_INCONS;
      printf("sanity check failed! (gsl_sf_hermite_prob_zero)\n");
    }
    res[0] = GSL_MIN(res[0],fabs(x-r.val));
    x = r.val;
  }
  gsl_test(sa, "gsl_sf_hermite_prob_zero(n, m, r)");

  n = 2121;
  x = -1.*n;
  res[0] = (double) n;
  sa = 0;
  for (m=1; m<=n/2; m++) {
    gsl_sf_hermite_zero_e(n, m, &r);
    if (x>=r.val) {
      sa += TEST_SF_INCONS;
      printf("sanity check failed! (gsl_sf_hermite_zero)\n");
    }
    res[0] = GSL_MIN(res[0],fabs(x-r.val));
    x = r.val;
  }
  gsl_test(sa, "gsl_sf_hermite_zero(n, m, r)");

  return s;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
 * mutation.h
 *
 *  Created on: 17.5.2017
 *      Author: heine
 */
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#ifndef MUTATION_H_
#define MUTATION_H_

//void set_mutationseed( int id );
void
mutation(
	 long N,
	 double *x,
	 int state_dim,
	 gsl_rng * rgen );

#endif /* MUTATION_H_ */

ext: h

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
/**
 *
 * @generated s Tue Jan  7 11:45:26 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <plasma.h>
#include <core_blas.h>
#include "auxiliary.h"

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

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

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

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

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

    normQ = LAPACKE_slansy_work(LAPACK_COL_MAJOR, 'i', 'u', minMN, Id, minMN, work);

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

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

    free(work); free(Id);

    return info_ortho;
}

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

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

    eps = LAPACKE_slamch_work('e');

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

    alpha=1.0;
    beta=0.0;

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

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

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

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

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, Residual, M, work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, A2, LDA, work);

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

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

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

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */

int s_check_LLTfactorization(int N, float *A1, float *A2, int LDA, int uplo)
{
    float Anorm, Rnorm;
    float alpha;
    int info_factorization;
    int i,j;
    float eps;

    eps = LAPACKE_slamch_work('e');

    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *L1       = (float *)malloc(N*N*sizeof(float));
    float *L2       = (float *)malloc(N*N*sizeof(float));
    float *work              = (float *)malloc(N*sizeof(float));

    memset((void*)L1, 0, N*N*sizeof(float));
    memset((void*)L2, 0, N*N*sizeof(float));

    alpha= 1.0;

    LAPACKE_slacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_strmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_strmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', N, N, Residual, N, work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

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

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*--------------------------------------------------------------
 * Check the gemm
 */
float s_check_gemm(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                   float alpha, float *A, int LDA,
                   float *B, int LDB,
                   float beta, float *Cplasma,
                   float *Cref, int LDC,
                   float *Cinitnorm, float *Cplasmanorm, float *Clapacknorm )
{
    float beta_const = -1.0;
    float Rnorm;
    float *work = (float *)malloc(max(K,max(M, N))* sizeof(float));

    *Cinitnorm   = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref,    LDC, work);
    *Cplasmanorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, Cplasma, LDC, work);

    cblas_sgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                (alpha), A, LDA, B, LDB, (beta), Cref, LDC);

    *Clapacknorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    cblas_saxpy(LDC * N, (beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the trsm
 */
float s_check_trsm(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                   int M, int NRHS, float alpha,
                   float *A, int LDA,
                   float *Bplasma, float *Bref, int LDB,
                   float *Binitnorm, float *Bplasmanorm, float *Blapacknorm )
{
    float beta_const = -1.0;
    float Rnorm;
    float *work = (float *)malloc(max(M, NRHS)* sizeof(float));
    /*float eps = LAPACKE_slamch_work('e');*/

    *Binitnorm   = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, Bref,    LDB, work);
    *Bplasmanorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bplasma, LDB, work);

    cblas_strsm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo,
                (CBLAS_TRANSPOSE)trans, (CBLAS_DIAG)diag, M, NRHS,
                (alpha), A, LDA, Bref, LDB);

    *Blapacknorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);

    cblas_saxpy(LDB * NRHS, (beta_const), Bplasma, 1, Bref, 1);

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);
    Rnorm = Rnorm / *Blapacknorm;
    /* max(M,NRHS) * eps);*/

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the solution
 */

float s_check_solution(int M, int N, int NRHS, float *A, int LDA,
                      float *B,  float *X, int LDB,
                      float *anorm, float *bnorm, float *xnorm )
{
/*     int info_solution; */
    float Rnorm = -1.00;
    float zone  =  1.0;
    float mzone = -1.0;
    float *work = (float *)malloc(max(M, N)* sizeof(float));

    *anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, N,    A, LDA, work);
    *xnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work);
    *bnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, (zone), A, LDA, X, LDB, (mzone), B, LDB);

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    free(work);

    return Rnorm;
}

ext: c


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include "Utils\buffer.h"
#include "MageSettings.h"           // for CameraIdentity
#include "ORBDescriptor.h"

#include "Memory\allocators.h"
#include "arcana\utils\serialization\serializable.h"

#include <opencv2\core\core.hpp>
#include <stdint.h>
#include <atomic>

#include <gsl\span>

namespace mage
{
    using ImageAllocator = memory::std_allocator<void*, block_splitting_allocation_strategy>;

    template<typename AllocT>
    struct ImageData : mira::serializable<ImageData<AllocT>>
    {
        ImageData(CameraIdentity cameraIdentity, size_t maxFeatures, float pyramidScale, size_t numLevels, float imageBorder, const AllocT& allocator)
            :   m_cameraIdentity{ cameraIdentity },
                m_featureCount{ 0 },
                m_pyramidScale { pyramidScale },
                m_numLevels { numLevels },
                m_imageBorder { imageBorder },
                m_published{ false },
                m_maxFeatures{ maxFeatures },
                m_keypointBuffer{ maxFeatures, allocator },
                m_descriptorBuffer{ maxFeatures, allocator }
        {}

        template<typename StreamT, typename = mira::is_stream_t<StreamT>>
        ImageData(CameraIdentity cameraIdentity, size_t maxFeatures, float pyramidScale, size_t numLevels, float imageBorder, const AllocT& allocator, StreamT& stream)
            : ImageData{ cameraIdentity, imageWidth, imageHeight, maxFeatures, pyramidScale, allocator }
        {
            deserialize(stream);
        }

        static constexpr auto members()
        {
            return declare_members(
                &ImageData::m_cameraIdentity,
                &ImageData::m_featureCount,
                &ImageData::m_pyramidScale,
                &ImageData::m_numLevels,
                &ImageData::m_imageBorder,
                &ImageData::m_published,
                &ImageData::m_maxFeatures,
                &ImageData::m_keypointBuffer,
                &ImageData::m_descriptorBuffer
            );
        }

        /*
            Inserts the keypoints into the image data, and returns how many
            items it copied. Because we only support a fixed size, we
            only copy items if we have space.
        */
        size_t Insert(gsl::span<const cv::KeyPoint> keypoints)
        {
            auto toCopy = std::min<int>(gsl::narrow_cast<int>(keypoints.size()), gsl::narrow_cast<int>(m_maxFeatures - m_featureCount));
            m_featureCount += m_keypointBuffer.copy(keypoints.data(), keypoints.data() + toCopy, m_featureCount);
            return toCopy;
        }

        size_t Insert(const cv::KeyPoint& kp)
        {
            if (m_featureCount == m_maxFeatures)
                return 0;

            m_keypointBuffer[m_featureCount] = kp;
            m_featureCount++;
            return 1;
        }

        size_t Insert(const cv::KeyPoint* beg, const cv::KeyPoint* end)
        {
            auto toCopy = std::min<int>(gsl::narrow_cast<int>(end - beg), gsl::narrow_cast<int>(m_maxFeatures - m_featureCount));
            m_featureCount += m_keypointBuffer.copy(beg, beg + toCopy, m_featureCount);
            return toCopy;
        }

        size_t GetMaxFeatures() const
        {
            return m_maxFeatures;
        }

        size_t GetFeatureCount() const
        {
            return m_featureCount;
        }

        float GetPyramidScale() const
        {
            return m_pyramidScale;
        }

        size_t GetNumLevels() const
        {
            return m_numLevels;
        }

        float GetImageBorder() const
        {
            return m_imageBorder;
        }

        CameraIdentity GetCameraIdentity() const
        {
            return m_cameraIdentity;
        }

        const cv::KeyPoint& GetKeypoint(size_t idx) const
        {
            return const_cast<ImageData*>(this)->GetKeypoint(idx);
        }

        cv::KeyPoint& GetKeypoint(size_t idx)
        {
            assert(0 <= idx && idx < m_featureCount);
            return m_keypointBuffer[idx];
        }

        gsl::span<const cv::KeyPoint> GetKeypoints() const
        {
            return { m_keypointBuffer.begin(), m_keypointBuffer.begin() + m_featureCount };
        }

        gsl::span<cv::KeyPoint> GetKeypoints()
        {
            return { m_keypointBuffer.begin(), m_keypointBuffer.begin() + m_featureCount };
        }

        const ORBDescriptor& GetDescriptor(size_t idx) const
        {
            return const_cast<ImageData*>(this)->GetDescriptor(idx);
        }

        ORBDescriptor& GetDescriptor(size_t idx)
        {
            assert(0 <= idx && idx < m_featureCount);
            return m_descriptorBuffer[idx];
        }

        gsl::span<const ORBDescriptor> GetDescriptors() const
        {
            return const_cast<ImageData*>(this)->GetDescriptors();
        }

        gsl::span<ORBDescriptor> GetDescriptors()
        {
            return { m_descriptorBuffer.begin(), m_descriptorBuffer.begin() + m_featureCount };
        }

        void SetDescriptors(gsl::span<const ORBDescriptor> values)
        {
            assert(m_featureCount == values.size());
            m_descriptorBuffer.copy(values.begin(), values.end());
        }

        void ZeroOutDescriptors()
        {
            m_descriptorBuffer.fill(0);
        }

        /*
            Returns whether or not this image was ever published to the map
        */
        bool IsPublished() const
        {
            return m_published;
        }

        void MarkAsPublished() const
        {
            m_published = true;
        }

        template<typename DerivedT>
        static size_t AllocationSizeInBytes(size_t maxFeatures)
        {
            return sizeof(DerivedT) + maxFeatures * sizeof(cv::KeyPoint) + maxFeatures * sizeof(ORBDescriptor);
        }

    private:
        template<typename T>
        using rebind_allocator = typename std::allocator_traits<AllocT>::template rebind_alloc<T>;

        size_t m_featureCount{};
        mutable std::atomic<bool> m_published{ false };

        const size_t m_maxFeatures{};
        const float m_pyramidScale{};
        const size_t m_numLevels{}; //TODO: change numlevels to be reasonable with image resolution
        const float m_imageBorder{};
        const CameraIdentity m_cameraIdentity{};

        buffer<cv::KeyPoint, rebind_allocator<cv::KeyPoint>> m_keypointBuffer{};
        buffer<ORBDescriptor, rebind_allocator<ORBDescriptor>> m_descriptorBuffer{};
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>

int
main(void)
{
  const size_t N = 500;                        /* length of time series */
  const size_t K = 51;                         /* window size */
  const double alpha[3] = { 0.5, 3.0, 10.0 };  /* alpha values */
  gsl_vector *x = gsl_vector_alloc(N);         /* input vector */
  gsl_vector *y1 = gsl_vector_alloc(N);        /* filtered output vector for alpha1 */
  gsl_vector *y2 = gsl_vector_alloc(N);        /* filtered output vector for alpha2 */
  gsl_vector *y3 = gsl_vector_alloc(N);        /* filtered output vector for alpha3 */
  gsl_vector *k1 = gsl_vector_alloc(K);        /* Gaussian kernel for alpha1 */
  gsl_vector *k2 = gsl_vector_alloc(K);        /* Gaussian kernel for alpha2 */
  gsl_vector *k3 = gsl_vector_alloc(K);        /* Gaussian kernel for alpha3 */
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  gsl_filter_gaussian_workspace *gauss_p = gsl_filter_gaussian_alloc(K);
  size_t i;
  double sum = 0.0;

  /* generate input signal */
  for (i = 0; i < N; ++i)
    {
      double ui = gsl_ran_gaussian(r, 1.0);
      sum += ui;
      gsl_vector_set(x, i, sum);
    }

  /* compute kernels without normalization */
  gsl_filter_gaussian_kernel(alpha[0], 0, 0, k1);
  gsl_filter_gaussian_kernel(alpha[1], 0, 0, k2);
  gsl_filter_gaussian_kernel(alpha[2], 0, 0, k3);

  /* apply filters */
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha[0], 0, x, y1, gauss_p);
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha[1], 0, x, y2, gauss_p);
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha[2], 0, x, y3, gauss_p);

  /* print kernels */
  for (i = 0; i < K; ++i)
    {
      double k1i = gsl_vector_get(k1, i);
      double k2i = gsl_vector_get(k2, i);
      double k3i = gsl_vector_get(k3, i);

      printf("%e %e %e\n", k1i, k2i, k3i);
    }

  printf("\n\n");

  /* print filter results */
  for (i = 0; i < N; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double y1i = gsl_vector_get(y1, i);
      double y2i = gsl_vector_get(y2, i);
      double y3i = gsl_vector_get(y3, i);

      printf("%.12e %.12e %.12e %.12e\n", xi, y1i, y2i, y3i);
    }

  gsl_vector_free(x);
  gsl_vector_free(y1);
  gsl_vector_free(y2);
  gsl_vector_free(y3);
  gsl_vector_free(k1);
  gsl_vector_free(k2);
  gsl_vector_free(k3);
  gsl_rng_free(r);
  gsl_filter_gaussian_free(gauss_p);

  return 0;
}


text:
#ifndef __COYOTE_H
#define __COYOTE_H

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "errorlist.h"
#include "maths.h"

/* Error codes */
#define coyote_base -2200
#define coyote_h        -1 + coyote_base
#define coyote_flat     -2 + coyote_base
#define coyote_range    -3 + coyote_base

/* =============================================== *
 * From other header files which have been removed *
 * Constants which are used also outside of Coyote *
 * (e.g. in cosmo.c) 				   *
 * =============================================== */

/* Coyote I */

/* Number of k values */
#define nsim 1995
const double ksim[nsim];

/* Coyote II */

#define fr_nsim 582
#define fr_rs   11
const double fr_ksim[fr_nsim];

/* v1: Functions from hubble.c and emu.c */
double getH0fromCMB(double omega_m, double omega_b, double w0_de, int physical);
void fill_xstar_wo_z(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de, double xstar[]);
double P_NL_coyote5(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de,
		    double h_100, double a, double k, error **err);
void emu(double *xstar, double *ystar, error **err);

/* v2: Functions from fr_emu.c */
double P_NL_coyote6(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de,
		    double h_100, double a, double k, double **ystar_allz, error **err);
void fr_check_range(const double *xstar, error **err);
void fr_fill_ystar_allz(double *ystar_allz, const double *xstar, error **err);
void fr_emu(const double *xstar, double *ystar, const double *ystar_allz, error **err);
void fill_xstar6_wo_z(double omega_m, double omega_b, double n_spec, double sigma_8, double w0_de, double h_100,
                      double xstar[]);

#endif

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
#ifndef ORIENTATION_CONTROLLER_VARIABLE_GAINS_H_
#define ORIENTATION_CONTROLLER_VARIABLE_GAINS_H_
#include <cmath>
#include <barrett/units.h>
#include <barrett/systems.h>
#include <barrett/products/product_manager.h>
#include <libconfig.h++>
#include <barrett/detail/ca_macro.h>
#include <barrett/cdlbt/calgrav.h>
#include <barrett/systems/abstract/system.h>
#include <barrett/systems/abstract/single_io.h>
#include <barrett/systems/kinematics_base.h>
#include <barrett/math/kinematics.h>
#include <stdio.h>
#include <iostream>
#include <string>
#include <cstdlib>  // For std::atexit()
#include <barrett/os.h>  // For btsleep()
#include <barrett/math.h>  // For barrett::math::saturate()
#include <Eigen/Dense>
#include <Eigen/Core>
#include <barrett/standard_main_function.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#define BARRETT_SMF_VALIDATE_ARGS
#include <barrett/standard_main_function.h>

using namespace barrett;
using namespace systems;
using systems::connect;
BARRETT_UNITS_FIXED_SIZE_TYPEDEFS;

template<size_t DOF>
class OrientationControllerVariableGains : public systems::System  , public KinematicsInput<DOF>
{
	BARRETT_UNITS_TEMPLATE_TYPEDEFS(DOF);

public:
	Input<Eigen::Quaterniond> ReferenceOrnInput;
	Input<Eigen::Quaterniond> FeedbackOrnInput;
	Input<cp_type> KpGains;
	Input<cp_type> KdGains;

public:
	Output<ct_type> CTOutput;

protected:
	typename Output<ct_type>::Value* ctOutputValue;

public:
	OrientationControllerVariableGains(
//    const units::CartesianPosition::type& ProportionalGains,
//    const units::CartesianPosition::type& DerivativeGains,
    const std::string& sysName = "OrientationControllerVariableGains") :
      System(sysName), KinematicsInput<DOF>(this), ReferenceOrnInput(this), FeedbackOrnInput(this),
      KpGains(this) , KdGains(this),
      CTOutput(this, &ctOutputValue)
//      , ProportionalGains(ProportionalGains), DerivativeGains(DerivativeGains)

{
}
virtual ~OrientationControllerVariableGains()
{
this->mandatoryCleanUp();
}

protected:
	cp_type ProportionalGains;
	cp_type DerivativeGains;

	Eigen::AngleAxisd error;
	ct_type ct;

	cp_type tempVect;
	cp_type tempVect2;
	cp_type errorVect;

	virtual void operate() {

		ProportionalGains = KpGains.getValue();
		DerivativeGains	  = KdGains.getValue();

//		error = this->referenceInput.getValue() * this->feedbackInput.getValue().inverse();  // I think it should be this way
		error = this->FeedbackOrnInput.getValue() * this->ReferenceOrnInput.getValue().inverse();  // but CD's math (that works, BTW) does it this way

		double angle = error.angle();
		// TODO(dc): I looked into Eigen's implementation and noticed that angle will always be between 0 and 2*pi. We should test for this so if Eigen changes, we notice.
		if (angle > M_PI) {
			angle -= 2.0*M_PI;
		}

		if (math::abs(angle) > 3.13) {	// a little dead-zone near the discontinuity at +/-180 degrees
			ct.setZero();
		} else {
		    tempVect = ProportionalGains ;				// copy the ProportiocalGains
		    errorVect = error.axis() * angle ;
		    gsl_vector_mul (tempVect.asGslType() , errorVect.asGslType()  ) ; 	// tempVect <- tempVect * errorVect
		    ct = this->ReferenceOrnInput.getValue().inverse() * ( tempVect );

		}

		tempVect2 = DerivativeGains;
		gsl_vector_mul ( tempVect2.asGslType() , this->kinInput.getValue().impl->tool_velocity_angular ) ;
		ct -= tempVect2 ;
//		ct -= DerivativeGains.dot( this->kinInput.getValue().impl->tool_velocity_angular ) ;
//		gsl_blas_daxpy( -kd, this->kinInput.getValue().impl->tool_velocity_angular, ct.asGslType());

		this->ctOutputValue->setData(&ct);
	}

private:
	DISALLOW_COPY_AND_ASSIGN(OrientationControllerVariableGains);

public:
	EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

#endif /*ORIENTATION_CONTROLLER_VARIABLE_GAINS_H_*/


text:
/*Python wrapper module for parameter space sampling;
the method used is the same as in
http://dan.iel.fm/posts/python-c-extensions/

The module is called _design and it defines the methods below (see docstrings)
*/

#include <gsl/gsl_matrix.h>

#include <Python.h>
#include <numpy/arrayobject.h>

#include "lenstoolsPy3.h"
#include "design.h"

#ifndef IS_PY3K
static struct module_state _state;
#endif

//Python module docstrings
static char module_docstring[] = "This module provides a python interface for sampling a N-dimensional parameter space";
static char diagonalCost_docstring[] = "Compute the cost function for a diagonal design";
static char cost_docstring[] = "Compute the cost function for a given";
static char sample_docstring[] = "Sample the N-dimensional parameter space and record cost function changes";

//Method declarations
static PyObject *_design_diagonalCost(PyObject *self,PyObject *args);
static PyObject *_design_cost(PyObject *self,PyObject *args);
static PyObject *_design_sample(PyObject *self,PyObject *args);

//_design method definitions
static PyMethodDef module_methods[] = {

	{"diagonalCost",_design_diagonalCost,METH_VARARGS,diagonalCost_docstring},
	{"cost",_design_cost,METH_VARARGS,cost_docstring},
	{"sample",_design_sample,METH_VARARGS,sample_docstring},
	{NULL,NULL,0,NULL}

} ;

//_design constructor

#ifdef IS_PY3K

static struct PyModuleDef moduledef = {

	PyModuleDef_HEAD_INIT,
	"_design",
	module_docstring,
	sizeof(struct module_state),
	module_methods,
	NULL,
	myextension_trasverse,
	myextension_clear,
	NULL

};

#define INITERROR return NULL

PyMODINIT_FUNC
PyInit__design(void)

#else

#define INITERROR return

void
init_design(void)
#endif

{
#ifdef IS_PY3K
	PyObject *m = PyModule_Create(&moduledef);
#else
	PyObject *m = Py_InitModule3("_design",module_methods,module_docstring);
#endif

	if(m==NULL)
		INITERROR;
	struct module_state *st = GETSTATE(m);

	st->error = PyErr_NewException("_design.Error", NULL, NULL);
    if (st->error == NULL) {
        Py_DECREF(m);
        INITERROR;
    }

	/*Load numpy functionality*/
	import_array();

	/*Return*/
#ifdef IS_PY3K
	return m;
#endif

}

//////////////////////////////////////////////
//////////////////////////////////////////////
//////////////////////////////////////////////

//////////////////////////////////////////////////
/*Function implementations using backend C code*/
/////////////////////////////////////////////////

//diagonalCost() implementation
static PyObject *_design_diagonalCost(PyObject *self,PyObject *args){

	int Npoints;
	double lambda,costValue;

	//Interpret arguments
	if(!PyArg_ParseTuple(args,"id",&Npoints,&lambda)){
		return NULL;
	}

	//Compute cost function for a diagonal design
	costValue = diagonalCost(Npoints,lambda);

	//Throw error if cost is not positive
	if(costValue<=0.0){
		PyErr_SetString(PyExc_RuntimeError,"Cost function is not positive!");
		return NULL;
	}

	//Build return value and return
	PyObject *ret = Py_BuildValue("d",costValue);
	return ret;

}

//cost() implementation
static PyObject *_design_cost(PyObject *self,PyObject *args){

	int i,j;
	double p,lambda,costValue;
	PyObject *data_obj;

	/*Parse the input tuple*/
	if(!PyArg_ParseTuple(args,"Odd",&data_obj,&p,&lambda)){
		return NULL;
	}

	/*Interpret the points in the design as a numpy array*/
	PyObject *data_array = PyArray_FROM_OTF(data_obj,NPY_DOUBLE,NPY_IN_ARRAY);
	if(data_array==NULL){
		return NULL;
	}

	/*Get a pointer to the array data*/
	double *data = (double *)PyArray_DATA(data_array);

	/*squeeze the dimensions of the parameter space and the number of points*/
	int Npoints = (int)PyArray_DIM(data_array,0);
	int Ndim = (int)PyArray_DIM(data_array,1);

	/*Wrap a gsl matrix object around the data_array*/
	gsl_matrix *m = gsl_matrix_alloc(Npoints,Ndim);

	/*Copy the elements into it*/
	for(i=0;i<Npoints;i++){
		for(j=0;j<Ndim;j++){
			gsl_matrix_set(m,i,j,data[i*Ndim+j]);
		}
	}

	/*Call the C backend to compute the cost function*/
	costValue = cost(m,Npoints,Ndim,p,lambda);

	/*Release the gsl resources*/
	gsl_matrix_free(m);

	/*Throw error if cost is not positive*/
	if(costValue<=0.0){
		PyErr_SetString(PyExc_RuntimeError,"Cost function is not positive!");
		Py_DECREF(data_array);
		return NULL;
	}

	/*Build the return value and return it*/
	PyObject *ret = Py_BuildValue("d",costValue);
	Py_DECREF(data_array);

	return ret;

}

//sampler() implementation
static PyObject *_design_sample(PyObject *self,PyObject *args){

	int i,j,maxIterations,seed;
	double p,lambda;
	PyObject *data_obj,*cost_obj;

	/*Parse the input tuple*/
	if(!PyArg_ParseTuple(args,"OddiiO",&data_obj,&p,&lambda,&maxIterations,&seed,&cost_obj)){
		return NULL;
	}

	/*Interpret the data_obj as an array*/
	PyObject *data_array = PyArray_FROM_OTF(data_obj,NPY_DOUBLE,NPY_IN_ARRAY);
	PyObject *cost_array = PyArray_FROM_OTF(cost_obj,NPY_DOUBLE,NPY_IN_ARRAY);

	if(data_array==NULL || cost_array==NULL){

		Py_XDECREF(data_array);
		Py_XDECREF(cost_array);

		return NULL;
	}

	/*Get a pointer to the array data*/
	double *data = (double *)PyArray_DATA(data_array);

	/*Get a pointer to the cost values*/
	double *cost_values = (double *)PyArray_DATA(cost_array);

	/*squeeze the dimensions of the parameter space and the number of points*/
	int Npoints = (int)PyArray_DIM(data_array,0);
	int Ndim = (int)PyArray_DIM(data_array,1);

	/*Wrap a gsl matrix object around the data_array*/
	gsl_matrix *m = gsl_matrix_alloc(Npoints,Ndim);

	/*Copy the elements into it*/
	for(i=0;i<Npoints;i++){
		for(j=0;j<Ndim;j++){
			gsl_matrix_set(m,i,j,data[i*Ndim+j]);
		}
	}

	/*Spread the points in the parameter space looking for the cost function minimum*/
	double deltaPerc = sample(Npoints,Ndim,p,lambda,seed,maxIterations,m,cost_values);

	/*Copy the points positions from the gsl matrix to the data array*/
	for(i=0;i<Npoints;i++){
		for(j=0;j<Ndim;j++){
			data[i*Ndim+j] = gsl_matrix_get(m,i,j);
		}
	}

	/*Release the gsl resources*/
	gsl_matrix_free(m);

	/*Release the other resources*/
	Py_DECREF(data_array);
	Py_DECREF(cost_array);

	/*Build the return value*/
	PyObject *ret = Py_BuildValue("d",deltaPerc);
	return ret;

}

max_stars_repo_name: asabyr/LensTools

max_issues_repo_name: asabyr/LensTools

max_forks_repo_name: asabyr/LensTools


text:
/******************************************************************************
 * Copyright 2018 The Apollo Authors. All Rights Reserved.
 *
 * 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.
 *****************************************************************************/

#pragma once

#include <cblas.h>
#include <cuda_runtime_api.h>

#include <boost/shared_ptr.hpp>
#include <fstream>
#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <vector>

#include "modules/perception/base/blob.h"
#include "modules/perception/base/image.h"

namespace apollo {
namespace perception {
namespace inference {

bool ResizeGPU(const base::Image8U &src,
               std::shared_ptr<apollo::perception::base::Blob<float>> dst,
               int stepwidth, int start_axis);

bool ResizeGPU(const apollo::perception::base::Blob<uint8_t> &src_gpu,
               std::shared_ptr<apollo::perception::base::Blob<float>> dst,
               int stepwidth, int start_axis, int mean_b, int mean_g,
               int mean_r, bool channel_axis, float scale);

bool ResizeGPU(const base::Image8U &src,
               std::shared_ptr<apollo::perception::base::Blob<float>> dst,
               int stepwidth, int start_axis, float mean_b, float mean_g,
               float mean_r, bool channel_axis, float scale);

}  // namespace inference
}  // namespace perception
}  // namespace apollo

max_stars_repo_name: ghdawn/apollo


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

/* Created: [GJ] Tue Apr 23 21:26:53 EDT 1996
 */
#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_errno.h>
/* #include "interpolation.h" */
#include <gsl/gsl_integration.h>
#if 0

static char * info_string = 0;

void set_integ_info(const char * mess)
{
  if(info_string != 0) free(info_string);
  info_string = 0;
  if(mess != 0) {
    info_string = (char *)malloc((strlen(mess)+1) * sizeof(char));
    strcpy(info_string, mess);
  }
}

/* Integration function for use with open_romberg(). */
#define FUNC(x) ((*func)(x))
static double midpnt(double (*func)(double), double a, double b, int n)
{
  double x,tnm,sum,del,ddel;
  static double s;
  static int it;
  int j;

  if (n == 1) {
    it=1;
    return (s=(b-a)*FUNC(0.5*(a+b)));
  }
  else {
    tnm=it;
    del=(b-a)/(3.0*tnm);
    ddel=del+del;
    x=a+0.5*del;
    sum=0.0;
    for (j=1;j<=it;j++) {
      sum += FUNC(x);
      x += ddel;
      sum += FUNC(x);
      x += del;
    }
    it *= 3;
    s=(s+(b-a)*sum/tnm)/3.0;
    return s;
  }
}
#undef FUNC

#define JMAX 10
#define JMAXP JMAX+1
#define K 5
double open_romberg(double(*func)(double), double a, double b, double eps)
{
  int j;
  double ss,dss,h[JMAXP+1],s[JMAXP+1];

  h[1]=1.0;
  for (j=1;j<=JMAX;j++) {
    s[j]=midpnt(func,a,b,j);
    if (j >= K) {
      /* local_polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss); */
      interp_poly(&h[j-K]+1,&s[j-K]+1,K,0.0,&ss,&dss);
      if (fabs(dss) < eps * fabs(ss)) return ss;
    }
    s[j+1]=s[j];
    h[j+1]=h[j]/9.0;
  }

  push_error("open_romberg: too many steps", Error_ConvFail_);
  push_generic_error("open_romberg:", info_string);

  return 0.;
}
#undef JMAX
#undef JMAXP
#undef K

double gauss_legendre_10(double (*func)(double), double a, double b)
{
  int j;
  static double x[] = {0.,
                         0.1488743389,
                         0.4333953941,
                         0.6794095682,
                         0.8650633666,
                         0.9739065285};
  static double w[] = {0.,
                         0.2955242247,
                         0.2692667193,
                         0.2190863625,
                         0.1494513491,
                         0.0666713443};

  double xm = 0.5 * (b + a);
  double xr = 0.5 * (b - a);

  double s = 0.;
  double dx;
  double result;

  for(j=1; j<=5; j++){
    dx = xr * x[j];
    s += w[j] * ((*func)(xm+dx) + (*func)(xm-dx));
  }
  result = s * xr;

  return result;
}

/************************************************************************
 *                                                                      *
 * Trapezoid rule.                                                      *
 *                                                                      *
 * The original trapzd() function from the Numerical Recipes worked     *
 * by side-effects; it tried to remember the value of the integral      *
 * at the current level of refinement. This is REAL BAD, because if     *
 * you try to use it to do a double integral you will get a surprise.   *
 * You cannot tell which "integral" it is remembering. This stems from  *
 * the stupid fact that there was only one, essentially global,         *
 * variable that was doing this memory job.                             *
 *                                                                      *
 * The solution is simple: pass the current refinement to the function  *
 * so that it can be remembered by an external environment, making it   *
 * easy to avoid confusion.                                             *
 *                                                                      *
 * So the new-method code-fragment for doing an integral looks like:    *
 *                                                                      *
 *       double answer;                                                 *
 *       for(j=1; j<=M+1; j++)                                          *
 *         trapezoid_rule(func, a, b, j, &answer);                      *
 *                                                                      *
 ************************************************************************/
void trapezoid_rule(double(*f)(double), double a, double b, int n, double *s)
{
  double x, tnm, sum, del;
  int it, j;

  if(n==1){
    *s = 0.5 * (b-a) * (f(b) + f(a));
  }
  else {
    for(it=1, j=1; j < n-1; j++)  it <<= 1;
    tnm = (double) it;
    del = (b-a) / tnm;
    x = a + 0.5 * del;

    for(sum=0., j=1; j<=it; j++, x+=del) { sum += f(x); }

    *s = 0.5 * (*s + del * sum);
  }
}

/************************************************************************
 *                                                                      *
 * Trapezoid rule.                                                      *
 * This version produces a tracing output.                              *
 *                                                                      *
 ************************************************************************/
#define FUNC(x) ((*func)(x))

void test_trapezoid_rule(double(*func)(double), double a, double b, int n,
                         double *s)
{
  double x, tnm, sum, del;
  int it, j;

  if(n==1){
    printf("t:  a= %g  b= %g   f(a)= %g  f(b)= %g\n",
           a, b, FUNC(a), FUNC(b));
  }

  if(n==1){
    *s = 0.5 * (b-a) * (FUNC(b) + FUNC(a));
    printf("s= %g\n", *s);
  }
  else {
    for(it=1, j=1; j < n-1; j++)  it <<= 1;
    tnm = (double) it;
    del = (b-a) / tnm;
    x = a + 0.5 * del;

    for(sum=0., j=1; j<=it; j++, x+=del) sum += FUNC(x);

    *s = 0.5 * (*s + del * sum);

    printf("sum= %g   tnm= %g  del= %g  s= %g\n", sum, tnm, del, *s);
  }
}
#undef FUNC

/* This is fixed to use the non-side-effecting version of the
 * trapezoidal rule, as implemented in trapezoid_rule() above.
 * See the discussion there for explanation of the original problem.
 */
#define JMAX 20
double gsl_integ_simpson(double (*func)(double), double a, double b, double eps)
{
  int j;
  double s, st, ost, os;

  ost = os = -1.e50;

  for(j=1; j<=JMAX; j++){

    trapezoid_rule(func, a, b, j, &st);
    s = (4.*st - ost) / 3.;

    if(fabs(s-os) < eps * fabs(os))
      return s;

    os = s;
    ost = st;
  }

  GSL_MESSAGE("simpson: too many steps");

  return 0.;
}
#undef JMAX

double gsl_integ_simpson_table(const double * x, const double * y, int n)
{
  int i;
  double result = 0.;
  for(i=0; i<n-1; i++) {
    result += 0.5 * (y[i+1] + y[i]) * (x[i+1] - x[i]);
  }
  return result;
}

/* apparatus for lorentzian variable change to remove pole */
static double (*dummy_f)(double);
static double dummy_x0;
static double dummy_w;
static inline double f_y(double y)
{ return dummy_f(dummy_x0 + dummy_w * tan(dummy_w * y)); }

double gsl_integ_lorenz(double (*f)(double),
                        double x0, double w,
                        double a, double b,
                        double eps)
{
  dummy_f = f;
  dummy_x0 = x0;
  dummy_w = fabs(w);

  if(dummy_w < 10.*sqrt(min_double)) {
    char buff[100];
    sprintf(buff,"lorenz_integ: width w= %g  too small", dummy_w);
    GSL_MESSAGE(buff);
    return 0.;
  }
  else {
    double lower_y = atan((a-x0)/dummy_w) / dummy_w;
    double upper_y = atan((b-x0)/dummy_w) / dummy_w;
    return gsl_integ_simpson(f_y, lower_y, upper_y, eps);
  }
}

#endif /* 0 */

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:
//
// Created by Harold on 2020/9/21.
//

#ifndef M_MATH_M_FILTER_H
#define M_MATH_M_FILTER_H

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

namespace M_MATH {
    class GaussianFilter {
    public:
        enum EndpointType {
            NO_PAD,
            ZERO_PAD,
            VALUE_PAD
        };
    public:
        static bool filter(double const* x,
                           size_t length,
                           size_t window_size,
                           double alpha,
                           EndpointType type,
                           size_t n_length,
                           double* n_x,
                           size_t order);
    private:
        static gsl_filter_end_t EndpointType2GSL(EndpointType);
    };

    gsl_filter_end_t GaussianFilter::EndpointType2GSL(GaussianFilter::EndpointType type) {
        switch (type) {
            case NO_PAD:
                return GSL_FILTER_END_TRUNCATE;
            case ZERO_PAD:
                return GSL_FILTER_END_PADZERO;
            case VALUE_PAD:
                return GSL_FILTER_END_PADVALUE;
        }
        // default is no pad
        // no padding is performed, and the windows are simply truncated as the end points are approached
        return GSL_FILTER_END_TRUNCATE;
    }

    // filter with nth order Gaussian filter, default order = 0
    bool GaussianFilter::filter(const double *x, size_t length, size_t window_size, double alpha,
                                GaussianFilter::EndpointType type, size_t n_length, double *n_x,
                                size_t order = 0) {
        // workspace
        gsl_filter_gaussian_workspace *gauss_p = gsl_filter_gaussian_alloc(window_size);
        if (gauss_p == nullptr)
            return false;
        // xx, yy
        gsl_vector *xx = gsl_vector_alloc(length);
        gsl_vector *yy = gsl_vector_alloc(length);
        for (auto i = 0; i < length; ++i) {
            gsl_vector_set(xx, i, x[i]);
        }

        /* print kernel to debug
        // kernel
        gsl_vector *k = gsl_vector_alloc(window_size);
        // compute kernel without normalization
        gsl_filter_gaussian_kernel(alpha, 0, 0, k);
        // print kernel
        printf("%s", "kernel: ");
        for (auto i = 0; i < window_size; ++i)
            printf("%e ", gsl_vector_get(k, i));
        printf("\n");
        gsl_vector_free(k);
         */

        // apply filter
        if (gsl_filter_gaussian(EndpointType2GSL(type), alpha, order, xx, yy, gauss_p) != 0) {
            gsl_vector_free(xx);
            gsl_vector_free(yy);
            gsl_filter_gaussian_free(gauss_p);
            return false;
        }

        for (auto i = 0; i < n_length; ++i) {
            n_x[i] = gsl_vector_get(yy, i);
        }

        // clean up
        gsl_vector_free(xx);
        gsl_vector_free(yy);
        gsl_filter_gaussian_free(gauss_p);

        return true;
    }
}

#endif //M_MATH_M_FILTER_H


text:
#include <cerrno>
#include <cstdio>
#include <cstdlib>
#include <sys/stat.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <time.h>

#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_math.h>
#include <math.h>       /* pow */
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>

// Common functions to all  .cpps

/* LCDM Hubble in 1/s */
double hubble(double om, double orad, double a){
	double a3 = pow(a,3);
	double a4 = a*a3;
	double ol = 1.-om-orad;
	return h0 * sqrt(om/a3+ orad/a4 + ol);
}

// Reheating time as a function of temperature --- Trh should be given in Kelvin
double timeofrh(double Trh){
	return sqrt(5./pow(M_PI,3)/gstar) * pow(10,lgmp)/pow(Trh*keltgev,2) *kgtgev / gevtds ;
}

/* Critical density in kg/m^3 */
double rhoc(double a){
	return 3.*pow(hubble(om, orad, a),2.)/(8.*M_PI*gnewton);
}

/* Omega_x for LCDM */
double omegalcdm(double ox, double a){
	return  ox/pow(a,3.) * pow(h0/hubble(om,orad,a),2.);
}

/* Log10 of Mass spectrum in units 1/kg */
// hi mass - Eq.2.5 of draft
//  Psi = 10^psihilg   as a function of Log10[M] = lgm and omega_pbh
double psihilg(double lgm, double opbh){
	double fpbh = opbh/oc;
	return -12.2764 + log10(fpbh) - 3./2. * lgm;
}

// low mass - Eq.2.2 of draft
//  Psi = 10^psilowlg  as a function of Log10[M] = lgm and peak mass
double psilowlg(double lgm, double peakm){
	double aofmf = (peakm+2.723)/(-0.6576);
	double expt = 2.85*(lgm-peakm);
	return 2.85*(aofmf + lgm) - pow(10,expt)/2.3026;
}

// full spectrum - polychromatic
// Psi = 10^psibroadlg as a function of Log10[M] = lgm and peak mass
double psibroadlg(double lgm, double peakm){
	//  if mass is less than planck mass, set to 0
	if(lgmp>lgm){
		 return -100.;
		}
	// if above peak mass, set to hi spectrum
	else if(lgm>peakm){
		return psihilg(lgm,oc);
		}
	//  if below peak mass, set to low spectrum
	else{
		return psilowlg(lgm,peakm);
	}
}

// calculate volume fration given number density
double epsilon(double lgni, double lgmass){
	double radius = 2. * gnewton * pow(10,lgmass)/pow(sofl,2);
	printf("%e \n", radius);
	return pow(10.,lgni) * 4. * M_PI * pow(radius,3) / 3.;
}

/* Entropy density */
// Temperature given in Kelvin
double entropy(double Trh,  double  a, double arh){
	return 2.*pow(M_PI,2)/45 * gstar * pow(Trh*keltgev*arh/a,3.);
}

// absolute maximum theoretical number density
// Some notes: A bh of mass 10^-5.96 kg has a volume of 1.76 x 10^-98 m^3  and so we can have number densities mathematically up to 10^98 x 0.74 where 0.74 is the close packing ratio (max ratio of cube to volume of packed spheres)
double maxn0(double  lgmass){
	double rsch = 2.*gnewton*pow(10.,lgmass)/pow(sofl,2); // schwarzschild radius in m
	double volume =4./3. * M_PI * pow(rsch,3);
	return 1./volume * 0.74;
}

max_stars_repo_name: nebblu/PBH

max_issues_repo_name: nebblu/PBH

max_forks_repo_name: nebblu/PBH


text:
// Copyright (c) 2013-2016 Anton Kozhevnikov, Thomas Schulthess
// 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.

/** \file sirius_internal.h
 *
 *  \brief Contains basic definitions and declarations.
 */

#ifndef __SIRIUS_INTERNAL_H__
#define __SIRIUS_INTERNAL_H__

#include <omp.h>
#include <signal.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf_bessel.h>
#include <fftw3.h>
#include <vector>
#include <complex>
#include <iostream>
#include <algorithm>
#include "config.h"
#include "communicator.hpp"
#include "runtime.h"
#include "sddk.hpp"
#include "utils.h"
#ifdef __MAGMA
#include "GPU/magma.hpp"
#endif
#include "constants.h"
#include "version.h"
#include "simulation_context.h"
#include "Beta_projectors/beta_projectors_base.h"

#ifdef __PLASMA
extern "C" void plasma_init(int num_cores);
#endif

#ifdef __LIBSCI_ACC
extern "C" void libsci_acc_init();
extern "C" void libsci_acc_finalize();
#endif

/// Namespace of the SIRIUS library.
namespace sirius {

    inline void initialize(bool call_mpi_init__ = true)
    {
        if (call_mpi_init__) {
            Communicator::initialize();
        }
        if (mpi_comm_world().rank() == 0) {
            printf("SIRIUS %i.%i, git hash: %s\n", major_version, minor_version, git_hash);
        }

        #ifdef __GPU
        if (acc::num_devices()) {
            acc::create_streams(omp_get_max_threads() + 1);
            cublas::create_stream_handles();
        }
        #endif
        #ifdef __MAGMA
        magma::init();
        #endif
        #ifdef __PLASMA
        plasma_init(omp_get_max_threads());
        #endif
        #ifdef __LIBSCI_ACC
        libsci_acc_init();
        #endif

        assert(sizeof(int) == 4);
        assert(sizeof(double) == 8);
    }

    inline void finalize(bool call_mpi_fin__ = true)
    {
        #ifdef __MAGMA
        magma::finalize();
        #endif
        #ifdef __LIBSCI_ACC
        libsci_acc_finalize();
        #endif

        Beta_projectors_base<1>::cleanup();
        Beta_projectors_base<3>::cleanup();
        Beta_projectors_base<9>::cleanup();

        #ifdef __GPU
        if (acc::num_devices()) {
            cublas::destroy_stream_handles();
            acc::destroy_streams();
            acc::reset();
        }
        #endif
        fftw_cleanup();

        json dict;
        dict["flat"] = sddk::timer::serialize_timers();
        dict["tree"] = sddk::timer::serialize_timers_tree();
        if (mpi_comm_world().rank() == 0) {
            std::ofstream ofs("timers.json", std::ofstream::out | std::ofstream::trunc);
            ofs << dict.dump(4);
        }

        //sddk::timer::print_tree();
        if (call_mpi_fin__) {
            Communicator::finalize();
        }

    }

     inline void terminate(int err_code__)
     {
        MPI_Abort(MPI_COMM_WORLD, err_code__);
     }
};

#endif // __SIRIUS_INTERNAL_H__

/**
\mainpage Welcome to SIRIUS
\section intro Introduction
SIRIUS is a domain-specific library for electronic structure calculations. It supports full-potential linearized
augmented plane wave (FP-LAPW) and pseudopotential plane wave (PP-PW) methods and is designed to work with codes
such as Exciting, Elk and Quantum ESPRESSO.
\section install Installation
First, you need to clone the source code:
\verbatim
git clone https://github.com/electronic-structure/SIRIUS.git
\endverbatim

Then you need to create a configuration \c json file where you specify your compiliers, compiler flags and libraries.
Examples of such configuration files can be found in the <tt>./platforms/</tt> folder. The following variables have to be
provided:
  - \c MPI_CXX -- the MPI wrapper for the C++11 compiler (this is the main compiler for the library)
  - \c MPI_CXX_OPT -- the C++ compiler options for the library
  - \c MPI_FC -- the MPI wrapper for the Fortran compiler (used to build ELPA and SIRIUS F90 interface)
  - \c MPI_FC_OPT -- Fortran compiler options
  - \c CC -- plain C compilers (used to build the external libraries)
  - \c CXX -- plain C++ compiler (used to build the external libraries)
  - \c FC -- plain Fortran compiler (used to build the external libraries)
  - \c FCCPP -- Fortran preprocessor (usually 'cpp', required by LibXC package)
  - \c SYSTEM_LIBS -- list of the libraries, necessary for the linking (typically BLAS/LAPACK/ScaLAPACK, libstdc++ and Fortran run-time)
  - \c install -- list of packages to download, configure and build

In addition, the following variables can also be specified:
  - \c CUDA_ROOT -- path to the CUDA toolkit (if you compile with GPU support)
  - \c NVCC -- name of the CUDA C++ compiler (usually \c nvcc)
  - \c NVCC_OPT -- CUDA compiler options
  - \c MAGMA_ROOT -- location of the compiled MAGMA library (if you compile with MAGMA support)

Below is an example of the configurtaion file for the Cray XC50 platform. Cray compiler wrappers for C/C++/Fortran are
ftn/cc/CC. The GNU compilers and MKL are used.
\verbatim
{
    "comment"     : "MPI C++ compiler and options",
    "MPI_CXX"     : "CC",
    "MPI_CXX_OPT" : "-std=c++11 -Wall -Wconversion -fopenmp -D__SCALAPACK -D__ELPA -D__GPU -D__MAGMA -I$(MKLROOT)/include/fftw/",

    "comment"     : "MPI Fortran compiler and oprions",
    "MPI_FC"      : "ftn",
    "MPI_FC_OPT"  : "-O3 -fopenmp -cpp",

    "comment"     : "plain C compler",
    "CC"          : "cc",

    "comment"     : "plain C++ compiler",
    "CXX"         : "CC",

    "comment"     : "plain Fortran compiler",
    "FC"          : "ftn",

    "comment"     : "Fortran preprocessor",
    "FCCPP"       : "cpp",

    "comment"     : "location of CUDA toolkit",
    "CUDA_ROOT"   : "$(CUDATOOLKIT_HOME)",

    "comment"     : "CUDA compiler and options",
    "NVCC"        : "nvcc",
    "NVCC_OPT"    : "-arch=sm_60 -m64 -DNDEBUG",

    "comment"     : "location of MAGMA library",
    "MAGMA_ROOT"  : "$(HOME)/src/daint/magma-2.2.0",

    "SYSTEM_LIBS" : "$(MKLROOT)/lib/intel64/libmkl_scalapack_lp64.a -Wl,--start-group  $(MKLROOT)/lib/intel64/libmkl_intel_lp64.a $(MKLROOT)/lib/intel64/libmkl_gnu_thread.a $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_blacs_intelmpi_lp64.a -Wl,--end-group -lpthread -lstdc++ -ldl",

    "install"     : ["spg", "gsl", "xc"]
}
\endverbatim
FFT library is part of the MKL. The corresponding C++ include directory is passed to the compiler: -I\$(MKLROOT)/include/fftw/.
The \c HDF5 library is installed as a module and handeled by the Cray wrappers. The remaining three libraries necessary for
SIRIUS (spglib, GSL, libXC) are not available and have to be installed.

Once the configuration \c json file is created, you can run
\verbatim
python configure.py path/to/config.json
\endverbatim

The Python script will download and configure external packages, specified in the <tt>"install"</tt> list and create
Makefile and make.inc files. That's it! The configuration is done and you can run
\verbatim
make
\endverbatim
*/

//! \page stdvarname Standard variable names
//!
//! Below is the list of standard names for some of the loop variables:
//!
//! l - index of orbital quantum number \n
//! m - index of azimutal quantum nuber \n
//! lm - combined index of (l,m) quantum numbers \n
//! ia - index of atom \n
//! ic - index of atom class \n
//! iat - index of atom type \n
//! ir - index of r-point \n
//! ig - index of G-vector \n
//! idxlo - index of local orbital \n
//! idxrf - index of radial function \n
//! xi - compbined index of lm and idxrf (product of angular and radial functions) \n
//! ik - index of k-point \n
//! itp - index of (theta, phi) spherical angles \n
//!
//! The _loc suffix is often added to the variables to indicate that they represent the local fraction of the elements
//! assigned to the given MPI rank.
//!

//! \page coding Coding style
//!
//! Below are some basic style rules that we follow:
//!     - Page width is approximately 120 characters. Screens are wide nowdays and 80 characters is an
//!       obsolete restriction. Going slightly over 120 characters is allowed if it is requird for the line continuity.
//!     - Identation: 4 spaces (no tabs)
//!     - Coments are inserted before the code with slash-star style starting with the lower case:
//!       \code{.cpp}
//!           /* call a very important function */
//!           do_something();
//!       \endcode
//!     - Spaces between most operators:
//!       \code{.cpp}
//!           if (i < 5) {
//!               j = 5;
//!           }
//!
//!           for (int k = 0; k < 3; k++)
//!
//!           int lm = l * l + l + m;
//!
//!           double d = std::abs(e);
//!
//!           int k = idx[3];
//!       \endcode
//!     - Spaces between function arguments:
//!       \code{.cpp}
//!           double d = some_func(a, b, c);
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           double d=some_func(a,b,c);
//!       \endcode
//!       or
//!       \code{.cpp}
//!           double d = some_func( a, b, c );
//!       \endcode
//!     - Spaces between template arguments, but not between <> brackets:
//!       \code{.cpp}
//!           std::vector<std::array<int, 2>> vec;
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           std::vector< std::array< int, 2 > > vec;
//!       \endcode
//!     - Curly braces for classes and functions start form the new line:
//!       \code{.cpp}
//!           class A
//!           {
//!               ....
//!           };
//!
//!           inline int num_points()
//!           {
//!               return num_points_;
//!           }
//!       \endcode
//!     - Curly braces for if-statements, for-loops, switch-case statements, etc. start at the end of the line:
//!       \code{.cpp}
//!           for (int i: {0, 1, 2}) {
//!               some_func(i);
//!           }
//!
//!           if (a == 0) {
//!               printf("a is zero");
//!           } else {
//!               printf("a is not zero");
//!           }
//!
//!           switch (i) {
//!               case 1: {
//!                   do_something();
//!                   break;
//!               case 2: {
//!                   do_something_else();
//!                   break;
//!               }
//!           }
//!       \endcode
//!     - Even single line 'if' statements and 'for' loops must have the curly brackes:
//!       \code{.cpp}
//!           if (i == 4) {
//!               some_variable = 5;
//!           }
//!
//!           for (int k = 0; k < 10; k++) {
//!               do_something(k);
//!           }
//!       \endcode
//!     - Reference and pointer symbols are part of type:
//!       \code{.cpp}
//!           std::vector<double>& vec = make_vector();
//!
//!           double* ptr = &vec[0];
//!
//!           auto& atom = unit_cell().atom(ia);
//!       \endcode
//!     - Const modifier follows the type declaration:
//!       \code{.cpp}
//!           std::vector<int> const& idx() const
//!           {
//!               return idx_;
//!           }
//!       \endcode
//!     - Names of class members end with underscore:
//!       \code{.cpp}
//!           class A
//!           {
//!               private:
//!                   int lmax_;
//!           };
//!       \endcode
//!     - Setter method starts from set_, getter method is a variable name itself:
//!       \code{.cpp}
//!           class A
//!           {
//!               private:
//!                   int lmax_;
//!               public:
//!                   int lmax() const
//!                   {
//!                       return lmax_;
//!                   }
//!                   void set_lmax(int lmax__)
//!                   {
//!                       lmax_ = lmax__;
//!                   }
//!           };
//!       \endcode
//!     - Single-line functions should not be flattened:
//!       \code{.cpp}
//!           struct A
//!           {
//!               int lmax() const
//!               {
//!                   return lmax_;
//!               }
//!           };
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           struct A
//!           {
//!               int lmax() const { return lmax_; }
//!           };
//!       \endcode
//!     - Header guards have a standard name: double underscore + file name in capital letters + double underscore
//!       \code{.cpp}
//!           #ifndef __SIRIUS_INTERNAL_H__
//!           #define __SIRIUS_INTERNAL_H__
//!           ...
//!           #endif // __SIRIUS_INTERNAL_H__
//!       \endcode
//! We use clang-format utility to enforce the basic formatting style. Please have a look at .clang-format config file
//! in the source root folder for the definitions.
//!
//! Class naming convention.
//!
//! Problem: all 'standard' naming conventions are not satisfactory. For example, we have a class
//! which does a DFT ground state. Following the common naming conventions it could be named like this:
//! DFTGroundState, DftGroundState, dft_ground_state. Last two are bad, because DFT (and not Dft or dft)
//! is a well recognized abbreviation. First one is band because capital G adds to DFT and we automaticaly
//! read DFTG round state.
//!
//! Solution: we can propose the following: DFTgroundState or DFT_ground_state. The first variant still
//! doens't look very good because one of the words is captalized (State) and one (ground) - is not. So we pick
//! the second variant: DFT_ground_state (by the way, this is close to the Bjarne Stroustrup's naiming convention,
//! where he uses first capital letter and underscores, for example class Io_obj).
//!
//! Some other examples:
//!     - class Ground_state (composed of two words)
//!     - class FFT_interface (composed of an abbreviation and a word)
//!     - class Interface_XC (composed of a word and abbreviation)
//!     - class Spline (single word)
//!
//! Exceptions are allowed if it makes sense. For example, low level utility classes like 'mdarray' (multi-dimentional
//! array) or 'pstdout' (parallel standard output) are named with small letters.
//!

/** \page fderiv Functional derivatives

    Definition:
    \f[
      \frac{dF[f+\epsilon \eta ]}{d \epsilon}\Bigg\rvert_{\epsilon = 0} := \int \frac{\delta F[f]}{\delta f(x')} \eta(x') dx'
    \f]
    Alternative definition is:
    \f[
      \frac{\delta F[f(x)]}{\delta f(x')} = \lim_{\epsilon \to 0} \frac{F[f(x) + \epsilon \delta(x-x')] - F[f(x)]}{\epsilon}
    \f]

*/