text:
//CELL operation: does output side of Jordan RNN layer.
//The inputs here are from the IN stage, so reduced to N driving input time-series in X.

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

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

int jordan_s (float *Y, const float *X, const float *U, float *Y1, const float *W, const float *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int jordan_d (double *Y, const double *X, const double *U, double *Y1, const double *W, const double *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int jordan_s (float *Y, const float *X, const float *U, float *Y1, const float *W, const float *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;

    float *H;
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in jordan_s: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        H[0] = 1.0f / (1.0f+expf(-X[0]-U[0]*Y1[0]));
        Y[0] = 1.0f / (1.0f+expf(-B[0]-W[0]*H[0]));
        for (size_t t=1; t<T; ++t)
        {
            H[0] = 1.0f / (1.0f+expf(-X[t]-U[0]*Y[t-1]));
            Y[t] = 1.0f / (1.0f+expf(-B[0]-W[0]*H[0]));
        }
    }
    else
    {
        if (dim==0u)
        {
            if (iscolmajor)
            {
                for (size_t n=0u; n<N; ++n) { cblas_scopy((int)T,&B[n],0,&Y[n],(int)N); }
                cblas_scopy((int)N,&X[0],1,H,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],1);
                for (size_t n=0u; n<N; ++n) { Y[n] = 1.0f/(1.0f+expf(-Y[n])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t*N],1,H,1);
                    cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,&Y[(t-1)*N],1,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                    cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t*N],1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0f/(1.0f+expf(-Y[t*N+n])); }
                }
            }
            else
            {
                for (size_t n=0u; n<N; ++n) { cblas_scopy((int)T,&B[n],0,&Y[n*T],1); }
                cblas_scopy((int)N,&X[0],(int)T,H,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],(int)T);
                for (size_t n=0u; n<N; ++n) { Y[n*T] = 1.0f/(1.0f+expf(-Y[n*T])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t],(int)T,H,1);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,&Y[t-1],(int)T,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t],(int)T);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0f/(1.0f+expf(-Y[t+n*T])); }
                }
            }
        }
        else if (dim==1u)
        {
            if (iscolmajor)
            {
                for (size_t n=0u; n<N; ++n) { cblas_scopy((int)T,&B[n],0,&Y[n*T],1); }
                cblas_scopy((int)N,&X[0],(int)T,H,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],(int)T);
                for (size_t n=0u; n<N; ++n) { Y[n*T] = 1.0f/(1.0f+expf(-Y[n*T])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t],(int)T,H,1);
                    cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,&Y[t-1],(int)T,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                    cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t],(int)T);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0f/(1.0f+expf(-Y[t+n*T])); }
                }
            }
            else
            {
                for (size_t n=0u; n<N; ++n) { cblas_scopy((int)T,&B[n],0,&Y[n],(int)N); }
                cblas_scopy((int)N,&X[0],1,H,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],1);
                for (size_t n=0u; n<N; ++n) { Y[n] = 1.0f/(1.0f+expf(-Y[n])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_scopy((int)N,&X[t*N],1,H,1);
                    cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,&Y[(t-1)*N],1,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0f/(1.0f+expf(-H[n])); }
                    cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t*N],1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0f/(1.0f+expf(-Y[t*N+n])); }
                }
            }
        }
        else
        {
            fprintf(stderr,"error in jordan_s: dim must be 0 or 1.\n"); return 1;
        }
    }

    return 0;
}

int jordan_d (double *Y, const double *X, const double *U, double *Y1, const double *W, const double *B, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;

    double *H;
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in jordan_d: problem with malloc. "); perror("malloc"); return 1; }

    if (N==1u)
    {
        H[0] = 1.0 / (1.0+exp(-X[0]-U[0]*Y1[0]));
        Y[0] = 1.0 / (1.0+exp(-B[0]-W[0]*H[0]));
        for (size_t t=1; t<T; ++t)
        {
            H[0] = 1.0 / (1.0+exp(-X[t]-U[0]*Y[t-1]));
            Y[t] = 1.0 / (1.0+exp(-B[0]-W[0]*H[0]));
        }
    }
    else
    {
        if (dim==0u)
        {
            if (iscolmajor)
            {
                for (size_t n=0u; n<N; ++n) { cblas_dcopy((int)T,&B[n],0,&Y[n],(int)N); }
                cblas_dcopy((int)N,&X[0],1,H,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],1);
                for (size_t n=0u; n<N; ++n) { Y[n] = 1.0/(1.0+exp(-Y[n])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t*N],1,H,1);
                    cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,&Y[(t-1)*N],1,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                    cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t*N],1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0/(1.0+exp(-Y[t*N+n])); }
                }
            }
            else
            {
                for (size_t n=0u; n<N; ++n) { cblas_dcopy((int)T,&B[n],0,&Y[n*T],1); }
                cblas_dcopy((int)N,&X[0],(int)T,H,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],(int)T);
                for (size_t n=0u; n<N; ++n) { Y[n*T] = 1.0/(1.0+exp(-Y[n*T])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t],(int)T,H,1);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,&Y[t-1],(int)T,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t],(int)T);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0/(1.0+exp(-Y[t+n*T])); }
                }
            }
        }
        else if (dim==1u)
        {
            if (iscolmajor)
            {
                for (size_t n=0u; n<N; ++n) { cblas_dcopy((int)T,&B[n],0,&Y[n*T],1); }
                cblas_dcopy((int)N,&X[0],(int)T,H,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],(int)T);
                for (size_t n=0u; n<N; ++n) { Y[n*T] = 1.0/(1.0+exp(-Y[n*T])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t],(int)T,H,1);
                    cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,&Y[t-1],(int)T,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                    cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t],(int)T);
                    for (size_t n=0u; n<N; ++n) { Y[t+n*T] = 1.0/(1.0+exp(-Y[t+n*T])); }
                }
            }
            else
            {
                for (size_t n=0u; n<N; ++n) { cblas_dcopy((int)T,&B[n],0,&Y[n],(int)N); }
                cblas_dcopy((int)N,&X[0],1,H,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,Y1,1,o,H,1);
                for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[0],1);
                for (size_t n=0u; n<N; ++n) { Y[n] = 1.0/(1.0+exp(-Y[n])); }
                for (size_t t=1; t<T; ++t)
                {
                    cblas_dcopy((int)N,&X[t*N],1,H,1);
                    cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,&Y[(t-1)*N],1,o,H,1);
                    for (size_t n=0u; n<N; ++n) { H[n] = 1.0/(1.0+exp(-H[n])); }
                    cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,W,(int)N,H,1,o,&Y[t*N],1);
                    for (size_t n=0u; n<N; ++n) { Y[t*N+n] = 1.0/(1.0+exp(-Y[t*N+n])); }
                }
            }
        }
        else
        {
            fprintf(stderr,"error in jordan_d: dim must be 0 or 1.\n"); return 1;
        }
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
/* Copyright 2015 The TensorFlow Authors. All Rights Reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
==============================================================================*/
// Portions Copyright (c) Microsoft Corporation

#pragma once

#include <functional>
#include <memory>
#include <string>
#include <unordered_map>
#include <vector>
#include <gsl/pointers>

#include "core/common/common.h"
#include "core/common/callback.h"
#include "core/platform/env_time.h"

#ifndef _WIN32
#include <sys/types.h>
#include <unistd.h>
#endif

namespace onnxruntime {

#ifdef _WIN32
using PIDType = unsigned long;
#else
using PIDType = pid_t;
#endif

/// \brief An interface used by the onnxruntime implementation to
/// access operating system functionality like the filesystem etc.
///
/// Callers may wish to provide a custom Env object to get fine grain
/// control.
///
/// All Env implementations are safe for concurrent access from
/// multiple threads without any external synchronization.
class Env {
 public:
  virtual ~Env() = default;

  /// \brief Returns a default environment suitable for the current operating
  /// system.
  ///
  /// Sophisticated users may wish to provide their own Env
  /// implementation instead of relying on this default environment.
  ///
  /// The result of Default() belongs to this library and must never be deleted.
  static const Env& Default();

  virtual int GetNumCpuCores() const = 0;

  /// \brief Returns the number of micro-seconds since the Unix epoch.
  virtual uint64_t NowMicros() const { return env_time_->NowMicros(); }

  /// \brief Returns the number of seconds since the Unix epoch.
  virtual uint64_t NowSeconds() const { return env_time_->NowSeconds(); }

  /// Sleeps/delays the thread for the prescribed number of micro-seconds.
  /// On Windows, it's the min time to sleep, not the actual one.
  virtual void SleepForMicroseconds(int64_t micros) const = 0;

#ifndef _WIN32
  /**
   *
   * \param file_path file_path must point to a regular file, which can't be a pipe/socket/...
   * \param[out] p  allocated buffer with the file data
   * \param[in] offset file offset. If offset>0, then len must also be >0.
   * \param[in, out] len length to read(or has read). If len==0, read the whole file.
   * @return
   */
  virtual common::Status ReadFileAsString(const char* file_path, off_t offset, void*& p, size_t& len,
      OrtCallback& deleter) const = 0;
#else
  virtual common::Status ReadFileAsString(const wchar_t* file_path, int64_t offset, void*& p, size_t& len,
                                          OrtCallback& deleter) const = 0;
#endif

#ifdef _WIN32
  //Mainly for use with protobuf library
  virtual common::Status FileOpenRd(const std::wstring& path, /*out*/ int& fd) const = 0;
  //Mainly for use with protobuf library
  virtual common::Status FileOpenWr(const std::wstring& path, /*out*/ int& fd) const = 0;
#endif
  //Mainly for use with protobuf library
  virtual common::Status FileOpenRd(const std::string& path, /*out*/ int& fd) const = 0;
  //Mainly for use with protobuf library
  virtual common::Status FileOpenWr(const std::string& path, /*out*/ int& fd) const = 0;
  //Mainly for use with protobuf library
  virtual common::Status FileClose(int fd) const = 0;
  //This functions is always successful. It can't fail.
  virtual PIDType GetSelfPid() const = 0;

  // \brief Load a dynamic library.
  //
  // Pass "library_filename" to a platform-specific mechanism for dynamically
  // loading a library.  The rules for determining the exact location of the
  // library are platform-specific and are not documented here.
  //
  // On success, returns a handle to the library in "*handle" and returns
  // OK from the function.
  // Otherwise returns nullptr in "*handle" and an error status from the
  // function.
  virtual common::Status LoadDynamicLibrary(const std::string& library_filename, void** handle) const = 0;

  virtual common::Status UnloadDynamicLibrary(void* handle) const = 0;

  // \brief Get a pointer to a symbol from a dynamic library.
  //
  // "handle" should be a pointer returned from a previous call to LoadDynamicLibrary.
  // On success, store a pointer to the located symbol in "*symbol" and return
  // OK from the function. Otherwise, returns nullptr in "*symbol" and an error
  // status from the function.
  virtual common::Status GetSymbolFromLibrary(void* handle, const std::string& symbol_name, void** symbol) const = 0;

  // \brief build the name of dynamic library.
  //
  // "name" should be name of the library.
  // "version" should be the version of the library or NULL
  // returns the name that LoadDynamicLibrary() can use
  virtual std::string FormatLibraryFileName(const std::string& name, const std::string& version) const = 0;

 protected:
  Env();

 private:
  ORT_DISALLOW_COPY_ASSIGNMENT_AND_MOVE(Env);
  EnvTime* env_time_ = EnvTime::Default();
};

}  // namespace onnxruntime

ext: h

max_stars_repo_name: dahburj/onnxruntime

max_issues_repo_name: dahburj/onnxruntime

max_forks_repo_name: dahburj/onnxruntime


text:
#pragma once

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

#include "defines.h"

// Pass matrices a, b
// Matrix c is already allocated and is the correct size
// c = [a b]'
int gsl_matrix_vstack(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b);

// Multiply matrices c = a*b
// Use gsl_blas_dgemm() instead:
//  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C);
int gsl_matrix_mul(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b);

// Helper function for matrix multiplication to multiply specified row by
// specified column
// Number of columns in a must be equal to number of rows in b
double gsl_matrix_rowcol_mul(const gsl_matrix *a, const gsl_matrix *b, int r, int c);

// Pass vectors a, b
// Vector c is already allocated and is the correct size
// c = [a b]'
int gsl_vector_vstack(gsl_vector *c, const gsl_vector *a, const gsl_vector *b);

double gsl_vector_infnorm(const gsl_vector *v);

void gsl_vector_print(const gsl_vector *v, char *name);

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
// An interface to GSL - hides all the gsl calls behind easier to use functions
// Jason Sanders
#ifndef GSLINTERFACE_H
#define GSLINTERFACE_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_permute.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_ellint.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_siman.h>
#include <stdlib.h>
#include <iostream>
//=================================================================================================
// RANDOM NUMBERS //
// random number generators - rand_uniform returns random numbers uniformly distributed in the
// interval [0,1], rand_gaussian returns gaussian distributed random numbers with sd = sigma

class rand_base{
	private:
		const gsl_rng_type * TYPE;
		unsigned long int seed;
    public:
       	gsl_rng * r;
     	rand_base(unsigned long int s){
     	// construct random number generator with seed s
     		seed = s;
     		gsl_rng_env_setup();
     	    TYPE = gsl_rng_default;
       		r = gsl_rng_alloc (TYPE);
       		gsl_rng_set(r, seed);
       		}
       	~rand_base(){gsl_rng_free (r);}
       	void reseed(unsigned long int newseed){
       	// give new seed
       		seed = newseed;
       		gsl_rng_set(r,newseed);
       		}
};

class rand_uniform:public rand_base{
	public:
		rand_uniform(unsigned long int SEED=0):rand_base(SEED){}
		~rand_uniform(){}
		// return uniformly distributed random numbers
		double nextnumber(){return gsl_rng_uniform (r);}
};

class rand_gaussian:public rand_base{
	public:
		double sigma;
		rand_gaussian(double s, unsigned long int SEED=0):rand_base(SEED){sigma = s;}
		~rand_gaussian(){}
		// return gaussian distributed random numbers
		double nextnumber(){return gsl_ran_gaussian (r,sigma);}
		void newsigma(double newsigma){sigma=newsigma;}
};

//=================================================================================================
// ROOT FINDING	  //
// finds root by Brent's method. Constructor initialises function and tolerances, findroot finds
// root in given interval. Function must be of form double(*func)(double,void*)

class root_find{
	private:
		int status;
       	const gsl_root_fsolver_type *T; gsl_root_fsolver *s;
       	gsl_function F;	double xlo,xhi,eps,root,tol;
       	int iter, max_iter;
    public:
    	root_find(double tol,int max_iter)
    		: tol(tol), max_iter(max_iter){
    		iter=0;eps=1.e-8;
       		T = gsl_root_fsolver_brent;
       		s = gsl_root_fsolver_alloc (T);
		}
		~root_find(){gsl_root_fsolver_free (s);}

    	void bracket(){
    		// crude bracketing routine. Expands interval till root found
			if(xhi<xlo){double xt = xlo; xlo=xhi; xhi=xt;}
			double diff = (xhi-xlo)/2.;
			while(GSL_FN_EVAL(&F,xlo)*GSL_FN_EVAL(&F,xhi)>0.){
				xlo-=diff;xhi+=diff;
			}
		}

     	double findroot(double(*func)(double,void *),double xlo1,double xhi1, void *p=NULL){
     		// finds root in interval [xlo1,xhi1]
    		F.function = func;F.params=p;xhi = xhi1; xlo = xlo1;
     	    if(GSL_FN_EVAL(&F,xlo)*GSL_FN_EVAL(&F,xhi)>0.){
     			bracket();
     		}
       		gsl_root_fsolver_set (s, &F, xlo, xhi);
			do{
				iter++;
				status = gsl_root_fsolver_iterate (s);
				root = gsl_root_fsolver_root (s);
				xlo = gsl_root_fsolver_x_lower (s);
				xhi = gsl_root_fsolver_x_upper (s);
				status = gsl_root_test_interval (xlo, xhi, tol, eps);
			  }
			while (status == GSL_CONTINUE && iter < max_iter);
       		return root;
		}

};

//=================================================================================================
// INTEGRATION //
// Simple 1d numerical integration using adaptive Gauss-Kronrod which can deal with singularities
// constructor takes function and tolerances and integrate integrates over specified region.
// integrand function must be of the form double(*func)(double,void*)
class integrator{
	private:
		gsl_integration_workspace *w;
		void *p;
		double result,err,eps;
		gsl_function F;
		size_t neval;
	public:
		integrator(double eps): eps(eps){
			w= gsl_integration_workspace_alloc (1000);
			F.params = &p;
       		}
		~integrator(){gsl_integration_workspace_free (w);}
		double integrate(double(*func)(double,void*),double xa, double xb){
		    F.function = func;
			gsl_integration_qags (&F, xa, xb, 0, eps, 1000,w, &result, &err);
			//gsl_integration_qng(&F, xa, xb, 0, eps, &result, &err, &neval);
			return result;
			}
		double error(){return err;}
};

inline double integrate(double(*func)(double,void*),double xa, double xb, double eps){
	double result,err; size_t neval;void *p;gsl_function F;F.function = func;F.params = &p;
	gsl_integration_qng(&F, xa, xb, 0, eps, &result, &err, &neval);
	return result;
}

class MCintegrator{
	private:
		gsl_monte_vegas_state *s;
		const gsl_rng_type *T;
 		gsl_rng *r;
 		size_t dim;
 	public:
 		MCintegrator(size_t Dim){
 			dim=Dim;
 			gsl_rng_env_setup ();
  			T = gsl_rng_default;
  			r = gsl_rng_alloc (T);
 			s=gsl_monte_vegas_alloc(Dim);
 		}
 		~MCintegrator(){
 			gsl_monte_vegas_free(s);
 			gsl_rng_free(r);
 		}
 		double integrate(double(*func)(double*,size_t,void*),double *xlow, double *xhigh,
 		size_t calls, double *err, int burnin=10000){
 			gsl_monte_function G = { func, dim, 0 }; double res;
 			if(burnin)gsl_monte_vegas_integrate(&G,xlow,xhigh,dim,burnin,r,s,&res,err);
 			gsl_monte_vegas_integrate(&G,xlow,xhigh,dim,calls,r,s,&res,err);
 			return res;
 		}
};

//=================================================================================================
// 1D INTERPOLATION //
// Interpolation using cubic splines

class interpolator{
	private:
		gsl_interp_accel *acc;
		gsl_spline *spline;
	public:
		interpolator(double *x, double *y, int n){
			acc = gsl_interp_accel_alloc();
			spline = gsl_spline_alloc(gsl_interp_cspline,n);
			gsl_spline_init (spline, x, y, n);
		}
		~interpolator(){
			gsl_spline_free (spline);
         	gsl_interp_accel_free (acc);
        }
        double interpolate(double xi){
        	return gsl_spline_eval (spline, xi, acc);
        }
        double derivative(double xi){
        	return gsl_spline_eval_deriv(spline, xi, acc);
        }
        void new_arrays(double *x, double *y,int n){
        	spline = gsl_spline_alloc(gsl_interp_cspline,n);
        	gsl_spline_init (spline, x, y, n);
        }
};

//=================================================================================================
// SORTING //
// sorting algorithm
// sort2 sorts first argument and then applies the sorted permutation to second list
class sorter{
	private:
		const gsl_rng_type * T;
       	gsl_rng * r;
    public:
    	sorter(){
    		gsl_rng_env_setup();
            T = gsl_rng_default;
     		r = gsl_rng_alloc (T);
     		}
     	~sorter(){gsl_rng_free (r);}
     	void sort(double *data, int n){
     		gsl_sort(data,1,n);
     	}
     	void sort2(double *data, int n, double *data2){
     		size_t p[n];
     		gsl_sort_index(p,data,1,n);
     		gsl_permute(p,data2,1,n);
     	}

};

//=================================================================================================
// ODE SOLVER //
// Simple ODE integrator using Runge-Kutta Dormand-Prince 8 adaptive stepping
// dy_i/dt = f_i(t) where int (*f)(double t, const double y, double f, void *params)

class ode{
	private:
		const gsl_odeiv2_step_type * T;
		gsl_odeiv2_step * s;
		gsl_odeiv2_control * c;
		gsl_odeiv2_evolve * e;
  		gsl_odeiv2_system sys;
		int N;double h,eps;
		double direction;
	public:
		ode(int (*derivs)(double,const double *,double*,void*),int N,double eps, void *params=NULL):N(N),eps(eps){
			sys.function=derivs;sys.jacobian=NULL;sys.dimension=N;sys.params=params;

			T = gsl_odeiv2_step_rk8pd;s = gsl_odeiv2_step_alloc (T, N);
			c = gsl_odeiv2_control_y_new (eps, 0.0);e = gsl_odeiv2_evolve_alloc (N);
			gsl_ieee_env_setup();
		}
		~ode(){
			gsl_odeiv2_evolve_free(e);
  			gsl_odeiv2_control_free(c);
  			gsl_odeiv2_step_free(s);
  		}
		void step(double tstart, double tfinish, double *y, double step){
			h = step;
			direction = (h>0?1.:-1.);
			while ((tfinish-tstart)*direction>0){
      			int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &tstart, tfinish, &h, y);
      			if (status != GSL_SUCCESS)break;
    		}
		}
};

//=================================================================================================
// MINIMISER //
// finds a minimum of a function of the form double(*func)(const gsl_vector *v, void *params)
// using a downhill simplex algorithm. Setup minimiser with initial guesses and required tolerance
// with constructor and then minimise with minimise().

class minimiser{
	private:
		const gsl_multimin_fminimizer_type *T; ;
		gsl_multimin_fminimizer *s;
		gsl_vector *ss, *x;
		gsl_multimin_function minex_func;
		size_t iter; int status,N_params; double size;
		double eps;
	public:
		minimiser(double(*func)(const gsl_vector *v, void *params),double *parameters,
		 int N, double *sizes, double eps):N_params(N), eps(eps){

			T = gsl_multimin_fminimizer_nmsimplex2rand;
			ss = gsl_vector_alloc (N_params);x = gsl_vector_alloc (N_params);
			for(int i=0;i<N_params;i++){
				gsl_vector_set (x, i, parameters[i]);gsl_vector_set(ss,i,sizes[i]);}

			minex_func.n = N; minex_func.f = func; minex_func.params = NULL;
			s = gsl_multimin_fminimizer_alloc (T, N_params);
			gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
			status = 0; iter = 0;

		}

		~minimiser(){
			gsl_vector_free(x);
			gsl_vector_free(ss);
			gsl_multimin_fminimizer_free (s);
		}

		double minimise(double *results,unsigned int maxiter,bool vocal){
			do
			  {
				iter++; status = gsl_multimin_fminimizer_iterate(s);
				if(status)break;
				size = gsl_multimin_fminimizer_size (s);
				status = gsl_multimin_test_size (size, eps);
				if(vocal){	std::cout<<iter<<" ";
							for(int i=0; i<N_params;i++)std::cout<<gsl_vector_get(s->x,i)<<" ";
							std::cout<<s->fval<<" "<<size<<std::endl;
							}
			}
			while (status == GSL_CONTINUE && iter < maxiter);
			for(int i=0;i<N_params;i++){results[i] = gsl_vector_get(s->x,i);}
			return s->fval;
		}
};

class minimiser1D{
	private:
		const gsl_min_fminimizer_type *T; ;
		gsl_min_fminimizer *s;
		size_t iter; int status; double size;
		double m, a, b, eps;
	public:
		minimiser1D(double(*func)(double, void *params), double m, double a, double b, double eps)
			:m(m), a(a), b(b), eps(eps){

			gsl_function F;F.function = func;F.params = 0;
			T = gsl_min_fminimizer_brent;
			s = gsl_min_fminimizer_alloc (T);
			gsl_min_fminimizer_set (s, &F, m, a, b);
			status = 0; iter = 0;
		}
		~minimiser1D(){
			gsl_min_fminimizer_free (s);
		}
		double minimise(unsigned int maxiter){
			do
			  {
				iter++;
				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);
				status = gsl_min_test_interval (a, b, eps, 0.0);
			}
			while (status == GSL_CONTINUE && iter < maxiter);
			return m;
		}
};
/*
double Distance(void *xp, void *yp){
       double x = *((double *) xp);
       double y = *((double *) yp);
       return fabs(x - y);
}

void Step(const gsl_rng * r, void *xp, double step_size){
    double old_x = *((double *) xp);
    double new_x;

    double u = gsl_rng_uniform(r);
    new_x = u * 2 * step_size - step_size + old_x;

    memcpy(xp, &new_x, sizeof(new_x));
}

void Print(void *xp){
    printf ("%12g", *((double *) xp));
}

class sim_anneal{
	private:
		const gsl_rng_type * T;
    	gsl_rng * r;
    	int N_TRIES, ITER_FIXED_T;
    	double STEP_SIZE, K, T_INITIAL, MU_T, T_MIN;
    	gsl_siman_params_t params;
    public:
    	sim_anneal(int N_TRIES, int ITER_FIXED_T, double STEP_SIZE):
    	N_TRIES(N_TRIES), ITER_FIXED_T(ITER_FIXED_T),STEP_SIZE(STEP_SIZE){
    		gsl_rng_env_setup();
            T = gsl_rng_default;
     	  	r = gsl_rng_alloc(T);
     	  	//params[0]=N_TRIES;params[1]=ITER_FIZED_T;params[2]=STEP_SIZE;
     	  	//K=1.; params[3]=K; T_INITIAL=0.008; params[4]=T_INITIAL;
     	  	//MU_T=1.003; params[5]=MU_T; T_MIN=2.0e-6; params[6]=T_MIN;
     	  	gsl_siman_params_t params
       = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
          K, T_INITIAL, MU_T, T_MIN};
    	}
    	~sim_anneal(){
    		gsl_rng_free(r);
    	}
    	double minimise(double(*func)(void *xp), double x){
    		double x_initial=x;
    		gsl_siman_solve(r, &x_initial, &func, Step, Distance, Print,
                       		NULL, NULL, NULL,
                       		sizeof(double), params);
    		return x_initial;
    	}
};*/

//=================================================================================================
// SPECIAL FUNCTIONS //
inline double erf(double x){return gsl_sf_erf (x);}
inline double erfc(double x){return gsl_sf_erfc (x);}
inline double besselI(double x, int n){return gsl_sf_bessel_In (n,x);}
inline double besselJ(double x, int n){return gsl_sf_bessel_Jn (n,x);}
inline double gamma(double x){return gsl_sf_gamma (x);}
inline double ellint_first(double phi, double k){ return gsl_sf_ellint_F(phi,k,1e-15);}
// F(\phi,k) = \int_0^\phi \d t \, \frac{1}{\sqrt{1-k^2\sin^2 t}}
inline double ellint_second(double phi, double k){ return gsl_sf_ellint_E(phi,k,1e-15);}
// E(\phi,k) = \int_0^\phi \d t \, \sqrt{1-k^2\sin^2 t}
inline double ellint_third(double phi, double k, double n){ return gsl_sf_ellint_P(phi,k,n,1e-15);}
// \Pi(\phi,k,n) = \int_0^\phi \d t \, \frac{1}{(1+n\sin^2 t)\sqrt{1-k^2\sin^2 t}}
#endif

max_stars_repo_name: jlsanders/genfunc

max_issues_repo_name: jlsanders/genfunc

max_forks_repo_name: jlsanders/genfunc


text:
/*
 *     Copyright 2021-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json.hpp>

// Used by most unit test files.
#include <folly/portability/GTest.h>

// Used throughout the codebase
#include <folly/SharedMutex.h>
#include <folly/Synchronized.h>

// MB-46844:
// Included by collections/vbucket_manifest.h, which in turn included
// by 50+ other files.
// Consider changing collections/vbucket_manifest.h to use pimpl for
// Manifest::map which would avoid the need to include F14Map.h.
#include <folly/container/F14Map.h>

#include <map>
#include <string>


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

#include <nlopt.h>

typedef struct {
     int N;
     double *x, *y; /* length N; */
} lorentzdata;

static double sqr(double x) { return x * x; }

static int count = 0;

static double lorentzerr(int n, const double *p, double *grad, void *data)
{
     lorentzdata *d = (lorentzdata *) data;
     int N = d->N;
     const double *xs = d->x;
     const double *ys = d->y;
     double val = 0;
     int i, j;

     for (i = 0; i < N; ++i) {
	  double x = xs[i], y = ys[i];
	  double lorsum = 0;

	  for (j = 0; j < n; j += 3) {
	       double A = p[j + 0];
	       double w = p[j + 1];
	       double G = p[j + 2];
	       double lor = A / (sqr(x - w) + G*G);

	       lorsum += lor;
	  }

	  val += sqr(y - lorsum);

	  if (grad)
	       for (j = 0; j < n; j += 3) {
		    double A = p[j + 0];
		    double w = p[j + 1];
		    double G = p[j + 2];
		    double deninv =  1.0 / (sqr(x - w) + G*G);

		    grad[j + 0] += -2 * (y - lorsum) * deninv;
		    grad[j + 1] += 4*A * (w - x) * (y - lorsum) * sqr(deninv);
		    grad[j + 2] += 4*A * G * (y - lorsum) * sqr(deninv);
	       }
     }
     ++count;
     // printf("%d: f(%g,%g,%g) = %g\n", count, p[0],p[1],p[2], val);
     return val;
}

extern double nlopt_urand(double a, double b);

int main(void)
{
     lorentzdata d;
     int i;
     double A = 1, w = 0, G = 1, noise=0.01;
     double lb[3] = {-HUGE_VAL,-HUGE_VAL,0};
     double ub[3] = {HUGE_VAL,HUGE_VAL,HUGE_VAL};
     double p[3] = {0,1,2}, minf;

     nlopt_srand_time();

     d.N = 200;
     d.x = (double *) malloc(sizeof(double) * d.N * 2);
     d.y = d.x + d.N;
     for (i = 0; i < d.N; ++i) {
	  d.x[i] = nlopt_urand(-0.5, 0.5) * 8*G + w;
	  d.y[i] = 2*noise*nlopt_urand(-0.5,0.5) + A / (sqr(d.x[i]-w) + G*G);
     }

     nlopt_minimize(NLOPT_LN_NEWUOA_BOUND, 3, lorentzerr, &d,
		    lb, ub, p, &minf,
		    -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);

     printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);

     count = 0;
     nlopt_minimize(NLOPT_LN_COBYLA, 3, lorentzerr, &d,
		    lb, ub, p, &minf,
		    -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);

     printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);

     count = 0;
     nlopt_minimize(NLOPT_LN_NELDERMEAD, 3, lorentzerr, &d,
		    lb, ub, p, &minf,
		    -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);

     printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);

     count = 0;
     nlopt_minimize(NLOPT_LN_SBPLX, 3, lorentzerr, &d,
		    lb, ub, p, &minf,
		    -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);

     printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);

     return 0;
}

max_stars_repo_name: daemonslayer/Notebook

max_issues_repo_name: daemonslayer/Notebook

max_forks_repo_name: daemonslayer/Notebook


text:
// Test program for comparison with brighten_float_image.

#include <assert.h>
#include <stdlib.h>

#include <glib.h>
#include <gsl/gsl_blas.h>

#include <asf_meta.h>

int
main (int argc, char **argv)
{
  assert (argc == 3);

  // Image to be brightened.
  GString *input_image_basename = g_string_new (argv[1]);

  // Basename to use for brightened version of image.
  GString *output_image_basename = g_string_new (argv[2]);

  // Form names for the input image data and metadata.
  GString *input_metadata_file = g_string_new (input_image_basename->str);
  g_string_append (input_metadata_file, ".meta");
  GString *input_data_file = g_string_new (input_image_basename->str);
  g_string_append (input_data_file, ".img");

  // Form names for the output image data and metadata.
  GString *output_metadata_file = g_string_new (output_image_basename->str);
  g_string_append (output_metadata_file, ".meta");
  GString *output_data_file = g_string_new (output_image_basename->str);
  g_string_append (output_data_file, ".img");

  // Read the metadata for the input image.
  meta_parameters *imd = meta_read (input_image_basename->str);

  // Shorthand for the dimension of the input image.
  size_t ixs = imd->general->sample_count;
  size_t iys = imd->general->line_count;

  // Set up gsl matrix for the input image.
  gsl_matrix_float *id = gsl_matrix_float_alloc (iys, ixs);

  // Set up gsl matrix for the output image.
  gsl_matrix_float *od = gsl_matrix_float_alloc (iys, ixs);

  // Do the brightening.
  size_t ii;
  // For each line.
  for ( ii = 0 ; ii < iys ; ii++ ) {
    size_t jj;
    // For each pixel of line.
    for ( jj = 0 ; jj < ixs ; jj++ ) {
      // Brighten pixel by factor of two and store in output image.
      float pixel_value = gsl_matrix_float_get (id, ii, jj);
      gsl_matrix_float_set (od, ii, jj, pixel_value * 2);
    }
  }

  // Store the output image data.
  FILE *od_stream = fopen (output_data_file->str, "wb");
  g_assert (od_stream != NULL);
  int return_code = gsl_matrix_float_fwrite (od_stream, od);
  g_assert (return_code == GSL_SUCCESS);
  return_code = fclose (od_stream);
  g_assert (return_code == 0);

  // Copy the input metadata to the output metadata.
  GString *system_command = g_string_new ("cp ");
  g_string_append_printf (system_command, "%s %s", input_metadata_file->str,
			  output_metadata_file->str);
  return_code = system (system_command->str);
  g_assert (return_code == 0);

  exit (EXIT_SUCCESS);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* matrices.c */

// 8.4.13 Example programs for matrices of gsl-ref.pdf
// GNU GSL GNU Scientific Library reference book

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

int main (void) {
  int i, j;
  gsl_matrix * m = gsl_matrix_alloc (10, 3);

  for (i = 0; i < 10; i++)
    for (j = 0; j < 3; j++)
      gsl_matrix_set (m, i, j, 0.23 + 100*i + j);

  for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */
    for (j = 0; j < 3; j++)
      printf("m(%d,%d) = %g\n", i, j,
	     gsl_matrix_get (m, i, j));

  gsl_matrix_free(m);

  return 0;
}

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


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

int
main (void)
{
  int i;
  /* coefficients of P(x) =  -1 + x^5  */
  double a[6] = { -1, 0, 0, 0, 0, 1 };
  double z[10];

  gsl_poly_complex_workspace * w
      = gsl_poly_complex_workspace_alloc (6);

  gsl_poly_complex_solve (a, 6, w, z);

  gsl_poly_complex_workspace_free (w);

  for (i = 0; i < 5; i++)
    {
      printf ("z%d = %+.18f %+.18f\n",
              i, z[2*i], z[2*i+1]);
    }

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* movstat/mvacc.c
 *
 * Moving window mean/variance accumulator - based on a modification
 * to Welford's algorithm, discussed here:
 *
 * https://stackoverflow.com/a/6664212
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>

typedef double ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  size_t n;      /* window size */
  size_t k;      /* number of samples currently in window */
  double mean;   /* current window mean */
  double M2;     /* current window M2 */
  ringbuf *rbuf; /* ring buffer storing current window */
} mvacc_state_t;

static size_t
mvacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(mvacc_state_t);
  size += ringbuf_size(n);

  return size;
}

static int
mvacc_init(const size_t n, void * vstate)
{
  mvacc_state_t * state = (mvacc_state_t *) vstate;

  state->n = n;
  state->k = 0;
  state->mean = 0.0;
  state->M2 = 0.0;

  state->rbuf = (ringbuf *) ((unsigned char *) vstate + sizeof(mvacc_state_t));
  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
mvacc_insert(const double x, void * vstate)
{
  mvacc_state_t * state = (mvacc_state_t *) vstate;

  if (ringbuf_is_full(state->rbuf))
    {
      /* remove oldest window element and add new one */
      double old = ringbuf_peek_back(state->rbuf);
      double prev_mean = state->mean;

      state->mean += (x - old) / (double) state->n;
      state->M2 += ((old - prev_mean) + (x - state->mean)) * (x - old);
    }
  else
    {
      double delta = x - state->mean;

      /*
       * Welford algorithm:
       *
       * mu_new = mu_old + (x - mu_old) / n
       * M2_new = M2_old + (x - mu_old) * (x - mu_new)
       */

      ++(state->k);
      state->mean += delta / (double) state->k;
      state->M2 += delta * (x - state->mean);
    }

  /* add new element to ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
mvacc_delete(void * vstate)
{
  mvacc_state_t * state = (mvacc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    {
      if (state->k > 1)
        {
          /*
           * mu_new = mu_old + (mu_old - x_old) / (n - 1)
           * M2_new = M2_old - (mu_old - x_old) * (mu_new - x_old)
           */

          double old = ringbuf_peek_back(state->rbuf);
          double prev_mean = state->mean;
          double delta = prev_mean - old;

          state->mean += delta / (state->k - 1.0);
          state->M2 -= delta * (state->mean - old);
        }
      else if (state->k == 1)
        {
          state->mean = 0.0;
          state->M2 = 0.0;
        }

      ringbuf_pop_back(state->rbuf);
      --(state->k);
    }

  return GSL_SUCCESS;
}

static int
mvacc_mean(void * params, double * result, const void * vstate)
{
  const mvacc_state_t * state = (const mvacc_state_t *) vstate;
  (void) params;
  *result = state->mean;
  return GSL_SUCCESS;
}

static int
mvacc_variance(void * params, double * result, const void * vstate)
{
  const mvacc_state_t * state = (const mvacc_state_t *) vstate;

  (void) params;

  if (state->k < 2)
    *result = 0.0;
  else
    *result = state->M2 / (state->k - 1.0);

  return GSL_SUCCESS;
}

static int
mvacc_sd(void * params, double * result, const void * vstate)
{
  double variance;
  int status = mvacc_variance(params, &variance, vstate);
  *result = sqrt(variance);
  return status;
}

static const gsl_movstat_accum mean_accum_type =
{
  mvacc_size,
  mvacc_init,
  mvacc_insert,
  mvacc_delete,
  mvacc_mean
};

const gsl_movstat_accum *gsl_movstat_accum_mean = &mean_accum_type;

static const gsl_movstat_accum variance_accum_type =
{
  mvacc_size,
  mvacc_init,
  mvacc_insert,
  mvacc_delete,
  mvacc_variance
};

const gsl_movstat_accum *gsl_movstat_accum_variance = &variance_accum_type;

static const gsl_movstat_accum sd_accum_type =
{
  mvacc_size,
  mvacc_init,
  mvacc_insert,
  mvacc_delete,
  mvacc_sd
};

const gsl_movstat_accum *gsl_movstat_accum_sd = &sd_accum_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <gsl/gsl_statistics.h>
#include <assert.h>

extern "C" {
#include <quadmath.h>
}

#ifndef N
#define N 64
#endif

#define K N/2

#ifndef FB
#define FB 64
#endif

#if(FB==64)
#define FT double
#elif(FB==80)
#define FT long double
#endif

int
main(int argc, char *argv[])
{
  assert(argc == 3);

  int i;
  char *inname = argv[1];
  char *outname = argv[2];
  FILE *infile = fopen(inname, "r");
  FILE *outfile = fopen(outname, "w");

  FT data[K], w[K];
  FT wsd;

  for (i = 0 ; i < K ; i++) {
    __float128 in_data;
    fread(&in_data, sizeof(__float128), 1, infile);
    data[i] = (float) in_data;
  }

  for (i = 0 ; i < K ; i++) {
    __float128 in_data;
    fread(&in_data, sizeof(__float128), 1, infile);
    w[i] = (float) in_data;
  }

  fclose(infile);

  // library call
  double wmean;
  #if(FB==64)
    wmean   = gsl_stats_wmean(w, 1, data, 1, K);
    wsd     = gsl_stats_wsd_with_fixed_mean(w, 1, data, 1, K, wmean);
  #elif(FB==80)
    wmean   = gsl_stats_long_double_wmean(w, 1, data, 1, K);
    wsd     = gsl_stats_long_double_wsd_with_fixed_mean(w, 1, data, 1, K, wmean);
  #else
    exit(-1);
  #endif

  //printf ("The sample wsd is %g\n", wsd);

  __float128 out_data;
  out_data = (__float128) wsd;
  fwrite(&out_data, sizeof(__float128), 1, outfile);
  fclose(outfile);

  return 0;
}

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


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

#ifndef _MCAS_COMMON_BYTE_SPAN_
#define _MCAS_COMMON_BYTE_SPAN_

#include <common/pointer_cast.h>
#include <cstddef>
#include <gsl/gsl_byte>
#include <gsl/span>
#include <sys/uio.h>

#ifndef MCAS_SPAN_USES_GSL
/* For minimal difference with older code; implement span as iovec (not gsl::span) */
#define MCAS_SPAN_USES_GSL 0
#endif
#ifndef MCAS_BYTE_USES_STD
/* For compilcation with C++14, use gsl::byte, not C++17 std::byte */
#define MCAS_BYTE_USES_STD 0
#endif

namespace common
{
#if MCAS_BYTE_USES_STD
	using byte = std::byte;
#else
	using byte = gsl::byte; /* can be std::byte in C++17 */
#endif
	template <typename T> using span = gsl::span<T>; /* can be std::span in C++20 */
	/* span of a const area. No equivalent in ::iovec, so always use span */
	using const_byte_span = span<const byte>;
	inline const_byte_span make_const_byte_span(const void *base, std::size_t len)
    {
      return const_byte_span(static_cast<const_byte_span::pointer>(base), len);
    }
}

namespace
{
	/* Accessors: Non-member in order to match ::iovec accessors */
	/* Start of area as a void *, for use with %p format and for conversion to an arbitrary type */
	inline const void *base(const common::const_byte_span &r) { return r.data(); }
    /* Length of area in bytes, for use in byte-based address calculations and comparisons */
	inline std::size_t size(const common::const_byte_span &r) { return r.size(); }
	/* Start of area as byte *, for use in byte-based address calculations and comparisons */
	inline const common::byte *data(const common::const_byte_span &r) { return r.data(); }
	/* End of area as byte *, for use in byte-based address calculations and comparisons */
	inline const common::byte *data_end(const common::const_byte_span &r) { return r.data() + r.size(); }
	/* End of are as a void *, for use with %p format */
	inline const void *end(const common::const_byte_span &r) { return ::data_end(r); }
}

namespace
{
	/* Same accessors as above, for ::iovec */
	inline void *base(const ::iovec &r) { return r.iov_base; }
	inline std::size_t size(const ::iovec &r) { return r.iov_len; }
	inline common::byte *data(const ::iovec &r) { return static_cast<common::byte *>(::base(r)); }
	inline common::byte *data_end(const ::iovec &r) { return ::data(r) + ::size(r); }
	inline void *end(const ::iovec &r) { return ::data_end(r); }
}

namespace common
{
	inline constexpr ::iovec make_iovec(void *base, std::size_t len) { return ::iovec{base, len}; }
}

#if MCAS_SPAN_USES_GSL
#include <common/pointer_cast.h>
namespace common
{
	using byte_span = span<byte>;
}

namespace
{
	/* Same accessors as above, for non-const byte span */
	inline void *base(const common::byte_span &r) { return r.data(); }
	inline std::size_t size(const common::byte_span &r) { return r.size(); }
	inline common::byte *data(const common::byte_span &r) { return r.data(); }
	inline common::byte *data_end(const common::byte_span &r) { return ::data(r) + ::size(r); }
	inline void *end(const common::byte_span &r) { return ::data_end(r); }
}

namespace common
{
	/* Construct a byte_span in syntax compatible with iovec (function, not constructor) */
	inline byte_span make_byte_span(void *base, std::size_t len) { return byte_span(common::pointer_cast<byte_span::value_type>(base), len); }
}
#else /* ! MCAS_SPAN_USES_GSL */
namespace common
{
	using byte_span = ::iovec;
	/* Construct an iovec in syntax compatible with span (no braces) */
	inline byte_span make_byte_span(void *base, std::size_t len) { return byte_span{base, len}; }
}
#endif

namespace common
{
	/* Converion from span of non-const to span of const */
	inline const_byte_span make_const_byte_span(const byte_span s) { return const_byte_span(::data(s), ::size(s)); }
}

#endif


text:
#ifndef ARRUS_CORE_DEVICES_TXRXPARAMETERS_H
#define ARRUS_CORE_DEVICES_TXRXPARAMETERS_H

#include <gsl/gsl>
#include <utility>
#include <ostream>

#include "arrus/core/api/common/Interval.h"
#include "arrus/core/api/common/Tuple.h"
#include "arrus/core/api/common/types.h"
#include "arrus/common/format.h"
#include "arrus/core/common/collections.h"
#include "arrus/core/api/ops/us4r/Pulse.h"

namespace arrus::devices {

class TxRxParameters {
public:
    static const TxRxParameters US4OEM_NOP;

    static TxRxParameters createRxNOPCopy(const TxRxParameters& op) {
        return TxRxParameters(
            op.txAperture,
            op.txDelays,
            op.txPulse,
            BitMask(op.rxAperture.size(), false),
            op.rxSampleRange,
            op.rxDecimationFactor,
            op.pri,
            op.rxPadding
        );
    }

    /**
     *
     * ** tx aperture, tx delays and rx aperture should have the same size
     * (tx delays is NOT limited to the tx aperture active elements -
     * the whole array must be provided).**
     *
     * @param txAperture
     * @param txDelays
     * @param txPulse
     * @param rxAperture
     * @param rxSampleRange [start, end) range of samples to acquire, starts from 0
     * @param rxDecimationFactor
     * @param pri
     * @param rxPadding how many 0-channels padd from the left and right
     */
    TxRxParameters(std::vector<bool> txAperture,
                   std::vector<float> txDelays,
                   const ops::us4r::Pulse &txPulse,
                   std::vector<bool> rxAperture,
                   Interval<uint32> rxSampleRange,
                   uint32 rxDecimationFactor, float pri,
                   Tuple<ChannelIdx> rxPadding = {0, 0})
        : txAperture(std::move(txAperture)), txDelays(std::move(txDelays)),
          txPulse(txPulse),
          rxAperture(std::move(rxAperture)), rxSampleRange(std::move(rxSampleRange)),
          rxDecimationFactor(rxDecimationFactor), pri(pri),
          rxPadding(std::move(rxPadding)){}

    [[nodiscard]] const std::vector<bool> &getTxAperture() const {
        return txAperture;
    }

    [[nodiscard]] const std::vector<float> &getTxDelays() const {
        return txDelays;
    }

    [[nodiscard]] const ops::us4r::Pulse &getTxPulse() const {
        return txPulse;
    }

    [[nodiscard]] const std::vector<bool> &getRxAperture() const {
        return rxAperture;
    }

    [[nodiscard]] const Interval<uint32> &getRxSampleRange() const {
        return rxSampleRange;
    }

    [[nodiscard]] uint32 getNumberOfSamples() const {
        return rxSampleRange.end() - rxSampleRange.start();
    }

    [[nodiscard]] int32 getRxDecimationFactor() const {
        return rxDecimationFactor;
    }

    [[nodiscard]] float getPri() const {
        return pri;
    }

    [[nodiscard]] const Tuple<ChannelIdx> &getRxPadding() const {
        return rxPadding;
    }

    [[nodiscard]] bool isNOP() const  {
        auto atLeastOneTxActive = ::arrus::reduce(
            std::begin(txAperture),
            std::end(txAperture),
            false, [](auto a, auto b) {return a | b;});
        auto atLeastOneRxActive = ::arrus::reduce(
            std::begin(rxAperture),
            std::end(rxAperture),
            false, [](auto a, auto b) {return a | b;});
        return !atLeastOneTxActive && !atLeastOneRxActive;
    }

    [[nodiscard]] bool isRxNOP() const {
        auto atLeastOneRxActive = ::arrus::reduce(
            std::begin(rxAperture),
            std::end(rxAperture),
            false, [](auto a, auto b) {return a | b;});
        return !atLeastOneRxActive;
    }

    friend std::ostream &
    operator<<(std::ostream &os, const TxRxParameters &parameters) {
        os << "Tx/Rx: ";
        os << "TX: ";
        os << "aperture: " << ::arrus::toString(parameters.getTxAperture())
           << ", delays: " << ::arrus::toString(parameters.getTxDelays())
           << ", center frequency: " << parameters.getTxPulse().getCenterFrequency()
           << ", n. periods: " << parameters.getTxPulse().getNPeriods()
           << ", inverse: " << parameters.getTxPulse().isInverse();
        os << "; RX: ";
        os << "aperture: " << ::arrus::toString(parameters.getRxAperture());
        os << "sample range: " << parameters.getRxSampleRange().start() << ", "
           << parameters.getRxSampleRange().end();
        os << ", fs divider: " << parameters.getRxDecimationFactor();
        os << std::endl;
        return os;
    }

    bool operator==(const TxRxParameters &rhs) const {
        return txAperture == rhs.txAperture &&
               txDelays == rhs.txDelays &&
               txPulse == rhs.txPulse &&
               rxAperture == rhs.rxAperture &&
               rxSampleRange == rhs.rxSampleRange &&
               rxDecimationFactor == rhs.rxDecimationFactor &&
               pri == rhs.pri;
    }

    bool operator!=(const TxRxParameters &rhs) const {
        return !(rhs == *this);
    }
private:
    ::std::vector<bool> txAperture;
    ::std::vector<float> txDelays;
    ::arrus::ops::us4r::Pulse txPulse;
    ::std::vector<bool> rxAperture;
    // TODO change to a simple pair
    Interval<uint32> rxSampleRange;
    int32 rxDecimationFactor;
    float pri;
    Tuple<ChannelIdx> rxPadding;
};

using TxRxParamsSequence = std::vector<TxRxParameters>;

/**
 * Returns the number of actual ops, that is, a the number of ops excluding RxNOPs.
 */
uint16 getNumberOfNoRxNOPs(const TxRxParamsSequence &seq);

}

#endif //ARRUS_CORE_DEVICES_TXRXPARAMETERS_H


text:
/*=========================================================================

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id$
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date$
  Version   : $Revision$
  Changes   : $Author$

Copyright (c) 1999-2014 and onwards, Imperial College London
All rights reserved.
See LICENSE for details

=========================================================================*/

#ifndef _IRTKVECTOR_H

#define _IRTKVECTOR_H

#ifdef USE_VXL
// VXL header file was loaded in irtkMatrix.h
#include <vnl/vnl_diag_matrix.h>
#else
#include <gsl/gsl_vector.h>
#endif

#include <iostream>
#include <cmath>

#include <irtkObject.h>

/**

  Vector class.

*/

class irtkVector : public irtkObject
{

protected:

  /// Number of rows
  int _rows;

  /// Data
  double *_vector;

public:

  /// Default constructor
  irtkVector();

  /// Constructor for given row dimensions
  irtkVector(int);

  /// Copy constructor
  irtkVector(const irtkVector &);

  /// Destructor
  ~irtkVector();

  /// Intialize matrix with number of rows
  void Initialize(int);

  //
  // Vector access functions
  //

  /// Returns number of rows
  int Rows() const;

  /// Puts vector value
  void   Put(int, double);

  /// Gets vector value
  double Get(int) const;

  //
  // Operators for vector access
  //

  /// Puts vector value
  double &operator()(int);

  /// Gets vector value
  double  operator()(int) const;

  //
  // Vector operators for doubles
  //

  /// Subtraction of a double
  irtkVector& operator-=(double);

  /// Addition of a double
  irtkVector& operator+=(double);

  /// Multiplication with a double
  irtkVector& operator*=(double);

  /// Division by a double
  irtkVector& operator/=(double);

  /// Return result of subtraction of a double
  irtkVector  operator- (double);

  /// Return result of addition of a double
  irtkVector  operator+ (double);

  /// Return result of multiplication with a double
  irtkVector  operator* (double);

  /// Return result of division by a double
  irtkVector  operator/ (double);

  //
  // Vector operators for vectors
  //

  /// Vector copy operator
  irtkVector& operator =(const irtkVector&);

  /// Vector subtraction operator
  irtkVector& operator-=(const irtkVector&);

  /// Vector addition operator
  irtkVector& operator+=(const irtkVector&);

  /// Vector componentwise multiplication operator (no scalar nor cross product)
  irtkVector& operator*=(const irtkVector&);

  /// Vector componentwise division operator
  irtkVector& operator/=(const irtkVector&);

  /// Return result for vector subtraction
  irtkVector  operator- (const irtkVector&);

  /// Return result for vector addition
  irtkVector  operator+ (const irtkVector&);

  /// Return result for componentwise vector multiplication (no scalar nor cross product)
  irtkVector  operator* (const irtkVector&);

  /// Return result for componentwise vector division
  irtkVector  operator/ (const irtkVector&);

  /// Comparison operator ==
  bool operator==(const irtkVector &);

#ifndef USE_STL
  /// Comparison operator != (if USE_STL is defined, negate == operator)
  bool operator!=(const irtkVector &);
#endif

  /// Comparison operator <
  bool operator<(const irtkVector &);

  /// Scalar/dot product
  double ScalarProduct(const irtkVector&);

  /// Vector/cross product
  irtkVector CrossProduct(const irtkVector&);

  /// Returns norm of a vector
  double Norm(void) const;

  /// Vector normalization
  void   Normalize(void);

  //
  // Vector in- and output
  //

  /// Interface to output stream
  friend std::ostream& operator<< (std::ostream&, const irtkVector&);

  /// Interface to input stream
  friend std::istream& operator>> (std::istream&, irtkVector&);

  /// Print vector
  void Print();

  /// Read vector from file
  void Read (char *);

  /// Write vector to file
  void Write(char *);

#ifdef USE_VXL
  template <class T> void Vector2Vnl(vnl_diag_matrix<T>*) const;

  template <class T> void Vnl2Vector(vnl_diag_matrix<T>*);

#else
  /// Conversion to GSL vector (memory must be allocated)
  void Vector2GSL(gsl_vector *) const;

  /// Conversion from GSL vector
  void GSL2Vector(gsl_vector *);
#endif

};

//
// Access operators
//

inline int irtkVector::Rows() const
{
  return _rows;
}

inline void irtkVector::Put(int rows, double vector)
{
  _vector[rows] = vector;
}

inline double irtkVector::Get(int rows) const
{
  return _vector[rows];
}

inline double& irtkVector::operator()(int rows)
{
  return _vector[rows];
}

inline double irtkVector::operator()(int rows) const
{
  return _vector[rows];
}

//
// Vector operators for doubles
//

inline irtkVector& irtkVector::operator-=(double x)
{
  int i;

  for (i = 0; i < _rows; i++) {
    _vector[i] -= x;
  }

  return *this;
}

inline irtkVector& irtkVector::operator+=(double x)
{
  int i;

  for (i = 0; i < _rows; i++) {
    _vector[i] += x;
  }

  return *this;
}

inline irtkVector& irtkVector::operator*=(double x)
{
  int i;

  for (i = 0; i < _rows; i++) {
    _vector[i] *= x;
  }

  return *this;
}

inline irtkVector& irtkVector::operator/=(double x)
{
  int i;

  for (i = 0; i < _rows; i++) {
    _vector[i] /= x;
  }

  return *this;
}

inline irtkVector irtkVector::operator- (double x)
{
  int i;

  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] - x;
  }
  return m;
}

inline irtkVector irtkVector::operator+ (double x)
{
  int i;

  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] + x;
  }
  return m;
}

inline irtkVector irtkVector::operator* (double x)
{
  int i;

  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] * x;
  }
  return m;
}

inline irtkVector irtkVector::operator/ (double x)
{
  int i;

  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] / x;
  }
  return m;
}

//
// Vector operators for vectors
//

inline irtkVector& irtkVector::operator =(const irtkVector& v)
{
  int i;

  // Copy size
  this->Initialize(v._rows);

  // Copy vector
  for (i = 0; i < _rows; i++) {
    _vector[i] = v._vector[i];
  }
  return *this;
}

inline irtkVector& irtkVector::operator-=(const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator-=: Size mismatch" << std::endl;
    exit(1);
  }
  for (i = 0; i < _rows; i++) {
    _vector[i] -= v._vector[i];
  }

  return *this;
}

inline irtkVector& irtkVector::operator+=(const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator+=: Size mismatch" << std::endl;
    exit(1);
  }
  for (i = 0; i < _rows; i++) {
    _vector[i] += v._vector[i];
  }

  return *this;
}

inline irtkVector& irtkVector::operator*=(const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator*=: Size mismatch" << std::endl;
    exit(1);
  }
  for (i = 0; i < _rows; i++) {
    _vector[i] *= v._vector[i];
  }

  return *this;
}

inline irtkVector& irtkVector::operator/=(const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator/=: Size mismatch" << std::endl;
    exit(1);
  }
  for (i = 0; i < _rows; i++) {
    _vector[i] /= v._vector[i];
  }

  return *this;
}

inline irtkVector irtkVector::operator- (const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator-: Size mismatch" << std::endl;
    exit(1);
  }
  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] - v._vector[i];
  }
  return m;
}

inline irtkVector irtkVector::operator+ (const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator+: Size mismatch" << std::endl;
    exit(1);
  }
  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] + v._vector[i];
  }
  return m;
}

inline irtkVector irtkVector::operator* (const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator*: Size mismatch" << std::endl;
    exit(1);
  }
  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] * v._vector[i];
  }
  return m;
}

inline irtkVector irtkVector::operator/ (const irtkVector& v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::operator/: Size mismatch" << std::endl;
    exit(1);
  }
  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = _vector[i] / v._vector[i];
  }
  return m;
}

//
// Comparison
//

inline bool  irtkVector::operator==(const irtkVector &v)
{
  if (_rows != v._rows) {
    return false;
  }
  for (int i = 0; i < _rows; i++) {
    if (_vector[i] != v._vector[i]) return false;
  }
  return true;

}

#ifndef USE_STL
inline bool  irtkVector::operator!=(const irtkVector &v)
{
  if (_rows != v._rows) {
    return true;
  }
  for (int i = 0; i < _rows; i++) {
    if (_vector[i] != v._vector[i]) return true;
  }
  return false;
}
#endif

inline bool  irtkVector::operator<(const irtkVector &v)
{
  if (_rows > v._rows) {
    return false;
  }
  for (int i = 0; i < _rows; i++) {
    if (_vector[i] >= v._vector[i]) return false;
  }
  return true;

}

//
// Vector products
//

inline double irtkVector::ScalarProduct(const irtkVector &v)
{
  int i;
  double scalar_product=0;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::ScalarProduct: Size mismatch" << std::endl;
    exit(1);
  }

  for (i = 0; i < _rows; i++) {
    scalar_product += _vector[i]*v._vector[i];
  }
  return scalar_product;
}

inline irtkVector irtkVector::CrossProduct(const irtkVector &v)
{
  int i;

  if (_rows != v._rows) {
    std::cerr << "irtkVector::CrossProduct: Size mismatch" << std::endl;
    exit(1);
  }

  irtkVector m(_rows);
  for (i = 0; i < _rows; i++) {
    m._vector[i] = (_vector[(i+1)%_rows]*v._vector[(i+2)%_rows]-
                    _vector[(i+2)%_rows]*v._vector[(i+1)%_rows]);
  }
  return m;
}

//
// Functions
//

inline double irtkVector::Norm(void) const
{
  double norm = 0;

  for (int i = 0; i < _rows; i++) {
    norm += _vector[i]*_vector[i];
  }
  return std::sqrt(norm);
}

inline void irtkVector::Normalize(void)
{
  double norm = Norm();

  if (norm != 0) {
    *this /= norm;
  }
}

#ifdef USE_VXL

template <class T>
inline void irtkVector::Vector2Vnl(vnl_diag_matrix<T>* m) const
{
  unsigned i;

  for (i = 0; i < (unsigned) _rows; i++) {
    (*m)(i) = (T) _vector[i];
  }
}

template <class T>
inline void irtkVector::Vnl2Vector(vnl_diag_matrix<T>* m)
{
  unsigned i;

  for (i = 0; i < (unsigned) _rows; i++) {
    _vector[i] = (float) (*m)(i);
  }
}

#endif

#endif

max_stars_repo_name: gordon-n-stevenson/fetalReconstruction

max_stars_repo_licenses: Zlib Unlicense Intel MIT

max_issues_repo_name: gordon-n-stevenson/fetalReconstruction

max_issues_repo_licenses: Zlib Unlicense Intel MIT

max_forks_repo_name: gordon-n-stevenson/fetalReconstruction

max_forks_repo_licenses: Zlib Unlicense Intel MIT