diff --git a/src/backend/lib/macros.h b/src/backend/lib/macros.h index 7d7dc4abe..3f6a2f53c 100644 --- a/src/backend/lib/macros.h +++ b/src/backend/lib/macros.h @@ -1,490 +1,492 @@ /*************************************************************************** File : macros.h Project : LabPlot Description : Various preprocessor macros -------------------------------------------------------------------- Copyright : (C) 2008 Tilman Benkert (thzs@gmx.net) Copyright : (C) 2013-2015 Alexander Semke (alexander.semke@web.de) Copyright : (C) 2016-2017 Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * 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 MACROS_H #define MACROS_H #include #include // C++ style warning (works on Windows) #include +#ifndef WARNING #define WARNING(x) std::cout << x << std::endl +#endif #ifndef NDEBUG #include #define QDEBUG(x) qDebug() << x // C++ style debugging (works on Windows) #define DEBUG(x) std::cout << x << std::endl #else #define QDEBUG(x) {} #define DEBUG(x) {} #endif #define ENUM_TO_STRING(class, enum, value) \ (class::staticMetaObject.enumerator(class::staticMetaObject.indexOfEnumerator(#enum)).valueToKey(value)) #define ENUM_COUNT(class, enum) \ (class::staticMetaObject.enumerator(class::staticMetaObject.indexOfEnumerator(#enum)).keyCount()) #define BASIC_ACCESSOR(type, var, method, Method) \ type method() const { return var; }; \ void set ## Method(const type value) { var = value; } #define CLASS_ACCESSOR(type, var, method, Method) \ type method() const { return var; }; \ void set ## Method(const type & value) { var = value; } #define BASIC_D_ACCESSOR_DECL(type, method, Method) \ type method() const; \ void set ## Method(const type value); #define BASIC_D_ACCESSOR_IMPL(classname, type, method, Method, var) \ void classname::set ## Method(const type value) \ { \ d->var = value; \ } \ type classname::method() const \ { \ return d->var; \ } #define BASIC_D_READER_IMPL(classname, type, method, var) \ type classname::method() const \ { \ return d->var; \ } #define BASIC_SHARED_D_READER_IMPL(classname, type, method, var) \ type classname::method() const \ { \ Q_D(const classname); \ return d->var; \ } #define CLASS_D_ACCESSOR_DECL(type, method, Method) \ type method() const; \ void set ## Method(const type & value); #define CLASS_D_ACCESSOR_IMPL(classname, type, method, Method, var) \ void classname::set ## Method(const type & value) \ { \ d->var = value; \ } \ type classname::method() const \ { \ return d->var; \ } #define CLASS_D_READER_IMPL(classname, type, method, var) \ type classname::method() const \ { \ return d->var; \ } #define CLASS_SHARED_D_READER_IMPL(classname, type, method, var) \ type classname::method() const \ { \ Q_D(const classname); \ return d->var; \ } #define POINTER_D_ACCESSOR_DECL(type, method, Method) \ type *method() const; \ void set ## Method(type *ptr); #define FLAG_D_ACCESSOR_DECL(Method) \ bool is ## Method() const; \ bool has ## Method() const; \ void set ## Method(const bool value=true); \ void enable ## Method(const bool value=true); #define FLAG_D_ACCESSOR_IMPL(classname, Method, var) \ void classname::set ## Method(const bool value) \ { \ d->var = value; \ } \ void classname::enable ## Method(const bool value) \ { \ d->var = value; \ } \ bool classname::is ## Method() const \ { \ return d->var; \ } \ bool classname::has ## Method() const \ { \ return d->var; \ } #define WAIT_CURSOR QApplication::setOverrideCursor(QCursor(Qt::WaitCursor)) #define RESET_CURSOR QApplication::restoreOverrideCursor() #define STD_SETTER_CMD_IMPL(class_name, cmd_name, value_type, field_name) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ }; #define STD_SETTER_CMD_IMPL_F(class_name, cmd_name, value_type, field_name, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void finalize() { m_target->finalize_method(); } \ }; // setter class with finalize() and signal emmiting. #define STD_SETTER_CMD_IMPL_S(class_name, cmd_name, value_type, field_name) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void finalize() { emit m_target->q->field_name##Changed(m_target->*m_field); } \ }; #define STD_SETTER_CMD_IMPL_F_S(class_name, cmd_name, value_type, field_name, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void finalize() { m_target->finalize_method(); emit m_target->q->field_name##Changed(m_target->*m_field); } \ }; // setter class with finalize() and signal emmiting for changing several properties in one single step (embedded in beginMacro/endMacro) #define STD_SETTER_CMD_IMPL_M_F_S(class_name, cmd_name, value_type, field_name, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardMacroSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardMacroSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void finalize() { m_target->finalize_method(); emit m_target->q->field_name##Changed(m_target->*m_field); } \ virtual void finalizeUndo() { emit m_target->q->field_name##Changed(m_target->*m_field); } \ }; #define STD_SETTER_CMD_IMPL_I(class_name, cmd_name, value_type, field_name, init_method) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void initialize() { m_target->init_method(); } \ }; #define STD_SETTER_CMD_IMPL_IF(class_name, cmd_name, value_type, field_name, init_method, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSetterCmd(target, &class_name::Private::field_name, newValue, description) {} \ virtual void initialize() { m_target->init_method(); } \ virtual void finalize() { m_target->finalize_method(); } \ }; #define STD_SWAP_METHOD_SETTER_CMD_IMPL(class_name, cmd_name, value_type, method_name) \ class class_name ## cmd_name ## Cmd: public StandardSwapMethodSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSwapMethodSetterCmd(target, &class_name::Private::method_name, newValue, description) {} \ }; #define STD_SWAP_METHOD_SETTER_CMD_IMPL_F(class_name, cmd_name, value_type, method_name, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardSwapMethodSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSwapMethodSetterCmd(target, &class_name::Private::method_name, newValue, description) {} \ virtual void finalize() { m_target->finalize_method(); } \ }; #define STD_SWAP_METHOD_SETTER_CMD_IMPL_I(class_name, cmd_name, value_type, method_name, init_method) \ class class_name ## cmd_name ## Cmd: public StandardSwapMethodSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSwapMethodSetterCmd(target, &class_name::Private::method_name, newValue, description) {} \ virtual void initialize() { m_target->init_method(); } \ }; #define STD_SWAP_METHOD_SETTER_CMD_IMPL_IF(class_name, cmd_name, value_type, method_name, init_method, finalize_method) \ class class_name ## cmd_name ## Cmd: public StandardSwapMethodSetterCmd { \ public: \ class_name ## cmd_name ## Cmd(class_name::Private *target, value_type newValue, const QString &description) \ : StandardSwapMethodSetterCmd(target, &class_name::Private::method_name, newValue, description) {} \ virtual void initialize() { m_target->init_method(); } \ virtual void finalize() { m_target->finalize_method(); } \ }; //xml-serialization/deserialization //QColor #define WRITE_QCOLOR(color) \ do { \ writer->writeAttribute( "color_r", QString::number(color.red()) ); \ writer->writeAttribute( "color_g", QString::number(color.green()) ); \ writer->writeAttribute( "color_b", QString::number(color.blue()) ); \ } while (0) #define READ_QCOLOR(color) \ do { \ str = attribs.value("color_r").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_r'")); \ else \ color.setRed( str.toInt() ); \ \ str = attribs.value("color_g").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_g'")); \ else \ color.setGreen( str.toInt() ); \ \ str = attribs.value("color_b").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_b'")); \ else \ color.setBlue( str.toInt() ); \ } while(0) //QPen #define WRITE_QPEN(pen) \ do { \ writer->writeAttribute( "style", QString::number(pen.style()) ); \ writer->writeAttribute( "color_r", QString::number(pen.color().red()) ); \ writer->writeAttribute( "color_g", QString::number(pen.color().green()) ); \ writer->writeAttribute( "color_b", QString::number(pen.color().blue()) ); \ writer->writeAttribute( "width", QString::number(pen.widthF()) ); \ } while (0) #define READ_QPEN(pen) \ do { \ str = attribs.value("style").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'style'")); \ else \ pen.setStyle( (Qt::PenStyle)str.toInt() ); \ \ QColor color; \ str = attribs.value("color_r").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_r'")); \ else \ color.setRed( str.toInt() ); \ \ str = attribs.value("color_g").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_g'")); \ else \ color.setGreen( str.toInt() ); \ \ str = attribs.value("color_b").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'color_b'")); \ else \ color.setBlue( str.toInt() ); \ \ pen.setColor(color); \ \ str = attribs.value("width").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'width'")); \ else \ pen.setWidthF( str.toDouble() ); \ } while(0) //QFont #define WRITE_QFONT(font) \ do { \ writer->writeAttribute( "fontFamily", font.family() ); \ writer->writeAttribute( "fontSize", QString::number(font.pixelSize()) ); \ writer->writeAttribute( "fontPointSize", QString::number(font.pointSize()));\ writer->writeAttribute( "fontWeight", QString::number(font.weight()) ); \ writer->writeAttribute( "fontItalic", QString::number(font.italic()) ); \ } while(0) #define READ_QFONT(font) \ do { \ str = attribs.value("fontFamily").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'fontFamily'")); \ else \ font.setFamily( str ); \ \ str = attribs.value("fontSize").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'fontSize'")); \ else { \ int size = str.toInt(); \ if (size != -1) \ font.setPixelSize(size); \ } \ \ str = attribs.value("fontPointSize").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'fontPointSize'")); \ else { \ int size = str.toInt(); \ if (size != -1) \ font.setPointSize(size); \ } \ \ str = attribs.value("fontWeight").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'fontWeight'")); \ else \ font.setWeight( str.toInt() ); \ \ str = attribs.value("fontItalic").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'fontItalic'")); \ else \ font.setItalic( str.toInt() ); \ } while(0) //QBrush #define WRITE_QBRUSH(brush) \ do { \ writer->writeAttribute("brush_style", QString::number(brush.style()) ); \ writer->writeAttribute("brush_color_r", QString::number(brush.color().red())); \ writer->writeAttribute("brush_color_g", QString::number(brush.color().green()));\ writer->writeAttribute("brush_color_b", QString::number(brush.color().blue())); \ } while(0) #define READ_QBRUSH(brush) \ do { \ str = attribs.value("brush_style").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'brush_style'")); \ else \ brush.setStyle( (Qt::BrushStyle)str.toInt() ); \ \ QColor color; \ str = attribs.value("brush_color_r").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'brush_color_r'")); \ else \ color.setRed( str.toInt() ); \ \ str = attribs.value("brush_color_g").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'brush_color_g'")); \ else \ color.setGreen( str.toInt() ); \ \ str = attribs.value("brush_color_b").toString(); \ if(str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg("'brush_color_b'")); \ else \ color.setBlue( str.toInt() ); \ \ brush.setColor(color); \ } while(0) //Column #define WRITE_COLUMN(column, columnName) \ do { \ if (column){ \ writer->writeAttribute( #columnName, column->path() ); \ } else { \ writer->writeAttribute( #columnName, "" ); \ } \ } while(0) //column names can be empty in case no columns were used before save //the actual pointers to the x- and y-columns are restored in Project::load() #define READ_COLUMN(columnName) \ do { \ str = attribs.value(#columnName).toString(); \ d->columnName ##Path = str; \ } while(0) #define READ_INT_VALUE(name, var, type) \ str = attribs.value(name).toString(); \ if (str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg(name)); \ else \ d->var = (type)str.toInt(); #define READ_DOUBLE_VALUE(name, var) \ str = attribs.value(name).toString(); \ if (str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg(name)); \ else \ d->var = str.toDouble(); #define READ_STRING_VALUE(name, var) \ str = attribs.value(name).toString(); \ if (str.isEmpty()) \ reader->raiseWarning(attributeWarning.arg(name)); \ else \ d->var = str; //used in Project::load() #define RESTORE_COLUMN_POINTER(obj, col, Col) \ do { \ if (!obj->col ##Path().isEmpty()) { \ foreach (Column* column, columns) { \ if (!column) continue; \ if (column->path() == obj->col ##Path()) { \ obj->set## Col(column); \ break; \ } \ } \ } \ } while(0) #define WRITE_PATH(obj, name) \ do { \ if (obj){ \ writer->writeAttribute( #name, obj->path() ); \ } else { \ writer->writeAttribute( #name, "" ); \ } \ } while(0) #define READ_PATH(name) \ do { \ str = attribs.value(#name).toString(); \ d->name ##Path = str; \ } while(0) #define RESTORE_POINTER(obj, name, Name, Type, list) \ do { \ if (!obj->name ##Path().isEmpty()) { \ foreach (AbstractAspect* aspect, list) { \ if (aspect->path() == obj->name ##Path()) { \ Type * a = dynamic_cast(aspect); \ if (!a) continue; \ obj->set## Name(a); \ break; \ } \ } \ } \ } while(0) #endif // MACROS_H diff --git a/src/backend/nsl/nsl_filter.c b/src/backend/nsl/nsl_filter.c index 3e9e1978f..93483fdf3 100644 --- a/src/backend/nsl/nsl_filter.c +++ b/src/backend/nsl/nsl_filter.c @@ -1,316 +1,318 @@ /*************************************************************************** File : nsl_filter.c Project : LabPlot Description : NSL Fourier filter functions -------------------------------------------------------------------- Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_filter.h" #include "nsl_common.h" #include "nsl_sf_poly.h" #include #include #include #ifdef HAVE_FFTW3 #include #endif const char* nsl_filter_type_name[] = { i18n("Low pass"), i18n("High pass"), i18n("Band pass"), i18n("Band reject") }; const char* nsl_filter_form_name[] = { i18n("Ideal"), i18n("Butterworth"), i18n("Chebyshev type I"), i18n("Chebyshev type II"), i18n("Legendre (Optimum L)"), i18n("Bessel (Thomson)") }; const char* nsl_filter_cutoff_unit_name[] = { i18n("Frequency"), i18n("Fraction"), i18n("Index") }; /* n - order, x = w/w0 */ double nsl_filter_gain_bessel(int n, double x) { #ifdef _MSC_VER COMPLEX z0 = {0.0, 0.0}; - COMPLEX z = {0.0, x}; + COMPLEX z; + z.re = 0.0; + z.im = x; double norm = cabs(nsl_sf_poly_reversed_bessel_theta(n, z)); COMPLEX value = nsl_sf_poly_reversed_bessel_theta(n, z0); return creal(value)/norm; #else return nsl_sf_poly_reversed_bessel_theta(n, 0)/cabs(nsl_sf_poly_reversed_bessel_theta(n, I*x)); #endif } int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, double cutindex, double bandwidth) { if (cutindex < 0) { printf("index for cutoff must be >= 0\n"); return -1; } if (cutindex > n/2+1) { printf("index for cutoff must be <= n/2+1\n"); return -1; } size_t i; double factor, centerindex = cutindex + bandwidth/2.; switch (type) { case nsl_filter_type_low_pass: switch (form) { case nsl_filter_form_ideal: - for (i = cutindex; i < n/2+1; i++) + for (i = (size_t)cutindex; i < n/2+1; i++) data[2*i] = data[2*i+1] = 0; break; case nsl_filter_form_butterworth: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1. + gsl_sf_pow_int(i/cutindex, 2*order)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_i: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1. + gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i/cutindex), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_ii: for (i = 1; i < n/2+1; i++) { /* i==0: factor=1 */ factor = 1./sqrt(1. + 1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, cutindex/i), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_legendre: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, i*i/(cutindex*cutindex) )); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_bessel: for (i = 0; i < n/2+1; i++) { factor = nsl_filter_gain_bessel(order, i/cutindex); data[2*i] *= factor; data[2*i+1] *= factor; } break; } break; case nsl_filter_type_high_pass: switch (form) { case nsl_filter_form_ideal: for (i = 0; i < cutindex; i++) data[2*i] = data[2*i+1] = 0; break; case nsl_filter_form_butterworth: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int(cutindex/i, 2*order)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_i: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, cutindex/i), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_ii: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i/cutindex), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_legendre: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, cutindex*cutindex/(i*i) )); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_bessel: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = nsl_filter_gain_bessel(order, cutindex/i); data[2*i] *= factor; data[2*i+1] *= factor; } break; } break; case nsl_filter_type_band_pass: switch (form) { case nsl_filter_form_ideal: - for (i = 0; i < cutindex; i++) + for (i = 0; i < (size_t)cutindex; i++) data[2*i] = data[2*i+1] = 0; for (i = cutindex + bandwidth; i < n/2+1; i++) data[2*i] = data[2*i+1] = 0; break; case nsl_filter_form_butterworth: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int((i*i - centerindex*centerindex)/i/bandwidth, 2*order)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_i: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, (i*i - centerindex*centerindex)/i/bandwidth), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_ii: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i*bandwidth/(i*i - centerindex*centerindex)), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_legendre: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, (i*i-2.*centerindex*centerindex+gsl_pow_4(centerindex)/(i*i))/gsl_pow_2(bandwidth) )); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_bessel: data[0]=data[1]=0; for (i = 1; i < n/2+1; i++) { factor = nsl_filter_gain_bessel(order, (i*i - centerindex*centerindex)/i/bandwidth); data[2*i] *= factor; data[2*i+1] *= factor; } break; } break; case nsl_filter_type_band_reject: switch (form) { case nsl_filter_form_ideal: - for (i = cutindex; i < cutindex + bandwidth; i++) + for (i = (size_t)cutindex; i < (size_t)(cutindex + bandwidth); i++) data[2*i] = data[2*i+1] = 0; break; case nsl_filter_form_butterworth: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int(i*bandwidth/(i*i - centerindex*centerindex), 2*order)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_i: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i*bandwidth/(i*i - centerindex*centerindex)), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_chebyshev_ii: for (i = 1; i < n/2+1; i++) { /* i==0: factor=1 */ factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, (i*i - centerindex*centerindex)/i/bandwidth), 2)); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_legendre: for (i = 0; i < n/2+1; i++) { factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, gsl_pow_2(i*bandwidth)/gsl_pow_2(i*i-centerindex*centerindex) )); data[2*i] *= factor; data[2*i+1] *= factor; } break; case nsl_filter_form_bessel: for (i = 0; i < n/2+1; i++) { if (i == centerindex) factor = 0; else factor = nsl_filter_gain_bessel(order, i*bandwidth/(i*i - centerindex*centerindex)); data[2*i] *= factor; data[2*i+1] *= factor; } break; } break; } return 0; } void print_fdata(double data[], size_t n) { size_t i; for(i=0; i < 2*(n/2+1); i++) printf("%g ", data[i]); puts(""); printf("real: "); for(i=0; i < n/2+1; i++) printf("%g ", data[2*i]); puts(""); printf("imag: "); for(i=0; i < n/2+1; i++) printf("%g ", data[2*i+1]); puts(""); } int nsl_filter_fourier(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, int cutindex, int bandwidth) { /* 1. transform */ double* fdata = (double*)malloc(2*n*sizeof(double)); /* contains re0,im0,re1,im1,re2,im2,... */ #ifdef HAVE_FFTW3 fftw_plan plan = fftw_plan_dft_r2c_1d(n, data, (fftw_complex *) fdata, FFTW_ESTIMATE); fftw_execute(plan); fftw_destroy_plan(plan); #else gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n); gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n); gsl_fft_real_transform(data, 1, n, real, work); gsl_fft_real_wavetable_free(real); gsl_fft_halfcomplex_unpack(data, fdata, 1, n); #endif /* 2. apply filter */ /*print_fdata(fdata, n);*/ int status = nsl_filter_apply(fdata, n, type, form, order, cutindex, bandwidth); /*print_fdata(fdata, n);*/ /* 3. back transform */ #ifdef HAVE_FFTW3 plan = fftw_plan_dft_c2r_1d(n, (fftw_complex *) fdata, data, FFTW_ESTIMATE); fftw_execute(plan); fftw_destroy_plan(plan); /* normalize*/ size_t i; for (i=0; i < n; i++) data[i] /= n; #else gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(n); gsl_fft_halfcomplex_inverse(data, 1, n, hc, work); gsl_fft_halfcomplex_wavetable_free(hc); gsl_fft_real_workspace_free (work); #endif free(fdata); return status; } diff --git a/src/backend/nsl/nsl_fit.c b/src/backend/nsl/nsl_fit.c index e400b404f..de9a4d3dd 100644 --- a/src/backend/nsl/nsl_fit.c +++ b/src/backend/nsl/nsl_fit.c @@ -1,647 +1,647 @@ /*************************************************************************** File : nsl_fit.c Project : LabPlot Description : NSL (non)linear fit functions -------------------------------------------------------------------- Copyright : (C) 2016-2017 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_fit.h" #include "nsl_common.h" #include #include #include #include #include const char* nsl_fit_model_category_name[] = {i18n("Basic functions"), i18n("Peak functions"), i18n("Growth (sigmoidal)"), i18n("Statistics (distributions)"), i18n("Custom")}; const char* nsl_fit_model_basic_name[] = {i18n("Polynomial"), i18n("Power"), i18n("Exponential"), i18n("Inverse exponential"), i18n("Fourier")}; const char* nsl_fit_model_basic_equation[] = {"c0 + c1*x", "a*x^b", "a*exp(b*x)", "a*(1-exp(b*x)) + c", "a0 + (a1*cos(w*x) + b1*sin(w*x))"}; const char* nsl_fit_model_basic_pic_name[] = {"polynom", "power", "exponential", "inv_exponential", "fourier"}; const char* nsl_fit_model_peak_name[] = {i18n("Gaussian (normal)"), i18n("Cauchy-Lorentz"), i18n("Hyperbolic secant (sech)"), i18n("Logistic (sech squared)")}; const char* nsl_fit_model_peak_equation[] = {"a/sqrt(2*pi)/s * exp(-((x-mu)/s)^2/2)", "a/pi * g/(g^2+(x-mu)^2)", "a/pi/s * sech((x-mu)/s)", "a/4/s * sech((x-mu)/2/s)**2"}; const char* nsl_fit_model_peak_pic_name[] = {"gaussian", "cauchy_lorentz", "sech", "logistic"}; const char* nsl_fit_model_growth_name[] = {i18n("Inverse tangent"), i18n("Hyperbolic tangent"), i18n("Algebraic sigmoid"), i18n("Logistic function"), i18n("Error function (erf)"), i18n("Hill"), i18n("Gompertz"), i18n("Gudermann (gd)")}; const char* nsl_fit_model_growth_equation[] = {"a * atan((x-mu)/s)", "a * tanh((x-mu)/s)", "a * (x-mu)/s/sqrt(1+((x-mu)/s)^2)", "a/(1+exp(-k*(x-mu)))", "a/2 * erf((x-mu)/s/sqrt(2))", "a * x^n/(s^n + x^n)", "a*exp(-b*exp(-c*x))", "a * asin(tanh((x-mu)/s))"}; const char* nsl_fit_model_growth_pic_name[] = {"atan", "tanh", "alg_sigmoid", "logistic_function", "erf", "hill", "gompertz", "gd"}; const char* nsl_fit_error_type_name[] = {"No", "Direct (col)", "Inverse (1/col)"}; const char* nsl_fit_weight_type_name[] = {"No", "Instrumental (1/col^2)", "Direct (col)", "Inverse (1/col)", "Statistical (1/data)", "Statistical (Fit)", "Relative (1/data^2)", "Relative (Fit)"}; /* see http://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf and https://lmfit.github.io/lmfit-py/bounds.html */ double nsl_fit_map_bound(double x, double min, double max) { if (max <= min) { printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); return DBL_MAX; } /* not bounded */ if (min == -DBL_MAX && max == DBL_MAX) return x; /* open bounds */ if (min == -DBL_MAX) return max + 1. - sqrt(x*x + 1.); if (max == DBL_MAX) return min - 1. + sqrt(x*x + 1.); return min + sin(x + 1.) * (max - min)/2.; /* alternative transformation for closed bounds return min + (max - min)/(1. + exp(-x)); */ } /* see http://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf and https://lmfit.github.io/lmfit-py/bounds.html */ double nsl_fit_map_unbound(double x, double min, double max) { if (max <= min) { printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); return DBL_MAX; } if (x < min || x > max) { printf("given value must be within bounds! Giving up.\n"); return -DBL_MAX; } /* not bounded */ if (min == -DBL_MAX && max == DBL_MAX) return x; /* open bounds */ if (min == -DBL_MAX) return sqrt(gsl_pow_2(max - x + 1.) - 1.); if (max == DBL_MAX) return sqrt(gsl_pow_2(x - min + 1.) - 1.); return asin(2. * (x - min)/(max - min) - 1.); /* alternative transformation for closed bounds return -log((max - x)/(x - min)); */ } /********************** parameter derivatives ******************/ /* basic */ double nsl_fit_model_polynomial_param_deriv(double x, int j, double weight) { return weight*pow(x, j); } double nsl_fit_model_power1_param_deriv(int param, double x, double a, double b, double weight) { if (param == 0) return weight*pow(x, b); if (param == 1) return weight*a*pow(x, b)*log(x); return 0; } double nsl_fit_model_power2_param_deriv(int param, double x, double b, double c, double weight) { if (param == 0) return weight; if (param == 1) return weight*pow(x, c); if (param == 2) return weight*b*pow(x, c)*log(x); return 0; } double nsl_fit_model_exponentialn_param_deriv(int param, double x, double *p, double weight) { if (param % 2 == 0) return weight*exp(p[param+1]*x); else return weight*p[param-1]*x*exp(p[param]*x); } double nsl_fit_model_inverse_exponential_param_deriv(int param, double x, double a, double b, double weight) { if (param == 0) return weight*(1. - exp(b*x)); if (param == 1) return -weight*a*x*exp(b*x); if (param == 2) return weight; return 0; } double nsl_fit_model_fourier_param_deriv(int param, int degree, double x, double w, double weight) { if (param == 0) return weight*cos(degree*w*x); if (param == 1) return weight*sin(degree*w*x); return 0; } /* peak */ double nsl_fit_model_gaussian_param_deriv(int param, double x, double s, double mu, double A, double weight) { double s2 = s*s, norm = weight/sqrt(2.*M_PI)/s, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); if (param == 0) return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; if (param == 1) return A * norm/s2 * (x-mu) * efactor; if (param == 2) return norm * efactor; return 0; } double nsl_fit_model_lorentz_param_deriv(int param, double x, double s, double t, double A, double weight) { double norm = weight/M_PI, denom = s*s+(x-t)*(x-t); if (param == 0) return A * norm * ((x-t)*(x-t) - s*s)/(denom*denom); if (param == 1) return A * norm * 2.*s*(x-t)/(denom*denom); if (param == 2) return norm * s/denom; return 0; } double nsl_fit_model_sech_param_deriv(int param, double x, double s, double mu, double A, double weight) { double y = (x-mu)/s, norm = weight/M_PI/s; if (param == 0) return A/s * norm * (y*tanh(y)-1.)/cosh(y); if (param == 1) return A/s * norm * tanh(y)/cosh(y); if (param == 2) return norm/cosh(y); return 0; } double nsl_fit_model_logistic_param_deriv(int param, double x, double s, double mu, double A, double weight) { double y = (x-mu)/2./s, norm = weight/4./s; if (param == 0) return A/s * norm * (2.*y*tanh(y)-1.)/cosh(y); if (param == 1) return A/s * norm * tanh(y)/cosh(y)/cosh(y); if (param == 2) return norm/cosh(y)/cosh(y); return 0; } /* growth */ double nsl_fit_model_atan_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight, y = (x-mu)/s; if (param == 0) return -A/s * norm * y/(1.+y*y); if (param == 1) return -A/s * norm * 1./(1+y*y); if (param == 2) return norm * atan(y); return 0; } double nsl_fit_model_tanh_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight, y = (x-mu)/s; if (param == 0) return -A/s * norm * y/cosh(y)/cosh(y); if (param == 1) return -A/s * norm * 1./cosh(y)/cosh(y); if (param == 2) return norm * tanh(y); return 0; } double nsl_fit_model_algebraic_sigmoid_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight, y = (x-mu)/s, y2 = y*y; if (param == 0) return -A/s * norm * y/pow(1.+y2, 1.5); if (param == 1) return -A/s * norm * 1./pow(1.+y2, 1.5); if (param == 2) return norm * y/sqrt(1.+y2); return 0; } double nsl_fit_model_sigmoid_param_deriv(int param, double x, double k, double mu, double A, double weight) { double norm = weight, y = k*(x-mu); if (param == 0) return A/k * norm * y*exp(-y)/gsl_pow_2(1. + exp(-y)); if (param == 1) return -A*k * norm * exp(-y)/gsl_pow_2(1. + exp(-y)); if (param == 2) return norm/(1. + exp(-y)); return 0; } double nsl_fit_model_erf_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight, y = (x-mu)/(sqrt(2.)*s); if (param == 0) return -A/sqrt(M_PI)/s * norm * y*exp(-y*y); if (param == 1) return -A/sqrt(2.*M_PI)/s * norm * exp(-y*y); if (param == 2) return norm/2. * gsl_sf_erf(y); return 0; } double nsl_fit_model_hill_param_deriv(int param, double x, double s, double n, double A, double weight) { double norm = weight, y = x/s; if (param == 0) return -A*n/s * norm * pow(y, n)/gsl_pow_2(1.+pow(y, n)); if (param == 1) return A * norm * log(y)*pow(y, n)/gsl_pow_2(1.+pow(y, n)); if (param == 2) return norm * pow(y, n)/(1.+pow(y, n)); return 0; } double nsl_fit_model_gompertz_param_deriv(int param, double x, double a, double b, double c, double weight) { if (param == 0) return weight*exp(-b*exp(-c*x)); if (param == 1) return -weight*a*exp(-c*x-b*exp(-c*x)); if (param == 2) return weight*a*b*x*exp(-c*x-b*exp(-c*x)); return 0; } double nsl_fit_model_gudermann_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight, y = (x-mu)/s; if (param == 0) return -A/s * norm * y/cosh(y); if (param == 1) return -A/s * norm * 1./cosh(y); if (param == 2) return -asin(tanh(y)); return 0; } /* distributions */ double nsl_fit_model_gaussian_tail_param_deriv(int param, double x, double s, double mu, double A, double a, double weight) { if (x < a) return 0; double s2 = s*s, N = erfc(a/s/M_SQRT2)/2., norm = weight/sqrt(2.*M_PI)/s/N, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); if (param == 0) return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; if (param == 1) return A * norm/s2 * (x-mu) * efactor; if (param == 2) return norm * efactor; if (param == 3) return A/norm/norm * efactor * exp(-a*a/(2.*s2)); return 0; } double nsl_fit_model_exponential_param_deriv(int param, double x, double l, double mu, double A, double weight) { if (x < mu) return 0; double y = l*(x-mu), efactor = exp(-y); if (param == 0) return weight * A * (1. - y) * efactor; if (param == 1) return weight * A * gsl_pow_2(l) * efactor; if (param == 2) return weight * l * efactor; return 0; } double nsl_fit_model_laplace_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight/(2.*s), y = fabs((x-mu)/s), efactor = exp(-y); if (param == 0) return A/s*norm * (y-1.) * efactor; if (param == 1) return A/(s*s)*norm * (x-mu)/y * efactor; if (param == 2) return norm * efactor; return 0; } double nsl_fit_model_exp_pow_param_deriv(int param, double x, double s, double mu, double b, double a, double weight) { double norm = weight/2./s/gsl_sf_gamma(1.+1./b), y = (x-mu)/s, efactor = exp(-pow(fabs(y), b)); if (param == 0) return norm * a/s * efactor * (b * y * pow(fabs(1./y), 1.-b) * GSL_SIGN(y) - 1.); if (param == 1) return norm * a*b/s * efactor * pow(fabs(y), b-1.) * GSL_SIGN(y); if (param == 2) return norm * a/b * gsl_sf_gamma(1.+1./b)/gsl_sf_gamma(1./b) * efactor * (gsl_sf_psi(1.+1./b) - gsl_pow_2(b) * pow(fabs(y), b) * log(fabs(y))); if (param == 3) return norm * efactor; return 0; } double nsl_fit_model_maxwell_param_deriv(int param, double x, double a, double c, double weight) { double a2 = a*a, a3 = a*a2, norm = weight*sqrt(2./M_PI)/a3, x2 = x*x, efactor = exp(-x2/2./a2); if (param == 0) return c * norm * x2*(x2-3.*a2)/a3 * efactor; if (param == 1) return norm * x2 * efactor; return 0; } double nsl_fit_model_poisson_param_deriv(int param, double x, double l, double A, double weight) { double norm = weight*pow(l, x)/gsl_sf_gamma(x+1.); if (param == 0) return A/l * norm *(x-l)*exp(-l); if (param == 1) return norm * exp(-l); return 0; } double nsl_fit_model_lognormal_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight/sqrt(2.*M_PI)/(x*s), y = log(x)-mu, efactor = exp(-(y/s)*(y/s)/2.); if (param == 0) return A * norm * (y*y - s*s) * efactor; if (param == 1) return A * norm * y/(s*s) * efactor; if (param == 2) return norm * efactor; return 0; } double nsl_fit_model_gamma_param_deriv(int param, double x, double t, double k, double A, double weight) { double factor = weight*pow(x, k-1.)/pow(t, k)/gsl_sf_gamma(k), efactor = exp(-x/t); if (param == 0) return A * factor/t * (x/t-k) * efactor; if (param == 1) return A * factor * (log(x/t) - gsl_sf_psi(k)) * efactor; if (param == 2) return factor * efactor; return 0; } double nsl_fit_model_flat_param_deriv(int param, double x, double a, double b, double A, double weight) { if (x < a || x > b) return 0; if (param == 0) return weight * A/gsl_pow_2(a-b); if (param == 1) return - weight * A/gsl_pow_2(a-b); if (param == 2) return weight/(b-a); return 0; } double nsl_fit_model_rayleigh_param_deriv(int param, double x, double s, double A, double weight) { double y=x/s, norm = weight*y/s, efactor = exp(-y*y/2.); if (param == 0) return A*y/(s*s) * (y*y-2.)*efactor; if (param == 1) return norm * efactor; return 0; } double nsl_fit_model_rayleigh_tail_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight*x/(s*s), y = (mu*mu - x*x)/2./(s*s); if (param == 0) return -2. * A * norm/s * (1. + y) * exp(y); if (param == 1) return A * mu * norm/(s*s) * exp(y); if (param == 2) return norm * exp(y); return 0; } double nsl_fit_model_levy_param_deriv(int param, double x, double g, double mu, double A, double weight) { double y=x-mu, norm = weight*sqrt(g/(2.*M_PI))/pow(y, 1.5), efactor = exp(-g/2./y); if (param == 0) return A/2.*norm/g/y * (y - g) * efactor; if (param == 1) return A/2.*norm/y/y * (3.*y - g) * efactor; if (param == 2) return norm * efactor; return 0; } double nsl_fit_model_landau_param_deriv(int param, double x, double weight) { if (param == 0) return weight * gsl_ran_landau_pdf(x); return 0; } double nsl_fit_model_chi_square_param_deriv(int param, double x, double n, double A, double weight) { double y=n/2., norm = weight*pow(x, y-1.)/pow(2., y)/gsl_sf_gamma(y), efactor = exp(-x/2.); if (param == 0) return A/2. * norm * (log(x/2.) - gsl_sf_psi(y)) * efactor; if (param == 1) return norm * efactor; return 0; } double nsl_fit_model_students_t_param_deriv(int param, double x, double n, double A, double weight) { if (param == 0) return weight * A * gsl_sf_gamma((n+1.)/2.)/2./pow(n, 1.5)/sqrt(M_PI)/gsl_sf_gamma(n/2.) * pow(1.+x*x/n, - (n+3.)/2.) * (x*x - 1. - (n+x*x)*log(1.+x*x/n) + (n+x*x)*(gsl_sf_psi((n+1.)/2.) - gsl_sf_psi(n/2.)) ) ; if (param == 1) return weight * gsl_ran_tdist_pdf(x, n); return 0; } double nsl_fit_model_fdist_param_deriv(int param, double x, double n1, double n2, double A, double weight) { double norm = weight * gsl_sf_gamma((n1+n2)/2.)/gsl_sf_gamma(n1/2.)/gsl_sf_gamma(n2/2.) * pow(n1, n1/2.) * pow(n2, n2/2.) * pow(x, n1/2.-1.); double y = n2+n1*x; if (param == 0) return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n2*(1.-x) + y*(log(n1) + log(x) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n1/2.))); if (param == 1) return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n1*(x-1.) + y*(log(n2) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n2/2.))); if (param == 2) return weight * gsl_ran_fdist_pdf(x, n1, n2); return 0; } double nsl_fit_model_beta_param_deriv(int param, double x, double a, double b, double A, double weight) { double norm = weight * A * gsl_sf_gamma(a+b)/gsl_sf_gamma(a)/gsl_sf_gamma(b) * pow(x, a-1.) * pow(1.-x, b-1.); if (param == 0) return norm * (log(x) - gsl_sf_psi(a) + gsl_sf_psi(a+b)); if (param == 1) return norm * (log(1.-x) - gsl_sf_psi(b) + gsl_sf_psi(a+b)); if (param == 2) return weight * gsl_ran_beta_pdf(x, a, b); return 0; } double nsl_fit_model_pareto_param_deriv(int param, double x, double a, double b, double A, double weight) { if (x < b) return 0; double norm = weight * A; if (param == 0) return norm * pow(b/x, a) * (1. + a * log(b/x))/x; if (param == 1) return norm * a*a * pow(b/x, a-1.)/x/x; if (param == 2) return weight * gsl_ran_pareto_pdf(x, a, b); return 0; } double nsl_fit_model_weibull_param_deriv(int param, double x, double k, double l, double mu, double A, double weight) { double y = (x-mu)/l, z = pow(y, k), efactor = exp(-z); if (param == 0) return weight * A/l * z/y*(k*log(y)*(1.-z) + 1.) * efactor; if (param == 1) return weight * A*k*k/l/l * z/y*(z-1.) * efactor; if (param == 2) return weight * A*k/l/l * z/y/y*(k*z + 1. - k) * efactor; if (param == 3) return weight * k/l * z/y * efactor; return 0; } double nsl_fit_model_frechet_param_deriv(int param, double x, double g, double mu, double s, double A, double weight) { double y = (x-mu)/s, efactor = exp(-pow(y, -g)); if (param == 0) return weight * A/s * pow(y, -2.*g-1.) * (g*log(y)*(1.-pow(y, g))+pow(y, g)) * efactor; if (param == 1) return A * weight * g/(s*s)*pow(y, -g-2.) * (g+1.-g*pow(y, -g)) * efactor; if (param == 2) return A * weight * gsl_pow_2(g/s)*pow(y, -2.*g-1.) * (pow(y, g)-1.) * efactor; if (param == 3) return g * weight/s * pow(y, -g-1.) * efactor; return 0; } double nsl_fit_model_gumbel1_param_deriv(int param, double x, double s, double b, double mu, double A, double weight) { double norm = weight/s, y = (x-mu)/s, efactor = exp(-y - b*exp(-y)); if (param == 0) return A/s * norm * (y - 1. - b*exp(-y)) * efactor; if (param == 1) return -A * norm * exp(-y) * efactor; if (param == 2) return A/s * norm * (1. - b*exp(-y)) * efactor; if (param == 3) return norm * efactor; return 0; } double nsl_fit_model_gumbel2_param_deriv(int param, double x, double a, double b, double mu, double A, double weight) { double y = x - mu, norm = A * weight * exp(-b * pow(y, -a)); if (param == 0) return norm * b * pow(y, -1. -2.*a) * (pow(y, a) -a*(pow(y, a)-b)*log(y)); if (param == 1) return norm * a * pow(y, -1. -2.*a) * (pow(y, a) - b); if (param == 2) return norm * a * b * pow(y, -2.*(a + 1.)) * ((1. + a)*pow(y, a) - a*b); if (param == 3) return weight * gsl_ran_gumbel2_pdf(y, a, b); return 0; } double nsl_fit_model_binomial_param_deriv(int param, double k, double p, double n, double A, double weight) { if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) return 0; k = round(k); n = round(n); - double norm = weight * gsl_sf_fact(n)/gsl_sf_fact(n-k)/gsl_sf_fact(k); + double norm = weight * gsl_sf_fact((unsigned int)n)/gsl_sf_fact((unsigned int)(n-k))/gsl_sf_fact((unsigned int)(k)); if (param == 0) return A * norm * pow(p, k-1.) * pow(1.-p, n-k-1.) * (k-n*p); if (param == 1) return A * norm * pow(p, k) * pow(1.-p, n-k) * (log(1.-p) + gsl_sf_psi(n+1.) - gsl_sf_psi(n-k+1.)); if (param == 2) - return weight * gsl_ran_binomial_pdf(k, p, n); + return weight * gsl_ran_binomial_pdf((unsigned int)k, p, (unsigned int)n); return 0; } double nsl_fit_model_negative_binomial_param_deriv(int param, double k, double p, double n, double A, double weight) { if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) return 0; double norm = A * weight * gsl_sf_gamma(n+k)/gsl_sf_gamma(k+1.)/gsl_sf_gamma(n); if (param == 0) return - norm * pow(p, n-1.) * pow(1.-p, k-1.) * (n*(p-1.) + k*p); if (param == 1) return norm * pow(p, n) * pow(1.-p, k) * (log(p) - gsl_sf_psi(n) + gsl_sf_psi(n+k)); if (param == 2) - return weight * gsl_ran_negative_binomial_pdf(k, p, n); + return weight * gsl_ran_negative_binomial_pdf((unsigned int)k, p, n); return 0; } double nsl_fit_model_pascal_param_deriv(int param, double k, double p, double n, double A, double weight) { return nsl_fit_model_negative_binomial_param_deriv(param, k, p, round(n), A, weight); } double nsl_fit_model_geometric_param_deriv(int param, double k, double p, double A, double weight) { if (param == 0) return A * weight * pow(1.-p, k-2.) * (1.-k*p); if (param == 1) - return weight * gsl_ran_geometric_pdf(k, p); + return weight * gsl_ran_geometric_pdf((unsigned int)k, p); return 0; } double nsl_fit_model_hypergeometric_param_deriv(int param, double k, double n1, double n2, double t, double A, double weight) { if (t > n1 + n2) return 0; - double norm = weight * gsl_ran_hypergeometric_pdf(k, n1, n2, t); + double norm = weight * gsl_ran_hypergeometric_pdf((unsigned int)k, (unsigned int)n1, (unsigned int)n2, (unsigned int)t); if (param == 0) return A * norm * (gsl_sf_psi(n1+1.) - gsl_sf_psi(n1-k+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); if (param == 1) return A * norm * (gsl_sf_psi(n2+1.) - gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); if (param == 2) return A * norm * (gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2-t+1.) - gsl_sf_psi(t-k+1.) + gsl_sf_psi(t+1.)); if (param == 3) return norm; return 0; } double nsl_fit_model_logarithmic_param_deriv(int param, double k, double p, double A, double weight) { if (param == 0) return A * weight * pow(1.-p, k-2.) * (1.-k*p); if (param == 1) - return weight * gsl_ran_logarithmic_pdf(k, p); + return weight * gsl_ran_logarithmic_pdf((unsigned int)k, p); return 0; } double nsl_fit_model_sech_dist_param_deriv(int param, double x, double s, double mu, double A, double weight) { double norm = weight/2./s, y = M_PI/2.*(x-mu)/s; if (param == 0) return -A/s * norm * (y*tanh(y)+1.)/cosh(y); if (param == 1) return A*M_PI/2./s * norm * tanh(y)/cosh(y); if (param == 2) return norm * 1./cosh(y); return 0; } diff --git a/src/backend/nsl/nsl_interp.c b/src/backend/nsl/nsl_interp.c index 3e151948e..a5d97170d 100644 --- a/src/backend/nsl/nsl_interp.c +++ b/src/backend/nsl/nsl_interp.c @@ -1,84 +1,84 @@ /*************************************************************************** File : nsl_interp.c Project : LabPlot Description : NSL interpolation functions -------------------------------------------------------------------- Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_interp.h" #include "nsl_common.h" const char* nsl_interp_type_name[] = { i18n("linear"), i18n("polynomial"), i18n("cubic spline (natural)"), i18n("cubic spline (periodic)"), i18n("Akima-spline (natural)"), i18n("Akima-spline (periodic)"), i18n("Steffen spline"), i18n("cosine"), i18n("exponential"), i18n("piecewise cubic Hermite (PCH)"), i18n("rational functions") }; const char* nsl_interp_pch_variant_name[] = { i18n("finite differences"), i18n("Catmull-Rom"), i18n("cardinal"), i18n("Kochanek-Bartels (TCB)")}; const char* nsl_interp_evaluate_name[] = { i18n("function"), i18n("derivative"), i18n("second derivative"), i18n("integral")}; int nsl_interp_ratint(double *x, double *y, int n, double xn, double *v, double *dv) { - int i,j,a=0,b=n-1; + int i,j,a = 0,b = n-1; while (b-a > 1) { /* find interval using bisection */ - j = floor((a+b)/2.); + j = (int)floor((a+b)/2.); if (x[j] > xn) b = j; else a = j; } int ns=a; /* nearest index */ if (fabs(xn-x[a]) > fabs(xn-x[b])) ns=b; if (xn == x[ns]) { /* exact point */ *v = y[ns]; *dv = 0; return 1; } double *c = (double*)malloc(n*sizeof(double)); double *d = (double*)malloc(n*sizeof(double)); for (i=0; i < n; i++) c[i] = d[i] = y[i]; *v = y[ns--]; double t, dd; int m; for (m=1; m < n; m++) { for (i=0; i < n-m; i++) { t = (x[i]-xn)*d[i]/(x[i+m]-xn); dd = t-c[i+1]; if (dd == 0.0) /* pole */ dd += DBL_MIN; dd = (c[i+1]-d[i])/dd; d[i] = c[i+1]*dd; c[i] = t*dd; } *dv = (2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]); *v += *dv; } free(c); free(d); return 0; } diff --git a/src/backend/nsl/nsl_sf_basic.h b/src/backend/nsl/nsl_sf_basic.h index d649b24a5..14cd60aee 100644 --- a/src/backend/nsl/nsl_sf_basic.h +++ b/src/backend/nsl/nsl_sf_basic.h @@ -1,186 +1,186 @@ /*************************************************************************** File : nsl_sf_basic.h Project : LabPlot Description : NSL special basic functions -------------------------------------------------------------------- Copyright : (C) 2017 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * 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 NSL_SF_BASIC_H #define NSL_SF_BASIC_H #include #include /* stdlib.h */ double nsl_sf_rand() { return rand(); } double nsl_sf_random() { return random(); } double nsl_sf_drand() { return random()/(double)RAND_MAX; } /* math.h */ double nsl_sf_sgn(double x) { #ifndef _WIN32 return copysign(1.0, x); #else if (x > 0) return 1; else if (x < 0) return -1; else return 0; #endif } double nsl_sf_theta(double x) { if (x >= 0) return 1; else return 0; } /* missing trig. functions */ double nsl_sf_sec(double x) { return 1./cos(x); } double nsl_sf_csc(double x) { return 1./sin(x); } double nsl_sf_cot(double x) { return 1./tan(x); } double nsl_sf_asec(double x) { return acos(1./x); } double nsl_sf_acsc(double x) { return asin(1./x); } double nsl_sf_acot(double x) { if (x > 0) return atan(1./x); else return atan(1./x) + M_PI; } double nsl_sf_sech(double x) { return 1./cosh(x); } double nsl_sf_csch(double x) { return 1./sinh(x); } double nsl_sf_coth(double x) { return 1./tanh(x); } double nsl_sf_asech(double x) { return gsl_acosh(1./x); } double nsl_sf_acsch(double x) { return gsl_asinh(1./x); } double nsl_sf_acoth(double x) { return gsl_atanh(1./x); } /* wrapper for GSL functions with integer parameters */ #define MODE GSL_PREC_DOUBLE /* mathematical functions */ -double nsl_sf_ldexp(double x, double expo) { return gsl_ldexp(x, round(expo)); } -double nsl_sf_powint(double x, double n) { return gsl_sf_pow_int(x, round(n)); } +double nsl_sf_ldexp(double x, double expo) { return gsl_ldexp(x, (int)round(expo)); } +double nsl_sf_powint(double x, double n) { return gsl_sf_pow_int(x, (int)round(n)); } /* Airy functions */ double nsl_sf_airy_Ai(double x) { return gsl_sf_airy_Ai(x, MODE); } double nsl_sf_airy_Bi(double x) { return gsl_sf_airy_Bi(x, MODE); } double nsl_sf_airy_Ais(double x) { return gsl_sf_airy_Ai_scaled(x, MODE); } double nsl_sf_airy_Bis(double x) { return gsl_sf_airy_Bi_scaled(x, MODE); } double nsl_sf_airy_Aid(double x) { return gsl_sf_airy_Ai_deriv(x, MODE); } double nsl_sf_airy_Bid(double x) { return gsl_sf_airy_Bi_deriv(x, MODE); } double nsl_sf_airy_Aids(double x) { return gsl_sf_airy_Ai_deriv_scaled(x, MODE); } double nsl_sf_airy_Bids(double x) { return gsl_sf_airy_Bi_deriv_scaled(x, MODE); } -double nsl_sf_airy_0_Ai(double s) { return gsl_sf_airy_zero_Ai(round(s)); } -double nsl_sf_airy_0_Bi(double s) { return gsl_sf_airy_zero_Bi(round(s)); } -double nsl_sf_airy_0_Aid(double s) { return gsl_sf_airy_zero_Ai_deriv(round(s)); } -double nsl_sf_airy_0_Bid(double s) { return gsl_sf_airy_zero_Bi_deriv(round(s)); } +double nsl_sf_airy_0_Ai(double s) { return gsl_sf_airy_zero_Ai((unsigned int)round(s)); } +double nsl_sf_airy_0_Bi(double s) { return gsl_sf_airy_zero_Bi((unsigned int)round(s)); } +double nsl_sf_airy_0_Aid(double s) { return gsl_sf_airy_zero_Ai_deriv((unsigned int)round(s)); } +double nsl_sf_airy_0_Bid(double s) { return gsl_sf_airy_zero_Bi_deriv((unsigned int)round(s)); } /* Bessel functions */ -double nsl_sf_bessel_Jn(double n, double x) { return gsl_sf_bessel_Jn(round(n), x); } -double nsl_sf_bessel_Yn(double n, double x) { return gsl_sf_bessel_Yn(round(n), x); } -double nsl_sf_bessel_In(double n, double x) { return gsl_sf_bessel_In(round(n), x); } -double nsl_sf_bessel_Ins(double n, double x) { return gsl_sf_bessel_In_scaled(round(n), x); } -double nsl_sf_bessel_Kn(double n, double x) { return gsl_sf_bessel_Kn(round(n), x); } -double nsl_sf_bessel_Kns(double n, double x) { return gsl_sf_bessel_Kn_scaled(round(n), x); } -double nsl_sf_bessel_jl(double l, double x) { return gsl_sf_bessel_jl(round(l), x); } -double nsl_sf_bessel_yl(double l, double x) { return gsl_sf_bessel_yl(round(l), x); } -double nsl_sf_bessel_ils(double l, double x) { return gsl_sf_bessel_il_scaled(round(l), x); } -double nsl_sf_bessel_kls(double l, double x) { return gsl_sf_bessel_kl_scaled(round(l), x); } -double nsl_sf_bessel_0_J0(double s) { return gsl_sf_bessel_zero_J0(round(s)); } -double nsl_sf_bessel_0_J1(double s) { return gsl_sf_bessel_zero_J1(round(s)); } -double nsl_sf_bessel_0_Jnu(double nu, double s) { return gsl_sf_bessel_zero_Jnu(nu, round(s)); } - -double nsl_sf_hydrogenicR(double n, double l, double z, double r) { return gsl_sf_hydrogenicR(round(n), round(l), z, r); } +double nsl_sf_bessel_Jn(double n, double x) { return gsl_sf_bessel_Jn((int)round(n), x); } +double nsl_sf_bessel_Yn(double n, double x) { return gsl_sf_bessel_Yn((int)round(n), x); } +double nsl_sf_bessel_In(double n, double x) { return gsl_sf_bessel_In((int)round(n), x); } +double nsl_sf_bessel_Ins(double n, double x) { return gsl_sf_bessel_In_scaled((int)round(n), x); } +double nsl_sf_bessel_Kn(double n, double x) { return gsl_sf_bessel_Kn((int)round(n), x); } +double nsl_sf_bessel_Kns(double n, double x) { return gsl_sf_bessel_Kn_scaled((int)round(n), x); } +double nsl_sf_bessel_jl(double l, double x) { return gsl_sf_bessel_jl((int)round(l), x); } +double nsl_sf_bessel_yl(double l, double x) { return gsl_sf_bessel_yl((int)round(l), x); } +double nsl_sf_bessel_ils(double l, double x) { return gsl_sf_bessel_il_scaled((int)round(l), x); } +double nsl_sf_bessel_kls(double l, double x) { return gsl_sf_bessel_kl_scaled((int)round(l), x); } +double nsl_sf_bessel_0_J0(double s) { return gsl_sf_bessel_zero_J0((unsigned int)round(s)); } +double nsl_sf_bessel_0_J1(double s) { return gsl_sf_bessel_zero_J1((unsigned int)round(s)); } +double nsl_sf_bessel_0_Jnu(double nu, double s) { return gsl_sf_bessel_zero_Jnu(nu, (unsigned int)round(s)); } + +double nsl_sf_hydrogenicR(double n, double l, double z, double r) { return gsl_sf_hydrogenicR((int)round(n), (int)round(l), z, r); } /* elliptic integrals */ double nsl_sf_ellint_Kc(double x) { return gsl_sf_ellint_Kcomp(x, MODE); } double nsl_sf_ellint_Ec(double x) { return gsl_sf_ellint_Ecomp(x, MODE); } double nsl_sf_ellint_Pc(double x, double n) { return gsl_sf_ellint_Pcomp(x, n, MODE); } double nsl_sf_ellint_F(double phi, double k) { return gsl_sf_ellint_F(phi, k, MODE); } double nsl_sf_ellint_E(double phi, double k) { return gsl_sf_ellint_E(phi, k, MODE); } double nsl_sf_ellint_P(double phi, double k, double n) { return gsl_sf_ellint_P(phi, k, n, MODE); } double nsl_sf_ellint_D(double phi, double k) { #if GSL_MAJOR_VERSION >= 2 return gsl_sf_ellint_D(phi,k,MODE); #else return gsl_sf_ellint_D(phi,k,0.0,MODE); #endif } double nsl_sf_ellint_RC(double x, double y) { return gsl_sf_ellint_RC(x, y, MODE); } double nsl_sf_ellint_RD(double x, double y, double z) { return gsl_sf_ellint_RD(x, y, z, MODE); } double nsl_sf_ellint_RF(double x, double y, double z) { return gsl_sf_ellint_RF(x, y, z, MODE); } double nsl_sf_ellint_RJ(double x, double y, double z, double p) { return gsl_sf_ellint_RJ(x, y, z, p, MODE); } -double nsl_sf_exprel_n(double n, double x) { return gsl_sf_exprel_n(round(n), x); } -double nsl_sf_fermi_dirac_int(double j, double x) { return gsl_sf_fermi_dirac_int(round(j), x); } +double nsl_sf_exprel_n(double n, double x) { return gsl_sf_exprel_n((int)round(n), x); } +double nsl_sf_fermi_dirac_int(double j, double x) { return gsl_sf_fermi_dirac_int((int)round(j), x); } /* Gamma */ -double nsl_sf_fact(double n) { return gsl_sf_fact(round(n)); } -double nsl_sf_doublefact(double n) { return gsl_sf_doublefact(round(n)); } -double nsl_sf_lnfact(double n) { return gsl_sf_lnfact(round(n)); } -double nsl_sf_lndoublefact(double n) { return gsl_sf_lndoublefact(round(n)); } -double nsl_sf_choose(double n, double m) { return gsl_sf_choose(round(n), round(m)); } -double nsl_sf_lnchoose(double n, double m) { return gsl_sf_lnchoose(round(n), round(m)); } -double nsl_sf_taylorcoeff(double n, double x) { return gsl_sf_taylorcoeff(round(n), x); } +double nsl_sf_fact(double n) { return gsl_sf_fact((unsigned int)round(n)); } +double nsl_sf_doublefact(double n) { return gsl_sf_doublefact((unsigned int)round(n)); } +double nsl_sf_lnfact(double n) { return gsl_sf_lnfact((unsigned int)round(n)); } +double nsl_sf_lndoublefact(double n) { return gsl_sf_lndoublefact((unsigned int)round(n)); } +double nsl_sf_choose(double n, double m) { return gsl_sf_choose((unsigned int)round(n), (unsigned int)round(m)); } +double nsl_sf_lnchoose(double n, double m) { return gsl_sf_lnchoose((unsigned int)round(n), (unsigned int)round(m)); } +double nsl_sf_taylorcoeff(double n, double x) { return gsl_sf_taylorcoeff((int)round(n), x); } -double nsl_sf_gegenpoly_n(double n, double l, double x) { return gsl_sf_gegenpoly_n(round(n), l, x); } +double nsl_sf_gegenpoly_n(double n, double l, double x) { return gsl_sf_gegenpoly_n((int)round(n), l, x); } #if (GSL_MAJOR_VERSION > 2) || (GSL_MAJOR_VERSION == 2) && (GSL_MINOR_VERSION >= 4) double nsl_sf_hermite_prob(double n, double x) { return gsl_sf_hermite_prob(round(n), x); } double nsl_sf_hermite_phys(double n, double x) { return gsl_sf_hermite_phys(round(n), x); } double nsl_sf_hermite_func(double n, double x) { return gsl_sf_hermite_func(round(n), x); } double nsl_sf_hermite_prob_der(double m, double n, double x) { return gsl_sf_hermite_prob_der(round(m), round(n), x); } double nsl_sf_hermite_phys_der(double m, double n, double x) { return gsl_sf_hermite_phys_der(round(m), round(n), x); } double nsl_sf_hermite_func_der(double m, double n, double x) { return gsl_sf_hermite_func_der(round(m), round(n), x); } #endif -double nsl_sf_hyperg_1F1i(double m, double n, double x) { return gsl_sf_hyperg_1F1_int(round(m), round(n), x); } -double nsl_sf_hyperg_Ui(double m, double n, double x) { return gsl_sf_hyperg_U_int(round(m), round(n), x); } -double nsl_sf_laguerre_n(double n, double a, double x) { return gsl_sf_laguerre_n(round(n), a, x); } +double nsl_sf_hyperg_1F1i(double m, double n, double x) { return gsl_sf_hyperg_1F1_int((int)round(m), (int)round(n), x); } +double nsl_sf_hyperg_Ui(double m, double n, double x) { return gsl_sf_hyperg_U_int((int)round(m), (int)round(n), x); } +double nsl_sf_laguerre_n(double n, double a, double x) { return gsl_sf_laguerre_n((int)round(n), a, x); } -double nsl_sf_legendre_Pl(double l, double x) { return gsl_sf_legendre_Pl(round(l), x); } -double nsl_sf_legendre_Ql(double l, double x) { return gsl_sf_legendre_Ql(round(l), x); } -double nsl_sf_legendre_Plm(double l, double m, double x) { return gsl_sf_legendre_Plm(round(l), round(m), x); } -double nsl_sf_legendre_sphPlm(double l, double m, double x) { return gsl_sf_legendre_sphPlm(round(l), round(m), x); } -double nsl_sf_conicalP_sphreg(double l, double L, double x) { return gsl_sf_conicalP_sph_reg(round(l), L, x); } -double nsl_sf_conicalP_cylreg(double m, double l, double x) { return gsl_sf_conicalP_sph_reg(round(m), l, x); } -double nsl_sf_legendre_H3d(double l, double L, double e) { return gsl_sf_legendre_H3d(round(l), L, e); } +double nsl_sf_legendre_Pl(double l, double x) { return gsl_sf_legendre_Pl((int)round(l), x); } +double nsl_sf_legendre_Ql(double l, double x) { return gsl_sf_legendre_Ql((int)round(l), x); } +double nsl_sf_legendre_Plm(double l, double m, double x) { return gsl_sf_legendre_Plm((int)round(l), (int)round(m), x); } +double nsl_sf_legendre_sphPlm(double l, double m, double x) { return gsl_sf_legendre_sphPlm((int)round(l), (int)round(m), x); } +double nsl_sf_conicalP_sphreg(double l, double L, double x) { return gsl_sf_conicalP_sph_reg((int)round(l), L, x); } +double nsl_sf_conicalP_cylreg(double m, double l, double x) { return gsl_sf_conicalP_sph_reg((int)round(m), l, x); } +double nsl_sf_legendre_H3d(double l, double L, double e) { return gsl_sf_legendre_H3d((int)round(l), L, e); } -double nsl_sf_psiint(double n) { return gsl_sf_psi_int(round(n)); } -double nsl_sf_psi1int(double n) { return gsl_sf_psi_1_int(round(n)); } -double nsl_sf_psin(double n, double x) { return gsl_sf_psi_n(round(n), x); } +double nsl_sf_psiint(double n) { return gsl_sf_psi_int((int)round(n)); } +double nsl_sf_psi1int(double n) { return gsl_sf_psi_1_int((int)round(n)); } +double nsl_sf_psin(double n, double x) { return gsl_sf_psi_n((int)round(n), x); } -double nsl_sf_zetaint(double n) { return gsl_sf_zeta_int(round(n)); } -double nsl_sf_zetam1int(double n) { return gsl_sf_zetam1_int(round(n)); } -double nsl_sf_etaint(double n) { return gsl_sf_eta_int(round(n)); } +double nsl_sf_zetaint(double n) { return gsl_sf_zeta_int((int)round(n)); } +double nsl_sf_zetam1int(double n) { return gsl_sf_zetam1_int((int)round(n)); } +double nsl_sf_etaint(double n) { return gsl_sf_eta_int((int)round(n)); } /* random number distributions */ -double nsl_sf_poisson(double k, double m) { return gsl_ran_poisson_pdf(round(k), m); } -double nsl_sf_bernoulli(double k, double p) { return gsl_ran_bernoulli_pdf(round(k), p); } -double nsl_sf_binomial(double k, double p, double n) { return gsl_ran_binomial_pdf(round(k), p, round(n)); } -double nsl_sf_negative_binomial(double k, double p, double n) { return gsl_ran_negative_binomial_pdf(round(k), p, n); } -double nsl_sf_pascal(double k, double p, double n) { return gsl_ran_pascal_pdf(round(k), p, round(n)); } -double nsl_sf_geometric(double k, double p) { return gsl_ran_geometric_pdf(round(k), p); } +double nsl_sf_poisson(double k, double m) { return gsl_ran_poisson_pdf((unsigned int)round(k), m); } +double nsl_sf_bernoulli(double k, double p) { return gsl_ran_bernoulli_pdf((unsigned int)round(k), p); } +double nsl_sf_binomial(double k, double p, double n) { return gsl_ran_binomial_pdf((unsigned int)round(k), p, (unsigned int)round(n)); } +double nsl_sf_negative_binomial(double k, double p, double n) { return gsl_ran_negative_binomial_pdf((unsigned int)round(k), p, n); } +double nsl_sf_pascal(double k, double p, double n) { return gsl_ran_pascal_pdf((unsigned int)round(k), p, (unsigned int)round(n)); } +double nsl_sf_geometric(double k, double p) { return gsl_ran_geometric_pdf((unsigned int)round(k), p); } double nsl_sf_hypergeometric(double k, double n1, double n2, double t) { - return gsl_ran_hypergeometric_pdf(round(k), round(n1), round(n2), round(t)); + return gsl_ran_hypergeometric_pdf((unsigned int)round(k), (unsigned int)round(n1), (unsigned int)round(n2), (unsigned int)round(t)); } -double nsl_sf_logarithmic(double k, double p) { return gsl_ran_logarithmic_pdf(round(k), p); } +double nsl_sf_logarithmic(double k, double p) { return gsl_ran_logarithmic_pdf((unsigned int)round(k), p); } #endif /* NSL_SF_BASIC_H */ diff --git a/src/backend/nsl/nsl_sf_poly.c b/src/backend/nsl/nsl_sf_poly.c index 2d9e5f5d1..5468840be 100644 --- a/src/backend/nsl/nsl_sf_poly.c +++ b/src/backend/nsl/nsl_sf_poly.c @@ -1,330 +1,338 @@ /*************************************************************************** File : nsl_sf_poly.c Project : LabPlot Description : NSL special polynomial functions -------------------------------------------------------------------- Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_sf_poly.h" #include #include #include /* see https://en.wikipedia.org/wiki/Chebyshev_polynomials */ double nsl_sf_poly_chebyshev_T(int n, double x) { if (fabs(x) <= 1) return cos(n * acos(x)); else if (x > 1) return cosh(n * gsl_acosh(x)); else return pow(-1., n)*cosh(n * gsl_acosh(-x)); } double nsl_sf_poly_chebyshev_U(int n, double x) { double sq = sqrt(x*x - 1); return (pow(x + sq, n + 1) - pow(x - sq, n + 1))/2./sq; } /* from http://www.crbond.com/papers/lopt.pdf */ double nsl_sf_poly_optimal_legendre_L(int n, double x) { if (n < 1 || n > 10) return 0.0; switch (n) { case 1: return x; case 2: return x*x; case 3: return (1. + (-3. + 3.*x)*x)*x; case 4: return (3. + (-8. + 6*x)*x)*x*x; case 5: return (1. + (-8. + (28. + (-40. + 20*x)*x)*x)*x)*x; case 6: return (6. + (-40. + (105. + (-120. + 50.*x)*x)*x)*x)*x*x; case 7: return (1. + (-15. + (105. + (-355. + (615. + (-525. + 175.*x)*x)*x)*x)*x)*x)*x; case 8: return (10. + (-120. + (615. + (-1624. + (2310. + (-1680. + 490.*x)*x)*x)*x)*x)*x)*x*x; case 9: return (1. + (-24. + (276. + (-1624. + (5376. + (-10416. + (11704. + (-7056. + 1764*x)*x)*x)*x)*x)*x)*x)*x)*x; case 10: return (15. + (-280. + (2310. + (-10416. + (27860. + (-45360. + (44100. + (23520. + 5292.*x)*x)*x)*x)*x)*x)*x)*x)*x*x; } return 0.0; } /* * https://en.wikipedia.org/wiki/Bessel_polynomials * using recursion */ COMPLEX nsl_sf_poly_bessel_y(int n, COMPLEX x) { #ifdef _MSC_VER if (n == 0) { COMPLEX z = {1.0, 0.0}; return z; } else if (n == 1) { - COMPLEX z = {creal(x) + 1.0, cimag(x)}; + COMPLEX z; + z.re = creal(x) + 1.0; + z.im = cimag(x); return z; } double factor = 2 * n - 1; COMPLEX z1 = nsl_sf_poly_bessel_y(n - 1, x); COMPLEX z2 = nsl_sf_poly_bessel_y(n - 2, x); - COMPLEX z = {factor * (creal(x)*creal(z1)-cimag(x)*cimag(z1)) + creal(z2), factor * (creal(x)*cimag(z1) - cimag(x)*creal(z1)) + cimag(z2)}; + COMPLEX z; + z.re = factor * (creal(x)*creal(z1)-cimag(x)*cimag(z1)) + creal(z2); + z.im = factor * (creal(x)*cimag(z1) - cimag(x)*creal(z1)) + cimag(z2); return z; #else if (n == 0) return 1.0; else if (n == 1) return 1.0 + x; return (2*n - 1) * x * nsl_sf_poly_bessel_y(n - 1, x) + nsl_sf_poly_bessel_y(n - 2, x); #endif } /* * https://en.wikipedia.org/wiki/Bessel_polynomials * using recursion */ COMPLEX nsl_sf_poly_reversed_bessel_theta(int n, COMPLEX x) { #ifdef _MSC_VER if (n == 0) { COMPLEX z = {1.0, 0.0}; return z; } else if (n == 1) { - COMPLEX z = {creal(x) + 1.0, cimag(x)}; + COMPLEX z; + z.re = creal(x) + 1.0; + z.im = cimag(x); return z; } double factor = 2 * n - 1; COMPLEX z1 = nsl_sf_poly_bessel_y(n - 1, x); COMPLEX z2 = nsl_sf_poly_bessel_y(n - 2, x); double rex2 = creal(x)*creal(x) - cimag(x)*cimag(x); double imx2 = 2.0*creal(x)*cimag(x); - COMPLEX z = {factor * creal(z1) + rex2*creal(z2) - imx2*cimag(z2) , factor * cimag(z1) + rex2*cimag(z2) + imx2*creal(z2)}; + COMPLEX z; + z.re = factor * creal(z1) + rex2*creal(z2) - imx2*cimag(z2); + z.im = factor * cimag(z1) + rex2*cimag(z2) + imx2*creal(z2); return z; #else if (n == 0) return 1.0; else if (n == 1) return 1.0 + x; return (2*n - 1) * nsl_sf_poly_reversed_bessel_theta(n - 1, x) + x * x * nsl_sf_poly_reversed_bessel_theta(n - 2, x); #endif } /***************** interpolating polynomials *************/ double nsl_sf_poly_interp_lagrange_0_int(double *x, double y) { /* rectangle rule */ return (x[1]-x[0])*y; } double nsl_sf_poly_interp_lagrange_1(double v, double *x, double *y) { return (y[0]*(x[1]-v) + y[1]*(v-x[0]))/(x[1]-x[0]); } double nsl_sf_poly_interp_lagrange_1_deriv(double *x, double *y) { return (y[0]-y[1])/(x[1]-x[0]); } double nsl_sf_poly_interp_lagrange_1_int(double *x, double *y) { /* trapezoid rule (2-point) */ return (x[1]-x[0])*(y[0]+y[1])/2.; } double nsl_sf_poly_interp_lagrange_1_absint(double *x, double *y) { double dx = x[1]-x[0]; if (y[0]*y[1] < 0) /* sign change */ return dx*( (fabs(y[0])-fabs(y[1]))/(fabs(y[1]/y[0])+1) + fabs(y[1]) )/2.; if (y[0] < 0 && y[1] < 0) return dx*(fabs(y[0])+fabs(y[1]))/2.; else return dx*(y[0]+y[1])/2.; } double nsl_sf_poly_interp_lagrange_2(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1]; double h12 = h1+h2; return y[0]*(v-x[1])*(v-x[2])/(h1*h12) + y[1]*(x[0]-v)*(v-x[2])/(h1*h2) + y[2]*(x[0]-v)*(x[1]-v)/(h12*h2); } double nsl_sf_poly_interp_lagrange_2_deriv(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1]; double h12 = h1+h2; return y[0]*(2.*v-x[1]-x[2])/(h1*h12) + y[1]*(x[0]-2.*v+x[2])/(h1*h2) + y[2]*(2.*v-x[0]-x[1])/(h12*h2); } double nsl_sf_poly_interp_lagrange_2_deriv2(double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1]; double h12 = h1+h2; return 2.*( y[0]/(h1*h12) - y[1]/(h1*h2) + y[2]/(h12*h2) ); } double nsl_sf_poly_interp_lagrange_2_int(double *x, double *y) { /* Simpson's 1/3 rule for non-uniform grid (3-point) */ double h1 = x[1]-x[0], h2 = x[2]-x[1]; double h12 = h1+h2, h1_2 = h1/h2; return h12/6.*( (2.-1./h1_2)*y[0] + (2.+h1_2+1./h1_2)*y[1] + (2.-h1_2)*y[2] ); } double nsl_sf_poly_interp_lagrange_3(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2]; double h12 = h1+h2, h23 = h2+h3, h13=h12+h3; return y[0]*(v-x[1])*(v-x[2])*(v-x[3])/(h1*h12*h13) + y[1]*(v-x[0])*(v-x[2])*(v-x[3])/(h1*h2*h23) - y[2]*(v-x[0])*(v-x[1])*(v-x[3])/(h12*h2*h3) + y[3]*(v-x[0])*(v-x[1])*(v-x[2])/(h13*h23*h3); } double nsl_sf_poly_interp_lagrange_3_deriv(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], S = x[0]+x[1]+x[2]+x[3]; double h12 = h1+h2, h23 = h2+h3, h13=h12+h3; return -y[0]*(3*v*v-2.*v*(S-x[0])+x[1]*x[2]+(S-x[0]-x[3])*x[3])/(h1*h12*h13) + y[1]*(3*v*v-2.*v*(S-x[1])+(x[0]+x[2])*x[3])/(h1*h2*h23) + y[2]*(3*v*v+x[0]*x[1]-2.*v*(S-x[2])+(x[0]+x[1])*x[3])/(h12*h2*h3) + y[3]*(3*v*v+x[0]*x[1]+(x[0]+x[1])*x[2]-2.*v*(S-x[3]))/(h13*h23*h3); } double nsl_sf_poly_interp_lagrange_3_deriv2(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], S = x[0]+x[1]+x[2]+x[3]; double h12 = h1+h2, h23 = h2+h3, h13=h12+h3; return 2.*( y[0]*(S-3.*v-x[0])/(h1*h12*h13) + y[1]*(3.*v-S+x[1])/(h1*h2*h23) + y[2]*(S-3.*v-x[2])/(h12*h2*h3) + y[3]*(3.*v-S+x[3])/(h13*h23*h3) ); } double nsl_sf_poly_interp_lagrange_3_deriv3(double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2]; double h12 = h1+h2, h23 = h2+h3, h13=h12+h3; return 6.*( -y[0]/(h1*h12*h13) + y[1]/(h1*h2*h23) - y[2]/(h12*h2*h3) +y[3]/(h13*h23*h3) ); } double nsl_sf_poly_interp_lagrange_3_int(double *x, double *y) { /* Simpson's 3/8 rule for non-uniform grid (4-point) */ double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2]; double h12 = h1+h2, h23 = h2+h3, h13 = h12+h3; return h13/12.*( (3.*x[0]*x[0]-4.*x[0]*(x[1]+x[2])+6.*x[1]*x[2]+2.*x[0]*x[3]-2.*(x[1]+x[2])*x[3]+x[3]*x[3])/(h1*h12)*y[0] + h13*h13*(h12-h3)/(h1*h2*h23)*y[1] + h13*h13*(h23-h1)/(h12*h2*h3)*y[2] + (x[0]*x[0]-2.*x[0]*(x[1]+x[2])+6.*x[1]*x[2]+2.*x[0]*x[3]-4.*(x[1]+x[2])*x[3]+3.*x[3]*x[3])/(h23*h3)*y[3] ); } /* 1/2 sum_{i != j}^n (x_i x_j) for n=4 */ double sum_of_2product_combinations_4_2(double a, double b, double c, double d) { return a*(b+c+d) + b*(c+d) + c*d; } /* 1/6 sum_{i != j, j != k, i != k}^n (x_i x_j x_k) for n=4 */ double sum_of_3product_combinations_4_6(double a, double b, double c, double d) { return a*(b*(c+d) + c*d) + b*c*d; } double nsl_sf_poly_interp_lagrange_4(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h13=h12+h3, h24 = h23+h4, h14 = h12+h34; return y[0]*(v-x[1])*(v-x[2])*(v-x[3])*(v-x[4])/(h1*h12*h13*h14) - y[1]*(v-x[0])*(v-x[2])*(v-x[3])*(v-x[4])/(h1*h2*h23*h24) + y[2]*(v-x[0])*(v-x[1])*(v-x[3])*(v-x[4])/(h12*h2*h3*h34) - y[3]*(v-x[0])*(v-x[1])*(v-x[2])*(v-x[4])/(h13*h23*h3*h4) + y[4]*(v-x[0])*(v-x[1])*(v-x[2])*(v-x[3])/(h14*h24*h34*h4); } double nsl_sf_poly_interp_lagrange_4_deriv(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], S = x[0]+x[1]+x[2]+x[3]+x[4]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h13=h12+h3, h24 = h23+h4, h14 = h12+h34; return y[0]*(4.*v*v*v - 3.*v*v*(S-x[0]) - sum_of_3product_combinations_4_6(x[1],x[2],x[3],x[4]) + 2.*v*sum_of_2product_combinations_4_2(x[1],x[2],x[3],x[4]) )/(h1*h12*h13*h14) - y[1]*(4.*v*v*v - 3.*v*v*(S-x[1]) - sum_of_3product_combinations_4_6(x[0],x[2],x[3],x[4]) + 2.*v*sum_of_2product_combinations_4_2(x[0],x[2],x[3],x[4]) )/(h1*h2*h23*h24) + y[2]*(4.*v*v*v - 3.*v*v*(S-x[2]) - sum_of_3product_combinations_4_6(x[0],x[1],x[3],x[4]) + 2.*v*sum_of_2product_combinations_4_2(x[0],x[1],x[3],x[4]) )/(h12*h2*h3*h34) - y[3]*(4.*v*v*v - 3.*v*v*(S-x[3]) - sum_of_3product_combinations_4_6(x[0],x[1],x[2],x[4]) + 2.*v*sum_of_2product_combinations_4_2(x[0],x[1],x[2],x[4]) )/(h13*h23*h3*h4) + y[4]*(4.*v*v*v - 3.*v*v*(S-x[4]) - sum_of_3product_combinations_4_6(x[0],x[1],x[2],x[3]) + 2.*v*sum_of_2product_combinations_4_2(x[0],x[1],x[2],x[3]) )/(h14*h24*h34*h4); } double nsl_sf_poly_interp_lagrange_4_deriv2(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], S = x[0]+x[1]+x[2]+x[3]+x[4]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h13=h12+h3, h24 = h23+h4, h14 = h12+h34; return 2.*( y[0]*(6.*v*v-3.*v*(S-x[0]) + sum_of_2product_combinations_4_2(x[1],x[2],x[3],x[4]))/(h1*h12*h13*h14) - y[1]*(6.*v*v-3.*v*(S-x[1]) + sum_of_2product_combinations_4_2(x[0],x[2],x[3],x[4]))/(h1*h2*h23*h24) + y[2]*(6.*v*v-3.*v*(S-x[2]) + sum_of_2product_combinations_4_2(x[0],x[1],x[3],x[4]))/(h12*h2*h3*h34) - y[3]*(6.*v*v-3.*v*(S-x[3]) + sum_of_2product_combinations_4_2(x[0],x[1],x[2],x[4]))/(h13*h23*h3*h4) + y[4]*(6.*v*v-3.*v*(S-x[4]) + sum_of_2product_combinations_4_2(x[0],x[1],x[2],x[3]))/(h14*h24*h34*h4) ); } double nsl_sf_poly_interp_lagrange_4_deriv3(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], S = x[0]+x[1]+x[2]+x[3]+x[4]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h13=h12+h3, h24 = h23+h4, h14 = h12+h34; return 6.*( y[0]*(4.*v-S+x[0])/(h1*h12*h13*h14) + y[1]*(S-4.*v-x[1])/(h1*h2*h23*h24) + y[2]*(4.*v-S+x[2])/(h12*h2*h3*h34) + y[3]*(S-4.*v-x[3])/(h13*h23*h3*h4) + y[4]*(4.*v-S+x[4])/(h14*h24*h34*h4) ); } double nsl_sf_poly_interp_lagrange_4_deriv4(double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h13=h12+h3, h24 = h23+h4, h14 = h12+h34; return 24.*( y[0]/(h1*h12*h13*h14) - y[1]/(h1*h2*h23*h24) + y[2]/(h12*h2*h3*h34)- y[3]/(h13*h23*h3*h4) + y[4]/(h14*h24*h34*h4) ); } /* 1/2 sum_{i != j}^n (x_i x_j) for n=6 */ double sum_of_product_combinations_6_2(double a, double b, double c, double d, double e, double f) { return a*(b+c+d+e+f) + b*(c+d+e+f) + c*(d+e+f) + d*(e+f) + e*f; } double nsl_sf_poly_interp_lagrange_6_deriv4(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], h5 = x[5]-x[4], h6 = x[6]-x[5]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h45 = h4+h5, h56 = h5+h6; double h13 = h12+h3, h24 = h23+h4, h35 = h34+h5, h46 = h45+h6, h14 = h13+h4, h25 = h24+h5, h36 = h35+h6; double h15 = h14+h5, h26 = h25 + h6, h16 = h15+h6, S = x[0]+x[1]+x[2]+x[3]+x[4]+x[5]+x[6]; return 24.*( y[0]*(15.*v*v-5.*v*(S-x[0]) + sum_of_product_combinations_6_2(x[1],x[2],x[3],x[4],x[5],x[6]))/(h1*h12*h13*h14*h15*h16) - y[1]*(15.*v*v-5.*v*(S-x[1]) + sum_of_product_combinations_6_2(x[0],x[2],x[3],x[4],x[5],x[6]))/(h1*h2*h23*h24*h25*h26) + y[2]*(15.*v*v-5.*v*(S-x[2]) + sum_of_product_combinations_6_2(x[0],x[1],x[3],x[4],x[5],x[6]))/(h12*h2*h3*h34*h35*h36) - y[3]*(15.*v*v-5.*v*(S-x[3]) + sum_of_product_combinations_6_2(x[0],x[1],x[2],x[4],x[5],x[6]))/(h13*h23*h3*h4*h45*h46) + y[4]*(15.*v*v-5.*v*(S-x[4]) + sum_of_product_combinations_6_2(x[0],x[1],x[2],x[3],x[5],x[6]))/(h14*h24*h34*h4*h5*h56) - y[5]*(15.*v*v-5.*v*(S-x[5]) + sum_of_product_combinations_6_2(x[0],x[1],x[2],x[3],x[4],x[6]))/(h15*h25*h35*h45*h5*h6) + y[6]*(15.*v*v-5.*v*(S-x[6]) + sum_of_product_combinations_6_2(x[0],x[1],x[2],x[3],x[4],x[5]))/(h16*h26*h36*h46*h56*h6) ); } double nsl_sf_poly_interp_lagrange_6_deriv5(double v, double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], h5 = x[5]-x[4], h6 = x[6]-x[5]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h45 = h4+h5, h56 = h5+h6; double h13 = h12+h3, h24 = h23+h4, h35 = h34+h5, h46 = h45+h6, h14 = h13+h4, h25 = h24+h5, h36 = h35+h6; double h15 = h14+h5, h26 = h25 + h6, h16 = h15+h6, S = x[0]+x[1]+x[2]+x[3]+x[4]+x[5]+x[6]; return 120.*( y[0]*(6.*v-S+x[0])/(h1*h12*h13*h14*h15*h16) - y[1]*(6.*v-S+x[1])/(h1*h2*h23*h24*h25*h26) + y[2]*(6.*v-S+x[2])/(h12*h2*h3*h34*h35*h36) - y[3]*(6.*v-S+x[3])/(h13*h23*h3*h4*h45*h46) + y[4]*(6.*v-S+x[4])/(h14*h24*h34*h4*h5*h56) - y[5]*(6.*v-S+x[5])/(h15*h25*h35*h45*h5*h6) + y[6]*(6.*v-S+x[6])/(h16*h26*h36*h46*h56*h6) ); } double nsl_sf_poly_interp_lagrange_6_deriv6(double *x, double *y) { double h1 = x[1]-x[0], h2 = x[2]-x[1], h3 = x[3]-x[2], h4 = x[4]-x[3], h5 = x[5]-x[4], h6 = x[6]-x[5]; double h12 = h1+h2, h23 = h2+h3, h34 = h3+h4, h45 = h4+h5, h56 = h5+h6; double h13 = h12+h3, h24 = h23+h4, h35 = h34+h5, h46 = h45+h6, h14 = h13+h4, h25 = h24+h5, h36 = h35+h6; double h15 = h14+h5, h26 = h25 + h6, h16 = h15+h6; return 720.*( y[0]/(h1*h12*h13*h14*h15*h16) - y[1]/(h1*h2*h23*h24*h25*h26) + y[2]/(h12*h2*h3*h34*h35*h36) - y[3]/(h13*h23*h3*h4*h45*h46) + y[4]/(h14*h24*h34*h4*h5*h56) - y[5]/(h15*h25*h35*h45*h5*h6) + y[6]/(h16*h26*h36*h46*h56*h6) ); } diff --git a/src/backend/nsl/nsl_smooth.c b/src/backend/nsl/nsl_smooth.c index b2ef2c1dd..4f5156f5f 100644 --- a/src/backend/nsl/nsl_smooth.c +++ b/src/backend/nsl/nsl_smooth.c @@ -1,544 +1,545 @@ /*************************************************************************** File : nsl_smooth.c Project : LabPlot Description : NSL smooth functions -------------------------------------------------------------------- Copyright : (C) 2010 by Knut Franke (knut.franke@gmx.de) Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) ***************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the Free Software * * Foundation, Inc., 51 Franklin Street, Fifth Floor, * * Boston, MA 02110-1301 USA * * * ***************************************************************************/ #include "nsl_smooth.h" #include "nsl_common.h" #include "nsl_sf_kernel.h" #include "nsl_stats.h" #include #include #include #include /* gsl_sf_choose */ const char* nsl_smooth_type_name[] = { i18n("moving average (central)"), i18n("moving average (lagged)"), i18n("percentile"), i18n("Savitzky-Golay") }; const char* nsl_smooth_pad_mode_name[] = { i18n("none"), i18n("interpolating"), i18n("mirror"), i18n("nearest"), i18n("constant"), i18n("periodic") }; const char* nsl_smooth_weight_type_name[] = { i18n("uniform (rectangular)"), i18n("triangular"), i18n("binomial"), i18n("parabolic (Epanechnikov)"), i18n("quartic (biweight)"), i18n("triweight"), i18n("tricube"), i18n("cosine") }; double nsl_smooth_pad_constant_lvalue = 0.0, nsl_smooth_pad_constant_rvalue = 0.0; int nsl_smooth_moving_average(double *data, unsigned int n, unsigned int points, nsl_smooth_weight_type weight, nsl_smooth_pad_mode mode) { unsigned int i,j; double *result = (double *)malloc(n*sizeof(double)); for (i=0; i(int)n-1) result[i] += w[j]*nsl_smooth_pad_constant_rvalue; else result[i] += w[j]*data[index]; break; case nsl_smooth_pad_periodic: if(index<0) index = n+index; else if(index>(int)n-1) index = index-n; result[i] += w[j]*data[index]; break; } } /*puts("");*/ free(w); } for (i=0; i(int)n-1) values[j] = nsl_smooth_pad_constant_rvalue; else values[j] = data[index]; break; case nsl_smooth_pad_periodic: if(index<0) index = n+index; else if(index>(int)n-1) index = index-n; /*printf(" %d",index);*/ values[j] = data[index]; break; } } /*puts("");*/ /*using type 4 as default */ result[i] = nsl_stats_quantile(values, 1, np, percentile, nsl_stats_quantile_type4); free(values); } for (i=0; i n) { printf("Tried to smooth over more points (points=%u) than given as input (%u).", points, n); return -1; } if (order < 1 || order > points-1) { printf("The polynomial order must be between 1 and %u (%u given).", points-1, order); return -2; } /* Savitzky-Golay coefficient matrix, y' = H y */ gsl_matrix *h = gsl_matrix_alloc(points, points); error = nsl_smooth_savgol_coeff(points, order, h); if (error) { printf("Internal error in Savitzky-Golay algorithm:\n%s",gsl_strerror(error)); gsl_matrix_free(h); return error; } double *result = (double *)malloc(n*sizeof(double)); for (i = 0; i < n; i++) result[i] = 0; /* left edge */ if(mode == nsl_smooth_pad_none) { for (i = 0; i < half; i++) { /*reduce points and order*/ unsigned int rpoints = 2*i+1, rorder = GSL_MIN(order, rpoints-GSL_MIN(rpoints, 2)); gsl_matrix *rh = gsl_matrix_alloc(rpoints, rpoints); error = nsl_smooth_savgol_coeff(rpoints, rorder, rh); if (error) { printf("Internal error in Savitzky-Golay algorithm:\n%s",gsl_strerror(error)); gsl_matrix_free(rh); free(result); return error; } for (k = 0; k < rpoints; k++) result[i] += gsl_matrix_get(rh, i, k) * data[k]; } } else { for (i = 0; i < half; i++) { for (k = 0; k < points; k++) switch(mode) { case nsl_smooth_pad_interp: result[i] += gsl_matrix_get(h, i, k) * data[k]; break; case nsl_smooth_pad_mirror: result[i] += gsl_matrix_get(h, half, k) * data[abs((int)(k+i-half))]; break; case nsl_smooth_pad_nearest: result[i] += gsl_matrix_get(h, half, k) * data[i+k-GSL_MIN(half,i+k)]; break; case nsl_smooth_pad_constant: if (k #include #include #include double nsl_stats_minimum(const double data[], const size_t n, size_t *index) { size_t i; double min = data[0]; if (index != NULL) *index = 0; for (i = 1; i < n; i++) { if (data[i] < min) { min = data[i]; if (index != NULL) *index = i; } } return min; } double nsl_stats_maximum(const double data[], const size_t n, size_t *index) { size_t i; double max = data[0]; if (index != NULL) *index = 0; for (i = 1; i < n; i++) { if (data[i] > max) { max = data[i]; if (index != NULL) *index = i; } } return max; } double nsl_stats_median(double data[], size_t stride, size_t n, nsl_stats_quantile_type type) { gsl_sort(data, stride, n); return nsl_stats_median_sorted(data,stride,n,type); } double nsl_stats_median_sorted(const double sorted_data[], size_t stride, size_t n, nsl_stats_quantile_type type) { return nsl_stats_quantile_sorted(sorted_data,stride,n,0.5,type); } double nsl_stats_median_from_sorted_data(const double sorted_data[], size_t stride, size_t n) { return nsl_stats_median_sorted(sorted_data,stride,n,nsl_stats_quantile_type7); } double nsl_stats_quantile(double data[], size_t stride, size_t n, double p, nsl_stats_quantile_type type) { gsl_sort(data, stride, n); return nsl_stats_quantile_sorted(data,stride,n,p,type); } double nsl_stats_quantile_sorted(const double d[], size_t stride, size_t n, double p, nsl_stats_quantile_type type) { switch(type) { case nsl_stats_quantile_type1: if (p == 0.0) return d[0]; else return d[((int)ceil(n*p)-1)*stride]; case nsl_stats_quantile_type2: if (p == 0.0) return d[0]; else if (p == 1.0) return d[(n-1)*stride]; else return (d[((int)ceil(n*p)-1)*stride]+d[((int)ceil(n*p+1)-1)*stride])/2.; case nsl_stats_quantile_type3: if(p <= 0.5/n) return d[0]; else #ifdef _WIN32 return d[((int)(n*p)-1)*stride]; #else return d[(lrint(n*p)-1)*stride]; #endif case nsl_stats_quantile_type4: if(p < 1./n) return d[0]; else if (p == 1.0) return d[(n-1)*stride]; else { - int i = floor(n*p); + int i = (int)floor(n*p); return d[(i-1)*stride]+(n*p-i)*(d[i*stride]-d[(i-1)*stride]); } case nsl_stats_quantile_type5: if(p < 0.5/n) return d[0]; else if (p >= (n-0.5)/n) return d[(n-1)*stride]; else { - int i = floor(n*p+0.5); + int i = (int)floor(n*p+0.5); return d[(i-1)*stride]+(n*p+0.5-i)*(d[i*stride]-d[(i-1)*stride]); } case nsl_stats_quantile_type6: if(p < 1./(n+1.)) return d[0]; else if (p > n/(n+1.)) return d[(n-1)*stride]; else { - int i = floor((n+1)*p); + int i = (int)floor((n+1)*p); return d[(i-1)*stride]+((n+1)*p-i)*(d[i*stride]-d[(i-1)*stride]); } case nsl_stats_quantile_type7: if (p == 1.0) return d[(n-1)*stride]; else { - int i = floor((n-1)*p+1); + int i = (int)floor((n-1)*p+1); return d[(i-1)*stride]+((n-1)*p+1-i)*(d[i*stride]-d[(i-1)*stride]); } case nsl_stats_quantile_type8: if (p < 2./3./(n+1./3.)) return d[0]; else if (p >= (n-1./3.)/(n+1./3.)) return d[(n-1)*stride]; else { - int i = floor((n+1./3.)*p+1./3.); + int i = (int)floor((n+1./3.)*p+1./3.); return d[(i-1)*stride]+((n+1./3.)*p+1./3.-i)*(d[i*stride]-d[(i-1)*stride]); } case nsl_stats_quantile_type9: if (p < 5./8./(n+1./4.)) return d[0]; else if (p >= (n-3./8.)/(n+1./4.)) return d[(n-1)*stride]; else { - int i = floor((n+1./4.)*p+3./8.); + int i = (int)floor((n+1./4.)*p+3./8.); return d[(i-1)*stride]+((n+1./4.)*p+3./8.-i)*(d[i*stride]-d[(i-1)*stride]); } } return 0; } double nsl_stats_quantile_from_sorted_data(const double sorted_data[], size_t stride, size_t n, double p) { return nsl_stats_quantile_sorted(sorted_data, stride, n, p, nsl_stats_quantile_type7); } /* R^2 and adj. R^2 */ double nsl_stats_rsquare(double sse, double sst) { return 1. - sse/sst; } double nsl_stats_rsquareAdj(double rsquare, size_t np, size_t dof) { size_t n = np + dof; return 1. - (1. - rsquare) * (n - 1.)/(dof - 1.); } /* t distribution */ double nsl_stats_tdist_t(double parameter, double error) { if (error > 0) return parameter/error; else return DBL_MAX; } double nsl_stats_tdist_p(double t, double dof) { double p = 2. * gsl_cdf_tdist_Q(fabs(t), dof); if (p < 1.e-9) p = 0; return p; } double nsl_stats_tdist_margin(double alpha, double dof, double error) { return gsl_cdf_tdist_Pinv(1. - alpha/2., dof) * error; } /* chi^2 distribution */ double nsl_stats_chisq_p(double t, double dof) { double p = gsl_cdf_chisq_Q(t, dof); if (p < 1.e-9) p = 0; return p; } /* F distribution */ double nsl_stats_fdist_F(double sst, double rms) { return sst/rms; } double nsl_stats_fdist_p(double F, size_t np, double dof) { - double p = gsl_cdf_fdist_Q(F, np, dof); + double p = gsl_cdf_fdist_Q(F, (double)np, dof); if (p < 1.e-9) p = 0; return p; } /* Akaike information criterion */ double nsl_stats_aic(double sse, size_t n, size_t np) { double aic = n * log(sse/n) + 2. * np; // standard formula if (n < 40 * np) // bias correction aic += 2. * np * (np + 1.)/(n - np - 1.); return aic; } /* Bayasian information criterion */ double nsl_stats_bic(double sse, size_t n, size_t np) { - return n * log(sse/n) + np * log(n); + return n * log(sse/(double)n) + np * log((double)n); }