text:
/*
 * test_shaw.c
 *
 * Test L-curve (Tikhonov) regression routines using Shaw
 * problem. See example 1.10 of
 *
 * [1] R.C. Aster, B. Borchers and C. H. Thurber,
 *     Parameter Estimation and Inverse Problems (2nd ed), 2012.
 */

#include <gsl/gsl_sf_trig.h>

/* alternate (and inefficient) method of computing G(lambda) */
static double
shaw_gcv_G(const double lambda, const gsl_matrix * X, const gsl_vector * y,
           gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;
  gsl_matrix * XTX = gsl_matrix_alloc(p, p);
  gsl_matrix * XI = gsl_matrix_alloc(p, n);
  gsl_matrix * XXI = gsl_matrix_alloc(n, n);
  gsl_vector * c = gsl_vector_alloc(p);
  gsl_vector_view d;
  double rnorm, snorm;
  double term1, term2, G;
  size_t i;

  /* compute regularized solution with this lambda */
  gsl_multifit_linear_solve(lambda, X, y, c, &rnorm, &snorm, work);

  /* compute X^T X */
  gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, X, 0.0, XTX);

  /* add lambda*I */
  d = gsl_matrix_diagonal(XTX);
  gsl_vector_add_constant(&d.vector, lambda * lambda);

  /* invert (X^T X + lambda*I) */
  gsl_linalg_cholesky_decomp1(XTX);
  gsl_linalg_cholesky_invert(XTX);
  gsl_matrix_transpose_tricpy('L', 0, XTX, XTX);

  /* XI = (X^T X + lambda*I)^{-1} X^T */
  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, XTX, X, 0.0, XI);

  /* XXI = X (X^T X + lambda*I)^{-1} X^T */
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, X, XI, 0.0, XXI);

  /* compute: term1 = Tr(I - X XI) */
  term1 = 0.0;
  for (i = 0; i < n; ++i)
    {
      double *Ai = gsl_matrix_ptr(XXI, i, i);
      term1 += 1.0 - (*Ai);
    }

  gsl_matrix_free(XTX);
  gsl_matrix_free(XI);
  gsl_matrix_free(XXI);
  gsl_vector_free(c);

  term2 = rnorm / term1;

  return term2 * term2;;
}

/* construct design matrix and rhs vector for Shaw problem */
static int
shaw_system(gsl_matrix * X, gsl_vector * y)
{
  int s = GSL_SUCCESS;
  const size_t n = X->size1;
  const size_t p = X->size2;
  const double dtheta = M_PI / (double) p;
  size_t i, j;
  gsl_vector *m = gsl_vector_alloc(p);

  /* build the design matrix */
  for (i = 0; i < n; ++i)
    {
      double si = (i + 0.5) * M_PI / n - M_PI / 2.0;
      double csi = cos(si);
      double sni = sin(si);

      for (j = 0; j < p; ++j)
        {
          double thetaj = (j + 0.5) * M_PI / p - M_PI / 2.0;
          double term1 = csi + cos(thetaj);
          double term2 = gsl_sf_sinc(sni + sin(thetaj));
          double Xij = term1 * term1 * term2 * term2 * dtheta;

          gsl_matrix_set(X, i, j, Xij);
        }
    }

  /* construct coefficient vector */
  {
    const double a1 = 2.0;
    const double a2 = 1.0;
    const double c1 = 6.0;
    const double c2 = 2.0;
    const double t1 = 0.8;
    const double t2 = -0.5;

    for (j = 0; j < p; ++j)
      {
        double tj = -M_PI / 2.0 + (j + 0.5) * dtheta;
        double mj = a1 * exp(-c1 * (tj - t1) * (tj - t1)) +
                    a2 * exp(-c2 * (tj - t2) * (tj - t2));
        gsl_vector_set(m, j, mj);
      }
  }

  /* construct rhs vector */
  gsl_blas_dgemv(CblasNoTrans, 1.0, X, m, 0.0, y);

  gsl_vector_free(m);

  return s;
}

static int
test_shaw_system_l(gsl_rng *rng_p, const size_t n, const size_t p,
                   const double lambda_expected,
                   gsl_vector *rhs)
{
  const size_t npoints = 1000; /* number of points on L-curve */
  const double tol1 = 1.0e-12;
  const double tol2 = 1.0e-9;
  const double tol3 = 1.0e-5;
  gsl_vector * reg_param = gsl_vector_alloc(npoints);
  gsl_vector * rho = gsl_vector_alloc(npoints);
  gsl_vector * eta = gsl_vector_alloc(npoints);

  gsl_matrix * X = gsl_matrix_alloc(n, p);
  gsl_matrix * cov = gsl_matrix_alloc(p, p);
  gsl_vector * c = gsl_vector_alloc(p);
  gsl_vector * ytmp = gsl_vector_alloc(n);
  gsl_vector * y;
  gsl_vector * r = gsl_vector_alloc(n);
  gsl_multifit_linear_workspace * work =
    gsl_multifit_linear_alloc (n, p);

  size_t reg_idx, i;
  double lambda, rnorm, snorm;

  /* build design matrix */
  shaw_system(X, ytmp);

  if (rhs)
    y = rhs;
  else
    {
      y = ytmp;

      /* add random noise to exact rhs vector */
      test_random_vector_noise(rng_p, y);
    }

  /* SVD decomposition */
  gsl_multifit_linear_svd(X, work);

  /* calculate L-curve */
  gsl_multifit_linear_lcurve(y, reg_param, rho, eta, work);

  /* test rho and eta vectors */
  for (i = 0; i < npoints; ++i)
    {
      double rhoi = gsl_vector_get(rho, i);
      double etai = gsl_vector_get(eta, i);
      double lami = gsl_vector_get(reg_param, i);

      /* solve regularized system and check for consistent rho/eta values */
      gsl_multifit_linear_solve(lami, X, y, c, &rnorm, &snorm, work);
      gsl_test_rel(rhoi, rnorm, tol3, "shaw rho n=%zu p=%zu lambda=%e",
                   n, p, lami);
      gsl_test_rel(etai, snorm, tol1, "shaw eta n=%zu p=%zu lambda=%e",
                   n, p, lami);
    }

  /* calculate corner of L-curve */
  gsl_multifit_linear_lcorner(rho, eta, &reg_idx);

  lambda = gsl_vector_get(reg_param, reg_idx);

  /* test against known lambda value if given */
  if (lambda_expected > 0.0)
    {
      gsl_test_rel(lambda, lambda_expected, tol1,
                   "shaw: n=%zu p=%zu L-curve corner lambda",
                   n, p);
    }

  /* compute regularized solution with optimal lambda */
  gsl_multifit_linear_solve(lambda, X, y, c, &rnorm, &snorm, work);

  /* compute residual norm ||y - X c|| */
  gsl_vector_memcpy(r, y);
  gsl_blas_dgemv(CblasNoTrans, 1.0, X, c, -1.0, r);

  /* test rnorm value */
  gsl_test_rel(rnorm, gsl_blas_dnrm2(r), tol2,
               "shaw: n=%zu p=%zu rnorm", n, p);

  /* test snorm value */
  gsl_test_rel(snorm, gsl_blas_dnrm2(c), tol2,
               "shaw: n=%zu p=%zu snorm", n, p);

  gsl_matrix_free(X);
  gsl_matrix_free(cov);
  gsl_vector_free(reg_param);
  gsl_vector_free(rho);
  gsl_vector_free(eta);
  gsl_vector_free(r);
  gsl_vector_free(c);
  gsl_vector_free(ytmp);
  gsl_multifit_linear_free(work);

  return 0;
} /* test_shaw_system_l() */

static int
test_shaw_system_gcv(gsl_rng *rng_p, const size_t n, const size_t p,
                     const double lambda_expected,
                     gsl_vector *rhs)
{
  const size_t npoints = 200; /* number of points on L-curve */
  const double tol1 = 1.0e-12;
  const double tol2 = 1.4e-10;
  const double tol3 = 1.0e-5;
  gsl_vector * reg_param = gsl_vector_alloc(npoints);
  gsl_vector * G = gsl_vector_alloc(npoints);

  gsl_matrix * X = gsl_matrix_alloc(n, p);
  gsl_matrix * cov = gsl_matrix_alloc(p, p);
  gsl_vector * c = gsl_vector_alloc(p);
  gsl_vector * ytmp = gsl_vector_alloc(n);
  gsl_vector * y;
  gsl_vector * r = gsl_vector_alloc(n);
  gsl_multifit_linear_workspace * work =
    gsl_multifit_linear_alloc (n, p);

  size_t reg_idx, i;
  double lambda, rnorm, snorm, G_lambda;

  /* build design matrix */
  shaw_system(X, ytmp);

  if (rhs)
    y = rhs;
  else
    {
      y = ytmp;

      /* add random noise to exact rhs vector */
      test_random_vector_noise(rng_p, y);
    }

  /* SVD decomposition */
  gsl_multifit_linear_svd(X, work);

  /* calculate GCV curve */
  gsl_multifit_linear_gcv(y, reg_param, G, &lambda, &G_lambda, work);

  /* test G vector */
  for (i = 0; i < npoints; ++i)
    {
      double lami = gsl_vector_get(reg_param, i);

      if (lami > 1.0e-5)
        {
          /* test unreliable for small lambda */
          double Gi = gsl_vector_get(G, i);
          double Gi_expected = shaw_gcv_G(lami, X, y, work);

          gsl_test_rel(Gi, Gi_expected, tol3, "shaw[%zu,%zu] gcv G i=%zu lambda=%e",
                       n, p, i, lami);
        }
    }

  /* test against known lambda value if given */
  if (lambda_expected > 0.0)
    {
      gsl_test_rel(lambda, lambda_expected, tol2,
                   "shaw gcv: n=%zu p=%zu lambda",
                   n, p);
    }

  /* compute regularized solution with optimal lambda */
  gsl_multifit_linear_solve(lambda, X, y, c, &rnorm, &snorm, work);

  /* compute residual norm ||y - X c|| */
  gsl_vector_memcpy(r, y);
  gsl_blas_dgemv(CblasNoTrans, 1.0, X, c, -1.0, r);

  /* test rnorm value */
  gsl_test_rel(rnorm, gsl_blas_dnrm2(r), tol2,
               "shaw gcv: n=%zu p=%zu rnorm", n, p);

  /* test snorm value */
  gsl_test_rel(snorm, gsl_blas_dnrm2(c), tol2,
               "shaw gcv: n=%zu p=%zu snorm", n, p);

  gsl_matrix_free(X);
  gsl_matrix_free(cov);
  gsl_vector_free(reg_param);
  gsl_vector_free(G);
  gsl_vector_free(r);
  gsl_vector_free(c);
  gsl_vector_free(ytmp);
  gsl_multifit_linear_free(work);

  return 0;
} /* test_shaw_system_gcv() */

void
test_shaw(void)
{
  gsl_rng * r = gsl_rng_alloc(gsl_rng_default);

  {
    double shaw20_y[] = {
      8.7547455124379323e-04, 5.4996835885761936e-04,
      1.7527999407005367e-06, 1.9552372913117047e-03,
      1.4411645433785081e-02, 5.2800013336393704e-02,
      1.3609152023257112e-01, 2.7203484587635818e-01,
      4.3752225136193390e-01, 5.7547667319875240e-01,
      6.2445052213539942e-01, 5.6252658286441348e-01,
      4.2322239923561566e-01, 2.6768469219560631e-01,
      1.4337901162734543e-01, 6.5614569346074361e-02,
      2.6013851831752945e-02, 9.2336933089481269e-03,
      3.2269066658993694e-03, 1.3999201459261811e-03 };
    gsl_vector_view rhs = gsl_vector_view_array(shaw20_y, 20);

    /* lambda and rhs values from [1] */
    test_shaw_system_l(r, 20, 20, 5.793190958069266e-06, &rhs.vector);

    test_shaw_system_gcv(r, 20, 20, 1.24921780949051038e-05, &rhs.vector);
  }

  {
    size_t n, p;

    for (n = 10; n <= 50; n += 2)
      {
        for (p = n - 6; p <= n; p += 2)
          test_shaw_system_l(r, n, p, -1.0, NULL);
      }
  }

  gsl_rng_free(r);
} /* test_shaw() */


text:
// lu2t.c
//
// test program for blocked LU decomposition
//
// Time-stamp: <09/08/31 00:41:24 makino>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef NOBLAS
#include <cblas.h>
#endif
#include <lu2tlib.h>
void timer_init();
double cpusec();
double wsec();
#define FTYPE double
void copymats( int n, double a[n+1][n], double a2[n+1][n])

{
    int i, j;
    for(j=0;j<n+1;j++)
	for(i=0;i<n;i++)
	    a2[j][i] = a[j][i];
}

void showresult(int n, double a[n+1][n], double x[])
{
    int i, j;
    double emax = 0;
    for(i=0;i<n;i++){
	int k;
	double b2=0;
	printf("%3d: ", i);
	//	for(j=0;j<n;j++) printf(" %10.3e", a[j][i]);
	for(j=0;j<n;j++) b2 += a[j][i] * x[j];
	double err = b2-a[n][i];
        emax =  (fabs(err) > emax) ? fabs(err):emax;
	printf(" %10.3e  %10.3e %10.3e %10.3e \n", x[i], a[n][i], b2, err);
    }
    printf("Emax= %10.3e\n", emax);
}

void readmat( int n, double a[n+1][n])
{
    int i, j;
    for(i=0;i<n;i++){
	for(j=0;j<n+1;j++) scanf("%le", &(a[j][i]));
    }
}

void printmat( int n, double a[n+1][n])

{
    int i, j;
    for(i=0;i<n;i++){
	printf("%3d: ", i);
	for(j=0;j<n+1;j++) printf(" %10.3e", a[j][i]);
	printf("\n");
    }
    printf("\n");
}

void backward_sub(int n,double a[n+1][n], double b[])
{
    int i,j,k;
    for (i=0;i<n;i++)b[i] = a[n][i];
    for(j=n-2;j>=0;j--)
	for(k=j+1;k<n;k++) b[j] -= b[k]*a[k][j];
}

void lu( int n, double a[n+1][n], double b[])
{
    int i, j, k;
    for(i=0;i<n-1;i++){
	// select pivot
	double amax = fabs(a[i][i]);
	int p=i;
	for(j=i+1;j<n;j++){
	    if (fabs(a[i][j]) > amax){
		amax = fabs(a[i][j]);
		p = j;
	    }
	}
	// exchange rows
	if (p != i){
	    for(j=i;j<n+1;j++){
		double tmp = a[j][p];
		a[j][p] = a[j][i];
		a[j][i]=tmp;
	    }
	}

	// normalize row i
	double ainv = 1.0/a[i][i];
	//	fprintf(stderr,"%d %e\n", i, ainv);
	for(k=i+1;k<n+1;k++) a[k][i]*= ainv;
	// subtract row i from all lower rows
	for(j=i+1;j<n;j++){
	    //	    fprintf(stderr,"j=%d \n",j);
	    for(k=i+1;k<n+1;k++) a[k][j] -= a[i][j] * a[k][i];
	}
    }
    printmat(n,a);

    a[n][n-1] /= a[n-1][n-1];
    backward_sub(n,a,b);
}

void lumcolumn( int n, double a[n+1][n], double b[], int m, int pv[],
		int recursive)
{
    int i;
    for(i=0;i<n;i+=m){
	if (recursive){
	    cm_column_decomposition_recursive(n, a+i, m, pv,i);
	}else{
	    cm_column_decomposition(n, a+i, m, pv,i);
	}
	cm_process_right_part(n,a+i,m,pv,i,n+1);
    }
    backward_sub(n,a,b);
}

main()
{
    int n;
    fprintf(stderr, "Enter n:");
    scanf("%d", &n);
    printf("N=%d\n", n);
#if 0
    double a[n+1][n];
    double b[n];
    double acopy[n+1][n];
    double bcopy[n];
    int pv[n];
#endif
    double (*a)[];
    double (*acopy)[];
    int * pv;
    double *b, *bcopy;
    a = (double(*)[]) malloc(sizeof(double)*n*(n+1));
    acopy = (double(*)[]) malloc(sizeof(double)*n*(n+1));
    b = (double*)malloc(sizeof(double)*n);
    bcopy = (double*)malloc(sizeof(double)*n);
    pv = (int*)malloc(sizeof(int)*n);
    readmat(n,a);
    copymats(n,a,acopy);
    //    printmat(n,a,b);
    //    lu2columnv2(n,a,b);
    //    lu2columnv2(n,a,b);
    //    lub(n,a,b,NBK);
    //    printmat(n,a,b);
    //    showresult(n,acopy, b, bcopy);
    //    copymats(n,acopy,bcopy,a, b);
    //    lu(n,a,b);
    timer_init();
    lumcolumn(n,a,b,16,pv,1);
    double ctime=cpusec();
    double wtime=wsec();
    showresult(n,acopy, b);
    printf("cpsec =  %g wsec=%g\n",   ctime, wtime);
    return 0;
}

ext: c


text:
#pragma once
#include "iflag.h"
#include "optioninfo.h"
#include "configaccess.h"
#include "format.h"
#include <gsl/gsl>
#include <cmdlime/customnames.h>
#include <memory>
#include <functional>
#include <utility>

namespace cmdlime::detail{

class Flag : public IFlag{
public:
    enum class Type{
        Normal,
        Exit
    };

    Flag(std::string name,
         std::string shortName,
         std::function<bool&()> flagGetter,
         Type type)
        : info_(std::move(name), std::move(shortName), {})
        , flagGetter_(std::move(flagGetter))
        , type_(type)
    {
    }

    OptionInfo& info() override
    {
        return info_;
    }

    const OptionInfo& info() const override
    {
        return info_;
    }

private:
    void set() override
    {
        flagGetter_() = true;
    }

    bool isSet() const override
    {
        return flagGetter_();
    }

    bool isExitFlag() const override
    {
        return type_ == Type::Exit;
    }

private:
    OptionInfo info_;
    std::function<bool&()> flagGetter_;
    Type type_;
};

template <typename TConfig>
class FlagCreator{
    using NameProvider = typename Format<ConfigAccess<TConfig>::format()>::nameProvider;

public:
    FlagCreator(TConfig& cfg,
                const std::string& varName,
                std::function<bool&()> flagGetter,
                Flag::Type flagType = Flag::Type::Normal)
        : cfg_(cfg)
    {
        Expects(!varName.empty());
        flag_ = std::make_unique<Flag>(NameProvider::name(varName),
                                       NameProvider::shortName(varName),
                                       std::move(flagGetter),
                                       flagType);
    }

    FlagCreator& operator<<(const std::string& info)
    {
        flag_->info().addDescription(info);
        return *this;
    }

    FlagCreator& operator<<(const Name& customName)
    {
        flag_->info().resetName(customName.value());
        return *this;
    }

    FlagCreator& operator<<(const ShortName& customName)
    {
        static_assert(Format<ConfigAccess<TConfig>::format()>::shortNamesEnabled,
                      "Current command line format doesn't support short names");
        flag_->info().resetShortName(customName.value());
        return *this;
    }

    FlagCreator& operator<<(const WithoutShortName&)
    {
        static_assert(Format<ConfigAccess<TConfig>::format()>::shortNamesEnabled,
                      "Current command line format doesn't support short names");
        flag_->info().resetShortName({});
        return *this;
    }

    operator bool()
    {
        ConfigAccess<TConfig>{cfg_}.addFlag(std::move(flag_));
        return false;
    }

private:
    std::unique_ptr<Flag> flag_;
    TConfig& cfg_;
};
}

ext: h

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef clip_c55b09bc_5fbd_463f_8d45_b12013b57f42_h
#define clip_c55b09bc_5fbd_463f_8d45_b12013b57f42_h

#include <gslib/tree.h>
#include <ariel/painterpath.h>

__ariel_begin__

enum clip_point_tag
{
    cpt_corner = 0x01,
    cpt_interpolate = 0x02,
};

enum clip_fill_type
{
    cft_even_odd = 0,
    cft_non_zero,
    cft_positive,
    cft_negative,
};

struct clip_vec2_hash
{
    size_t operator()(const vec2& pt) const { return hash_bytes((const byte*)&pt, sizeof(pt)); }
};

typedef _treenode_wrapper<painter_path> clip_result_wrapper;
typedef tree<painter_path, clip_result_wrapper> clip_result;
typedef typename clip_result::iterator clip_result_iter;
typedef typename clip_result::const_iterator clip_result_const_iter;
typedef unordered_map<vec2, uint, clip_vec2_hash> clip_point_attr;

ariel_export extern void clip_remap_points(painter_linestrips& output, clip_point_attr& attrmap, const painter_path& input, uint attr_selector, float step_len = -1.f);
ariel_export extern void clip_offset(painter_linestrips& lss, const painter_linestrips& input, float offset);
ariel_export extern void clip_stroker(painter_linestrips& lss, const painter_linestrips& input, float offset);
ariel_export extern void clip_simplify(painter_linestrips& lss, const painter_linestrip& input, clip_fill_type ft = cft_even_odd);
ariel_export extern void clip_simplify(painter_linestrips& lss, const painter_linestrips& input, clip_fill_type ft = cft_even_odd);
ariel_export extern void clip_simplify(clip_result& output, const painter_path& input);
ariel_export extern void clip_union(clip_result& output, const painter_path& subjects, const painter_path& clips);
ariel_export extern void clip_intersect(clip_result& output, const painter_path& subjects, const painter_path& clips);
ariel_export extern void clip_substract(clip_result& output, const painter_path& subjects, const painter_path& clips);
ariel_export extern void clip_exclude(clip_result& output, const painter_path& subjects, const painter_path& clips);
ariel_export extern void clip_convert(painter_path& path, const clip_result& result);

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
FILE DETAILS
description: header file needed for matrix operations in LCP solver, and some macros for the main
LCP function
project: MPT 3.0
filename: lcp_matrix.h
author: Colin N. Jones, Automatic Control Laboratory, ETH Zurich, 2006

LICENSE:
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.

REVISION HISTORY:
date: Nov, 2010
revised by: Martin Herceg, Automatic Control Laboratory, ETH Zurich, 2010
details: Included BLAS and LAPACK libraries which are standard in MATLAB. Replaced int
declarations with ptrdiff_t to be 32/64 compatible. Added declaration of LCP_options structure.
*/

#ifndef LCP_MATRIX__H
#define LCP_MATRIX__H

/**** MODULES ****/
/* Standard headers */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stddef.h> /* to have ptrdiff_t */
#include <math.h>
#include <time.h>

/* BLAS and LAPACK headers are in $MATLABROOT/extern/include
and are on the default search path when compiling with standard options */
#include <blas.h>
#include <lapack.h>

/**** MACROS ****/

/* Select an element from a column-major matrix */
#define C_SEL(mat,r,c) ((mat)->A[(r)+(c)*((mat)->rows_alloc)])

/* Select an element from a row-major matrix */
#define R_SEL(mat,r,c) ((mat)->A[(c)+(r)*((mat)->cols_alloc)])

/* Get a matrix */
#define MAT(mat) ((mat).A)

/* Get a matrix pointer */
#define pMAT(mat) ((mat)->A)

/* Get a pointer to a given column in a matrix */
#define COL(mat,n) &([C_SEL(mat,0,n)])

/* Get number of rows of a given matrix */
#define Matrix_Rows(pA) ((pA)->rows)

/* Get number of columns of a given matrix */
#define Matrix_Cols(pA) ((pA)->cols)

/* Dimension of an index set */
#define Index_Length(pI) ((pI)->len)

/* Assign value "x" to position "i" in the index set "I" */
#define Index_Set(pI,i,x) (((pI)->I[(i)]) = (x))

/* Obtain value of the index set "I" at the position "i" */
#define Index_Get(pI,i) ((pI)->I[(i)])

/**** TYPEDEFS AND STRUCTURES ****/

/* Matrices */
/* If a matrix is column-major, then the form is: */
/*      [M(0) M(rows_alloc)   M(2rows_alloc)   ... ] */
/*      [M(1) M(rows_alloc+1) M(2rows_alloc+1) ... ] */
/*      [...  ...             ....             ... ] */
typedef struct LCP_Matrix {
	double *A;
	ptrdiff_t rows,cols;
	ptrdiff_t rows_alloc, cols_alloc;
	char *name;
} T_Matrix, *PT_Matrix;

/* Index sets for basis defintion */
typedef struct LCP_IndexSet {
	/* C-based indexing used throughout. i.e. 0 is the first element */
	ptrdiff_t *I;
	ptrdiff_t len, len_alloc;
	char *name;
} T_Index, *PT_Index;

/* The basis */
typedef struct LCP_Basis {
	PT_Index pW;  /* Slacks in the basis */
	PT_Index pZ;  /* Variables in the basis */
	PT_Index pWc; /* Complement of W */

	PT_Index pP;  /* Work index (used to list positive variables) */

	/* B      = [I G;0 X]; */
	/* inv(B) = [I -G*iX; 0 iX] */
	/* G = M_W,Z */
	/* X = M_Wc,X */

	PT_Matrix pG, pX, pF, pLX, pUX, pUt;  /* F is used as factorization of X on exit from dgesv routine,
						                         Ut is full, column major upper triangular matrix */
	double *y,*z,*w, *v, *r;    /* work vectors */
	ptrdiff_t *p; /* work vector for dgesv, dgetrs */
	double *s; /* work vector for dgels */
} T_Basis, *PT_Basis;

/* options structure */
typedef struct LCP_options {
	double zerotol; /* zero tolerance*/
	double lextol; /* lexicographic tolerance - a small treshold to determine if values are equal */
	ptrdiff_t maxpiv; /* maximum number of pivots */
	ptrdiff_t nstepf; /* after n steps refactorize the basis */
	int clock;  /* print computation time */
	int verbose; /* verbose mode */
	int routine; /* routine from LAPACK library to solve linear equations in Basis_solve*/
	double timelimit; /* time limit for iterations in seconds */
	int normalize; /* normalize matrices M, q */
    double normalizethres; /* normalization threshold */
} T_Options;

/**** FUNCTIONS ****/

/* Basis functions */
PT_Basis Basis_Init(ptrdiff_t m);
void Basis_Free(PT_Basis pB);
void Basis_Print(PT_Basis pB);
ptrdiff_t  Basis_Pivot(PT_Basis pB, PT_Matrix pA, ptrdiff_t enter, ptrdiff_t leave);
ptrdiff_t Basis_Solve(PT_Basis pB, double *q, double *x, T_Options options);
void Reinitialize_Basis(ptrdiff_t m, PT_Basis pB);

void LU_Expand(PT_Matrix L, PT_Matrix U, double *y, double *z, double *w);
void LU_Replace_Col(PT_Matrix L, PT_Matrix U, ptrdiff_t kcol, double *y, double *z, double *w);
void LU_Replace_Row(PT_Matrix L, PT_Matrix U, ptrdiff_t krow, double *y, double *z, double *w);
void LU_Shrink(PT_Matrix L, PT_Matrix U, ptrdiff_t krow, ptrdiff_t kcol, double *y, double *z, double *w);
void LU_Solve0(PT_Matrix pL, PT_Matrix pU, double *vec, double *x);
void LU_Solve1(PT_Matrix pL, PT_Matrix pU, double *vec, double *x, ptrdiff_t *info);
void LU_Refactorize(PT_Basis pB);

/* Matrix functions */
PT_Matrix Matrix_Init(ptrdiff_t m, ptrdiff_t n, const char *name);
PT_Matrix Matrix_Init_NoAlloc(ptrdiff_t m, ptrdiff_t n, const char *name);
void Matrix_Free(PT_Matrix pMat);
void Matrix_Free_NoAlloc(PT_Matrix pMat);
void Matrix_Print_col(PT_Matrix pMat);
void Matrix_Print_row(PT_Matrix pMat);
void Matrix_Print_utril_row(PT_Matrix pMat);
void Matrix_Print_size(PT_Matrix pMat);
void Matrix_AddRow(PT_Matrix pA, double *row);
void Matrix_AddCol(PT_Matrix pA, double *col);
void Matrix_RemoveRow(PT_Matrix pA, ptrdiff_t row);
void Matrix_RemoveCol(PT_Matrix pA, ptrdiff_t col);
void Matrix_GetRow(PT_Matrix pA, double *row, ptrdiff_t row_index, PT_Index col_index);
void Matrix_GetCol(PT_Matrix pA, double *col, PT_Index row_index, ptrdiff_t col_index);
void Matrix_ReplaceCol(PT_Matrix pA, ptrdiff_t kcol, double *vec);
void Matrix_ReplaceRow(PT_Matrix pA, ptrdiff_t krow, double *vec);
void Matrix_Copy(PT_Matrix pA, PT_Matrix pB, double *w);
void Matrix_Triu(PT_Matrix pF, PT_Matrix pU, double *w);
void Matrix_Uzeros(PT_Matrix pF);
void Matrix_Full2RowTriangular(PT_Matrix pF, PT_Matrix  pU, double *w);
void Matrix_Transpose(PT_Matrix pA, PT_Matrix pB, double *w);
void Vector_Print_raw(double *vec, ptrdiff_t m);

/* Index set functions */
PT_Index Index_Init(ptrdiff_t len, const char *name);
void Index_Free(PT_Index pIndex);
void Index_Print(PT_Index pI);
void Index_Remove(PT_Index pI, ptrdiff_t leave);
void Index_Add(PT_Index pI, ptrdiff_t num);
void Index_Replace(PT_Index pI, ptrdiff_t index, ptrdiff_t newval);
ptrdiff_t Index_FindElement(PT_Index pI, ptrdiff_t element);

#endif

max_stars_repo_name: adrilow/bary-polytopes

max_issues_repo_name: adrilow/bary-polytopes

max_forks_repo_name: adrilow/bary-polytopes


text:
/**
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_spline.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_PET.h"
#include "inout_aper.h"
//#include "trace_conf.h"
//#include "aper_conf.h"
//#include "spc_sex.h"
//#include "disp_conf.h"
//#include "spc_wl_calib.h"
//#include "spce_binning.h"
//#include "spc_spc.h"
//#include "spce_pgp.h"
//#include "spce_output.h"
//#include "spc_FITScards.h"
//#include "spc_flatfield.h"
#include "fringe_conf.h"
#include "fringe_utils.h"

#define AXE_IMAGE_PATH "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

int
main (int argc, char *argv[])
{
  char *opt;
  char aper_file[MAXCHAR];
  char aper_file_path[MAXCHAR];

  char conf_file[MAXCHAR];
  char conf_file_path[MAXCHAR];

  char fconf_file[MAXCHAR];
  char fconf_file_path[MAXCHAR];

  char grism_file[MAXCHAR];
  char grism_file_path[MAXCHAR];

  char FF_file[MAXCHAR];
  //char FF_file_path[MAXCHAR];

  char OPET_file[MAXCHAR];
  char OPET_file_path[MAXCHAR];

  char BPET_file[MAXCHAR];
  char BPET_file_path[MAXCHAR];

  aperture_conf *conf=NULL;
  fringe_conf   *fconf=NULL;

  object **oblist;
  observation *obs=NULL;

  ap_pixel *OPET=NULL;
  ap_pixel *BPET=NULL;

  fitsfile *OPET_ptr=NULL;
  fitsfile *BPET_ptr=NULL;

  beam act_beam;

  int index;
  int i;
  int f_status = 0;
  int bckmode = 0;
  int oaperID, obeamID;
  int baperID, bbeamID;

  if ((argc < 4) || (opt = get_online_option ("help", argc, argv)))
    {
      fprintf (stdout,
	       "aXe_FRINGECORR Version %s:\n",RELEASE);
      exit (1);
    }

  fprintf (stdout, "aXe_FRINGECORR: Starting...\n\n");

  // copy the first parameter to the grism filename
  index = 0;
  strcpy (grism_file, argv[++index]);
  build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);

  // copy the next parameter to the aXe configuration file
  strcpy (conf_file, argv[++index]);
  build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);

  // copy the next parameter to the fringe configuration file
  strcpy (fconf_file, argv[++index]);
  build_path (AXE_CONFIG_PATH, fconf_file, fconf_file_path);

  // Determine if we are using the special bck mode
  // In this mode, file names are handled diferently
  if ((opt = get_online_option("bck", argc, argv)))
    bckmode = 1;

  // load the aXe configuration file
  conf = get_aperture_descriptor (conf_file_path);

  // load the fringe configuration file
  fconf = load_fringe_conf(fconf_file_path);

  // Determine where the various extensions are in the FITS file
  get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

  // Build aperture file name
  replace_file_extension (grism_file, aper_file, ".fits",
			  ".OAF", conf->science_numext);
  build_path (AXE_OUTPUT_PATH, aper_file, aper_file_path);

  // Build object PET file name
  replace_file_extension (grism_file, OPET_file, ".fits",
			  ".PET.fits", conf->science_numext);
  build_path (AXE_OUTPUT_PATH, OPET_file, OPET_file_path);

  // build the background PET file name
  // if necessary
  if (bckmode)
    {
      replace_file_extension (grism_file, BPET_file, ".fits",
			      ".BCK.PET.fits", conf->science_numext);
      build_path (AXE_OUTPUT_PATH, BPET_file, BPET_file_path);
  }

  // Give a short feedback on the filenames
  // of the various input
  fprintf (stdout, "aXe_FRINGECORR: Grism image file name:          %s\n",
	   grism_file_path);
  fprintf (stdout, "aXe_FRINGECORR: Aperture file name:             %s\n",
	   aper_file_path);
  fprintf (stdout, "aXe_FRINGECORR: aXe configuration file name:    %s\n",
	   conf_file_path);
  fprintf (stdout, "aXe_FRINGECORR: Fringe configuration file name: %s\n",
	   fconf_file_path);
  fprintf (stdout, "aXe_FRINGECORR: Object PET file name:           %s\n",
	   OPET_file_path);
  if (bckmode)
    fprintf (stdout, "aXe_FRINGECORR: Background PET file name:       %s\n",
	     BPET_file_path);

  // check whether all necessary information
  // is in place, provide defaults
  check_fringe_conf(fconf);

  fprintf (stdout, "aXe_FRINGECORR: Loading object aperture list...");
  fflush(stdout);
  oblist = file_to_object_list_seq (aper_file_path, obs);
  fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));

  // Open the OPET file for reading/writing
  fits_open_file (&OPET_ptr, OPET_file_path, READWRITE, &f_status);
  if (f_status)
    {
      ffrprt (stdout, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_FRINGECORR: Could not open file: %s\n",
		   OPET_file_path);
    }

  // if necessary, open the background PET
  if (bckmode)
    {
      fits_open_file (&BPET_ptr, BPET_file_path, READWRITE, &f_status);
      if (f_status)
	{
	  ffrprt (stdout, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_FRINGECORR: Could not open file: %s\n",
		       BPET_file_path);
	}
    }

  i = 0;
  if ((oblist!=NULL) && (strcmp(FF_file,"None")))
    {
      while (1)
	{
	  // Get the PET for this object
	  OPET = get_ALL_from_next_in_PET(OPET_ptr, &oaperID, &obeamID);
	  if ((oaperID==-1) && (obeamID==-1))
	    break;

	  // Get the PET for this object
	  if (bckmode)
	    {
	      BPET = get_ALL_from_next_in_PET(BPET_ptr, &baperID, &bbeamID);
	      if ((baperID != oaperID) ||  (bbeamID != obeamID))
		{
		  fprintf(stderr, "%i, %i, %i, %i\n", oaperID, baperID, obeamID, bbeamID);
		  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			       "aXe_FRINGECORR: The beam order in the\n"
			       "object PET %s and the background PET\n"
			       "%s does not coincide!\n",
			       OPET_file_path, BPET_file_path);
		}
	    }

	  // check whether there is content in the PET;
	  // skip it if empty
	  if (OPET==NULL)
	    continue;

	  // refer about the actual beam in process
	  fprintf (stdout, "aXe_FRINGECORR: BEAM %d%c",oaperID, BEAM(obeamID));

	  // get the beam information from the OAF
	  act_beam = find_beam_in_object_list(oblist, oaperID, obeamID);
	  fprintf (stdout, " Done. \n");

	  fringe_correct_PET(fconf, act_beam, OPET, BPET);

	  // update PET table with the new FF info
	  {
	    char ID[60];
	    sprintf (ID, "%d%c", oaperID, BEAM (act_beam.ID));
	    add_ALL_to_PET (OPET, ID, OPET_ptr,1);

	    if (bckmode)
	      add_ALL_to_PET (BPET, ID, BPET_ptr,1);
            }

	  // free the memory of the object pixels
	  if (OPET!=NULL)
	    {
	      free(OPET);
	      OPET=NULL;
	    }

	  // free the memory of the background pixels
	  if (BPET!=NULL)
	    {
	      free(BPET);
	      BPET=NULL;
	    }
	}
    }

  // free the object list
  if (oblist!=NULL)
    free_oblist (oblist);

  // close the object PET fits-file
  fits_close_file (OPET_ptr, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "aXe_PETFF: " "Error closing PET: %s\n",
		   OPET_file_path);
    }

  // if necessary, close the object PET fits-file
  if (bckmode)
    {
      fits_close_file (BPET_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_PETFF: " "Error closing PET: %s\n",
		       BPET_file_path);
	}
    }

  fprintf (stdout, "aXe_FRINGECORR: Done...\n\n");
  exit (0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#ifndef __SOLVERS_H_
#define __SOLVERS_H_

#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_multifit_nlin.h>

const gsl_multiroot_fdfsolver_type * get_multiroot_fdfsolver_type(int index);

int
solve_multiroot(size_t dim, double * x, void * params, double atol,
		int multiroot_fdfsolver_type_idx, int itermax, int print_,
		int store_intermediate, double * intermediate, int * iter,
		int * nfev_, int * njev_, int * nfjev_);

void
print_multiroot_state (int iter, gsl_multiroot_fdfsolver * s, size_t dim);

int
solve_multifit(size_t ne, size_t nx, double * x, void * params,
	       double atol, double rtol, int multifit_fdfsolver_type_idx,
	       int itermax, int criterion, int print_, int store_intermediate,
	       double * intermediate, int * iter, int * nfev_, int * njev_, int * nfjev_);

void
print_multifit_state (int iter, gsl_multifit_fdfsolver * s, size_t nx, size_t ne);

void
c_func(size_t nx, double * x, double * params, double * out);

void
c_jac(size_t nx, double * x, double * params, double * out);

#endif

max_stars_repo_name: bjodah/symneqsys

max_issues_repo_name: bjodah/symneqsys

max_forks_repo_name: bjodah/symneqsys


text:
#include "gs2crmod_ext.h"
#include <math.h>
#include <string.h>
#include <gsl/gsl_spline.h>

VALUE gs2crmod_tensor_complexes_field_gsl_tensor_complex_2(VALUE self, VALUE options)
{
	VALUE field, field_complex, field_narray;
	VALUE field_complex_narray;
	VALUE cgsl_tensor_complex;
	VALUE ccomplex;
	VALUE shape;
	int *c_shape;
	int i, j, k;

	ccomplex = RGET_CLASS_TOP("Complex");
	field = RFCALL_11("field_gsl_tensor", options);
	cgsl_tensor_complex = RGET_CLASS(cgsl, "TensorComplex");
	shape = rb_funcall(field, rb_intern("shape"), 0);
	/*rb_p(shape);*/

	CR_INT_ARY_R2C_STACK(shape, c_shape);
	field_complex = rb_funcall3(cgsl_tensor_complex, rb_intern("alloc"), 3, RARRAY_PTR(shape));
	field_narray = RFCALL_10_ON(field, "narray");
	field_complex_narray = RFCALL_10_ON(field_complex, "narray");
	/*rb_p(field_complex);*/
	for (i=0;i<c_shape[0];i++)
	  for (j=0;j<c_shape[1];j++)
			for (k=0;k<c_shape[2];k++)
			{

				/*printf("%d %d %d\n", i, j, k);*/
				VALUE cmplex, re, im;
				re = rb_funcall(
								field_narray,
								rb_intern("[]"),
								4,
								INT2FIX(0),
								INT2FIX(k), INT2FIX(j), INT2FIX(i)
							);
				/*printf("%d %d %d re\n", i, j, k);*/
				im = rb_funcall(
								field_narray,
								rb_intern("[]"),
								4,
								INT2FIX(1),
								INT2FIX(k), INT2FIX(j), INT2FIX(i)
							);
				/*printf("%d %d %d im\n", i, j, k);*/
				cmplex =
					 	rb_funcall(
							ccomplex,
							rb_intern("rectangular"),
							2, re, im
						);
				/*printf("%d %d %d again\n", i, j, k);*/
				rb_funcall(
						field_complex_narray,
						rb_intern("[]="),
						4,
						INT2FIX(k), INT2FIX(j), INT2FIX(i),
						cmplex
					);

			}
	/*rb_p(field);*/

		/*field_complex = */
	return field_complex;
}

VALUE gs2crmod_tensor_field_gsl_tensor(VALUE self, VALUE options)
{
	VALUE cgsl_tensor;
	VALUE field, field_narray;
	VALUE field_real_space, field_real_space_narray;
	VALUE field_real_space_new, field_real_space_new_narray;
	VALUE shape, shape_real_space, shape_real_space_new;
	VALUE workspacex, workspacey;
	VALUE ccomplex, cgsl_complex;
	VALUE re, im, cmplex;
	VALUE *ary_ptr;
	VALUE range_x, range_y;
	int shape0;
	double* c_field, *c_field_real_space;
	int* c_shape, *c_shape_real_space;
	int * c_shape_real_space_new;
	int i,j,k,m;
	int kint,km,kold;
	double frac, interp_value;

	ccomplex = RGET_CLASS_TOP("Complex");
	cgsl_complex = RGET_CLASS(cgsl, "Complex");
	cgsl_tensor = RGET_CLASS(cgsl, "Tensor");
	/*cgsl_vector = RGET_CLASS(cgsl, "Vector");*/

	Check_Type(options, T_HASH);
	if(RTEST(CR_HKS(options, "field_netcdf"))){
		field = CR_HKS(options, "field_netcdf");
		CR_CHECK_CLASS(field, cgsl_tensor);
	}
	else {
	 field = RFCALL_11("field_gsl_tensor", options);
	}
	field_narray = RFCALL_10_ON(field, "narray");
	shape = RFCALL_10_ON(field, "shape");

	workspacex = RFCALL_11_ON(cgsl_vector_complex, "alloc", RARRAY_PTR(shape)[1]);
	shape0 = NUM2INT(RARRAY_PTR(shape)[0]);
	workspacey = RFCALL_11_ON(cgsl_vector, "alloc", INT2NUM(shape0 * 2 - 2 + shape0%2));

	field_real_space = rb_funcall(
			cgsl_tensor,
			rb_intern("alloc"),
			3,
			RFCALL_10_ON(workspacey, "size"),
			RARRAY_PTR(shape)[1],
			RARRAY_PTR(shape)[2]
			);
	field_real_space_narray = RFCALL_10_ON(field_real_space, "narray");
	shape_real_space = RFCALL_10_ON(field_real_space, "shape");

	CR_INT_ARY_R2C_STACK(shape, c_shape);
	CR_INT_ARY_R2C_STACK(RFCALL_10_ON(field_real_space, "shape"), c_shape_real_space);

	c_field = ALLOC_N(double, c_shape[0] * c_shape[1] * c_shape[2] * c_shape[3]);

  /*printf("Allocated stuff\n");*/
	for (j=0; j<c_shape[2]; j++) /*theta loop*/
	{
		for (i=0; i<c_shape[0]; i++) /*ky loop*/
		{
			for (k=0; k<c_shape[1]; k++) /*kx loop*/
			{
				/*First copy the field onto the
				 * x workspace, Fourier transform
				 * the workspace and copy it onto
				 * the c array*/
				re = rb_funcall(
								field_narray,
								rb_intern("[]"),
								4,
								INT2FIX(0),
								INT2FIX(j), INT2FIX(k), INT2FIX(i)
							);
				/*printf("%d %d %d re\n", i, j, k);*/
				im = rb_funcall(
								field_narray,
								rb_intern("[]"),
								4,
								INT2FIX(1),
								INT2FIX(j), INT2FIX(k), INT2FIX(i)
							);
				/*printf("%d %d %d im\n", i, j, k);*/
				cmplex =
					 	rb_funcall(
							cgsl_complex,
							rb_intern("alloc"),
							2, re, im
						);
				/*printf("%d %d %d again\n", i, j, k);*/
				/*printf("Made complex\n");*/
				rb_funcall(
					workspacex,
					rb_intern("[]="),
					2,
					INT2FIX(k),
					cmplex);
				/*printf("Added complex\n");*/

			}
			/*printf("Made complex vector\n");*/
			workspacex = RFCALL_10_ON(workspacex, "backward");
			/*printf("Done x transform\n");*/
			for (k=0;k<c_shape[1];k++){
				cmplex = RFCALL_11_ON(workspacex,"[]",INT2FIX(k));
				c_field[
					j*c_shape[0]*c_shape[1]*2 +
					i*c_shape[1]*2 +
					k*2
					] = NUM2DBL(RFCALL_10_ON(cmplex,"real"));
				c_field[
					j*c_shape[0]*c_shape[1]*2 +
					i*c_shape[1]*2 +
					k*2 + 1
					] = NUM2DBL(RFCALL_10_ON(cmplex,"imag"));

			}
		}
		/*printf("Done x\n");*/
		/* Now copy onto the y workspace,
		 * Fourier transform and copy onto
		 * the second c array*/
		for (k=0;k<c_shape[1];k++)
		{
		  m=0;
		  for (i=0;i<c_shape[0];i++)
			{
				rb_funcall(
						workspacey,
						rb_intern("[]="),
						2,
						INT2FIX(m),
						rb_float_new(c_field[
								j*c_shape[0]*c_shape[1]*2 +
								i*c_shape[1]*2 +
								k*2
								])
						);
				m++;
				/* We are converting a complex array
				 * to a half complex array in
				 * preparation for transforming
				 * it to real space, and so there
				 * are one or two elements we leave
				 * unfilled.*/
				if (i==0 || (c_shape[0]%2==0 && i == c_shape[0]/2 + 1)) continue;
				rb_funcall(
						workspacey,
						rb_intern("[]="),
						2,
						INT2FIX(m),
						rb_float_new(c_field[
								j*c_shape[0]*c_shape[1]*2 +
								i*c_shape[1]*2 +
								k*2 + 1
								])
						);
				m++;
			}
			/*printf("Made y vector\n");*/
			/*rb_p(workspacey);*/
			/*rb_p(RFCALL_10_ON(workspacey, "class"));*/
				workspacey = RFCALL_10_ON(workspacey, "backward");
				/*printf("Done y transform\n");*/
			for (i=0;i<FIX2INT(RFCALL_10_ON(workspacey, "size"));i++)
			{
				rb_funcall(
						field_real_space_narray,
						rb_intern("[]="), 4,
						INT2FIX(j), INT2FIX(k), INT2FIX(i),
						RFCALL_11_ON(workspacey, "[]", INT2FIX(i))
						);
			}
		}
	}
	/*printf("HERE!\n");*/
  range_x = CR_RANGE_INC(
				(RTEST(CR_HKS(options, "ymin")) ?
					FIX2INT(CR_HKS(options, "ymin")) :
			 		0),
				(RTEST(CR_HKS(options, "ymax")) ?
					FIX2INT(CR_HKS(options, "ymax")) :
					c_shape_real_space[0]-1)
				);
	/*rb_p(range_x);*/
	range_y = CR_RANGE_INC(
				(RTEST(CR_HKS(options, "xmin")) ?
					FIX2INT(CR_HKS(options, "xmin")) :
			 		0),
				(RTEST(CR_HKS(options, "xmax")) ?
					FIX2INT(CR_HKS(options, "xmax")) :
					c_shape_real_space[1]-1)
				);
	/*printf("Made ranges\n");*/
	/*rb_p(rb_funcall(field_real_space, rb_intern("[]"), 3, INT2FIX(0), INT2FIX(2), INT2FIX(3)));*/
	field_real_space = rb_funcall(
			field_real_space,
			rb_intern("[]"), 3,
			range_x,
			range_y,
			Qtrue);
	/*rb_p(rb_funcall(field_real_space, rb_intern("[]"), 3, INT2FIX(0), INT2FIX(2), INT2FIX(3)));*/
	/*printf("SHORTENED!\n");*/
	if (RTEST(CR_HKS(options, "interpolate_theta")))
	{
		shape_real_space_new = RFCALL_10_ON(shape_real_space, "dup");
		kint = NUM2INT(CR_HKS(options, "interpolate_theta"));
		rb_funcall(
				shape_real_space_new,
				rb_intern("[]="), 2, INT2FIX(-1),
				INT2FIX((c_shape_real_space[2] - 1)*kint+1)
				);
		CR_INT_ARY_R2C_STACK(shape_real_space_new, c_shape_real_space_new);
		ary_ptr = RARRAY_PTR(shape_real_space_new);
		field_real_space_new = rb_funcall(
				cgsl_tensor,
				rb_intern("float"),
				3,
				ary_ptr[0], ary_ptr[1], ary_ptr[2]);
		field_real_space_new_narray = RFCALL_10_ON(field_real_space_new, "narray");
    /*printf("Allocated");*/
		/*rb_p(shape_real_space_new);*/
		for (i=0;i<c_shape_real_space_new[0];i++)
		for (j=0;j<c_shape_real_space_new[1];j++)
		{
			/*printf("i %d j %d k %d\n", i, j, c_shape_real_space_new[2]-1); */
			rb_funcall(
					field_real_space_new_narray,
					rb_intern("[]="),
					4,INT2FIX(c_shape_real_space_new[2]-1),
					INT2FIX(j),INT2FIX(i),
					rb_funcall(
						field_real_space_narray,
						rb_intern("[]"), 3,
						INT2FIX(c_shape_real_space[2]-1),
						INT2FIX(j),INT2FIX(i)
						)
					);
			for (k=0;k<c_shape_real_space_new[2]-1;k++)
			{
				/*printf("i %d j %d k %d\n", i, j, k); */
				km = k%kint;
				frac = (double)km/(double)kint;
				kold = (k-km)/kint;
				interp_value = NUM2DBL(
						rb_funcall(
							field_real_space_narray,
							rb_intern("[]"),3,
							INT2FIX(kold),INT2FIX(j),INT2FIX(i)
							)
						)*(1.0-frac) + NUM2DBL(
						rb_funcall(
							field_real_space_narray,
							rb_intern("[]"),3,
							INT2FIX(kold+1),INT2FIX(j),INT2FIX(i)
							)
						) * frac;
				rb_funcall(
						field_real_space_new_narray,
						rb_intern("[]="), 4,
						INT2FIX(k),INT2FIX(j),INT2FIX(i),
						rb_float_new(interp_value)
						);
				/*if (i==0 && j==2 && k==3){*/
					/*printf("frac %f\n", frac);*/
				/*}*/
			}
		}

		field_real_space = field_real_space_new;
	}
  /*printf("Finished field_real_space_gsl_tensor\n");*/
	/*rb_p(shape_real_space_new);*/
	return field_real_space;
}

void Init_gs2crmod_ext()
{
	/*printf("HERE!!!");*/
	ccode_runner =  RGET_CLASS_TOP("CodeRunner");
	ccode_runner_gs2 =  rb_define_class_under(ccode_runner, "Gs2",
			RGET_CLASS(
				RGET_CLASS(ccode_runner, "Run"),
				"FortranNamelist"
				)
			);

	ccode_runner_gs2_gsl_tensor_complexes = rb_define_module_under(ccode_runner_gs2, "GSLComplexTensors");
	rb_include_module(ccode_runner_gs2, ccode_runner_gs2_gsl_tensor_complexes);

	ccode_runner_gs2_gsl_tensors = rb_define_module_under(ccode_runner_gs2, "GSLTensors");
	rb_include_module(ccode_runner_gs2, ccode_runner_gs2_gsl_tensors);

	cgsl = RGET_CLASS_TOP("GSL");
	cgsl_vector = RGET_CLASS(cgsl, "Vector");
	cgsl_vector_complex = RGET_CLASS(cgsl_vector, "Complex");

	rb_define_method(ccode_runner_gs2_gsl_tensor_complexes, "field_gsl_tensor_complex_2", gs2crmod_tensor_complexes_field_gsl_tensor_complex_2, 1);
	rb_define_method(ccode_runner_gs2_gsl_tensors, "field_real_space_gsl_tensor", gs2crmod_tensor_field_gsl_tensor, 1);

	/*rb_define_method(ccode_runner_ext, "hello_world", code_runner_ext_hello_world, 0);*/
}


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>

#include <math.h>
#include <options/options.h>

#include "../csm/csm_all.h"

struct {
	double interval;
	int seed;
}p ;

LDP ld_resample(LDP ld);

gsl_rng * rng;

int main(int argc, const char ** argv) {
	sm_set_program_name(argv[0]);

	struct option* ops = options_allocate(3);
	options_double(ops, "interval", &p.interval, sqrt(2.0), " 1 = no resampling");

	if(!options_parse_args(ops, argc, argv)) {
		options_print_help(ops, stderr);
		return -1;
	}

	int count = -1;
	LDP ld;
	while( (ld = ld_read_smart(stdin))) {
		count++;

		if(!ld_valid_fields(ld))  {
			sm_error("Invalid laser data (#%d in file)\n", count);
			continue;
		}

/*		if(count & 1) {*/
			LDP ld2 = ld_resample(ld);
			ld_write_as_json(ld2, stdout);
			ld_free(ld2);
			ld_free(ld);
/*		} else {
			ld_write_as_json(ld, stdout);
			ld_free(ld);
		}*/

		count++;
	}

	return 0;
}

LDP ld_resample(LDP ld) {
	/* FIXME magic number */
	int n = (int) (floor(ld->nrays / p.interval));

	LDP ld2 = ld_alloc_new(n);
	int k;
	for(k=0;k<n;k++) {
		double index = k * p.interval;
		int i = (int) floor(index);
		int j = i + 1;
		double a = 1 - (index - i);

		if(  (j>= ld->nrays) || !ld->valid[i] || !ld->valid[j]) {
			ld2->valid[k] = 0;
			ld2->readings[k] = NAN;
			ld2->alpha_valid[k] = 0;
			ld2->alpha[k] = NAN;
			ld2->theta[k] = ld->theta[i];
		} else {
			ld2->theta[k] = a * ld->theta[i] + (1-a) * ld->theta[j];

			if(is_nan(ld2->theta[k])) {
				sm_debug("Hey, k=%d theta[%d]=%f theta[%d]=%f\n",
					k,i,ld->theta[i],j,ld->theta[j]);
			}

			ld2->readings[k] = a * ld->readings[i] + (1-a) * ld->readings[j];
			ld2->valid[k] = 1;
		}

/*		sm_debug("k=%d index=%f i=%d a=%f valid %d reading %f\n", k,index,i,a,ld2->valid[k],ld2->readings[k]);*/
	}

	ld2->min_theta = ld2->theta[0];
	ld2->max_theta = ld2->theta[n-1];
	ld2->tv = ld->tv;

	copy_d(ld->odometry, 3, ld2->odometry);
	copy_d(ld->estimate, 3, ld2->estimate);

	return ld2;
}

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


text:
/* SO3blas.c */

// SO(3) Lie group implemented with GNU GSL
// cf. wikipedia, "Rotation_group_SO(3)#Infinitesimal_rotations"

// EY note : 20160519 Bizarrely, matrix multiplication isn't implemented
// in a straightforward manner with the gsl_matrix class;
// see blas (which I don't know what it is yet) and dgemm for the libraries

/*
  Compiling advice:
  gcc matrixio.c -lgsl -lgslcblas

 */

#include <stdio.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

// Function prototype
/* int commute(gsl_matrix_view, gsl_matrix_view); */

void commute1(void);

// global variables

double l_x[] = {  0,  0,  0,
		    0,  0, -1,
		    0,  1,  0};

double l_y[] = {  0,  0,  1,
		    0,  0,  0,
		   -1,  0, 0};

double l_z[] = {  0, -1,  0,
		    1,  0,  0,
		    0,  0,  0};
// cannot put gsl_matrix view commands in global; I need to check why:
// error: initializer element is not constant

int main(void)
{
  gsl_matrix_view L_x = gsl_matrix_view_array(l_x,3,3);
  gsl_matrix_view L_y = gsl_matrix_view_array(l_y,3,3);
  gsl_matrix_view L_z = gsl_matrix_view_array(l_z,3,3);

  printf("[ %g, %g, %g]\n", l_x[0], l_x[1], l_x[2] );
  printf("[ %g, %g, %g]\n", l_x[3], l_x[4], l_x[5] );
  printf("[ %g, %g, %g]\n\n", l_x[6], l_x[7], l_x[8] );

  printf("[ %g, %g, %g]\n", l_y[0], l_y[1], l_y[2] );
  printf("[ %g, %g, %g]\n", l_y[3], l_y[4], l_y[5] );
  printf("[ %g, %g, %g]\n\n", l_y[6], l_y[7], l_y[8] );

  printf("[ %g, %g, %g]\n", l_z[0], l_z[1], l_z[2] );
  printf("[ %g, %g, %g]\n", l_z[3], l_z[4], l_z[5] );
  printf("[ %g, %g, %g]\n\n", l_z[6], l_z[7], l_z[8] );

  double xy[] = { 0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00 };

  gsl_matrix_view XY = gsl_matrix_view_array(xy,3,3);

  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,
		 &L_x.matrix,&L_z.matrix,
		 0.0, &XY.matrix);

  printf("[ %g, %g, %g ]\n", xy[0], xy[1], xy[2]);
  printf("[ %g, %g, %g ]\n", xy[3], xy[4], xy[5]);
  printf("[ %g, %g, %g ]\n\n", xy[6], xy[7], xy[8]);

  commute1();

  return 0;
}

void commute1(void) {

  gsl_matrix_view L_x = gsl_matrix_view_array(l_x,3,3);
  gsl_matrix_view L_y = gsl_matrix_view_array(l_y,3,3);
  gsl_matrix_view L_z = gsl_matrix_view_array(l_z,3,3);

  double xy[] = { 0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00 };

  gsl_matrix_view XY = gsl_matrix_view_array(xy,3,3);

  double yx[] = { 0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00 };

  gsl_matrix_view YX = gsl_matrix_view_array(yx,3,3);

  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,
		 &L_x.matrix,&L_y.matrix,
		 0.0, &XY.matrix);

  gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,
		 &L_y.matrix,&L_x.matrix,
		 0.0, &YX.matrix);

   printf("[ %g, %g, %g ]\n", yx[0], yx[1], yx[2]);
   printf("[ %g, %g, %g ]\n", yx[3], yx[4], yx[5]);
   printf("[ %g, %g, %g ]\n", yx[6], yx[7], yx[8]);

}

/*

int commute(gsl_matrix_view X, gsl_matrix_view Y)
{
  double xy[] = { 0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00 };

  double yx[] = { 0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00,
		  0.00, 0.00, 0.00 };

  double bracket[] = { 0.00, 0.00, 0.00,
		       0.00, 0.00, 0.00,
		       0.00, 0.00, 0.00 };

  gsl_matrix_view XY = gsl_matrix_view_array(xy,3,3);
  gsl_matrix_view YX = gsl_matrix_view_array(yx,3,3);
  gsl_matrix_view Bracket = gsl_matrix_view_array(bracket,3,3);

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,1.0,
		 &X.matrix,&Y.matrix,
		 0.0, &XY.matrix);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,1.0,
		 &X.matrix,&Y.matrix,
		 0.0, &XY.matrix);
  gls_blas_dgemm(CblasNoTrans, CblasNoTrans,1.0,
		 &XY.matrix, &YX.matrix,
		 0.0, &Bracket.matrix);

  return 1;
}

*/

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


text:
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "cconfigspace_internal.h"
#include "distribution_internal.h"

struct _ccs_distribution_uniform_data_s {
	_ccs_distribution_common_data_t common_data;
	ccs_numeric_t                   lower;
	ccs_numeric_t                   upper;
	ccs_scale_type_t                scale_type;
	ccs_numeric_t                   quantization;
	ccs_numeric_type_t              internal_type;
	ccs_numeric_t                   internal_lower;
        ccs_numeric_t                   internal_upper;
	int                             quantize;
};
typedef struct _ccs_distribution_uniform_data_s _ccs_distribution_uniform_data_t;

static ccs_result_t
_ccs_distribution_del(ccs_object_t o) {
	(void)o;
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_uniform_get_bounds(_ccs_distribution_data_t *data,
                                     ccs_interval_t           *interval_ret);

static ccs_result_t
_ccs_distribution_uniform_samples(_ccs_distribution_data_t *data,
                                  ccs_rng_t                 rng,
                                  size_t                    num_values,
                                  ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_uniform_strided_samples(_ccs_distribution_data_t *data,
                                          ccs_rng_t                 rng,
                                          size_t                    num_values,
                                          size_t                    stride,
                                          ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_uniform_soa_samples(_ccs_distribution_data_t  *data,
                                      ccs_rng_t                  rng,
                                      size_t                     num_values,
                                      ccs_numeric_t            **values);

static _ccs_distribution_ops_t _ccs_distribution_uniform_ops = {
	{ &_ccs_distribution_del },
	&_ccs_distribution_uniform_samples,
	&_ccs_distribution_uniform_get_bounds,
	&_ccs_distribution_uniform_strided_samples,
	&_ccs_distribution_uniform_soa_samples
 };

static ccs_result_t
_ccs_distribution_uniform_get_bounds(_ccs_distribution_data_t *data,
                                     ccs_interval_t           *interval_ret) {
	_ccs_distribution_uniform_data_t *d = (_ccs_distribution_uniform_data_t *)data;
        ccs_numeric_t          l;
	ccs_bool_t             li;
	ccs_numeric_t          u;
	ccs_bool_t             ui;

	l = d->lower;
	li = CCS_TRUE;
	u = d->upper;
        ui = CCS_FALSE;
	interval_ret->type = d->common_data.data_types[0];
	interval_ret->lower = l;
	interval_ret->upper = u;
	interval_ret->lower_included = li;
	interval_ret->upper_included = ui;
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_uniform_strided_samples(_ccs_distribution_data_t *data,
                                          ccs_rng_t                 rng,
                                          size_t                    num_values,
                                          size_t                    stride,
                                          ccs_numeric_t            *values) {
	_ccs_distribution_uniform_data_t *d = (_ccs_distribution_uniform_data_t *)data;
	size_t i;
	const ccs_numeric_type_t  data_type      = d->common_data.data_types[0];
	const ccs_scale_type_t    scale_type     = d->scale_type;
	const ccs_numeric_t       quantization   = d->quantization;
	const ccs_numeric_t       lower          = d->lower;
	const ccs_numeric_t       internal_lower = d->internal_lower;
	const ccs_numeric_t       internal_upper = d->internal_upper;
	const int                 quantize       = d->quantize;
	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));

	if (data_type == CCS_NUM_FLOAT) {
		for (i = 0; i < num_values; i++) {
			values[i*stride].f = gsl_ran_flat(grng, internal_lower.f, internal_upper.f);
		}
		if (scale_type == CCS_LOGARITHMIC) {
			for (i = 0; i < num_values; i++)
				values[i*stride].f = exp(values[i*stride].f);
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i*stride].f = floor((values[i*stride].f - lower.f)/quantization.f) * quantization.f + lower.f;
		} else
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i*stride].f = floor(values[i*stride].f) * quantization.f + lower.f;
			else
				for (i = 0; i < num_values; i++)
					values[i*stride].f += lower.f;
	} else {
		if (scale_type == CCS_LOGARITHMIC) {
			for (i = 0; i < num_values; i++) {
				values[i*stride].i = floor(exp(gsl_ran_flat(grng, internal_lower.f, internal_upper.f)));
			}
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i*stride].i = ((values[i*stride].i - lower.i)/quantization.i) * quantization.i + lower.i;
		} else {
			for (i = 0; i < num_values; i++) {
				values[i*stride].i = gsl_rng_uniform_int(grng, internal_upper.i);
			}
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i*stride].i = values[i*stride].i * quantization.i + lower.i;
			else
				for (i = 0; i < num_values; i++)
					values[i*stride].i += lower.i;
		}
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_uniform_samples(_ccs_distribution_data_t *data,
                                  ccs_rng_t                 rng,
                                  size_t                    num_values,
                                  ccs_numeric_t            *values) {
	_ccs_distribution_uniform_data_t *d = (_ccs_distribution_uniform_data_t *)data;
	size_t i;
	const ccs_numeric_type_t  data_type      = d->common_data.data_types[0];
	const ccs_scale_type_t    scale_type     = d->scale_type;
	const ccs_numeric_t       quantization   = d->quantization;
	const ccs_numeric_t       lower          = d->lower;
	const ccs_numeric_t       internal_lower = d->internal_lower;
	const ccs_numeric_t       internal_upper = d->internal_upper;
	const int                 quantize       = d->quantize;
	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));

	if (data_type == CCS_NUM_FLOAT) {
		for (i = 0; i < num_values; i++) {
			values[i].f = gsl_ran_flat(grng, internal_lower.f, internal_upper.f);
		}
		if (scale_type == CCS_LOGARITHMIC) {
			for (i = 0; i < num_values; i++)
				values[i].f = exp(values[i].f);
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i].f = floor((values[i].f - lower.f)/quantization.f) * quantization.f + lower.f;
		} else
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i].f = floor(values[i].f) * quantization.f + lower.f;
			else
				for (i = 0; i < num_values; i++)
					values[i].f += lower.f;
	} else {
		if (scale_type == CCS_LOGARITHMIC) {
			for (i = 0; i < num_values; i++) {
				values[i].i = floor(exp(gsl_ran_flat(grng, internal_lower.f, internal_upper.f)));
			}
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i].i = ((values[i].i - lower.i)/quantization.i) * quantization.i + lower.i;
		} else {
			for (i = 0; i < num_values; i++) {
				values[i].i = gsl_rng_uniform_int(grng, internal_upper.i);
			}
			if (quantize)
				for (i = 0; i < num_values; i++)
					values[i].i = values[i].i * quantization.i + lower.i;
			else
				for (i = 0; i < num_values; i++)
					values[i].i += lower.i;
		}
	}
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_uniform_soa_samples(_ccs_distribution_data_t  *data,
                                      ccs_rng_t                  rng,
                                      size_t                     num_values,
                                      ccs_numeric_t            **values) {
	if (*values)
		return _ccs_distribution_uniform_samples(data, rng, num_values, *values);
	return CCS_SUCCESS;
}

ccs_result_t
ccs_create_uniform_distribution(ccs_numeric_type_t  data_type,
                                ccs_numeric_t       lower,
                                ccs_numeric_t       upper,
                                ccs_scale_type_t    scale_type,
                                ccs_numeric_t       quantization,
                                ccs_distribution_t *distribution_ret) {
	CCS_CHECK_PTR(distribution_ret);
	if (data_type != CCS_NUM_FLOAT && data_type != CCS_NUM_INTEGER)
		return -CCS_INVALID_TYPE;
	if (scale_type != CCS_LINEAR && scale_type != CCS_LOGARITHMIC)
		return -CCS_INVALID_SCALE;
	if (data_type == CCS_NUM_INTEGER && (
		lower.i >= upper.i ||
		(scale_type == CCS_LOGARITHMIC && lower.i <= 0) ||
		quantization.i < 0 ||
		quantization.i > upper.i - lower.i ) )
		return -CCS_INVALID_VALUE;
	if (data_type == CCS_NUM_FLOAT && (
		lower.f >= upper.f ||
		(scale_type == CCS_LOGARITHMIC && lower.f <= 0.0) ||
		quantization.f < 0.0 ||
		quantization.f > upper.f - lower.f ) )
		return -CCS_INVALID_VALUE;
	uintptr_t mem = (uintptr_t)calloc(1, sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_uniform_data_t) + sizeof(ccs_numeric_type_t));

	if (!mem)
		return -CCS_OUT_OF_MEMORY;
	ccs_distribution_t distrib = (ccs_distribution_t)mem;
	_ccs_object_init(&(distrib->obj), CCS_DISTRIBUTION, (_ccs_object_ops_t *)&_ccs_distribution_uniform_ops);
        _ccs_distribution_uniform_data_t * distrib_data = (_ccs_distribution_uniform_data_t *)(mem + sizeof(struct _ccs_distribution_s));
	distrib_data->common_data.data_types    = (ccs_numeric_type_t *)(mem + sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_uniform_data_t));
	distrib_data->common_data.type          = CCS_UNIFORM;
	distrib_data->common_data.dimension     = 1;
	distrib_data->common_data.data_types[0] = data_type;
	distrib_data->scale_type                = scale_type;
	distrib_data->quantization              = quantization;
	distrib_data->lower                     = lower;
	distrib_data->upper                     = upper;

	if (data_type == CCS_NUM_FLOAT) {
		if (quantization.f != 0.0)
			distrib_data->quantize = 1;
		if (scale_type == CCS_LOGARITHMIC) {
			distrib_data->internal_lower.f = log(lower.f);
			distrib_data->internal_upper.f = log(upper.f);
		} else {
			distrib_data->internal_lower.f = 0.0;
			distrib_data->internal_upper.f = upper.f - lower.f;
			if (distrib_data->quantize)
				distrib_data->internal_upper.f /= quantization.f;
		}
	} else {
		if (quantization.i != 0)
			distrib_data->quantize = 1;
		if (scale_type == CCS_LOGARITHMIC) {
			distrib_data->internal_lower.f = log(lower.i);
			distrib_data->internal_upper.f = log(upper.i);
		} else {
			distrib_data->internal_lower.i = 0;
			distrib_data->internal_upper.i = upper.i - lower.i;
			if (quantization.i != 0)
				distrib_data->internal_upper.i /= quantization.i;
		}
	}
	distrib->data = (_ccs_distribution_data_t *)distrib_data;
	*distribution_ret = distrib;
	return CCS_SUCCESS;
}

ccs_result_t
ccs_uniform_distribution_get_parameters(ccs_distribution_t  distribution,
                                        ccs_numeric_t      *lower_ret,
                                        ccs_numeric_t      *upper_ret,
                                        ccs_scale_type_t   *scale_type_ret,
                                        ccs_numeric_t      *quantization_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_UNIFORM);
	if (!lower_ret && !upper_ret && !scale_type_ret && !quantization_ret)
		return -CCS_INVALID_VALUE;
	_ccs_distribution_uniform_data_t * data = (_ccs_distribution_uniform_data_t *)distribution->data;

	if (lower_ret)
		*lower_ret = data->lower;
	if (upper_ret)
		*upper_ret = data->upper;
	if (scale_type_ret)
		*scale_type_ret = data->scale_type;
	if (quantization_ret)
		*quantization_ret = data->quantization;
	return CCS_SUCCESS;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#ifdef PEDANTIC_MEMORY_HANDLER

#define MAXBLOCKS 5000

#ifdef PEDANTIC_MEMORY_CEILING
static size_t TotBytes;
static void *Base;
#endif

static unsigned long Nblocks;
static double Highmark;

static void **Table;
static size_t *BlockSize;

void mymalloc_init(void)
{
#ifdef PEDANTIC_MEMORY_CEILING
  size_t n;
#endif

  BlockSize = (size_t *) malloc(MAXBLOCKS * sizeof(size_t));
  Table = (void **) malloc(MAXBLOCKS * sizeof(void *));

#ifdef PEDANTIC_MEMORY_CEILING
  n = PEDANTIC_MEMORY_CEILING * 1024.0 * 1024.0;

  if(!(Base = malloc(n)))
    {
      printf("Failed to allocate memory for `Base' (%d Mbytes).\n", (int) PEDANTIC_MEMORY_CEILING);
      endrun(122);
    }

  TotBytes = FreeBytes = n;
#endif

  AllocatedBytes = 0;
  Nblocks = 0;
  Highmark = 0;
}

void *mymalloc(size_t n)
{
  if((n % 8) > 0)
    n = (n / 8 + 1) * 8;

  if(n < 8)
    n = 8;

  if(Nblocks >= MAXBLOCKS)
    {
      printf("Task=%d: No blocks left in mymalloc().\n", ThisTask);
      endrun(813);
    }

#ifdef PEDANTIC_MEMORY_CEILING
  if(n > FreeBytes)
    {
      printf("Task=%d: Not enough memory in mymalloc(n=%g MB).  FreeBytes=%g MB\n",
	     ThisTask, n / (1024.0 * 1024.0), FreeBytes / (1024.0 * 1024.0));
      endrun(812);
    }
  Table[Nblocks] = Base + (TotBytes - FreeBytes);
  FreeBytes -= n;
#else
  Table[Nblocks] = malloc(n);
  if(!(Table[Nblocks]))
    {
      printf("failed to allocate %g MB of memory. (presently allocated=%g MB)\n",
	     n / (1024.0 * 1024.0), AllocatedBytes / (1024.0 * 1024.0));
      endrun(18);
    }
#endif

  AllocatedBytes += n;
  BlockSize[Nblocks] = n;

  Nblocks += 1;

  /*
     if(AllocatedBytes / (1024.0 * 1024.0) > Highmark)
     {
     Highmark = AllocatedBytes / (1024.0 * 1024.0);
     printf("Task=%d:   new highmark=%g MB\n",
     ThisTask, Highmark);
     fflush(stdout);
     }
   */

  return Table[Nblocks - 1];
}

void *mymalloc_msg(size_t n, char *message)
{
  if((n % 8) > 0)
    n = (n / 8 + 1) * 8;

  if(n < 8)
    n = 8;

  if(Nblocks >= MAXBLOCKS)
    {
      printf("Task=%d: No blocks left in mymalloc().\n", ThisTask);
      endrun(813);
    }

#ifdef PEDANTIC_MEMORY_CEILING
  if(n > FreeBytes)
    {
      printf
	("Task=%d: Not enough memory in for allocating (n=%g MB) for block='%s'.  FreeBytes=%g MB  AllocatedBytes=%g MB.\n",
	 ThisTask, n / (1024.0 * 1024.0), message, FreeBytes / (1024.0 * 1024.0),
	 AllocatedBytes / (1024.0 * 1024.0));
      endrun(812);
    }
  Table[Nblocks] = Base + (TotBytes - FreeBytes);
  FreeBytes -= n;
#else
  Table[Nblocks] = malloc(n);
  if(!(Table[Nblocks]))
    {
      printf("failed to allocate %g MB of memory for block='%s'. (presently allocated=%g MB)\n",
	     n / (1024.0 * 1024.0), message, AllocatedBytes / (1024.0 * 1024.0));
      endrun(18);
    }
#endif

  AllocatedBytes += n;
  BlockSize[Nblocks] = n;

  Nblocks += 1;

  return Table[Nblocks - 1];
}

void myfree(void *p)
{
  if(Nblocks == 0)
    endrun(76878);

  if(p != Table[Nblocks - 1])
    {
      printf("Task=%d: Wrong call of myfree() - not the last allocated block!\n", ThisTask);
      fflush(stdout);
      endrun(814);
    }

  Nblocks -= 1;
  AllocatedBytes -= BlockSize[Nblocks];
#ifdef PEDANTIC_MEMORY_CEILING
  FreeBytes += BlockSize[Nblocks];
#else
  free(p);
#endif
}

void myfree_msg(void *p, char *msg)
{
  if(Nblocks == 0)
    endrun(76878);

  if(p != Table[Nblocks - 1])
    {
      printf("Task=%d: Wrong call of myfree() - '%s' not the last allocated block!\n", ThisTask, msg);
      fflush(stdout);
      endrun(8141);
    }

  Nblocks -= 1;
  AllocatedBytes -= BlockSize[Nblocks];
#ifdef PEDANTIC_MEMORY_CEILING
  FreeBytes += BlockSize[Nblocks];
#else
  free(p);
#endif
}

void *myrealloc(void *p, size_t n)
{
  if((n % 8) > 0)
    n = (n / 8 + 1) * 8;

  if(n < 8)
    n = 8;

  if(Nblocks == 0)
    endrun(76879);

  if(p != Table[Nblocks - 1])
    {
      printf("Task=%d: Wrong call of myrealloc() - not the last allocated block!\n", ThisTask);
      fflush(stdout);
      endrun(815);
    }

  AllocatedBytes -= BlockSize[Nblocks - 1];
#ifdef PEDANTIC_MEMORY_CEILING
  FreeBytes += BlockSize[Nblocks - 1];
#endif

#ifdef PEDANTIC_MEMORY_CEILING
  if(n > FreeBytes)
    {
      printf("Task=%d: Not enough memory in myremalloc(n=%g MB). previous=%g FreeBytes=%g MB\n",
	     ThisTask, n / (1024.0 * 1024.0), BlockSize[Nblocks - 1] / (1024.0 * 1024.0),
	     FreeBytes / (1024.0 * 1024.0));
      endrun(812);
    }
  Table[Nblocks - 1] = Base + (TotBytes - FreeBytes);
  FreeBytes -= n;
#else
  Table[Nblocks - 1] = realloc(Table[Nblocks - 1], n);
  if(!(Table[Nblocks - 1]))
    {
      printf("failed to reallocate %g MB of memory. previous=%g FreeBytes=%g MB\n",
	     n / (1024.0 * 1024.0), BlockSize[Nblocks - 1] / (1024.0 * 1024.0),
	     FreeBytes / (1024.0 * 1024.0));
      endrun(18123);
    }
#endif

  AllocatedBytes += n;
  BlockSize[Nblocks - 1] = n;

  return Table[Nblocks - 1];
}

#else

void mymalloc_init(void)
{
  AllocatedBytes = 0;
}

void *mymalloc(size_t n)
{
  void *ptr;

  if(n < 8)
    n = 8;

  ptr = malloc(n);

  if(!(ptr))
    {
      printf("failed to allocate %g MB of memory.\n", n / (1024.0 * 1024.0));
      endrun(14);
    }

  return ptr;
}

void *mymalloc_msg(size_t n, char *message)
{
  void *ptr;

  if(n < 8)
    n = 8;

  ptr = malloc(n);

  if(!(ptr))
    {
      printf("failed to allocate %g MB of memory when trying to allocate '%s'.\n",
	     n / (1024.0 * 1024.0), message);
      endrun(14);
    }

  return ptr;
}

void *myrealloc(void *p, size_t n)
{
  void *ptr;

  ptr = realloc(p, n);

  if(!(ptr))
    {
      printf("failed to re-allocate %g MB of memory.\n", n / (1024.0 * 1024.0));
      endrun(15);
    }

  return ptr;
}

void myfree(void *p)
{
  free(p);
}

void myfree_msg(void *p, char *msg)
{
  free(p);
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/* fft/errs.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.
 */

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

#include <gsl/gsl_complex.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_dft_complex.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#include <gsl/gsl_test.h>
#include "compare.h"
#include "complex_internal.h"

#include "urand.c"

int verbose = 0;

size_t tests = 0;
size_t passed = 0;
size_t failed = 0;

int
main (int argc, char *argv[])
{
  int status, factor_sum;
  size_t i, start, end, n;
  double *complex_data, *complex_tmp;
  double rms, total;

  gsl_fft_complex_wavetable * cw;

  if (argc == 2)
    {
      start = strtol (argv[1], NULL, 0);
      end = start + 1;
    }
  else
    {
      start = 1 ;
      end = 1000 ;
    }

  for (n = start; n < end; n++)
    {

      complex_data = (double *) malloc (n * 2 * sizeof (double));
      complex_tmp = (double *) malloc (n * 2 * sizeof (double));

      cw = gsl_fft_complex_wavetable_alloc (n);
      status = gsl_fft_complex_init (n, cw);
      status = gsl_fft_complex_generate (n, cw);

      for (i = 0; i < n; i++)
        {
          REAL(complex_data,1,i) = urand();
          IMAG(complex_data,1,i) = urand();
        }

      memcpy (complex_tmp, complex_data, n * 2 * sizeof (double));
      gsl_fft_complex_forward (complex_data, 1, n, cw);
      gsl_fft_complex_inverse (complex_data, 1, n, cw);

      total = 0.0;
      for (i = 0; i < n; i++)
        {
          double dr = REAL(complex_data,1,i) - REAL(complex_tmp,1,i);
          double di = IMAG(complex_data,1,i) - IMAG(complex_tmp,1,i);
          total += dr * dr + di * di;
        }

      rms = sqrt (total / n);

      factor_sum = 0;
      for (i = 0; i < cw->nf; i++)
        {
          int j = cw->factor[i];
          factor_sum += j;
        }

      printf ("n = %d factor_sum = %d rms = %e\n", n, factor_sum, rms);

      free (complex_data);
      free (complex_tmp);

    }

  return 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:
/* deriv/test.c
 *
 * Copyright (C) 2000 David Morrison
 *
 * 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.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

double
f1 (double x, void *params)
{
  return exp (x);
}

double
df1 (double x, void *params)
{
  return exp (x);
}

double
f2 (double x, void *params)
{
  if (x >= 0.0)
    {
      return x * sqrt (x);
    }
  else
    {
      return 0.0;
    }
}

double
df2 (double x, void *params)
{
  if (x >= 0.0)
    {
      return 1.5 * sqrt (x);
    }
  else
    {
      return 0.0;
    }
}

double
f3 (double x, void *params)
{
  if (x != 0.0)
    {
      return sin (1 / x);
    }
  else
    {
      return 0.0;
    }
}

double
df3 (double x, void *params)
{
  if (x != 0.0)
    {
      return -cos (1 / x) / (x * x);
    }
  else
    {
      return 0.0;
    }
}

double
f4 (double x, void *params)
{
  return exp (-x * x);
}

double
df4 (double x, void *params)
{
  return -2.0 * x * exp (-x * x);
}

double
f5 (double x, void *params)
{
  return x * x;
}

double
df5 (double x, void *params)
{
  return 2.0 * x;
}

double
f6 (double x, void *params)
{
  return 1.0 / x;
}

double
df6 (double x, void *params)
{
  return -1.0 / (x * x);
}

typedef int (deriv_fn) (const gsl_function * f, double x, double h, double * res, double *abserr);

void
test (deriv_fn * deriv, gsl_function * f, gsl_function * df, double x,
      const char * desc)
{
  double result, abserr;
  double expected = GSL_FN_EVAL (df, x);
  (*deriv) (f, x, 1e-4, &result, &abserr);

  gsl_test_abs (result, expected, GSL_MIN(1e-4,fabs(expected)) + GSL_DBL_EPSILON, desc);

  if (abserr < fabs(result-expected))
    {
      gsl_test_factor (abserr, fabs(result-expected), 2, "%s error estimate", desc);
    }
  else if (result == expected || expected == 0.0)
    {
      gsl_test_abs (abserr, 0.0, 1e-6, "%s abserr", desc);
    }
  else
    {
      double d = fabs(result - expected);
      gsl_test_abs (abserr, fabs(result-expected), 1e6*d, "%s abserr", desc);
    }
}

int
main ()
{
  gsl_function F1, DF1, F2, DF2, F3, DF3, F4, DF4, F5, DF5, F6, DF6;

  gsl_ieee_env_setup ();

  F1.function = &f1;
  DF1.function = &df1;

  F2.function = &f2;
  DF2.function = &df2;

  F3.function = &f3;
  DF3.function = &df3;

  F4.function = &f4;
  DF4.function = &df4;

  F5.function = &f5;
  DF5.function = &df5;

  F6.function = &f6;
  DF6.function = &df6;

  test (&gsl_deriv_central, &F1, &DF1, 1.0, "exp(x), x=1, central deriv");
  test (&gsl_deriv_forward, &F1, &DF1, 1.0, "exp(x), x=1, forward deriv");
  test (&gsl_deriv_backward, &F1, &DF1, 1.0, "exp(x), x=1, backward deriv");

  test (&gsl_deriv_central, &F2, &DF2, 0.1, "x^(3/2), x=0.1, central deriv");
  test (&gsl_deriv_forward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, forward deriv");
  test (&gsl_deriv_backward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, backward deriv");

  test (&gsl_deriv_central, &F3, &DF3, 0.45, "sin(1/x), x=0.45, central deriv");
  test (&gsl_deriv_forward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, forward deriv");
  test (&gsl_deriv_backward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, backward deriv");

  test (&gsl_deriv_central, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, central deriv");
  test (&gsl_deriv_forward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, forward deriv");
  test (&gsl_deriv_backward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, backward deriv");

  test (&gsl_deriv_central, &F5, &DF5, 0.0, "x^2, x=0, central deriv");
  test (&gsl_deriv_forward, &F5, &DF5, 0.0, "x^2, x=0, forward deriv");
  test (&gsl_deriv_backward, &F5, &DF5, 0.0, "x^2, x=0, backward deriv");

  test (&gsl_deriv_central, &F6, &DF6, 10.0, "1/x, x=10, central deriv");
  test (&gsl_deriv_forward, &F6, &DF6, 10.0, "1/x, x=10, forward deriv");
  test (&gsl_deriv_backward, &F6, &DF6, 10.0, "1/x, x=10, backward deriv");

  exit (gsl_test_summary ());
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*! \file
  Compute 3D shape features and moments.

The following functions compute several shape features, including
central moments, center of gravity, and volume size.

\par Author:
Gabriele Lohmann, MPI-CBS
*/

/* From the Vista library: */
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <via/via.h>

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

/* From the standard C libaray: */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define SQR(x) ((x)*(x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

static double power(double a, int k)
{
  register double res;
  register int i;

  res = 1.0;
  for (i=0; i<k; i++)
    res *= a;

  return res;
}

/*!
\fn  void VolumeCentroid(Volume v,double mean[3])
\brief Compute center of gravity of a volume.
\param v  input volume.
\param *mean output array.
The center of gravity is returned in mean[0],mean[1],mean[2].
Sufficient memory must have been allocated for <mean> prior to the call.
*/
void VolumeCentroid(Volume v, double mean[3])
{
  double npixels;
  long c,i,ce;
  VTrack t;

  mean[0] = 0;
  mean[1] = 0;
  mean[2] = 0;
  npixels = 0;

  for (i=0; i<v->nbuckets; i++) {
    for (t = v->bucket[i].first; t != NULL; t = t->next)  {
      mean[0] += (double) (t->band * t->length);
      mean[1] += (double) (t->row * t->length);

      ce = t->col + t->length;
      for (c = t->col; c<ce; c++)
	mean[2] += (double) c;

      npixels += t->length;
    }
  }
  if (npixels > 0) {
    mean[0] /= npixels;
    mean[1] /= npixels;
    mean[2] /= npixels;
  }
}

/*!
\fn double VolumeMoment(Volume v,double mean[3],int m0,int m1,int m2)
\brief Compute central moments of a volume.
\param v  input volume
\param *mean input array containing
the center of gravity as mean[0],mean[1],mean[2]  (slice,row,column).
If mean is NULL, the center of gravity is taken to be (0,0,0).
\param m0  first index of moment
\param m1  second index of moment
\param m2  third index of moment
*/
double VolumeMoment(Volume v,double mean[3],long m0,long m1,long m2)
{
  double res=0;
  double a = 1.0,sum = 0;
  long i,ca,ce,ci;
  double b,r,c;
  VTrack t;
  double g0,g1,g2;

  g0 = mean[0];
  g1 = mean[1];
  g2 = mean[2];

  res = 0;
  for (i=0; i<v->nbuckets; i++) {
    for (t = v->bucket[i].first; t != NULL; t = t->next)  {
      b  = (double) t->band - g0;
      r  = (double) t->row  - g1;
      a  = power(b,m0) * power(r,m1);
      ca = t->col;
      ce = ca + t->length;
      sum = 0;
      for (ci=ca; ci<ce; ci++) {
	c = (double)(ci - g2);
	sum += power(c,m2);
      }
      res += sum * a;
    }
  }
  return res;
}

/*!
\fn  void VBinCentroid(VImage src,double mean[3])
\brief Compute central moments of a binary raster image.
\param src  input image (bit repn)
\param *mean output array.
The center of gravity is returned in mean[0],mean[1],mean[2].
Sufficient memory must have been allocated for <mean> prior to the call.
*/
void VBinCentroid(VImage src,double mean[3])
{
  double npixels;
  long b,r,c;

  mean[0] = 0;
  mean[1] = 0;
  mean[2] = 0;
  npixels = 0;
  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	if (VPixel(src,b,r,c,VBit) > 0) {
	  mean[0] += b;
	  mean[1] += r;
	  mean[2] += c;
	  npixels++;
	}
      }
    }
  }
  if (npixels > 0) {
    mean[0] /= npixels;
    mean[1] /= npixels;
    mean[2] /= npixels;
  }
}

/*!
\fn  double VBinMoment(VImage src,double mean[3],int m0, int m1, int m2)
\brief Compute central moments of a binary raster image.
\param src  input image (bit repn)
\param *mean input array containing
the center of gravity as mean[0],mean[1],mean[2]  (slice,row,column).
If mean is NULL, the center of gravity is taken to be (0,0,0);
\param m0  first index of moment
\param m1  second index of moment
\param m2  third index of moment
*/
double VBinMoment(VImage src, double mean[3], long m0, long m1, long m2)
{
  double res;
  long b,r,c;
  double g0,g1,g2;

  g0 = mean[0];
  g1 = mean[1];
  g2 = mean[2];

  res = 0;
  for (b=0; b<VImageNBands(src); b++) {
    for (r=0; r<VImageNRows(src); r++) {
      for (c=0; c<VImageNColumns(src); c++) {
	if (VPixel(src,b,r,c,VBit) > 0)
	  res += power((double)(b - g0),m0)
	    * power((double)(r - g1),m1)
	    * power((double)(c - g2),m2);
      }
    }
  }
  return res;
}

/*!
  \fn long VolumeSize(Volume v)
  \brief compute volume size, output the number of voxels of a single volume
  \param v   a single volume
*/
long VolumeSize(Volume v)
{
  long i,isize;
  VTrack t;

  isize = 0;
  for (i=0; i<VolumeNBuckets(v); i++) {
    for (t = VFirstTrack(v,i); VTrackExists(t); t = VNextTrack(t))
      isize += VTrackLength(t);
  }
  return isize;
}

/*!
  \fn long VBinSize(VImage src)
  \brief count number of foreground voxels.
  \param src  input image (bit repn)
*/
long VBinSize(VImage src)
{
  long n;
  long i;
  VBit *bin_pp;

  if (VPixelRepn(src) != VBitRepn) VError(" input image must be bit repn");

  n = 0;
  bin_pp = (VBit *) VPixelPtr(src,0,0,0);
  for (i=0; i<VImageNPixels(src); i++) {
    if (*bin_pp > 0) n++;
    bin_pp++;
  }
  return n;
}

/*!
  \fn float VolumeDir(Volume vol,float *e,float x[3])
  \brief computer principal direction of a volume from its interia matrix
  \param vol   input volume
  \param e     output largest eigenvalue
  \param x     output first eigenvector
*/
float VolumeDir(Volume vol,float *e,float x[3])
{
  gsl_matrix *a=NULL;
  static gsl_matrix *evec=NULL;
  static gsl_vector *eval=NULL;
  static gsl_eigen_symmv_workspace *workspace=NULL;
  double m020,m002,m200,m110,m101,m011;
  double norm,angle;
  double mean[3];
  float tiny=1.0e-5;

  mean[0] = mean[1] = mean[2] = 0;
  VolumeCentroid(vol,mean);

  m020 = VolumeMoment(vol,mean,0,2,0);
  m002 = VolumeMoment(vol,mean,0,0,2);
  m200 = VolumeMoment(vol,mean,2,0,0);

  m110 = VolumeMoment(vol,mean,1,1,0);
  m101 = VolumeMoment(vol,mean,1,0,1);
  m011 = VolumeMoment(vol,mean,0,1,1);

  /* inertia matrix */
  if (a == NULL) {
    a = gsl_matrix_calloc(3,3);
    workspace = gsl_eigen_symmv_alloc(3);
  }
  gsl_matrix_set(a,0,0,m020 + m002);
  gsl_matrix_set(a,1,1,m200 + m002);
  gsl_matrix_set(a,2,2,m200 + m020);

  gsl_matrix_set(a,0,1,-m110);
  gsl_matrix_set(a,1,0,-m110);

  gsl_matrix_set(a,0,2,-m101);
  gsl_matrix_set(a,2,0,-m101);

  gsl_matrix_set(a,1,2,-m011);
  gsl_matrix_set(a,2,1,-m011);

  gsl_eigen_symmv(a,eval,evec,workspace);
  gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_DESC);

  x[0] = gsl_matrix_get(evec,0,0);
  x[1] = gsl_matrix_get(evec,1,0);
  x[2] = gsl_matrix_get(evec,2,0);

  *e = gsl_vector_get(eval,0);

  angle = 0;
  norm = sqrt((double)(SQR(x[0]) + SQR(x[1]) + SQR(x[2])));
  if (norm > tiny) {
    angle = ABS(x[1])/ norm;
  }

  return angle;
}

/*!
  \fn void VolumeEigen(Volume vol,gsl_vector *eval,gsl_matrix *evec)
  \brief computer principal directions of a volume from its interia matrix
  \param vol    input volume
  \param eval   output eigenvalues
  \param evec   output matrix of eigenvectors (columns)
*/
void VolumeEigen(Volume vol,gsl_vector *eval,gsl_matrix *evec)
{
  static gsl_matrix *a=NULL;
  static gsl_eigen_symmv_workspace *workspace=NULL;
  double m020,m002,m200,m110,m101,m011;
  double mean[3];

  mean[0] = mean[1] = mean[2] = 0;
  VolumeCentroid(vol,mean);

  m020 = VolumeMoment(vol,mean,0,2,0);
  m002 = VolumeMoment(vol,mean,0,0,2);
  m200 = VolumeMoment(vol,mean,2,0,0);

  m110 = VolumeMoment(vol,mean,1,1,0);
  m101 = VolumeMoment(vol,mean,1,0,1);
  m011 = VolumeMoment(vol,mean,0,1,1);

  /* inertia matrix */
  if (a == NULL) {
    a = gsl_matrix_calloc(3,3);
    workspace = gsl_eigen_symmv_alloc(3);
  }

  gsl_matrix_set(a,0,0,m020 + m002);
  gsl_matrix_set(a,1,1,m200 + m002);
  gsl_matrix_set(a,2,2,m200 + m020);

  gsl_matrix_set(a,0,1,-m110);
  gsl_matrix_set(a,1,0,-m110);

  gsl_matrix_set(a,0,2,-m101);
  gsl_matrix_set(a,2,0,-m101);

  gsl_matrix_set(a,1,2,-m011);
  gsl_matrix_set(a,2,1,-m011);

  gsl_eigen_symmv(a,eval,evec,workspace);
  gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_DESC);

  return;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#pragma once

// v1.2 By Ook
//

#include <concepts>
#include <gsl/gsl>

namespace EnumConcepts
{
	template <class T>
	concept IsEnum = std::is_enum_v<T>;

	template<class _Ty>
	struct is_Numeric
		: std::integral_constant<bool,
		(std::is_arithmetic_v<_Ty> || std::is_same_v<_Ty, std::byte> || std::is_enum_v<_Ty>)
		&& !std::is_same_v<_Ty, wchar_t>
		&& !std::is_same_v<_Ty, char>
		&& !std::is_pointer_v<_Ty>
		//&& !std::is_enum_v<_Ty>
		>
	{
	};

	// determine whether T is a Number type (excluding char / wchar), same as is_Numeric<T>::value
	template <typename T>
	constexpr bool is_Numeric_v = is_Numeric<T>::value;

	template <class T>
	concept IsNumeric = is_Numeric_v<T>;
}

//template <EnumConcepts::IsEnum E>
//constexpr E& operator |=(E& eStyle, const E& dStyle) noexcept
//{
//	if constexpr (sizeof(E) <= sizeof(size_t))
//		eStyle = static_cast<E>(static_cast<size_t>(eStyle) | static_cast<size_t>(dStyle));
//	else
//		eStyle = static_cast<E>(static_cast<uint64_t>(eStyle) | static_cast<uint64_t>(dStyle));
//	return eStyle;
//}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E& operator |=(E& eStyle, const T& dStyle) noexcept
{
	if constexpr(sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
		eStyle = static_cast<E>(static_cast<size_t>(eStyle) | static_cast<size_t>(dStyle));
	else
		eStyle = static_cast<E>(static_cast<uint64_t>(eStyle) | static_cast<uint64_t>(dStyle));
	return eStyle;
}

//template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
//constexpr T& operator |=(T& dStyle, const E& eStyle) noexcept
//{
//	//if constexpr (sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
//	//	dStyle = static_cast<T>(static_cast<size_t>(eStyle) | static_cast<size_t>(dStyle));
//	//else
//	//	dStyle = static_cast<T>(static_cast<uint64_t>(eStyle) | static_cast<uint64_t>(dStyle));
//	//return dStyle;
//
//	dStyle |= static_cast<T>(eStyle);
//
//	return dStyle;
//}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E operator |(const E& eStyle, const T& dStyle) noexcept
{
	if constexpr (sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
		return static_cast<E>(static_cast<size_t>(eStyle) | static_cast<size_t>(dStyle));
	else
		return static_cast<E>(static_cast<uint64_t>(eStyle) | static_cast<uint64_t>(dStyle));
}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E& operator &=(E& eStyle, const T& dStyle) noexcept
{
	if constexpr (sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
		eStyle = static_cast<E>(static_cast<size_t>(eStyle) & static_cast<size_t>(dStyle));
	else
		eStyle = static_cast<E>(static_cast<uint64_t>(eStyle) & static_cast<uint64_t>(dStyle));
	return eStyle;
}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E operator &(const E& eStyle, const T& dStyle) noexcept
{
	if constexpr (sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
		return static_cast<E>(static_cast<size_t>(eStyle) & static_cast<size_t>(dStyle));
	else
		return static_cast<E>(static_cast<uint64_t>(eStyle) & static_cast<uint64_t>(dStyle));
}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E operator ^(const E& eStyle, const T& dStyle) noexcept
{
	if constexpr (sizeof(E) <= sizeof(size_t) && sizeof(T) <= sizeof(size_t))
		return static_cast<E>(static_cast<size_t>(eStyle) ^ static_cast<size_t>(dStyle));
	else
		return static_cast<E>(static_cast<uint64_t>(eStyle) ^ static_cast<uint64_t>(dStyle));
}

template <EnumConcepts::IsEnum E>
constexpr E operator ~(const E& eStyle) noexcept
{
	if constexpr (sizeof(E) <= sizeof(size_t))
		return static_cast<E>(~static_cast<size_t>(eStyle));
	else
		return static_cast<E>(~static_cast<uint64_t>(eStyle));
}

template <EnumConcepts::IsNumeric T, EnumConcepts::IsEnum E>
constexpr E to_Enum(T t) noexcept
{
	return static_cast<E>(t);
}

max_stars_repo_name: twig/dcxdll

max_issues_repo_name: twig/dcxdll

max_forks_repo_name: twig/dcxdll