Common.h

Go to the documentation of this file.

#ifndef __COMMON_DEFINITIONS_H
#define __COMMON_DEFINITIONS_H

#include <list>
#include <valarray>
#include <limits>
#include <sstream>
#include <iostream>
#include <iterator>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include "gsl/gsl_poly.h"
#include <cassert>

#include "../Core/SharedLibraryExportMacros.h"
#include "Error.h"
#include "AstronomicalConstants.h"
#include "IncludeEigen.h"


namespace Core {

    const double NaN = std::numeric_limits<double>::quiet_NaN();

    const double Inf = std::numeric_limits<double>::infinity();

    enum EvolModeType {
        LOCKED_SURFACE_SPIN,

        BINARY,

        SINGLE,

        TABULATION
    };

    template<typename ITERATOR_TYPE>
        std::valarray<double> valarray_from_iterator(ITERATOR_TYPE values,
                                                     unsigned size)
        {
            std::valarray<double> result(size);
            double *dest=&result[0];
            for(; size>0; --size) {
                *dest = *values;
                ++dest;
                ++values;
            }
            return result;
        }

    LIB_LOCAL inline std::valarray<double> list_to_valarray(
        const std::list<double> &inlist
    ) {return valarray_from_iterator(inlist.begin(), inlist.size());}

    LIB_LOCAL std::valarray<double> solve_cubic(double c0,
                                                double c1,
                                                double c2,
                                                double c3);

    LIB_LOCAL double estimate_zerocrossing(
        double x0,

        double y0,

        double x1,

        double y1,

        double dy0=NaN,

        double dy1=NaN);

    LIB_LOCAL double quadratic_zerocrossing(double x0, double y0,
                                            double x1, double y1,
                                            double x2, double y2,
                                            double require_range_low = NaN,
                                            double require_range_high = NaN);

    LIB_LOCAL double cubic_zerocrossing(double x0, double y0,
                                        double x1, double y1,
                                        double x2, double y2,
                                        double x3, double y3,
                                        double require_range_low = NaN,
                                        double require_range_high = NaN);


    LIB_LOCAL double estimate_extremum(
        double x0,

        double y0,

        double x1,

        double y1,

        double dy0,

        double dy1,

        double *extremum_y = NULL);


    LIB_LOCAL double quadratic_extremum(double x0, double y0,
                                        double x1, double y1,
                                        double x2, double y2,
                                        double *extremum_y = NULL);

    LIB_LOCAL double cubic_extremum(double x0, double y0,
                                    double x1, double y1,
                                    double x2, double y2,
                                    double x3, double y3,
                                    double *extremum_y = NULL,
                                    double require_range_low = NaN,
                                    double require_range_high = NaN);

    LIB_LOCAL gsl_vector *cubic_coefficients(double x0, double y0,
                                             double x1, double y1,
                                             double x2, double y2,
                                             double x3, double y3);

    template<class ITERATOR>
        gsl_vector *polynomial_coefficients(ITERATOR x_i,
                                            ITERATOR y_i,
                                            size_t num_points)
        {
            Eigen::VectorXd y_vec(num_points);
            Eigen::MatrixXd xpowers(num_points, num_points);
            for(size_t i=0; i<num_points; i++) {
                double x=*x_i, xpow=1.0;
                y_vec(i)=*y_i;
                for(size_t pow=0; pow<num_points; pow++) {
                    xpowers(i,pow)=xpow;
                    xpow*=x;
                }
                x_i++; y_i++;
            }
            static Eigen::JacobiSVD<Eigen::MatrixXd,
                Eigen::FullPivHouseholderQRPreconditioner>
                    svd(num_points, num_points,
                        Eigen::ComputeFullU | Eigen::ComputeFullV);
            svd.compute(xpowers);
            Eigen::VectorXd solution=svd.solve(y_vec);
            gsl_vector *result=gsl_vector_alloc(num_points);
            for(size_t i=0; i<num_points; ++i)
                gsl_vector_set(result, i, solution(i));
            return result;

            /*  gsl_vector *y_vec=gsl_vector_alloc(num_points);
                gsl_matrix *xpowers=gsl_matrix_alloc(num_points, num_points);
                for(size_t i=0; i<num_points; i++) {
                double x=*x_i, xpow=1.0;
                gsl_vector_set(y_vec, i, *y_i);
                for(size_t pow=0; pow<num_points; pow++) {
                gsl_matrix_set(xpowers, i, pow, xpow);
                xpow*=x;
                }
                x_i++; y_i++;
                }
                gsl_matrix *xpowers_LU=gsl_matrix_alloc(num_points, num_points);
                gsl_matrix_memcpy(xpowers_LU, xpowers);

                gsl_permutation *permutation=gsl_permutation_alloc(num_points);
                gsl_vector *coefficients=gsl_vector_alloc(num_points);
                int permutation_sign;
                gsl_linalg_LU_decomp(xpowers_LU, permutation, &permutation_sign);
                gsl_linalg_LU_solve(xpowers_LU, permutation, y_vec, coefficients);
                gsl_vector *residuals=gsl_vector_alloc(num_points);
                gsl_linalg_LU_refine(xpowers, xpowers_LU, permutation, y_vec,
                coefficients, residuals);
                gsl_permutation_free(permutation);
                gsl_vector_free(y_vec);
                gsl_vector_free(residuals);
                gsl_matrix_free(xpowers);
                gsl_matrix_free(xpowers_LU);
                return coefficients;*/
        }

    LIB_LOCAL std::ostream &operator<<(std::ostream &os,
                                       const Core::EvolModeType &evol_mode);

} //End Core namespace.

namespace std {
    LIB_LOCAL std::ostream &operator<<(std::ostream &os,
                                       const std::valarray<double> &arr);

    LIB_LOCAL std::ostream &operator<<(std::ostream &os,
                                       const std::list<double> &l);
}

#endif