text:
#ifndef __CNN_FC_H__
#define __CNN_FC_H__

#include <cblas.h>

#include "cnn_types.h"

#ifdef CNN_WITH_CUDA
#include <cublas_v2.h>
#include <cuda_runtime.h>

#include "cnn_init.h"
#endif

static inline void cnn_forward_fc(union CNN_LAYER* layerRef,
                                  struct CNN_CONFIG* cfgRef, int layerIndex)
{
#ifdef CNN_WITH_CUDA
    float alpha = 1.0;
    float beta = 0.0;

    struct CNN_LAYER_FC* layerPtr = &layerRef[layerIndex].fc;

    struct CNN_MAT* wPtr = &layerPtr->weight;
    struct CNN_MAT* outData = &layerPtr->outMat.data;
    struct CNN_MAT* preOutData = &layerRef[layerIndex - 1].outMat.data;

    // Weight matrix multiplication
    cnn_assert_cublas(cublasSgemm(                       //
        cnnInit.blasHandle, CUBLAS_OP_N, CUBLAS_OP_N,    //
        outData->cols, cfgRef->batch, preOutData->cols,  //
        &alpha,                                          //
        wPtr->mat, wPtr->cols,                           //
        preOutData->mat, preOutData->cols,               //
        &beta,                                           //
        outData->mat, outData->cols));

    // Add bias
    beta = 1.0;
    cnn_assert_cudnn(cudnnAddTensor(cnnInit.cudnnHandle,                    //
                                    &alpha,                                 //
                                    layerPtr->biasTen, layerPtr->bias.mat,  //
                                    &beta,                                  //
                                    layerPtr->dstTen, outData->mat));
#else
    // Weight matrix multiplication
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                layerRef[layerIndex - 1].outMat.data.rows,
                layerRef[layerIndex].outMat.data.cols,
                layerRef[layerIndex - 1].outMat.data.cols, 1.0,
                layerRef[layerIndex - 1].outMat.data.mat,
                layerRef[layerIndex - 1].outMat.data.cols,
                layerRef[layerIndex].fc.weight.mat,
                layerRef[layerIndex].fc.weight.cols, 0.0,
                layerRef[layerIndex].outMat.data.mat,
                layerRef[layerIndex].outMat.data.cols);

    // Add bias
    for (int j = 0; j < cfgRef->batch; j++)
    {
        int dstIndex = j * layerRef[layerIndex].outMat.data.cols;
        cblas_saxpy(layerRef[layerIndex].fc.bias.cols, 1.0,
                    layerRef[layerIndex].fc.bias.mat, 1,
                    &layerRef[layerIndex].outMat.data.mat[dstIndex], 1);
    }
#endif
}

static inline void cnn_backward_fc(union CNN_LAYER* layerRef,
                                   struct CNN_CONFIG* cfgRef, int layerIndex)
{
#ifdef CNN_WITH_CUDA
    float alpha = 1.0;
    float beta = 1.0;

    struct CNN_LAYER_FC* layerPtr = &layerRef[layerIndex].fc;

    struct CNN_MAT* wPtr = &layerPtr->weight;
    struct CNN_MAT* outData = &layerPtr->outMat.data;
    struct CNN_MAT* preOutData = &layerRef[layerIndex - 1].outMat.data;

    // Sum weight gradient matrix
    cnn_assert_cublas(                                             //
        cublasSgemm(cnnInit.blasHandle, CUBLAS_OP_N, CUBLAS_OP_T,  //
                    wPtr->cols, wPtr->rows, cfgRef->batch,         //
                    &alpha,                                        //
                    outData->grad, outData->cols,                  //
                    preOutData->mat, preOutData->cols,             //
                    &beta,                                         //
                    wPtr->grad, wPtr->cols));

    // Sum bias gradient matrix
    if (cfgRef->batch > 1)
    {
        cnn_assert_cudnn(
            cudnnReduceTensor(cnnInit.cudnnHandle, layerPtr->reduDesc,
                              layerPtr->indData, layerPtr->indSize,  //
                              cnnInit.wsData, cnnInit.wsSize,        //
                              &alpha,                                //
                              layerPtr->dstTen, outData->grad,       //
                              &beta,                                 //
                              layerPtr->biasTen, layerPtr->bias.grad));
    }
    else
    {
        cnn_assert_cublas(                                        //
            cublasSaxpy(cnnInit.blasHandle, layerPtr->bias.cols,  //
                        &alpha,                                   //
                        outData->grad, 1,                         //
                        layerPtr->bias.grad, 1));
    }
#else
    // Sum weight gradient matrix
    cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
                layerRef[layerIndex].fc.weight.rows,
                layerRef[layerIndex].fc.weight.cols, cfgRef->batch, 1.0,
                layerRef[layerIndex - 1].outMat.data.mat,
                layerRef[layerIndex - 1].outMat.data.cols,
                layerRef[layerIndex].outMat.data.grad,
                layerRef[layerIndex].outMat.data.cols, 1.0,
                layerRef[layerIndex].fc.weight.grad,
                layerRef[layerIndex].fc.weight.cols);

    // Sum bias gradient matrix
    for (int j = 0; j < cfgRef->batch; j++)
    {
        int srcShift = j * layerRef[layerIndex].outMat.data.cols;
        cblas_saxpy(layerRef[layerIndex].fc.bias.cols, 1.0,
                    &layerRef[layerIndex].outMat.data.grad[srcShift], 1,
                    layerRef[layerIndex].fc.bias.grad, 1);
    }
#endif

    // Find layer gradient
    if (layerIndex > 1)
    {
#ifdef CNN_WITH_CUDA
        beta = 0.0;
        cnn_assert_cublas(                                             //
            cublasSgemm(cnnInit.blasHandle, CUBLAS_OP_T, CUBLAS_OP_N,  //
                        wPtr->rows, cfgRef->batch, wPtr->cols,         //
                        &alpha,                                        //
                        wPtr->mat, wPtr->cols,                         //
                        outData->grad, outData->cols,                  //
                        &beta,                                         //
                        preOutData->grad, preOutData->cols));
#else
        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
                    layerRef[layerIndex - 1].outMat.data.rows,
                    layerRef[layerIndex - 1].outMat.data.cols,
                    layerRef[layerIndex].outMat.data.cols, 1.0,
                    layerRef[layerIndex].outMat.data.grad,
                    layerRef[layerIndex].outMat.data.cols,
                    layerRef[layerIndex].fc.weight.mat,
                    layerRef[layerIndex].fc.weight.cols, 0.0,
                    layerRef[layerIndex - 1].outMat.data.grad,
                    layerRef[layerIndex - 1].outMat.data.cols);
#endif
    }
}

#endif

ext: h

max_stars_repo_name: jamesljlster/cnn

max_issues_repo_name: jamesljlster/cnn

max_forks_repo_name: jamesljlster/cnn


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "DirectionalLight.h"

namespace Library
{
	class ProxyModel;
}

namespace Rendering
{
	class DiffuseLightingMaterial;

	class DiffuseLightingDemo final : public Library::DrawableGameComponent
	{
	public:
		DiffuseLightingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		DiffuseLightingDemo(const DiffuseLightingDemo&) = delete;
		DiffuseLightingDemo(DiffuseLightingDemo&&) = default;
		DiffuseLightingDemo& operator=(const DiffuseLightingDemo&) = default;
		DiffuseLightingDemo& operator=(DiffuseLightingDemo&&) = default;
		~DiffuseLightingDemo();

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);
		void ToggleAnimation();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float DirectionalLightIntensity() const;
		void SetDirectionalLightIntensity(float intensity);

		const DirectX::XMFLOAT3& LightDirection() const;
		void RotateDirectionalLight(DirectX::XMFLOAT2 amount);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<DiffuseLightingMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		Library::DirectionalLight mDirectionalLight;
		std::unique_ptr<Library::ProxyModel> mProxyModel;
		float mModelRotationAngle{ 0.0f };
		bool mAnimationEnabled{ false };
		bool mUpdateMaterial{ true };
	};
}


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

#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

/*
 * Spline the linear power spectrum for mu-Sigma MG cosmologies.
 * @param cosmo Cosmological parameters
 ^ @param psp The linear power spectrum to spline.
 * @param status, integer indicating the status
 */
void ccl_rescale_musigma_s8(ccl_cosmology* cosmo, ccl_f2d_t *psp,
                            int mg_rescale, int* status) {

  int do_mg = mg_rescale && (fabs(cosmo->params.mu_0) > 1e-14);
  int do_s8 = isfinite(cosmo->params.sigma8) && (!isfinite(cosmo->params.A_s));

  if (!do_mg && !do_s8)
    return;

  size_t na;
  double *aa;
  double rescale_extra_musig = 1;
  double *rescale_factor = NULL;

  if (*status == 0) {
    // Get array of scale factors
    if(psp->fa != NULL) {
      na = psp->fa->size;
      aa = psp->fa->x;
    }
    else if(psp->fka != NULL) {
      na = psp->fka->interp_object.ysize;
      aa = psp->fka->yarr;
    }
    else {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,
             "ccl_musigma.c: ccl_rescale_musigma_s8(): "
             "input pk2d has no splines.\n");
    }
  }

  // Alloc rescale_factor
  if(*status == 0) {
    rescale_factor = malloc(na * sizeof(double));
    if(rescale_factor == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
             "ccl_musigma.c: ccl_rescale_musigma_s8(): memory allocation\n");
    }
  }

  // set to 1 by default
  if(*status == 0) {
    for(int i=0; i<na; i++)
      rescale_factor[i]=1;
  }

  // If scale-independent mu / Sigma modified gravity is in use
  // and mu ! = 0 : get the unnormalized growth factor in MG and for
  // corresponding GR case, to rescale CLASS power spectrum
  if(do_mg) {
    // Set up another cosmology which is exactly the same as the
    // current one but with mu_0 and Sigma_0=0, for scaling P(k)

    // Get a list of the three neutrino masses already calculated
    double norm_pk;
    double mnu_list[3] = {0, 0, 0};
    ccl_parameters params_GR;
    params_GR.m_nu = NULL;
    params_GR.z_mgrowth = NULL;
    params_GR.df_mgrowth = NULL;
    ccl_cosmology * cosmo_GR = NULL;

    for (int i=0; i< cosmo->params.N_nu_mass; i=i+1)
      mnu_list[i] = cosmo->params.m_nu[i];

    if (isfinite(cosmo->params.A_s)) {
      norm_pk = cosmo->params.A_s;
    }
    else if (isfinite(cosmo->params.sigma8)) {
      norm_pk = cosmo->params.sigma8;
    }
    else {
      *status = CCL_ERROR_PARAMETERS;
      ccl_cosmology_set_status_message(cosmo,
             "ccl_musigma.c: ccl_rescale_musigma_s8(): "
             "neither A_s nor sigma8 defined.\n");
    }

    if (*status == 0) {
      params_GR = ccl_parameters_create(
        cosmo->params.Omega_c, cosmo->params.Omega_b, cosmo->params.Omega_k,
        cosmo->params.Neff, mnu_list, cosmo->params.N_nu_mass,
        cosmo->params.w0, cosmo->params.wa, cosmo->params.h,
        norm_pk, cosmo->params.n_s,
        cosmo->params.bcm_log10Mc, cosmo->params.bcm_etab,
        cosmo->params.bcm_ks, 0., 0., cosmo->params.nz_mgrowth,
        cosmo->params.z_mgrowth, cosmo->params.df_mgrowth, status);

      if (*status) {
        *status = CCL_ERROR_PARAMETERS;
        ccl_cosmology_set_status_message(cosmo,
               "ccl_musigma.c: ccl_rescale_musigma_s8(): "
               "could not make MG params.\n");
      }
    }

    if (*status == 0) {
      cosmo_GR = ccl_cosmology_create(params_GR, cosmo->config);

      if (cosmo_GR == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(cosmo,
               "ccl_musigma.c: ccl_rescale_musigma_s8(): "
               "error initializing cosmology\n");
      }
    }

    if (*status == 0) {
      ccl_cosmology_compute_growth(cosmo_GR, status);

      if (*status) {
        *status = CCL_ERROR_PARAMETERS;
        ccl_cosmology_set_status_message(cosmo,
               "ccl_musigma.c: ccl_rescale_musigma_s8(): "
               "could not init GR growth.\n");
      }
    }

    // Populate rescale_factor
    if (*status == 0) {
      for (int i=0; i<na; i++) {
        double D_mu = ccl_growth_factor_unnorm(cosmo, aa[i], status);
        double D_GR = ccl_growth_factor_unnorm(cosmo_GR, aa[i], status);
        double renorm = D_mu/D_GR;
        rescale_factor[i] *= renorm*renorm;
      }
      rescale_extra_musig = rescale_factor[na-1];

      if (*status) {
        *status = CCL_ERROR_PARAMETERS;
        ccl_cosmology_set_status_message(cosmo,
               "ccl_musigma.c: ccl_rescale_musigma_s8(): "
               "could not make MG and GR growth.\n");
      }
    }

    ccl_parameters_free(&params_GR);
    ccl_cosmology_free(cosmo_GR);
  }

  if(do_s8) {
    if (*status == 0) {
      double renorm = cosmo->params.sigma8/ccl_sigma8(cosmo, psp, status);
      renorm *= renorm;
      renorm /= rescale_extra_musig;

      for (int i=0; i<na; i++)
        rescale_factor[i] *= renorm;
    }
  }

  if (*status == 0) {
    size_t nk=0;
    if (psp->fka != NULL)
      nk=psp->fka->interp_object.xsize;

    // Rescale
    for (int i=0; i<na; i++) {
      if (psp->fa != NULL) {
        if(psp->is_log)
          psp->fa->y[i] += log(rescale_factor[i]);
        else
          psp->fa->y[i] *= rescale_factor[i];
      }
      else {
        for(int j=0; j<nk; j++) {
          if(psp->is_log)
            psp->fka->zarr[i*nk+j] += log(rescale_factor[i]);
          else
            psp->fka->zarr[i*nk+j] *= rescale_factor[i];
        }
      }
    }

    //Respline
    if(psp->fa != NULL) {
      gsl_spline *fa = gsl_spline_alloc(gsl_interp_cspline,
                                        psp->fa->size);
      if(fa == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(cosmo,
                                         "ccl_musigma.c: ccl_rescale_musigma_s8(): "
                                         "memory allocation\n");
      }
      if(*status==0) {
        int spstatus = gsl_spline_init(fa, psp->fa->x,
                                       psp->fa->y, psp->fa->size);
        if(spstatus) {
          *status = CCL_ERROR_MEMORY;
          ccl_cosmology_set_status_message(cosmo,
                                           "ccl_musigma.c: ccl_rescale_musigma_s8(): "
                                           "Error initializing spline\n");
        }
      }
      if(*status==0) {
        gsl_spline_free(psp->fa);
        psp->fa=fa;
      }
    }
    else {
      gsl_spline2d *fka = gsl_spline2d_alloc(gsl_interp2d_bicubic,
                                             psp->fka->interp_object.xsize,
                                             psp->fka->interp_object.ysize);
      if(fka == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(cosmo,
                                         "ccl_musigma.c: ccl_rescale_musigma_s8(): "
                                         "memory allocation\n");
      }
      if(*status==0) {
        int spstatus = gsl_spline2d_init(fka, psp->fka->xarr,
                                         psp->fka->yarr, psp->fka->zarr,
                                         psp->fka->interp_object.xsize,
                                         psp->fka->interp_object.ysize);
        if(spstatus) {
          *status = CCL_ERROR_MEMORY;
          ccl_cosmology_set_status_message(cosmo,
                                           "ccl_musigma.c: ccl_rescale_musigma_s8(): "
                                           "Error initializing spline\n");
        }
      }
      if(*status==0) {
        gsl_spline2d_free(psp->fka);
        psp->fka=fka;
      }
    }
  }

  free(rescale_factor);
}

ext: c

max_stars_repo_name: borisbolliet/CCL

max_issues_repo_name: borisbolliet/CCL

max_forks_repo_name: borisbolliet/CCL


text:
#include "matrix_support.h"
#include <cblas.h>

int main(int argc, char *argv[])
{
    /* enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102}; */
    /* enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113}; */
    /* enum CBLAS_UPLO {CblasUpper=121, CblasLower=122}; */
    /* enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132}; */
    /* enum CBLAS_SIDE {CblasLeft=141, CblasRight=142}; */

    CBLAS_ORDER order = CblasRowMajor;
    CBLAS_TRANSPOSE transA = CblasNoTrans;
    CBLAS_TRANSPOSE transB = CblasNoTrans;
    int m = 2;
    int k = 3;
    int n = 4;
    double alpha = 1;
    double beta = 1;
    double *A = create_matrix(m, k);
    double *B = create_matrix(k, n);
    double *C = create_matrix(m, n);
    double *D = create_matrix(n, k);
    double *x = create_vector(k);
    double *d = create_vector(m);
    printf("A\n");
    init_matrix_vector(A, m * k, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0);
    print_matrix(A, m, k);
    printf("\nB\n");
    init_matrix_vector(B, k * n, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
            10.0, 11.0);
    print_matrix(B, k, n);
    printf("\nD\n");
    init_matrix_vector(D, n * k, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
            10.0, 11.0);
    print_matrix(D, n, k);
    printf("\nx\n");
    init_matrix_vector(x, k * 1, -2.0, 2.0, 10.0);
    print_vector(x, k);

    /* void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, */
    /*                  const enum CBLAS_TRANSPOSE TransB, const int M, const int N, */
    /*                  const int K, const double alpha, const double *A, */
    /*                  const int lda, const double *B, const int ldb, */
    /*                  const double beta, double *C, const int ldc); */

    cblas_dgemm(order,
            transA,
            transB,
            m,
            n,
            k,
            alpha,
            A,
            k,
            B,
            n,
            beta,
            C,
            n);
    printf("\nC = A.dot(B)\n");
    print_matrix(C, m, n);

    printf("\nd = A.dot(x)\n");

    /*
     * dasum takes the sum of the absolute values
     */

    /**
     * ddot dot product between two vectors
     */

    /*
     * dnrm2 Eucledian norm of a vector
     */

    /*
     * dscal scales a vector by a constant
     */

    /*
     * daxpy constant times a vector plus a vector. uses unrolled loops for
     * increments equal to one.
     */

    /* void cblas_dgemv(const enum CBLAS_ORDER order, */
    /*              const enum CBLAS_TRANSPOSE TransA, const int M, const int N, */
    /*              const double alpha, const double *A, const int lda, */
    /*              const double *X, const int incX, const double beta, */
    /*              double *Y, const int incY); */
    cblas_dgemv(order,
            transA, // not transposed
            m,   // rows of a
            k,   // cols of a
            alpha, // scalar to multiply A with
            A,   //matrix A
            k, // lda: leading dimension of A
            x, // vector x
            1, // incx
            beta, // beta: scalar to multiply y with (0 having y as output)
            d, // vector d (Y in formula)
            1); // incy

    print_vector(d, m);
    return 0;
}


text:
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>

#include <gsl/gsl_fft_complex_float.h>

#define BASE_FLOAT
#include "templates_on.h"
#include "c_pass.h"
#include "c_init.c"
#include "c_main.c"
#include "c_pass_2.c"
#include "c_pass_3.c"
#include "c_pass_4.c"
#include "c_pass_5.c"
#include "c_pass_6.c"
#include "c_pass_7.c"
#include "c_pass_n.c"
#include "bitreverse.c"
#include "c_radix2.c"
#include "templates_off.h"
#undef  BASE_FLOAT

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:
#ifndef KGS_NULLSPACE_H
#define KGS_NULLSPACE_H

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

#include "math/SVD.h"
//#include "math/QR.h"

class Molecule; //Forward declaration

/**
 * Computes, stores, and maintains the nullspace of a gsl_matrix.
 * Supports projections of any gradient represented as a gsl_vector onto the null-space
 * and also analyzing and updating rigidity. Note that the Nullspace object only
 * accurately reflects the matrix (SVD) if the UpdateFromMatrix function has been called.
 * This is not done on construction.
 *
 * An SVD-pointer is passed to the Nullspace on construction. This pointer is kept through
 * the life-time of the Nullspace object (but ownership is not assumed) and if the contents
 * of the matrix is changed, the Nullspace can be update to reflect the change using the
 * abstract UpdateFromMatrix function.
 */
class Nullspace {
 protected:
  Nullspace(gsl_matrix* M);

 public:

  virtual ~Nullspace();

  /** Projects a vector on the nullspace */
  void projectOnNullSpace(gsl_vector *to_project, gsl_vector *after_project) const;

  /** Analyzes which dihedrals and hydrogen bonds are rigidified by constraints */
  void performRigidityAnalysis(gsl_matrix *HBondJacobian,gsl_matrix *DBondJacobian,gsl_matrix *HydrophobicBondJacobian);

  /** Update the Nullspace (and underlying SVD/QR) to reflect an updated state of the matrix */
  virtual void updateFromMatrix() = 0;

  /** Return the nullspace size */
  int getNullspaceSize() const { return m_nullspaceSize; }

  /** Return number of cycle DOFs */
  int getNumDOFs() const { return n; }

  /** Return number of dihedral DOFs that were rigidified in the last call to RigidityAnalysis. */
  int getNumRigidDihedrals() const { return numRigidDihedrals; }

  /** Return number of h-bond dihedrals that were rigidified in the last call to RigidityAnalysis. */
  int getNumRigidHBonds() const { return numRigidHBonds; }
  int getNumRigidDBonds() const { return numRigidDBonds; }
  int getNumRigidHydrophobicBonds() const { return numRigidHydrophobicBonds;}
  /** Return the underlying matrix. */
  gsl_matrix* getMatrix() const{ return m_matrix; }

  /** Return the basis of the nullspace as columns of a matrix */
  gsl_matrix *getBasis() const;

  /**
   * Returns true iff the angle specified by the argument is rigidified.
   * This result is only accurate if UpdateFromMatrix and RigidityAnalysis have both
   * been called.
   */
  bool isCovBondRigid(int angle_id){ return fabs(gsl_vector_get(m_rigidCovBonds, angle_id)-1.0)<0.001; }

  bool isHBondRigid(int bond_id){ return fabs(gsl_vector_get(m_rigidHBonds, bond_id)-1.0)<0.001; }

  bool isDBondRigid(int bond_id){ return fabs(gsl_vector_get(m_rigidDBonds, bond_id)-1.0)<0.001; }

  bool isHydrophobicBondRigid(int bond_id){ return fabs(gsl_vector_get(m_rigidHydrophobicBonds,bond_id)-1.0)<0.001; }

 protected:
  int m_nullspaceSize;         ///< Size of nullspace (rank of jacobian)
  int m, n;                    ///< Dimensions of underlying matrix (jacobian)
  gsl_matrix* m_matrix;        ///< Original matrix. Only set if using QR decomp (as it needs to be transposed)

  gsl_matrix* m_nullspaceBasis;///< Basis of the nullspace

 private:
  gsl_vector* m_rigidCovBonds; ///< Binary vector indicating which m_dofs are rigid
  gsl_vector* m_rigidHBonds;   ///< Binary vector indicating which h-bonds are rigid
  gsl_vector* m_rigidDBonds;   ///< Binary vector indicating which d-bonds are rigid
  gsl_vector* m_rigidHydrophobicBonds;  ///< Binary vector indicating which hydrophobic-bonds are rigid
  int numCoordinatedDihedrals; ///< Coordinated dihedrals
  int numRigidDihedrals;       ///< Rigid dihedrals
  int numRigidHBonds;///< Rigid hydrogen bonds
  int numRigidDBonds;///< Rigid hydrogen bonds
  int numRigidHydrophobicBonds; ///< Rigid hydrophobic bonds
  /// These values have to be chosen according to the numerical analysis
//  static constexpr double SINGVAL_TOL = 1.0e-12; //0.000000000001; // only generic 10^-12
  static constexpr double RIGID_TOL =   1.0e-9; //0.0000000001; //most molecules work between 1e-4 and 1e-10, exceptions only between 1e-8 and 1e-9

//  friend class Configuration;
};

#endif //KGS_nullptrSPACE_H

max_stars_repo_name: ExcitedStates/KGS

max_issues_repo_name: ExcitedStates/KGS

max_forks_repo_name: ExcitedStates/KGS


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

#ifndef __CCPM_INTERFACES_H__
#define __CCPM_INTERFACES_H__

#include <common/byte_span.h>
#include <common/types.h>
#include <gsl/pointers>
#include <cstddef>
#include <functional>
#include <vector>

namespace ccpm
{
	struct persister
	{
		using byte_span = common::byte_span;
		virtual void persist(byte_span span) = 0;
	protected:
		~persister() {}
	};

/*
 * Ownership callback can resolve the ambiguity about area
 *        ownership which occurs during phase (b) of allocate() and free()
 *        Intended to be used by IHeap::reconstitute when a crash during
 *        an IHeap::allocate or IHeap::free call has left area ownership
 *        in doubt.
 *
 * In order to support ownership resolution, it is expected that a caller to
 * allocate or free will
 *   (1) store and persist the initial value of ptr
 *   (2) call allocate or free
 *
 * It is expected that IHeap::allocate() will
 *   (1) determine an appropriate new value for ptr
 *   (2) persist a note indicating that the allocation status of the new value
 *       of ptr is indeterminate.
 *   (3) store and persist the new value of ptr
 *   (4) persist an invalidation of the note written in step 2.
 *
 * It is expected that IHeap::free() will
 *   (1) persist a note indicating that the allocation status of the value
 *       of ptr is indeterminate.
 *   (2) store and persist nullptr as the new value of ptr
 *   (3) return the area located by the original value of ptr to free status
 *   (4) persist an invalidation of the note written in step 1.
 */

using ownership_callback_type = std::function<bool(const void * ptr)>;
using ownership_callback_t = ownership_callback_type;

inline bool accept_all(const void *) { return true; }
/*
 * Allocators can expand to more than one coarse-grained region of
 * memory.
 */
struct region_vector_t : public std::vector<common::byte_span>
{
	using base = std::vector<common::byte_span>;
  explicit region_vector_t(void * ptr_, std::size_t size_)
    : region_vector_t(common::make_byte_span(static_cast<common::byte *>(ptr_), size_))
  {}
  explicit region_vector_t(const value_type &v) {
    push_back(v);
  }
  region_vector_t() {
  }
	const base &cbase() const { return *this; }
};

using region_span = gsl::span<common::byte_span>;

enum class Type_id : int64_t
  {
    None        = 0,
    Fixed_array = 0xF0,
  };

/**
 * Heap allocator (variable sized allocations)
 */
class IHeap
{
public:
  virtual ~IHeap() {}

  /* Reconstitute/initialize slab from existing memory
   *
   * @param regions Pointer/length regions of contiguous memory
   * @param resolver Access to an object which can resolve the ambiguity
   *        over area ownership which occurs during a phase of allocate()
   *        and free() calls. If the callee of the ownership_callback_t
   *        function might own the area located by the arument, the callee
   *        must return true. If the callee does not own the area, it should
   *        return false.
   * @param force_init If true, force re-setting to empty
   *
   *
   * @return : True if memory was reset to empty
   **/
  virtual bool reconstitute(const region_span regions,
                            ownership_callback_t resolver = [] (const void *) -> bool { return true; },
                            const bool force_init = false) = 0;

  /* Allocate memory
   *
   * @param ptr in/out: nullptr -> pointer to newly allocated memory
   *                    On a successful allocation the allocator will
   *                    first write and then persist ptr. Allocation
   *                    has three phases:
   *                    (a) ptr is not yet written; the allocator unambiguously
   *                        owns the free "area" the address of which it will
   *                        later write into ptr.
   *                    (b) allocator has written ptr, but has not yet persisted
   *                        ptr. Caller must be able to tell the allocator (see
   *                        reconstitute) whether caller has accepted ownership
   *                        of the allocated area. Acceptance is implicit.
   *                        Visibility of a non-null value in ptr, which happens
   *                        upon write in the non-crash case and upon persist in
   *                        the crash+reconstitute case, consititutes acceptance.
   *                    (c) allocator has persisted the ptr. Caller will know upon
   *                        successful return (normal case) or discovery of a
   *                        non-null value (crash+reconstitute case) that it owns
   *                        the area.
   *
   *                    The pointer shall be altered iff the function returns S_OK.
   *
   * @param size Size of memory to allocate in bytes
   * @param alignment Alignment for memory
   *
   * @return : S_OK, E_INVAL, E_BAD_PARAM, E_EMPTY
   **/
  virtual status_t allocate(void * & ptr,
                            std::size_t bytes,
                            std::size_t alignment) = 0;

  /* Free previously allocated memory.
   * @param ptr in/out: Pointer to memory to free, which must have been previously
   *                    allocated with the same size and alignment.
   *                    On a successful free the allocator will first nullify ptr
   *                    and then persist ptr. Free has has three phases, similar
   *                    to allocation:
   *                    (a) ptr is not yet written; the caller unambiguously
   *                        owns the "area".
   *                    (b) allocator has nullified ptr, but has not yet persisted
   *                        ptr. Caller must be able to tell the allocator (see
   *                        reconstitute) whether caller has relinquished ownership
   *                        of the allocated area. Relinquishment is implicit.
   *                        Visibility of a null value in ptr, which happens upon
   *                        write in the non-crash case and upon persist in the
   *                        crash+reconstitute case, consititutes relinquishment.
   *                    (c) allocator has persisted the (null-valued) ptr. Caller
   *                        will know upon successful return (normal case) or
   *                        discovery of null value (crash+reconstitute case) that
   *                        the allocator now owns the area.
   *
   *                    The pointer shall be altered iff the function returns S_OK.
   *
   * @param size Size of memory to free in bytes
   *
   * @return : S_OK, E_INVAL
   **/
  virtual status_t free(void * & ptr,
                        std::size_t bytes = 0) = 0;

  /* Return remaining space in bytes
   * @param out_size Size remaining in bytes
   *
   * @return S_OK or E_NOT_IMPL
   **/
  virtual status_t remaining(std::size_t& out_size) const = 0;

  /* Return a vector of all regions
   *
   * @return vector containing the initial regions and all subsequently added
   * regions, if the implementation keeps track of regions. Otherwise, an
   * empty region_vector_t.
   **/
  virtual region_vector_t get_regions() const = 0;
};

/**
 * Heap allocator (variable sized allocations)
 */
class IHeap_expandable : public IHeap
{
public:
  /* Add an additional regions to the heap
   * @param regions Pointer/length regions of contiguous memory
   **/
  virtual void add_regions(const region_span regions) = 0;

  /** Test whether an address is in any heap region
   * @param addr the address to test
   *
   * @erturn true iff the address is in some heap region
   */
  virtual bool includes(const void *ptr) const = 0;
};

class ILog
{
public:
  virtual ~ILog() = default;

  /*
   * Record old value of a region to the log.
   */
  virtual void add(void *begin, std::size_t size) = 0;
  /*
   * Record an allocation to the log.
   */
  virtual void allocated(void *&p, std::size_t size) = 0;
  /*
   * Record a free to the log.
   */
  virtual void freed(void *&p, std::size_t size) = 0;

  /*
   * commit all previous add/allocated/freed commands
   */
  virtual void commit() = 0;
  /*
   * Restore all data areas added after initialization (or the most recent clear)
   * to the values present at their last add command.
   */
  virtual void rollback() = 0;

};
}
#endif //  __CCPM_INTERFACES_H__

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
/* block/gsl_block_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_BLOCK_UINT_H__
#define __GSL_BLOCK_UINT_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>

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

__BEGIN_DECLS

struct gsl_block_uint_struct
{
  size_t size;
  unsigned int *data;
};

typedef struct gsl_block_uint_struct gsl_block_uint;

gsl_block_uint *gsl_block_uint_alloc (const size_t n);
gsl_block_uint *gsl_block_uint_calloc (const size_t n);
void gsl_block_uint_free (gsl_block_uint * b);

int gsl_block_uint_fread (FILE * stream, gsl_block_uint * b);
int gsl_block_uint_fwrite (FILE * stream, const gsl_block_uint * b);
int gsl_block_uint_fscanf (FILE * stream, gsl_block_uint * b);
int gsl_block_uint_fprintf (FILE * stream, const gsl_block_uint * b, const char *format);

int gsl_block_uint_raw_fread (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
int gsl_block_uint_raw_fwrite (FILE * stream, const unsigned int * b, const size_t n, const size_t stride);
int gsl_block_uint_raw_fscanf (FILE * stream, unsigned int * b, const size_t n, const size_t stride);
int gsl_block_uint_raw_fprintf (FILE * stream, const unsigned int * b, const size_t n, const size_t stride, const char *format);

size_t gsl_block_uint_size (const gsl_block_uint * b);
unsigned int * gsl_block_uint_data (const gsl_block_uint * b);

__END_DECLS

#endif /* __GSL_BLOCK_UINT_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_dcopy (const int N, const double *X, const int incX, double *Y,
	     const int incY)
{
#define BASE double
#include "source_copy_r.h"
#undef BASE
}


text:
#pragma once
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <cfloat>

#include <algorithm>
#include <vector>
#include <memory>
#include <stdexcept>
#include <map>

#include <hip/hip_runtime.h>
#include <hiprand_kernel.h>
#include <gsl/gsl_assert.h>
#include <nonstd/byte.hpp>
#include <nonstd/expected.hpp>
#include <nonstd/optional.hpp>
#include <nonstd/span.hpp>
#include <nonstd/ring_span.hpp>
#include <nonstd/string_view.hpp>
#include <nonstd/variant.hpp>
#include <nonstd/value_ptr.hpp>

#include "kernel_config.h"
#include "logger.h"
#include "math.h"

using nonstd::variant;
using nonstd::value_ptr;
using nonstd::make_value;
using nonstd::ring_span;
using nonstd::expected;
using nonstd::make_unexpected;
using nonstd::optional;
using nonstd::nullopt;
using nonstd::make_optional;
using nonstd::byte;
using nonstd::to_integer;
using nonstd::to_byte;
using nonstd::to_string;
using nonstd::string_view;
using nonstd::span;
using nonstd::make_span;
using bytearray = std::vector<byte>;

#define HIP_ASSERT(x) (Ensures((x)==hipSuccess))

struct hip_free_deleter {
    void operator()(void *p) const noexcept { hipFree(p); }
};

template<typename T>
using host_unique_ptr = std::unique_ptr<T, hip_free_deleter>;

template<class T> struct _host_unique_if {
    typedef void _single_object;
};

template<class T> struct _host_unique_if<T[]> {
    typedef host_unique_ptr<T[]> _unknown_bound;
};

template<class T, std::size_t N> struct _host_unique_if<T[N]> {
    typedef void _known_bound;
};

template<class T, class... Args>
typename _host_unique_if<T>::_single_object
hip_alloc_managed_unique(Args&&... args) = delete;

template<class T>
typename _host_unique_if<T>::_unknown_bound
hip_alloc_managed_unique(std::size_t n) {
    typedef typename std::remove_extent<T>::type U;
    T *tmp;
    HIP_ASSERT(hipMallocManaged(&tmp, n * sizeof(U)));
    return host_unique_ptr<T>(reinterpret_cast<U*>(tmp));
}

template<class T, class... Args>
typename _host_unique_if<T>::_known_bound
hip_alloc_managed_unique(Args&&...) = delete;

max_stars_repo_name: sargarass/Raytracer

max_issues_repo_name: sargarass/Raytracer

max_forks_repo_name: sargarass/Raytracer


text:
#ifndef _SPC_FITTING_H
#define _SPC_FITTING_H

#include <math.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>
#include "aXe_errors.h"

extern void
comp_vector_average (const double *xs, double *ys,
		       double *ws, double *yi, const int n, const int final);
extern void
comp_vector_median (const double *xs, double *ys,
		      double *ws, double *yi, const int n, const int final);

extern void
comp_vector_linear (const double *xs, double *ys,
		      double *ws, double *yi, const int n, const int final);

extern void
comp_vector_polyN (const int m, const double *xs, double *ys,
		    double *ws, double *yi, const int n, const int final);

extern void
fill_const_value(double *ys, double *ws, double *yi, const int n,
		 double cval, double stdev, const int final);

extern void
det_vector_average (const double *xs, double *ys,
		    double *ws, const int n, double *avg, double *std);
extern void
det_vector_median (const double *xs, double *ys,
		    double *ws, const int n, double *med, double *std);

extern gsl_vector *
det_vector_linear(const double *xs, double *ys, double *ws,
		  const int n, const int weight);

extern gsl_vector *
det_vector_poly_N (int m, const double *const xs, double *const ys,
		   double *const ws, const int n, gsl_vector *c,
		   gsl_matrix *cov);

extern void
fill_linear_interp(const double *const xs, double *const ys,
		   double *const ws, double *yi, const int n,
		   gsl_vector *interp, const int final);

extern void
fill_polyN_interp(const double *const xs, double *const ys,
		  double *const ws, double *yi, const int n,
		  gsl_vector *coeffs, gsl_matrix *cov, gsl_vector *interp,
		  const int final);
#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#ifndef VolViz_Types_h
#define VolViz_Types_h

#include <Eigen/Core>
#include <Eigen/Geometry>
#include <gsl/gsl>
#include <phys/units/quantity.hpp>

#include <array>

namespace VolViz {

using namespace gsl;

using Vector3f = Eigen::Vector3f;

using Matrix3 = Eigen::Matrix3f;
using Matrix4 = Eigen::Matrix4f;

/// Position in 3D euclidean space
using Position = Eigen::Vector3f;
/// Position in homogenous coordinates
using PositionH = Eigen::Vector4f;

/// Position in 2D space
using Position2 = Eigen::Vector2f;

using Size2 = Eigen::Matrix<std::size_t, 2, 1>;
using Size3 = Eigen::Matrix<std::size_t, 3, 1>;

using Size3f = Eigen::Vector3f;

/// An 1D range
template <class T> struct Range {
  T min, max;
  constexpr auto length() const noexcept { return max - min; }
};

/// Normalized RGB color
using Color = Eigen::Vector3f;

namespace Colors {
inline auto Black() noexcept { return Color::Zero(); }
inline auto White() noexcept { return Color::Ones(); }
inline auto Red() noexcept { return Color::UnitX(); }
inline auto Green() noexcept { return Color::UnitY(); }
inline auto Blue() noexcept { return Color::UnitZ(); }
inline auto Yellow() noexcept { return Red() + Green(); }
inline auto Magenta() noexcept { return Red() + Blue(); }
inline auto Cyan() noexcept { return Blue() + Green(); }
}

/// 6-DOF orientation, represented as a quaternion
using Orientation = Eigen::Quaternionf;

struct DepthRange {
  float near{1.f}, far{-1.f};
};

enum class Axis { X, Y, Z };

using Scale = float;

using Length = phys::units::quantity<phys::units::length_d>;
using Angle = double;
using PhysicalPosition = Eigen::Matrix<Length, 3, 1>;

using phys::units::meter;
using phys::units::rad;
using phys::units::degree_angle;

using phys::units::centi;
using phys::units::milli;
using phys::units::micro;
using phys::units::nano;

using phys::units::abs;

namespace literals = phys::units::literals;

using VoxelSize = std::array<Length, 3>;

/// @brief Conveniance wrapper for span<>
/// @{
template <class T>
inline constexpr auto as_span(T *ptr,
                              typename span<T>::index_type size) noexcept {
  return span<T>(ptr, size);
}

template <
    class Container,
    class = std::enable_if_t<
        !details::is_span<Container>::value &&
        !details::is_std_array<Container>::value &&
        std::is_convertible<typename Container::pointer,
                            decltype(std::declval<Container>().data())>::value>>
inline constexpr auto as_span(Container &c) {
  return span<typename Container::value_type>(c);
}

template <
    class Container,
    class = std::enable_if_t<
        !details::is_span<Container>::value &&
        !details::is_std_array<Container>::value &&
        std::is_convertible<typename Container::pointer,
                            decltype(std::declval<Container>().data())>::value>>
inline constexpr auto as_span(Container const &c) {
  return span<typename Container::value_type const>(c);
}
///@}

} // namespace VolViz

#if __cplusplus < 201703L

namespace std {

template<class T>
constexpr const T& clamp( const T& v, const T& lo, const T& hi ) {
    return clamp( v, lo, hi, std::less<>() );
}

template<class T, class Compare>
constexpr const T& clamp( const T& v, const T& lo, const T& hi, Compare comp ) {
      return assert( !comp(hi, lo) ),
             comp(v, lo) ? lo : comp(hi, v) ? hi : v;
  }

}

#endif

#endif // VolViz_Types_h

max_stars_repo_name: ithron/VolViz

max_issues_repo_name: ithron/VolViz

max_forks_repo_name: ithron/VolViz


text:
/* min/test.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 <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#include "test.h"
#include "min.h"

/* stopping parameters */

const double EPSABS = 0.001 ;
const double EPSREL = 0.001 ;

const unsigned int MAX_ITERATIONS = 100;

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

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

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

  const gsl_min_fminimizer_type * fminimizer[4] ;
  const gsl_min_fminimizer_type ** T;

  gsl_ieee_env_setup ();

  fminimizer[0] = gsl_min_fminimizer_goldensection;
  fminimizer[1] = gsl_min_fminimizer_brent;
  fminimizer[2] = gsl_min_fminimizer_quad_golden;
  fminimizer[3] = 0;

  F_cos = create_function (f_cos) ;
  F_func1 = create_function (func1) ;
  F_func2 = create_function (func2) ;
  F_func3 = create_function (func3) ;
  F_func4 = create_function (func4) ;

  gsl_set_error_handler (&my_error_handler);

  for (T = fminimizer ; *T != 0 ; T++)
    {
      test_f (*T, "cos(x) [0 (3) 6]", &F_cos, 0.0, 3.0, 6.0, M_PI);
      test_f (*T, "x^4 - 1 [-3 (-1) 17]", &F_func1, -3.0, -1.0, 17.0, 0.0);
      test_f (*T, "sqrt(|x|) [-2 (-1) 1.5]", &F_func2, -2.0, -1.0, 1.5, 0.0);
      test_f (*T, "func3(x) [-2 (3) 4]", &F_func3, -2.0, 3.0, 4.0, 1.0);
      test_f (*T, "func4(x) [0 (0.782) 1]", &F_func4, 0, 0.782, 1.0, 0.8);

      test_f_e (*T, "invalid range check [4, 0]", &F_cos, 4.0, 3.0, 0.0, M_PI);
      test_f_e (*T, "invalid range check [1, 1]", &F_cos, 1.0, 1.0, 1.0, M_PI);
      test_f_e (*T, "invalid range check [-1, 1]", &F_cos, -1.0, 0.0, 1.0, M_PI);
    }
  test_bracket("cos(x) [1,2]",&F_cos,1.0,2.0,15);
  test_bracket("sqrt(|x|) [-1,0]",&F_func2,-1.0,0.0,15);
  test_bracket("sqrt(|x|) [-1,-0.6]",&F_func2,-1.0,-0.6,15);
  test_bracket("sqrt(|x|) [-1,1]",&F_func2,-1.0,1.0,15);

  exit (gsl_test_summary ());
}

void
test_f (const gsl_min_fminimizer_type * T,
        const char * description, gsl_function *f,
        double lower_bound, double middle, double upper_bound,
        double correct_minimum)
{
  int status;
  size_t iterations = 0;
  double m, a, b;
  double x_lower, x_upper;
  gsl_min_fminimizer * s;

  x_lower = lower_bound;
  x_upper = upper_bound;

  s = gsl_min_fminimizer_alloc (T) ;
  gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;

  do
    {
      iterations++ ;

      status = gsl_min_fminimizer_iterate (s);

      m = gsl_min_fminimizer_x_minimum(s);
      a = gsl_min_fminimizer_x_lower(s);
      b = gsl_min_fminimizer_x_upper(s);

#ifdef DEBUG
      printf("%.12f %.18f %.12f %.18f %.12f %.18f status=%d\n",
             a, GSL_FN_EVAL(f, a), m, GSL_FN_EVAL(f, m), b, GSL_FN_EVAL(f, b), status);
#endif

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

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

      if (status) break ;

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

  gsl_test (status, "%s, %s (%g obs vs %g expected) ",
            gsl_min_fminimizer_name(s), description,
            gsl_min_fminimizer_x_minimum(s), correct_minimum);

  /* check the validity of the returned result */

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

  gsl_min_fminimizer_free (s);

}

void
test_f_e (const gsl_min_fminimizer_type * T,
          const char * description, gsl_function *f,
          double lower_bound, double middle, double upper_bound,
          double correct_minimum)
{
  int status;
  size_t iterations = 0;
  double x_lower, x_upper;
  double a, b;
  gsl_min_fminimizer * s;

  x_lower = lower_bound;
  x_upper = upper_bound;

  s = gsl_min_fminimizer_alloc (T) ;
  status = gsl_min_fminimizer_set (s, f, middle, x_lower, x_upper) ;

  if (status != GSL_SUCCESS)
    {
      gsl_min_fminimizer_free (s) ;
      gsl_test (status == GSL_SUCCESS, "%s, %s", T->name, description);
      return ;
    }

  do
    {
      iterations++ ;
      gsl_min_fminimizer_iterate (s);
      a = gsl_min_fminimizer_x_lower(s);
      b = gsl_min_fminimizer_x_upper(s);

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

  gsl_test (!status, "%s, %s", gsl_min_fminimizer_name(s), description,
            gsl_min_fminimizer_x_minimum(s) - correct_minimum);

  gsl_min_fminimizer_free (s);
}

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

int
test_bracket (const char * description,gsl_function *f,double lower_bound,
              double upper_bound, unsigned int max)
{
  int status;
  double x_lower, x_upper;
  double f_upper,f_lower,f_minimum;
  double x_minimum;

  x_lower=lower_bound;
  x_upper=upper_bound;
  SAFE_FUNC_CALL (f,x_lower,&f_lower);
  SAFE_FUNC_CALL (f,x_upper,&f_upper);
  status=gsl_min_find_bracket(f,&x_minimum,&f_minimum,&x_lower,&f_lower,&x_upper,&f_upper,max);
  gsl_test (status,"%s, interval: [%g,%g], values: (%g,%g), minimum at: %g, value: %g",
            description,x_lower,x_upper,f_lower,f_upper,x_minimum,f_minimum);
  return status;
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#include "userFunc.h"
#include <math.h>
#include <gsl/gsl_sf_bessel.h>

/**************************************************************************/
/* This defines userFunc routines for the Pringle (1981) ring test        */
/**************************************************************************/

void
userEOS(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	void *params,
	double *gamma, double *delta) {
  fprintf(stderr,
	  "Warning: userEOS function called but not implemented!\n");
  return;
}

void
userAlpha(const double t, const double dt, const grid *grd,
	  const double *col, const double *pres, const double *eInt,
	  const double *gamma, const double *delta,
	  void *params,
	  double *alpha) {
  /* alpha = nu col/pres vphi / r */

  int i;
  double nu = *((double *) params);

  for (i=0; i<grd->nr; i++)
    alpha[i] = nu * col[i]/pres[i] * grd->vphi_g[i+1] / grd->r_g[i+1];

}

void
userMassSrc(const double t, const double dt, const grid *grd,
	    const double *col, const double *pres, const double *eInt,
	    const double *gamma, const double *delta,
	    void *params,
	    double *massSrc) {
  fprintf(stderr,
	  "Warning: userMassSrc function called but not implemented!\n");
  return;
}

void
userIntEnSrc(const double t, const double dt, const grid *grd,
	     const double *col, const double *pres, const double *eInt,
	     const double *gamma, const double *delta,
	     void *params,
	     double *intEnSrc) {
  fprintf(stderr,
	  "Warning: userIntEnSrc function called but not implemented!\n");
  return;
}

void
userIBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type ibc_pres, const enth_bc_type ibc_enth,
	void *params,
	double *ibc_pres_val, double *ibc_enth_val) {
  double nu = ((double *) params)[0];
  double r0 = ((double *) params)[1];
  double m0 = ((double *) params)[2];
  double colMin = ((double *) params)[3];
  double pOverCol = ((double *) params)[4];
  double x = grd->r_g[0]/r0;
  double tau = 12.0*nu*t/SQR(r0) + SMALL;
  double sigma0 = m0/(M_PI*SQR(r0));
  double sigma;

  sigma = sigma0/(pow(x, 0.25)*tau) * exp(-SQR(x-1.0)/tau) *
    gsl_sf_bessel_Inu_scaled(0.25, 2*x/tau);
  sigma = (sigma < colMin) ? colMin : sigma;
  *ibc_pres_val = -3.0*M_PI*grd->r_g[0]*nu*grd->vphi_g[0]*sigma;
  *ibc_enth_val = gamma[0]/(gamma[0]-1)*pOverCol;
}

void
userOBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type obc_pres, const enth_bc_type obc_enth,
	void *params,
	double *obc_pres_val, double *obc_enth_val) {
  double nu = ((double *) params)[0];
  double r0 = ((double *) params)[1];
  double m0 = ((double *) params)[2];
  double colMin = ((double *) params)[3];
  double pOverCol = ((double *) params)[4];
  double x = grd->r_g[grd->nr+1]/r0;
  double tau = 12.0*nu*t/SQR(r0) + SMALL;
  double sigma0 = m0/(M_PI*SQR(r0));
  double sigma;

  sigma = sigma0/(pow(x, 0.25)*tau) * exp(-SQR(x-1.0)/tau) *
    gsl_sf_bessel_Inu_scaled(0.25, 2*x/tau);
  sigma = (sigma < colMin) ? colMin : sigma;
  *obc_pres_val = -3.0*M_PI*grd->r_g[grd->nr+1]*nu*
    grd->vphi_g[grd->nr+1]*sigma;
  *obc_enth_val = gamma[0]/(gamma[0]-1)*pOverCol;
}

void
userPreTimestep(const double t, const double dt,
		const grid *grd, double *col, double *pres,
		double *eInt, double *mBnd, double *eBnd,
		double *mSrc, double *eSrc,
		void *params, const unsigned long nUserOut,
		double *userOut) {
  fprintf(stderr,
	  "Warning: userPreTimestep function called but not implemented!\n");
  return;
}

void
userPostTimestep(const double t, const double dt,
		 const grid *grd, double *col, double *pres,
		 double *eInt, double *mBnd, double *eBnd,
		 double *mSrc, double *eSrc,
		 void *params, const unsigned long nUserOut,
		 double *userOut) {
  fprintf(stderr,
	  "Warning: userPostTimestep function called but not implemented!\n");
  return;
}

void
userCheckRead(
	      FILE *fp, grid *grd, const unsigned long nOut,
	      double *tOut, double *colOut,
	      double *presOut, double *eIntOut, double *mBndOut,
	      double *eBndOut, double *mSrcOut, double *eSrcOut,
	      const unsigned long nUserOut, double *userOut,
	      void *params
	      ) {
  fprintf(stderr,
	  "Warning: userCheckRead function called but not implemented!\n");
  return;
}

void
userCheckWrite(
	      FILE *fp,
	      const grid *grd, const unsigned long nOut,
	      const double *tOut, const double *colOut,
	      const double *presOut, const double *eIntOut,
	      const double *mBndOut, const double *eBndOut,
	      const double *mSrcOut, const double *eSrcOut,
	      const unsigned long nUserOut, const double *userOut,
	      const void *params
	      ) {
  fprintf(stderr,
	  "Warning: userCheckWrite function called but not implemented!\n");
  return;
}

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
#ifndef SOURCE_ROBUSTNESSMAIN_H_
#define SOURCE_ROBUSTNESSMAIN_H_

// *** source/robustness_main.h ***
// Author: Kevin Wolz, date: 08/2018
//
// Numerical sampling of the Robustness distribution (GP and WPL).

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <stdlib.h>
#include <sys/stat.h>
#include <ctime>
#include <time.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "utils.h"
#include "lin_gauss.h"
#include "rob_analyt.h"
#include "print_to_file.h"
#include "parse_from_init_file.h"

double RobustnessWeakPriorLimit(double chi_squared_1,     // subset 1
                                double chi_squared_2,     // subset 2
                                double chi_squared_tot,   // total data
                                double det_fisher_1,      // determinants
                                double det_fisher_2,
                                double det_fisher_tot,
                                double det_prior);

double RobustnessGeneralPrior(double chi_squared_1,
                              double chi_squared_2,
                              double chi_squared_tot,
                              gsl_matrix* fisher_matrix_1,
                              gsl_matrix* fisher_matrix_2,
                              gsl_matrix* prior_matrix,
                              gsl_vector* param_values_1, // fiducial parameters
                              gsl_vector* param_values_2,
                              gsl_vector* param_values_tot,
                              gsl_vector* prior_param_values,// prior parameters
                              unsigned int pardim);      // number of parameters

int NumericalSamplingOfRobustness(unsigned int repetitions, // number of samples
                                  unsigned int datdim,  // number of data points
                                  bool is_weak_prior_limit,
                                  LinGauss model1,          // subset 1
                                  LinGauss model2,          // subset 2
                                  LinGauss model_tot,       // total data
                                  std::string& filepath);   // path of data logs

#endif // SOURCE_ROBUSTNESSMAIN_H_

max_stars_repo_name: kevguitar/Robustness

max_issues_repo_name: kevguitar/Robustness

max_forks_repo_name: kevguitar/Robustness


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 string_81448b68_30d7_4d64_a67a_cecef7c90139_h
#define string_81448b68_30d7_4d64_a67a_cecef7c90139_h

#include <string>
#include <assert.h>
#include <gslib/type.h>

__gslib_begin__

template<class _element>
class _std_string:
    public std::basic_string<_element, std::char_traits<_element>, std::allocator<_element> >
{
public:
    typedef _element protoch;
    typedef std::char_traits<protoch> prototr;
    typedef std::allocator<protoch> protoalloc;
    typedef std::basic_string<protoch, prototr, protoalloc> protoref;
    friend gs_export int _vsprintf(_std_string<char>& str, const char* fmt, va_list ap);
    friend gs_export int _vsprintf(_std_string<wchar>& str, const wchar* fmt, va_list ap);
};

template<class _element>
struct _string_tool
{
    typedef _element element;
};

template<>
struct _string_tool<char>
{
    static void copy(char* des, size_t size, const char* src) { strcpy_s(des, size, src); }
    static int length(const char* str) { return (int)strlen(str); }
    static int compare(const char* s1, const char* s2) { return (int)strcmp(s1, s2); }
    static int compare_cl(const char* s1, const char* s2) { return (int)_stricmp(s1, s2); }
    static int compare(const char* s1, const char* s2, int cnt) { return (int)strncmp(s1, s2, cnt); }
    static int compare_cl(const char* s1, const char* s2, int cnt) { return (int)_strnicmp(s1, s2, cnt); }
    static int vprintf(char* des, int size, const char* fmt, va_list ap) { return _vsprintf_s_l(des, size, fmt, 0, ap); }
    static int printf(char* des, int size, const char* fmt, ...)
    {
        va_list ptr;
        va_start(ptr, fmt);
        return _vsprintf_s_l(des, size, fmt, 0, ptr);
    }
    static int vsscanf(const char* src, const char* fmt, va_list ap) { return vsscanf_s(src, fmt, ap); }
    static int sscanf(const char* src, const char* fmt, ...)
    {
        va_list ptr;
        va_start(ptr, fmt);
        return vsscanf_s(src, fmt, ptr);
    }
    static int ctlprintf(const char* src, ...)
    {
        va_list ptr;
        va_start(ptr, src);
        return ::vprintf(src, ptr);
    }
    static int to_int(const char* src) { return atoi(src); }
    static int to_int(const char* src, int radix) { return strtol(src, 0, radix); }
    static uint to_uint(const char* src, int radix) { return strtoul(src, 0, radix); }
    static int64 to_int64(const char* src) { return _atoi64(src); }
    static int64 to_int64(const char* src, int radix) { return _strtoi64(src, 0, radix); }
    static uint64 to_uint64(const char* src, int radix) { return _strtoui64(src, 0, radix); }
    static real to_real(const char* src) { return atof(src); }
    static void from_int(int i, char* des, int size, int radix) { _itoa_s(i, des, size, radix); }
    static void from_int64(int64 i, char* des, int size, int radix) { _i64toa_s(i, des, size, radix); }
    static void from_real(char* des, int size, real d) { from_real(des, size, d, size-2); }
    static void from_real(char* des, int size, real d, int precis) { _gcvt_s(des, size, d, precis); }
    static void to_lower(char* str, int size) { _strlwr_s(str, size); }
    static void to_upper(char* str, int size) { _strupr_s(str, size); }
    static const char* find(const char* src, const char* tar) { return strstr(src, tar); }
    static int test(const char* src, const char* cpset) { return (int)strcspn(src, cpset); }
    static const char* _test(const char* src, const char* cpset) { return strpbrk(src, cpset); }
};

template<>
struct _string_tool<wchar>
{
    static void copy(wchar* des, size_t size, const wchar* src) { wcscpy_s(des, size, src); }
    static int length(const wchar* str) { return (int)wcslen(str); }
    static int compare(const wchar* s1, const wchar* s2) { return (int)wcscmp(s1, s2); }
    static int compare_cl(const wchar* s1, const wchar* s2) { return (int)_wcsicmp(s1, s2); }
    static int compare(const wchar* s1, const wchar* s2, int cnt) { return (int)wcsncmp(s1, s2, cnt); }
    static int compare_cl(const wchar* s1, const wchar* s2, int cnt) { return (int)_wcsnicmp(s1, s2, cnt); }
    static int vprintf(wchar* des, int size, const wchar* fmt, va_list ap) { return _vswprintf_s_l(des, size, fmt, 0, ap); }
    static int printf(wchar* des, int size, const wchar* fmt, ...)
    {
        va_list ptr;
        va_start(ptr, fmt);
        return _vswprintf_s_l(des, size, fmt, 0, ptr);
    }
    static int vsscanf(const wchar* src, const wchar* fmt, va_list ap) { return vswscanf_s(src, fmt, ap); }
    static int sscanf(const wchar* src, const wchar* fmt, ...)
    {
        va_list ptr;
        va_start(ptr, fmt);
        return vswscanf_s(src, fmt, ptr);
    }
    static int ctlprintf(const wchar* src, ...)
    {
        va_list ptr;
        va_start(ptr, src);
        return ::vwprintf(src, ptr);
    }
    static int to_int(const wchar* src) { return _wtoi(src); }
    static int to_int(const wchar* src, int radix) { return wcstol(src, 0, radix); }
    static uint to_uint(const wchar* src, int radix) { return wcstoul(src, 0, radix); }
    static int64 to_int64(const wchar* src) { return _wtoi64(src); }
    static int64 to_int64(const wchar* src, int radix) { return _wcstoi64(src, 0, radix); }
    static uint64 to_uint64(const wchar* src, int radix) { return _wcstoui64(src, 0, radix); }
    static real to_real(const wchar* src) { return _wtof(src); }
    static void from_int(int i, wchar* des, int size, int radix) { _itow_s(i, des, size, radix); }
    static void from_int64(int64 i, wchar* des, int size, int radix) { _i64tow_s(i, des, size, radix); }
    static void from_real(wchar* des, int size, real d) { from_real(des, size, d, size-2); }
    static void from_real(wchar* des, int size, real d, int precis)
    {
        _gcvt_s((char*)des, size, d, precis);
        int len = (int)strlen((char*)des); /* call it in asm wasn't available in /o2 */
        __asm
        {
            mov  ecx, len;
            mov  esi, des;
            add  esi, ecx;
            mov  edi, esi;
            add  edi, ecx;
            mov  word ptr[edi], 0;
            std;
        rptag:
            movsb;
            mov  byte ptr[edi], 0;
            dec  edi;
            loop rptag;
            cld;
        }
    }
    static void to_lower(wchar* str, int size) { _wcslwr_s(str, size); }
    static void to_upper(wchar* str, int size) { _wcsupr_s(str, size); }
    static const wchar* find(const wchar* src, const wchar* tar) { return wcsstr(src, tar); }
    static int test(const wchar* src, const wchar* cpset) { return (int)wcsspn(src, cpset); }
    static const wchar* _test(const wchar* src, const wchar* cpset) { return wcspbrk(src, cpset); }
};

gs_export extern const char* get_mbcs_elem(const char* str, uint& c);
gs_export extern const char* get_mbcs_elem(const char* str, uint& c, const char* end);
gs_export extern int convert_to_wide(wchar out[], int size, const char* str, int len = -1);
gs_export extern int convert_to_byte(char out[], int size, const wchar* str, int len = -1);
gs_export extern int convert_utf8_to_wide(wchar out[], int size, const char* str, int len = -1);
gs_export extern int convert_wide_to_utf8(char out[], int size, const wchar* str, int len = -1);

struct _string_caseful {};
struct _string_caseless {};

template<class _element,
    class _cpcase = _string_caseful>
class _string:
    public _std_string<_element>
{
public:
    typedef _element element;
    typedef _std_string<element> inheritref;
    typedef _string_tool<element> _strtool;
    typedef typename inheritref::protoref protoref;
    typedef _string<_element, _cpcase> myref;

    template<class _cpcase> static int comparefunc(const _string* s, const element* str);
    template<class _cpcase> static  int comparefunc(const _string* s, const element* str, int len);
    template<> static int comparefunc<_string_caseful>(const _string* s, const element* str) { return _strtool::compare(s->c_str(), str); }
    template<> static int comparefunc<_string_caseless>(const _string* s, const element* str) { return _strtool::compare_cl(s->c_str(), str); }
    template<> static int comparefunc<_string_caseful>(const _string* s, const element* str, int len) { return _strtool::compare(s->c_str(), str, len); }
    template<> static int comparefunc<_string_caseless>(const _string* s, const element* str, int len) { return _strtool::compare_cl(s->c_str(), str, len); }

protected:
#ifdef _MSC_VER
#if (_MSC_VER >= 1920)
    element* _stl_rawstr() { return const_cast<element*>(c_str()); }
    const element* _stl_rawstr() const { return c_str(); }
#elif (_MSC_VER >= 1914)
    element* _stl_rawstr() { return this->_Get_data()._Myptr(); }
    const element* _stl_rawstr() const { return this->_Get_data()._Myptr(); }
#else
    element* _stl_rawstr() { return this->_Myptr(); }
    const element* _stl_rawstr() const { return this->_Myptr(); }
#endif
#if (_MSC_VER >= 1920)
    int _stl_cap() const { return (int)capacity(); }
    void _stl_fix() { resize(_strtool::length(c_str())); }
#elif (_MSC_VER >= 1914)
    int& _stl_cap() { return (int&)this->_Get_data()._Myres; }
    int& _stl_size() { return (int&)this->_Get_data()._Mysize; }
    void _stl_fix() { this->_stl_size() = _strtool::length(this->c_str()); }
#elif ((_MSC_VER >= 1900) && (_MSC_VER < 1910))
    int& _stl_cap() { return (int&)this->_Myres(); }
    int& _stl_size() { return (int&)this->_Mysize(); }
    void _stl_fix() { this->_Mysize() = (size_t)_strtool::length(this->c_str()); }
#else
    int& _stl_cap() { return (int&)this->_Myres; }
    int& _stl_size() { return (int&)this->_Mysize; }
    void _stl_fix() { this->_Mysize = (size_t)_strtool::length(this->c_str()); }
#endif
#if (_MSC_VER >= 1920)
    void _stl_eos(int pos) { this->resize(pos); }
#else
    void _stl_eos(int pos) { this->_Eos(pos); }
#endif
#endif      /* endif _MSC_VER */

private:
    template<class _element>
    void convert_from(const wchar* str);
    template<class _element>
    void convert_from(const char* str);
    template<>
    void convert_from<wchar>(const wchar* str) { assign(str); }
    template<>
    void convert_from<char>(const wchar* str)
    {
        int len = convert_to_byte(0, 0, str);
        resize(len);
        convert_to_byte(_stl_rawstr(), len, str);
        _stl_fix();
    }
    template<>
    void convert_from<wchar>(const char* str)
    {
        int len = convert_to_wide(0, 0, str);
        resize(len);
        convert_to_wide(_stl_rawstr(), len, str);
        _stl_fix();
    }
    template<>
    void convert_from<char>(const char* str) { assign(str); }
    template<class _element>
    void convert_from(const char* str, int len);
    template<class _element>
    void convert_from(const wchar* str, int len);
    template<>
    void convert_from<wchar>(const wchar* str, int len) { assign(str, len); }
    template<>
    void convert_from<wchar>(const char* str, int len)
    {
        int l = convert_to_wide(0, 0, str, len);
        resize(l + 1);
        convert_to_wide(_stl_rawstr(), l, str, len);
        _stl_eos(l);
    }
    template<>
    void convert_from<char>(const wchar* str, int len)
    {
        int l = convert_to_byte(0, 0, str, len);
        resize(l + 1);
        convert_to_byte(_stl_rawstr(), l, str, len);
        _stl_eos(l);
    }
    template<>
    void convert_from<char>(const char* str, int len) { assign(str, len); }

public:
    _string() {}
    _string(const element* str) { if(str) this->assign(str); }
    _string(const element* str, int len) { if(str) this->assign(str, len); }
    _string(element c, int ctr) { this->assign(ctr, c); }
    void destroy() { clear(); }
    int length() const { return (int)this->size(); }
    element& front() { return this->at(0); }
    const element& front() const { return this->at(0); }
    element& back() { return this->at(this->size() - 1); }
    const element& back() const { return this->at(this->size() - 1); }
    element pop_back()
    {
        element e = back();
        resize(this->size() - 1);
        return e;
    }
    int test(const element* cpset) const
    {
        if(this->empty())
            return (int)this->npos;
        int pos = _strtool::test(this->c_str(), cpset);
        if(pos == length())
            return (int)this->npos;
        return pos;
    }
    int test(int off, const element* cpset) const
    {
        if(this->empty() || off >= length())
            return (int)this->npos;
        int pos = _strtool::test(this->c_str() + off, cpset);
        if(pos == length())
            return (int)this->npos;
        return pos;
    }
    const element* _test(const element* cpset) const
    {
        if(this->empty())
            return nullptr;
        return _strtool::_test(this->c_str(), cpset);
    }
    const element* _test(int off, const element* cpset) const
    {
        if(this->empty() || off >= length())
            return nullptr;
        return _strtool::_test(this->c_str() + off, cpset);
    }
    void format(const element* fmt, ...)
    {
        va_list ptr;
        va_start(ptr, fmt);
        _vsprintf(*this, fmt, ptr);
    }
    void formatv(const element* fmt, va_list ptr)
    {
        _vsprintf(*this, fmt, ptr);
    }
    /*
    void format(int size, const element* fmt, ...)
    {
        _stl_grow(size, false);
        va_list ptr;
        va_start(ptr, fmt);
        _stl_size() = _strtool::vprintf(_stl_rawstr(), _stl_cap(), fmt, ptr);
    }
    */
    _string& to_lower()
    {
        if(this->empty())
            return *this;
        _strtool::to_lower(_stl_rawstr(), _stl_cap());
        return *this;
    }
    _string& to_upper()
    {
        if(this->empty())
            return *this;
        _strtool::to_upper(_stl_rawstr(), _stl_cap());
        return *this;
    }
    int to_int() const { return _strtool::to_int(_stl_rawstr()); }
    int to_int(int radix) const { return _strtool::to_int(_stl_rawstr(), radix); }
    int64 to_int64() const { return _strtool::to_int64(_stl_rawstr()); }
    int64 to_int64(int radix) const { return _strtool::to_int64(_stl_rawstr(), radix); }
    real to_real() const { return _strtool::to_real(_stl_rawstr()); }
    _string& from(const char* str)
    {
        convert_from<element>(str);
        return *this;
    }
    _string& from(const wchar* str)
    {
        convert_from<element>(str);
        return *this;
    }
    _string& from(const char* str, int len)
    {
        convert_from<element>(str, len);
        return *this;
    }
    _string& from(const wchar* str, int len)
    {
        convert_from<element>(str, len);
        return *this;
    }
    _string& from_int(int i, int radix = 10)
    {
        resize(12);
        _strtool::from_int(i, _stl_rawstr(), _stl_cap(), radix);
        _stl_fix();
        return *this;
    }
    _string& from_int64(int64 i, int radix = 10)
    {
        resize(24);
        _strtool::from_int64(i, _stl_rawstr(), _stl_cap(), radix);
        _stl_fix();
        return *this;
    }
    _string& from_real(real d, int precis = 22)
    {
        resize(precis + 2);
        _strtool::from_real(_stl_rawstr(), _stl_cap(), d, precis);
        _stl_fix();
        return *this;
    }
    bool starts_with(const element* str, int len) const
    {
        assert(str);
        int s = length();
        if(len > s)
            return false;
        return compare(str, len) == 0;
    }
    bool ends_with(const element* str, int len) const
    {
        assert(str);
        int s = length();
        if(len > s)
            return false;
        return compare(str + s - len, len) == 0;
    }
    bool starts_with(const element* str) const { return starts_with(str, _strtool::length(str)); }
    bool starts_with(const myref& str) const { return starts_with(str.c_str(), str.length()); }
    bool starts_with(const myref& str, int len) const { return starts_with(str.c_str(), len); }
    bool ends_with(const element* str) const { return ends_with(str, _strtool::length(str)); }
    bool ends_with(const myref& str) const { return ends_with(str.c_str(), str.length()); }
    bool ends_with(const myref& str, int len) const { return ends_with(str.c_str(), len); }
    int compare(const element* str) const { return _strtool::compare(_stl_rawstr(), str); }
    int compare(const element* str, int len) const { return _strtool::compare(_stl_rawstr(), str, len); }
    int compare(int off, const element* str, int len) const { return _strtool::compare(_stl_rawstr() + off, str, len); }
    int compare_cl(const element* str) const { return _strtool::compare_cl(_stl_rawstr(), str); }
    int compare_cl(const element* str, int len) const { return _strtool::compare_cl(_stl_rawstr(), str, len); }
    int compare_cl(int off, const element* str, int len) const { return _strtool::compare_cl(_stl_rawstr() + off, str, len); }
    bool greater(const element* str) const { return comparefunc<_cpcase>(this, str) > 0; }
    bool greater(const element* str, int len) const { return comparefunc<_cpcase>(this, str, len) > 0; }
    bool equal(const element* str) const { return comparefunc<_cpcase>(this, str) == 0; }
    bool equal(const element* str, int len) const { return comparefunc<_cpcase>(this, str, len) == 0; }
    bool less(const element* str) const { return comparefunc<_cpcase>(this, str) < 0; }
    bool less(const element* str, int len) const { return comparefunc<_cpcase>(this, str, len) < 0; }
    bool operator < (const element* str) const { return less(str); }
    bool operator < (const _string* str) const { return less(str->c_str(), str->length()); }
    bool operator < (const _string& str) const { return less(str.c_str(), str.length()); }
    bool operator == (const element* str) const { return equal(str); }
    bool operator == (const _string* str) const { return equal(str->c_str(), str->length()); }
    bool operator == (const _string& str) const { return equal(str.c_str(), str.length()); }
    bool operator != (const element* str) const { return !equal(str); }
    bool operator != (const _string* str) const { return !equal(str->c_str(), str->length()); }
    bool operator != (const _string& str) const { return !equal(str.c_str(), str.length()); }
    bool operator > (const element* str) const { return greater(str); }
    bool operator > (const _string* str) const { return greater(str->c_str(), str->length()); }
    bool operator > (const _string& str) const { return greater(str.c_str(), str.length()); }
};

typedef _string_tool<gchar> strtool;

/*
 * string compare tool for short const string, notify that the s1 increase:
 * string symmetrical check
 */
template<class _element>
class strsymcheck
{
public:
    typedef _element mchar;
    typedef strsymcheck<mchar> mref;

public:
    template<int _len>
    static bool check(const mchar* s1, const mchar* s2) { return *s1 ++ == *s2 ++ && mref::check<_len-1>(s1, s2); }
    template<>
    static bool check<0>(const mchar* s1, const mchar* s2) { return true; }
    template<>
    static bool check<1>(const mchar* s1, const mchar* s2) { return *s1 ++ == *s2 ++; }
    template<int _len>
    static bool check_run(const mchar*& s1, const mchar* s2)
    {
        if(check<_len>(s1, s2)) {
            s1 += _len;
            return true;
        }
        return false;
    }
};

#ifdef _cstrcmp
#undef _cstrcmp
#endif
#define _cstrcmp(sr, cmp)    strsymcheck<gchar>::check<_cststrlen(cmp)>((sr), (cmp))

#define _cstrncmp(sr, len, cmp) \
    ((len >= _cststrlen(cmp)) && strsymcheck<gchar>::check<_cststrlen(cmp)>((sr), (cmp)))

#ifdef _cstrcmprun
#undef _cstrcmprun
#endif
#define _cstrcmprun(sr, cmp) strsymcheck<gchar>::check_run<_cststrlen(cmp)>((sr), (cmp))

#define mkstr(str1, str2) (((string(str1) += str2)).c_str())
#define mkstr3(str1, str2, str3) (((string(str1) += str2) += str3).c_str())
#define mkstr4(str1, str2, str3, str4) ((((string(str1) += str2) += str3) += str4).c_str())

template<int _cmpstr> inline
bool strequ(const char* str) { return *(const int*)str == _cmpstr; }
template<char _cmpch> inline
bool chequ(const char* str) { return *str == _cmpch; }
template<char _c1, char _c2>
struct chless { static const bool result = _c1 < _c2; };
template<char _c1, char _c2, bool _less = chless<_c1, _c2>::result>
struct chmin {};
template<char _c1, char _c2>
struct chmin<_c1, _c2, true> { static const char value = _c1; };
template<char _c1, char _c2>
struct chmin<_c1, _c2, false> { static const char value = _c2; };
template<char _c1, char _c2, bool _less = chless<_c1, _c2>::result>
struct chmax {};
template<char _c1, char _c2>
struct chmax<_c1, _c2, true> { static const char value = _c2; };
template<char _c1, char _c2>
struct chmax<_c1, _c2, false> { static const char value = _c1; };
template<int _len> inline
char strhit(const char* str, const char* cmp) { return *str == *cmp++ ? *str : strhit<_len + 1>(str, cmp); }
template<> inline
char strhit<1>(const char* str, const char* cmp) { return *str == *cmp ? *str : 0; }
template<> inline
char strhit<0>(const char* str, const char* cmp) { return *str; }
template<char _from, char _to>
char strhit(const char* str)
{
    return (*str >= chmin<_from, _to>::value &&
        *str <= chmax<_from, _to>::value) ?
        *str : 0;
}

/* used by the tool chain. */
inline int conv_quadstr(const char* str)
{
    int n = 0;
    if(str[0] == 0)
        return n;
    n |= str[0];
    if(str[1] == 0)
        return n;
    n |= ((uint)str[1]) << 8;
    if(str[2] == 0)
        return n;
    n |= ((uint)str[2]) << 16;
    if(str[3] == 0)
        return n;
    n |= ((uint)str[3]) << 24;
    return n;
}
inline char _1of_quadstr(int quadstr) { return (char)(quadstr & 0xff); }
inline char _2of_quadstr(int quadstr) { return (char)((quadstr & 0xff00) >> 8); }
inline char _3of_quadstr(int quadstr) { return (char)((quadstr & 0xff0000) >> 16); }
inline char _4of_quadstr(int quadstr) { return (char)((quadstr & 0xff000000) >> 24); }
inline int len_of_quadstr(int quadstr)
{
    if(!(quadstr & 0xff))
        return 0;
    if(!(quadstr & 0xff00))
        return 1;
    if(!(quadstr & 0xff0000))
        return 2;
    return (quadstr & 0xff000000) ? 4 : 3;
}
inline int utf8_unit_len(const char* str)
{
    if(str[0] == 0)
        return 0;
    if(!(str[0] & 0x80))
        return 1;
    if((str[0] & 0xe0) == 0xc0)
        return 2;
    if((str[0] & 0xf0) == 0xe0)
        return 3;
    if((str[0] & 0xf8) == 0xf0)
        return 4;
    if((str[0] & 0xfc) == 0xf8)
        return 5;
    if((str[0] & 0xfe) == 0xfc)
        return 6;
    return -1;
}
inline bool is_utf8_head(const char* str)
{
    assert(str);
    return (str[0] & 0xc0) != 0x80;
}
inline const char* to_utf8_head(const char* str)
{
    assert(str);
    for( ; str[0] && !is_utf8_head(str); str --);
    return str;
}
inline const char* next_utf8_unit(const char* str)
{
    assert(str);
    if(str[0] == 0)
        return str;
    for(str ++; !is_utf8_head(str); str ++);
    return str;
}
inline const char* prev_utf8_unit(const char* str)
{
    assert(str);
    return to_utf8_head(-- str);
}
inline bool utf8_unit_valid(const char* str)
{
    assert(str);
    return !is_utf8_head(str) ? false :
        next_utf8_unit(str) - str == utf8_unit_len(str);
}
inline int utf8_len(const char* str)
{
    assert(str);
    int len = 0;
    for( ; str[0]; str = next_utf8_unit(str), len ++);
    return len;
}

class string:
    public _string<gchar>
{
public:
    typedef std::basic_string<gchar> stdstr;
    typedef _string<gchar> superref;
    typedef superref::iterator iterator;
    typedef superref::const_iterator const_iterator;
    typedef superref::protoref protoref;

public:
    using superref::npos;
    using superref::begin;
    using superref::c_str;
    using superref::capacity;
    using superref::compare;
    using superref::copy;
    using superref::empty;
    using superref::find;
    using superref::find_first_not_of;
    using superref::find_first_of;
    using superref::find_last_not_of;
    using superref::find_last_of;
    using superref::max_size;
    using superref::operator+=;
    using superref::operator=;
    using superref::rbegin;
    using superref::reserve;
    using superref::rfind;
    using superref::substr;
    using superref::test;
    using superref::_test;
    using superref::to_lower;
    using superref::to_upper;
    using superref::to_int;
    using superref::to_int64;
    using superref::to_real;
    using superref::compare_cl;
    using superref::greater;
    using superref::less;
    using superref::equal;
    using superref::operator<;
    using superref::operator>;
    using superref::operator!=;
    using superref::operator==;
    using superref::front;
    using superref::back;
    using superref::assign;
    using superref::append;
    using superref::clear;
    using superref::erase;
    using superref::insert;
    using superref::push_back;
    using superref::replace;
    using superref::resize;
    using superref::swap;
    using superref::at;
    using superref::operator[];
    using superref::format;
    using superref::from_int;
    using superref::from_int64;
    using superref::from_real;
    using superref::length;
    using superref::pop_back;

public:
    string() {}
    string(const gchar* str): superref(str) {}
    string(const gchar* str, int len): superref(str, len) {}
    string(gchar ch, int cnt): superref(ch, cnt) {}
    string(protoref& rhs) { assign(rhs.c_str(), (int)rhs.size()); }
    int size() const { return (int)superref::size(); }
    bool operator<(const string& str) const { return less(str.c_str(), str.length()); }
    bool operator>(const string& str) const { return greater(str.c_str(), str.length()); }
    bool operator==(const string& str) const { return equal(str.c_str(), str.length()); }
    bool operator!=(const string& str) const { return !equal(str.c_str(), str.length()); }
};

inline bool is_mbcs_half(const char* str)
{
    if(!str || !(str[0] & 0x80))
        return false;
    int i = 1;
    for( ; str[i] && (str[i] & 0x80); i ++);
    return i % 2 != 0;
}

inline const gchar* get_next_char(const gchar* str, uint& c)
{
#ifdef _UNICODE
    if(!str || !str[0]) { c = 0; return 0; }
    c = (uint)str[0];
    return ++ str;
#else
    return get_mbcs_elem(str, c);
#endif
}

inline const gchar* get_next_char(const gchar* str, uint& c, const gchar* end)
{
#ifdef _UNICODE
    if(!str || str == end) { c = 0; return 0; }
    c = (uint)str[0];
    return ++ str;
#else
    return get_mbcs_elem(str, c, end);
#endif
}

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

#ifndef __GSL_VECTOR_INT_H__
#define __GSL_VECTOR_INT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_int.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  int *data;
  gsl_block_int *block;
  int owner;
}
gsl_vector_int;

typedef struct
{
  gsl_vector_int vector;
} _gsl_vector_int_view;

typedef _gsl_vector_int_view gsl_vector_int_view;

typedef struct
{
  gsl_vector_int vector;
} _gsl_vector_int_const_view;

typedef const _gsl_vector_int_const_view gsl_vector_int_const_view;

/* Allocation */

GSL_EXPORT gsl_vector_int *gsl_vector_int_alloc (const size_t n);
GSL_EXPORT gsl_vector_int *gsl_vector_int_calloc (const size_t n);

GSL_EXPORT gsl_vector_int *gsl_vector_int_alloc_from_block (gsl_block_int * b,
                                                                const size_t offset,
                                                                const size_t n,
                                                                const size_t stride);

GSL_EXPORT gsl_vector_int *gsl_vector_int_alloc_from_vector (gsl_vector_int * v,
                                                                 const size_t offset,
                                                                 const size_t n,
                                                                 const size_t stride);

GSL_EXPORT void gsl_vector_int_free (gsl_vector_int * v);

/* Views */

GSL_EXPORT
_gsl_vector_int_view
gsl_vector_int_view_array (int *v, size_t n);

GSL_EXPORT
_gsl_vector_int_view
gsl_vector_int_view_array_with_stride (int *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_vector_int_const_view_array (const int *v, size_t n);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_vector_int_const_view_array_with_stride (const int *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_int_view
gsl_vector_int_subvector (gsl_vector_int *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_int_view
gsl_vector_int_subvector_with_stride (gsl_vector_int *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_vector_int_const_subvector (const gsl_vector_int *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_int_const_view
gsl_vector_int_const_subvector_with_stride (const gsl_vector_int *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT int gsl_vector_int_get (const gsl_vector_int * v, const size_t i);
GSL_EXPORT void gsl_vector_int_set (gsl_vector_int * v, const size_t i, int x);

GSL_EXPORT int *gsl_vector_int_ptr (gsl_vector_int * v, const size_t i);
GSL_EXPORT const int *gsl_vector_int_const_ptr (const gsl_vector_int * v, const size_t i);

GSL_EXPORT void gsl_vector_int_set_zero (gsl_vector_int * v);
GSL_EXPORT void gsl_vector_int_set_all (gsl_vector_int * v, int x);
GSL_EXPORT int gsl_vector_int_set_basis (gsl_vector_int * v, size_t i);

GSL_EXPORT int gsl_vector_int_fread (FILE * stream, gsl_vector_int * v);
GSL_EXPORT int gsl_vector_int_fwrite (FILE * stream, const gsl_vector_int * v);
GSL_EXPORT int gsl_vector_int_fscanf (FILE * stream, gsl_vector_int * v);
GSL_EXPORT int gsl_vector_int_fprintf (FILE * stream, const gsl_vector_int * v,
                                         const char *format);

GSL_EXPORT int gsl_vector_int_memcpy (gsl_vector_int * dest, const gsl_vector_int * src);

GSL_EXPORT int gsl_vector_int_reverse (gsl_vector_int * v);

GSL_EXPORT int gsl_vector_int_swap (gsl_vector_int * v, gsl_vector_int * w);
GSL_EXPORT int gsl_vector_int_swap_elements (gsl_vector_int * v, const size_t i, const size_t j);

GSL_EXPORT int gsl_vector_int_max (const gsl_vector_int * v);
GSL_EXPORT int gsl_vector_int_min (const gsl_vector_int * v);
GSL_EXPORT void gsl_vector_int_minmax (const gsl_vector_int * v, int * min_out, int * max_out);

GSL_EXPORT size_t gsl_vector_int_max_index (const gsl_vector_int * v);
GSL_EXPORT size_t gsl_vector_int_min_index (const gsl_vector_int * v);
GSL_EXPORT void gsl_vector_int_minmax_index (const gsl_vector_int * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_int_add (gsl_vector_int * a, const gsl_vector_int * b);
GSL_EXPORT int gsl_vector_int_sub (gsl_vector_int * a, const gsl_vector_int * b);
GSL_EXPORT int gsl_vector_int_mul (gsl_vector_int * a, const gsl_vector_int * b);
GSL_EXPORT int gsl_vector_int_div (gsl_vector_int * a, const gsl_vector_int * b);
GSL_EXPORT int gsl_vector_int_scale (gsl_vector_int * a, const double x);
GSL_EXPORT int gsl_vector_int_add_constant (gsl_vector_int * a, const double x);

GSL_EXPORT int gsl_vector_int_isnull (const gsl_vector_int * v);

#ifdef HAVE_INLINE

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

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

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

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

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_INT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws