Homec4science
Diffusion LibMultiScale (features/pybind) Browse

libmultiscale/src/math/math_tools.cc
b05d5f2aba78features/pybind

Last Change
Raw File

math_tools.cc
View Options

/**
* @file math_tools.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
*
* @date Thu Jun 27 14:49:39 2013
*
* @brief This provides access to a set of math functions for various usage
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#include "lm_common.hh"
#include <cmath>
#include <cstdlib>
#include <vector>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
namespace MathTools {
Real ipow(Real x,int y) {
int i;
Real res=1.0;
for (i=y;i>0;--i){
res *= x;
}
return res;
};
/* -------------------------------------------------------------------------- */
Real zeta7[20],zeta13[20],zeta12[20],zeta6[20];
/* -------------------------------------------------------------------------- */
Real zeta(UInt puis,UInt ordre){
Real res=0;
UInt i,pos;
pos = ordre;
if (ordre > 19) pos = 19;
if(puis==7) {
if (zeta7[pos]==0.0) {
for (i = 1; i <= ordre ; ++i) {
zeta7[pos] += 1/ipow(i,puis);
}
}
return zeta7[pos];
}
if(puis==13){
if (zeta13[pos]==0.0){
for (i=1;i <= ordre ; ++i){
zeta13[pos] += 1/ipow(i,puis);
}
}
return zeta13[pos];
}
if(puis==12){
if (zeta12[pos]==0.0){
for (i=1;i <= ordre ; ++i){
zeta12[pos] += 1/ipow(i,puis);
}
}
return zeta12[pos];
}
if(puis==6){
if (zeta6[pos]==0.0){
for (i=1;i <= ordre ; ++i){
zeta6[pos] += 1/ipow(i,puis);
}
}
return zeta6[pos];
}
for (i=1;i <= ordre ; ++i){
res += 1/ipow(i,puis);
}
return res;
}
/* ------------------------------------------------------------------------ */
static long seed;
/* -------------------------------------------------------------------------- */
void setSeed(long s) {
seed = s;
srand48(seed);
}
/* ------------------------------------------------------------------------- */
/* A uniform distribution (between a and b) */
Real uniformRandom(Real a, Real b){
return drand48() * (b - a) + a;
}
/* -------------------------------------------------------------------------- */
/* A uniform distribution (between 0 and 1) */
Real uniform() {
return drand48();
}
/* -------------------------------------------------------------------------- */
/* A normal (Gaussian) distribution (with mean mu and stdev sig) */
Real gaussian() {
//static int iset=0;
Real fac, rsq, v1, v2;
/* If selec=0 we use zero mean and variance 1 */
/* else we use the specified mean and variance */
do {
v1=(2.0*uniform()-1.0);
v2=(2.0*uniform()-1.0);
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
return (v2*fac);
}
/* -------------------------------------------------------------------------- */
//! function to compute standard deviation of a set of values stored in an array
Real stdev(Real *x, Real mu, UInt n) {
UInt i;
Real std=0.0;
Real tmp=0.0;
for (i=0; i<n; i++) {
tmp = (*(x+i)-mu);
std += (tmp*tmp);
}
return sqrt(std/n);
}
Real stdev(std::vector<Real> & x, Real mu) {
return stdev(&x[0], mu,x.size());
}
/* ------------------------------------------------------------------------ */
//! function to compute the kurtosis of a set of values stored in an array
Real kurtosis(Real *x, Real mu, Real std, int n) {
int i;
Real kur=0.0;
Real tmp=0.0;
for (i=0; i<n; i++) {
tmp = (*(x+i)-mu);
kur += (tmp*tmp*tmp*tmp);
}
kur /= (n*std*std*std*std);
return kur-3;
}
/* -------------------------------------------------------------------------- */
//! function to compute the skewness of a set of values stored in an array
Real skewness(Real *x, Real mu, Real std, int n) {
int i;
Real skw=0.0;
Real tmp=0.0;
for (i=0; i<n; i++) {
tmp = (*(x+i)-mu);
skw += (tmp*tmp*tmp);
}
skw /= (n*std*std*std);
return skw;
}
/* -------------------------------------------------------------------------- */
//! function to compute the average of a set of values stored in an array
Real average(Real *x, UInt n) {
UInt i;
Real avg=0.0;
for (i=0; i<n; i++) {
avg += *(x+i);
}
return avg/n;
}
Real average(std::vector<Real> & x) {
return average(&x[0],x.size());
}
/* -------------------------------------------------------------------------- */
}
__END_LIBMULTISCALE__
c4science · Help