From 6ca2e6ea3c549c8673e2e78368ed653584730add Mon Sep 17 00:00:00 2001 From: jkriege2 Date: Mon, 18 Mar 2024 12:10:38 +0100 Subject: [PATCH] Speed-Improvemend for 2D-KDE calculation by reordering loops and allowing for inlined function calls --- lib/jkqtmath/jkqtpstatkde.cpp | 54 -------------- lib/jkqtmath/jkqtpstatkde.h | 136 ++++++++++++++++++++++++++++------ 2 files changed, 113 insertions(+), 77 deletions(-) diff --git a/lib/jkqtmath/jkqtpstatkde.cpp b/lib/jkqtmath/jkqtpstatkde.cpp index 8fd36e85fe..cbaa2ebf83 100644 --- a/lib/jkqtmath/jkqtpstatkde.cpp +++ b/lib/jkqtmath/jkqtpstatkde.cpp @@ -22,64 +22,10 @@ #include "jkqtmath/jkqtpstatkde.h" -double jkqtpstatKernel1DGaussian(double t) { - return exp(-0.5*t*t)/JKQTPSTATISTICS_SQRT_2PI; -} - - -double jkqtpstatKernel1DCauchy(double t) { - return 1.0/(JKQTPSTATISTICS_PI*(1.0+t*t)); -} -double jkqtpstatKernel1DPicard(double t) { - return exp(-0.5*fabs(t))/2.0; -} - - -double jkqtpstatKernel1DEpanechnikov(double t) { - return (fabs(t)<1.0)?(0.75*(1.0-t*t)):0.0; -} - - -double jkqtpstatKernel1DUniform(double t) { - return (fabs(t)<=1.0)?0.5:0.0; -} - - -double jkqtpstatKernel1DTriangle(double t) { - return (fabs(t)<=1.0)?(1.0-fabs(t)):0.0; -} -double jkqtpstatKernel1DQuartic(double t) { - return (fabs(t)<=1.0)?(15.0/16.0*jkqtp_sqr(1.0-t*t)):0.0; -} - -double jkqtpstatKernel1DTriweight(double t) { - return (fabs(t)<1.0)?(35.0/32.0*jkqtp_cube(1.0-t*t)):0.0; -} - - - -double jkqtpstatKernel1DTricube(double t) { - return (fabs(t)<1.0)?(70.0/81.0*jkqtp_cube(1.0-jkqtp_cube(fabs(t)))):0.0; -} - - -double jkqtpstatKernel1DCosine(double t) { - return (fabs(t)<1.0)?(JKQTPSTATISTICS_PI/4.0*cos(t*JKQTPSTATISTICS_PI/2.0)):0.0; -} - - -double jkqtpstatKernel2DGaussian(double tx, double ty) -{ - return exp(-0.5*(tx*tx+ty*ty))/(2.0*JKQTPSTATISTICS_PI); -} - -double jkqtpstatKernel2DUniform(double tx, double ty) { - return (fabs(tx)<1.0 && fabs(ty)<=1.0)?0.25:0.0; -} diff --git a/lib/jkqtmath/jkqtpstatkde.h b/lib/jkqtmath/jkqtpstatkde.h index 48644758b6..b011000565 100644 --- a/lib/jkqtmath/jkqtpstatkde.h +++ b/lib/jkqtmath/jkqtpstatkde.h @@ -49,64 +49,84 @@ \f[ k(t):=\frac{1}{\sqrt{2\pi}}\exp \left(-\frac{1}{2}t^2\right) \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DGaussian(double t); +inline double jkqtpstatKernel1DGaussian(double t) { + return exp(-0.5*t*t)/JKQTPSTATISTICS_SQRT_2PI; +} /*! \brief a 1D Cauchy kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t):=\frac{1}{\pi(1+t^2)} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DCauchy(double t); +inline double jkqtpstatKernel1DCauchy(double t) { + return 1.0/(JKQTPSTATISTICS_PI*(1.0+t*t)); +} /*! \brief a 1D Picard kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t):=\frac{1}{2}\exp(-|t|) \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DPicard(double t); +inline double jkqtpstatKernel1DPicard(double t) { + return exp(-0.5*fabs(t))/2.0; +} /*! \brief a 1D Epanechnikov kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}\frac{3}{4} ( 1- t^2 ), & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DEpanechnikov(double t); +inline double jkqtpstatKernel1DEpanechnikov(double t) { + return (fabs(t)<1.0)?(0.75*(1.0-t*t)):0.0; +} /*! \brief a 1D uniform kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}0.5, & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DUniform(double t); +inline double jkqtpstatKernel1DUniform(double t) { + return (fabs(t)<=1.0)?0.5:0.0; +} /*! \brief a 1D Epanechnikov kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}1-|t|, & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DTriangle(double t); +inline double jkqtpstatKernel1DTriangle(double t) { + return (fabs(t)<=1.0)?(1.0-fabs(t)):0.0; +} /*! \brief a 1D quartic kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}\frac{15}{16}(1-t^2)^2, & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DQuartic(double t); +inline double jkqtpstatKernel1DQuartic(double t) { + return (fabs(t)<=1.0)?(15.0/16.0*jkqtp_sqr(1.0-t*t)):0.0; +} /*! \brief a 1D triweight kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}\frac{35}{32}(1-t^2)^3, & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DTriweight(double t); +inline double jkqtpstatKernel1DTriweight(double t) { + return (fabs(t)<1.0)?(35.0/32.0*jkqtp_cube(1.0-t*t)):0.0; +} /*! \brief a 1D tricube kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}\frac{70}{81}(1-|t|^3)^3, & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DTricube(double t); +inline double jkqtpstatKernel1DTricube(double t) { + return (fabs(t)<1.0)?(70.0/81.0*jkqtp_cube(1.0-jkqtp_cube(fabs(t)))):0.0; +} /*! \brief a 1D cosine kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_1dkde_kernels \f[ k(t) :=\begin{cases}\frac{\pi}{4}\cos\left(\frac{\pi}{2}t\right), & \text{if }t\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel1DCosine(double t); +inline double jkqtpstatKernel1DCosine(double t) { + return (fabs(t)<1.0)?(JKQTPSTATISTICS_PI/4.0*cos(t*JKQTPSTATISTICS_PI/2.0)):0.0; +} @@ -125,14 +145,19 @@ jkqtmath_LIB_EXPORT double jkqtpstatKernel1DCosine(double t); \f[ k(t_x, t_y):=\frac{1}{2\pi}\exp \left(-\frac{t_x^2+t_y^2}{2}\right) \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel2DGaussian(double tx, double ty); +inline double jkqtpstatKernel2DGaussian(double tx, double ty) +{ + return exp(-0.5*(tx*tx+ty*ty))/(2.0*JKQTPSTATISTICS_PI); +} /*! \brief a 2D Uniform kernel function, e.g. for Kernel Density Estimation \ingroup jkqtptools_math_statistics_2dkde_kernels \f[ k(t_x, t_y):=\begin{cases}\frac{1}{4}, & \text{if }t_x,t_y\in [-1;1]\\0, & \text{else}\end{cases} \f] */ -jkqtmath_LIB_EXPORT double jkqtpstatKernel2DUniform(double tx, double ty); +inline double jkqtpstatKernel2DUniform(double tx, double ty) { + return (fabs(tx)<1.0 && fabs(ty)<=1.0)?0.25:0.0; +} @@ -473,6 +498,45 @@ inline double jkqtpstatEvaluateKernelSum2D(double x, double y, InputItX firstX, +/*! \brief stores a part-result of the Kernel Density Estimator (KDE) at a given position + \ingroup jkqtptools_math_statistics_2dkde + \internal + + Iterates over all datapoints in firstX..lastX and firstY..lastY, for each one determines whether the floating-point values are valid + and the calls fStore, which is a functor that is supposed to iterate over the KDE 2D array and over calls sum up the KDE. + + This is internally used by jkqtpstatKDE2D(). The construction is not as straight-forward, as iterating over all positions in the 2D-KJDE + and the calling jkqtpstatEvaluateKernelSum2D(), bt by reordering the loops it is significatyl faster (~20%) + + \tparam InputItX standard iterator type of \a firstX and \a lastX. + \tparam InputItY standard iterator type of \a firstY and \a lastY. + \tparam FF functor + \param x where to evaluate the kernel sum, x-coordinate + \param y where to evaluate the kernel sum, y-coordinate + \param firstX iterator pointing to the first x-position item in the dataset to use \f$ X_1 \f$ + \param lastX iterator pointing behind the last x-position item in the dataset to use \f$ X_N \f$ + \param firstY iterator pointing to the first y-position item in the dataset to use \f$ Y_1 \f$ + \param lastY iterator pointing behind the last y-position item in the dataset to use \f$ Y_N \f$ + \param kernel the kernel function to use (e.g. jkqtpstatKernel1DGaussian() ) + \param bandwidthX x-bandwidth used for the KDE + \param bandwidthY y-bandwidth used for the KDE + +*/ +template +inline int jkqtpstatStoreKernelSum2D(FF&& fStore, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, const std::function& kernel, double bandwidthX, double bandwidthY) { + auto itX=firstX; + auto itY=firstY; + int cnt=0; + for (; (itX!=lastX)&&(itY!=lastY); ++itX, ++itY) { + const double vx=jkqtp_todouble(*itX); + const double vy=jkqtp_todouble(*itY); + if (JKQTPIsOKFloat(vx) && JKQTPIsOKFloat(vy)) { + fStore(vx,vy, kernel,bandwidthX, bandwidthY); + cnt++; + } + } + return cnt; +} /*! \brief calculate an autoranged 2-dimensional Kernel Density Estimation (KDE) from the given data range \a firstX / \a firstY ... \a lastY / \a lastY \ingroup jkqtptools_math_statistics_2dkde @@ -504,18 +568,44 @@ inline void jkqtpstatKDE2D(InputItX firstX, InputItX lastX, InputItY firstY, Inp const double binwx=fabs(xmax-xmin)/static_cast(xbins); const double binwy=fabs(ymax-ymin)/static_cast(ybins); - double y=ymin; - auto itOut=histogramImgOut; - for (size_t by=0; by& kernel, double bandwidthX, double bandwidthY){ + double y=ymin; + auto itOut=histogramImgOut; + + + for (size_t by=0; by