\begin{figure}[h]
\centering
\includegraphics[angle=0,width=8cm]{Chapters/chapter11/gregory1_plot1.pdf}
-\caption{Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data from \cite{Gregory1982}. Cubic spline fails to preserve monotonicity of the data.}
+\caption[Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data]{Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data from \cite{Gregory1982}. Cubic spline fails to preserve monotonicity of the data.}
\label{ch11:fig1}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[angle=00,width=8cm]{Chapters/chapter11/gregory1_plot2_b.pdf}
-\caption{Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data from \cite{Gregory1982} with non-negative prescribed slopes. Despite non-negative slopes, Hermite cubic spline is not monotone.}
+\caption[Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data]{Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data from \cite{Gregory1982} with non-negative prescribed slopes. Despite non-negative slopes, Hermite cubic spline is not monotone.}
\label{ch11:fig2}
\end{figure}
At the spline evaluation stage we need to compute $s(z_k)$ for a sequence of query values ${z_k}, k=1,\ldots,K$. For each $z_k$ we locate the interval $[t_i,t_{i+1}]$ containing $z_k$, using bisection algorithm presented in Figure \ref{ch11:algeval}, and then apply the appropriate coefficients of the quadratic function. This is also done in parallel.
The bisection algorithm could be implemented using texture memory (to cache the array \texttt{z}), but this is not shown in Figure \ref{ch11:algeval}.
-\begin{figure}[!hp]
-\renewcommand{\baselinestretch}{1}
- \begin{alltt}
-\begin{center}
-\begin{minipage}{13cm}\small
-
-template<typename Tx, typename Ty>
-\_\_global\_\_ void CalculateCoefficientsKnots( Tx *u, Ty *v, double *b, double *c,
- double *t, double *alpha, double *beta, double *gamma, int N )
-\{
- int tid = threadIdx.x + blockIdx.x * blockDim.x;
- int s = tid*2;
- while(tid<=(N-2))
- \{
- // decide whether an additional knot is necessary
- if(fabs(c[tid]+c[tid+1]- 2*b[tid])<=0.1e-5) // tolerance
- \{ //no additional knot
- h[s]=h[s+1]=u[tid];
- alpha[s]=alpha[s+1]=v[tid];
- beta[s]=beta[s+1]=c[tid];
- gamma[s]=gamma[s+1]=(c[tid+1]-c[tid])/(2*(fmax(1e-10,u[tid+1]-u[tid])));
- \} else \{ //adding a knot
- h[s]=u[tid];
- //determine the position of the knot
- if((c[tid+1] - b[tid])*(c[tid] - b[tid])<0)
- h[s+1]=u[tid+1] + (c[tid] - b[tid])*(fmax(1e-10,u[tid+1]-u[tid]))/
- fmax(1e-10,(c[tid+1] - c[tid]));
- else
- h[s+1]=0.5*(u[tid+1] + u[tid]);
- //calculate coefficients
- double dtemp = (2*b[tid] - c[tid+1])+((c[tid+1] - c[tid])*(h[s+1] - u[tid]))/
- fmax(1e-10,(u[tid+1] - u[tid]));
- alpha[s]=v[tid]; beta[s]=c[tid];
- gamma[s]=(dtemp - c[tid])/(2*fmax(1e-10,(h[s+1] - u[tid])));
- alpha[s+1]=v[tid] + c[tid]*(h[s+1] - u[tid]) +
- (dtemp - c[tid])*(h[s+1] - u[tid])/2;
- gamma[s+1]=(c[tid+1] - dtemp)/(2*fmax(1e-10,(u[tid+1] - h[s+1])));
- beta[s+1]=dtemp;
- \}
- tid += blockDim.x * gridDim.x; s = tid*2;
- \}
- \_\_syncthreads();
- // Select a single thread to perform the last operation
- if((threadIdx.x ) == 0) \{
- s = (N-1) * 2; h[s]=u[N-1];
- \}
- \_\_syncthreads();
-\}
-\end{minipage}
-\end{center}
-\end{alltt}
-\caption{Implementation of the kernel for calcuating spline knots and coefficients. Function fmax is used to avoid division by zero for data with coinciding abscissae.}
-\label{ch11:algcoef}
-\renewcommand{\baselinestretch}{2}
-\end{figure}
-
-\begin{figure}[!hp]
-\renewcommand{\baselinestretch}{1}
- \begin{alltt}
-\begin{center}
-\begin{minipage}{13cm}\small
-
-template<typename Tx, typename Ty>
-\_\_global\_\_ void CalculateBeta(Tx *u, Ty *v, double *b, int N)
-\{
- int tid = threadIdx.x + blockIdx.x * blockDim.x;
- while(tid<=(N-2)) \{
- b[tid]=(v[tid+1]-v[tid])/fmax(1e-20,double(u[tid+1]-u[tid]));
- tid += blockDim.x * gridDim.x;
- \}
- \_\_syncthreads();
-\}
-\_\_global\_\_ void CalculateDGeneral( double *b, double *c, int N)
-\{
- int tid = threadIdx.x + blockIdx.x * blockDim.x;
- while(tid<=(N-2)) \{
- if((b[tid-1]*b[tid])<=0) c[tid]=0;
- else c[tid]=(2*b[tid-1]*b[tid])/(b[tid-1]+b[tid]);
- \}
- tid += blockDim.x * gridDim.x;
- \}
- \_\_syncthreads();
-\}
-\_\_global\_\_ void CalculateD( double *b, double *c, int N )
-\{
- if((b[0]*(2*b[0]-c[1]))<=0) c[0]=0;
- else c[0]=2*b[0] - c[1];
- if((b[N-2]*(2*b[N-2]-c[N-2]))<=0) c[N-1]=0;
- else c[N-1]=2*b[N-2] - c[N-2];
- \_\_syncthreads();
-\}
-template<typename Tx, typename Ty>
-int BuildMonotonSpline(Tx *d_X, Ty *d_Y, int N,
- double *t, double *alpha, double *beta, double *gamma)
-\{
- int T = (N-1)*2+1; // length of the output array
- double *b, *c; // temp variables
- cudaMalloc( (void**)&b, 1*N*sizeof(double) );
- cudaMalloc( (void**)&c, 2*N*sizeof(double) );
- int threads=256;
- int blocks = (N-1)/threads + 1;
- CalculateBeta<<<blocks,threads>>>(d_X,d_Y,b,N);
- CalculateDGeneral<<<blocks,threads>>>(b,c,N);
- CalculateD<<<1,1>>>(b,c,NN); // calculate d_1 and d_N
- CalculateCoefficientsKnots<<<blocks,threads>>>(d_X,
- d_Y,b,c,h,alpha,beta,gamma,N);
- cudaFree(b); cudaFree(c);
- return T;
-\}
-\end{minipage}
-\end{center}
-\end{alltt}
-\caption{Calculation of monotone spline knots and coefficients.}
-\label{ch11:algcoef1}
-\renewcommand{\baselinestretch}{2}
-\end{figure}
-
-\begin{figure}[!hp]
-\renewcommand{\baselinestretch}{1}
- \begin{alltt}
-\begin{center}
-\begin{minipage}{13cm}\small
-
-template<typename T>
-\_\_device\_\_ void Bisection\_device(T z, T* t, int mi, int ma, int* l)
-\{
- int i; ma--;
- while(1) \{
- i=(mi+ma)/2;
- if(z >= t[i]) mi=i+1;
- else ma=i;
- if(mi>=ma) break;
- \}
- *l = mi-1;
-\}
-
-/* Kernel to evaluates monotone spline for a sequence of query points
- residing in the array z of size m
-*/
-template<typename Tx, typename Ty>
-\_\_global\_\_ void d\_MonSplineValue(Tx* z, int K, double* t,
- double * alpha, double * beta, double * gamma, int T, Ty *value)
-\{
- int tid = threadIdx.x + blockIdx.x * blockDim.x;
- int mi=0, ma=T, i=0;
- Ty r;
- while(tid<K)
- \{
- Bisection\_device(z[tid], t, mi, ma, &i);
- r= z[tid]-t[i];
- r= alpha[i] + r*(beta[i] + gamma[i]*r);
- value[tid]=r;
- tid += blockDim.x * gridDim.x;
- \}
- \_\_syncthreads();
-\}
-
-template<typename Tx, typename Ty>
-void MonotoneSplineValue(Tx *z, int K, double* t,
- double * alpha, double * beta, double * gamma, int T, Ty* result)
-\{
- int blocks,threads=256;
- blocks=(K-1)/threads+1;
- d\_MonSplineValue<<<blocks,threads>>>(z,K,t,alpha,beta,gamma,T,result);
-\}
-\end{minipage}
-\end{center}
-\end{alltt}
-\caption{Implementation of the spline evaluation algorithm for GPU.}
-\label{ch11:algeval}
-\renewcommand{\baselinestretch}{2}
-\end{figure}
-
-
+\lstinputlisting[label=ch11:algcoef,caption=Implementation of the kernel for calcuating spline knots and coefficients. Function fmax is used to avoid division by zero for data with coinciding abscissae.]{Chapters/chapter11/code1.cu}
+
+
+%% \begin{figure}[!hp]
+%% \renewcommand{\baselinestretch}{1}
+%% \begin{alltt}
+%% \begin{center}
+%% \begin{minipage}{13cm}\small
+
+%% template<typename Tx, typename Ty>
+%% \_\_global\_\_ void CalculateCoefficientsKnots( Tx *u, Ty *v, double *b, double *c,
+%% double *t, double *alpha, double *beta, double *gamma, int N )
+%% \{
+%% int tid = threadIdx.x + blockIdx.x * blockDim.x;
+%% int s = tid*2;
+%% while(tid<=(N-2))
+%% \{
+%% // decide whether an additional knot is necessary
+%% if(fabs(c[tid]+c[tid+1]- 2*b[tid])<=0.1e-5) // tolerance
+%% \{ //no additional knot
+%% h[s]=h[s+1]=u[tid];
+%% alpha[s]=alpha[s+1]=v[tid];
+%% beta[s]=beta[s+1]=c[tid];
+%% gamma[s]=gamma[s+1]=(c[tid+1]-c[tid])/(2*(fmax(1e-10,u[tid+1]-u[tid])));
+%% \} else \{ //adding a knot
+%% h[s]=u[tid];
+%% //determine the position of the knot
+%% if((c[tid+1] - b[tid])*(c[tid] - b[tid])<0)
+%% h[s+1]=u[tid+1] + (c[tid] - b[tid])*(fmax(1e-10,u[tid+1]-u[tid]))/
+%% fmax(1e-10,(c[tid+1] - c[tid]));
+%% else
+%% h[s+1]=0.5*(u[tid+1] + u[tid]);
+%% //calculate coefficients
+%% double dtemp = (2*b[tid] - c[tid+1])+((c[tid+1] - c[tid])*(h[s+1] - u[tid]))/
+%% fmax(1e-10,(u[tid+1] - u[tid]));
+%% alpha[s]=v[tid]; beta[s]=c[tid];
+%% gamma[s]=(dtemp - c[tid])/(2*fmax(1e-10,(h[s+1] - u[tid])));
+%% alpha[s+1]=v[tid] + c[tid]*(h[s+1] - u[tid]) +
+%% (dtemp - c[tid])*(h[s+1] - u[tid])/2;
+%% gamma[s+1]=(c[tid+1] - dtemp)/(2*fmax(1e-10,(u[tid+1] - h[s+1])));
+%% beta[s+1]=dtemp;
+%% \}
+%% tid += blockDim.x * gridDim.x; s = tid*2;
+%% \}
+%% \_\_syncthreads();
+%% // Select a single thread to perform the last operation
+%% if((threadIdx.x ) == 0) \{
+%% s = (N-1) * 2; h[s]=u[N-1];
+%% \}
+%% \_\_syncthreads();
+%% \}
+%% \end{minipage}
+%% \end{center}
+%% \end{alltt}
+%% \caption{Implementation of the kernel for calcuating spline knots and coefficients. Function fmax is used to avoid division by zero for data with coinciding abscissae.}
+%% \label{ch11:algcoef}
+%% \renewcommand{\baselinestretch}{2}
+%% \end{figure}
+
+
+\lstinputlisting[label=ch11:algcoef1,caption=Calculation of monotone spline knots and coefficients.]{Chapters/chapter11/code2.cu}
+
+%% \begin{figure}[!hp]
+%% \renewcommand{\baselinestretch}{1}
+%% \begin{alltt}
+%% \begin{center}
+%% \begin{minipage}{13cm}\small
+
+%% template<typename Tx, typename Ty>
+%% \_\_global\_\_ void CalculateBeta(Tx *u, Ty *v, double *b, int N)
+%% \{
+%% int tid = threadIdx.x + blockIdx.x * blockDim.x;
+%% while(tid<=(N-2)) \{
+%% b[tid]=(v[tid+1]-v[tid])/fmax(1e-20,double(u[tid+1]-u[tid]));
+%% tid += blockDim.x * gridDim.x;
+%% \}
+%% \_\_syncthreads();
+%% \}
+%% \_\_global\_\_ void CalculateDGeneral( double *b, double *c, int N)
+%% \{
+%% int tid = threadIdx.x + blockIdx.x * blockDim.x;
+%% while(tid<=(N-2)) \{
+%% if((b[tid-1]*b[tid])<=0) c[tid]=0;
+%% else c[tid]=(2*b[tid-1]*b[tid])/(b[tid-1]+b[tid]);
+%% \}
+%% tid += blockDim.x * gridDim.x;
+%% \}
+%% \_\_syncthreads();
+%% \}
+%% \_\_global\_\_ void CalculateD( double *b, double *c, int N )
+%% \{
+%% if((b[0]*(2*b[0]-c[1]))<=0) c[0]=0;
+%% else c[0]=2*b[0] - c[1];
+%% if((b[N-2]*(2*b[N-2]-c[N-2]))<=0) c[N-1]=0;
+%% else c[N-1]=2*b[N-2] - c[N-2];
+%% \_\_syncthreads();
+%% \}
+%% template<typename Tx, typename Ty>
+%% int BuildMonotonSpline(Tx *d_X, Ty *d_Y, int N,
+%% double *t, double *alpha, double *beta, double *gamma)
+%% \{
+%% int T = (N-1)*2+1; // length of the output array
+%% double *b, *c; // temp variables
+%% cudaMalloc( (void**)&b, 1*N*sizeof(double) );
+%% cudaMalloc( (void**)&c, 2*N*sizeof(double) );
+%% int threads=256;
+%% int blocks = (N-1)/threads + 1;
+%% CalculateBeta<<<blocks,threads>>>(d_X,d_Y,b,N);
+%% CalculateDGeneral<<<blocks,threads>>>(b,c,N);
+%% CalculateD<<<1,1>>>(b,c,NN); // calculate d_1 and d_N
+%% CalculateCoefficientsKnots<<<blocks,threads>>>(d_X,
+%% d_Y,b,c,h,alpha,beta,gamma,N);
+%% cudaFree(b); cudaFree(c);
+%% return T;
+%% \}
+%% \end{minipage}
+%% \end{center}
+%% \end{alltt}
+%% \caption{Calculation of monotone spline knots and coefficients.}
+%% \label{ch11:algcoef1}
+%% \renewcommand{\baselinestretch}{2}
+%% \end{figure}
+
+%% \begin{figure}[!hp]
+%% \renewcommand{\baselinestretch}{1}
+%% \begin{alltt}
+%% \begin{center}
+%% \begin{minipage}{13cm}\small
+
+%% template<typename T>
+%% \_\_device\_\_ void Bisection\_device(T z, T* t, int mi, int ma, int* l)
+%% \{
+%% int i; ma--;
+%% while(1) \{
+%% i=(mi+ma)/2;
+%% if(z >= t[i]) mi=i+1;
+%% else ma=i;
+%% if(mi>=ma) break;
+%% \}
+%% *l = mi-1;
+%% \}
+
+%% /* Kernel to evaluates monotone spline for a sequence of query points
+%% residing in the array z of size m
+%% */
+%% template<typename Tx, typename Ty>
+%% \_\_global\_\_ void d\_MonSplineValue(Tx* z, int K, double* t,
+%% double * alpha, double * beta, double * gamma, int T, Ty *value)
+%% \{
+%% int tid = threadIdx.x + blockIdx.x * blockDim.x;
+%% int mi=0, ma=T, i=0;
+%% Ty r;
+%% while(tid<K)
+%% \{
+%% Bisection\_device(z[tid], t, mi, ma, &i);
+%% r= z[tid]-t[i];
+%% r= alpha[i] + r*(beta[i] + gamma[i]*r);
+%% value[tid]=r;
+%% tid += blockDim.x * gridDim.x;
+%% \}
+%% \_\_syncthreads();
+%% \}
+
+%% template<typename Tx, typename Ty>
+%% void MonotoneSplineValue(Tx *z, int K, double* t,
+%% double * alpha, double * beta, double * gamma, int T, Ty* result)
+%% \{
+%% int blocks,threads=256;
+%% blocks=(K-1)/threads+1;
+%% d\_MonSplineValue<<<blocks,threads>>>(z,K,t,alpha,beta,gamma,T,result);
+%% \}
+%% \end{minipage}
+%% \end{center}
+%% \end{alltt}
+%% \caption{Implementation of the spline evaluation algorithm for GPU.}
+%% \label{ch11:algeval}
+%% \renewcommand{\baselinestretch}{2}
+%% \end{figure}
+
+\lstinputlisting[label=ch11:algeval,caption=Implementation of the spline evaluation algorithm for GPU.]{Chapters/chapter11/code3.cu}
\subsection{Monotone Hermite splines}
$$
with $\hat y(k,l)$ being the unrestricted maximum likelihood estimator of $y_k\ldots,y_l$. For quadratic cost function $\hat y(k,l)$ is the mean, as in PAV and MLS algorithms, for the absolute deviations it becomes the median, and for other cost functions an M-estimator of location. The MLS algorithm can be applied to such isotone regression problems with very little modification, while linear in time algorithm may not be available. Our parallel MLS algorithm will be valuable in such cases.
-%\renewcommand{\baselinestretch}{1}
+
+
+%% %\renewcommand{\baselinestretch}{1}
\begin{table}[!h]
\begin{center}
\caption{The average CPU time (sec) of the serial PAVA, MLS and parallel MLS algorithms. } \label{ch11:table1}
\end{tabular}
\end{center}
\end{table}
-%\renewcommand{\baselinestretch}{2}
-
-
-\begin{figure}[!hp]
- \begin{alltt}
-\begin{center}
-\begin{minipage}{13cm}\small
-template<typename Tx>
-__device__ Tx Aver(Tx z,int i,int j, Tx *z) \{return (z-z[j+1])/(j-i+1);\}
-
-template<typename Tx>
-__global__ void monotonizekernel(Tx *y, Tx *z, Tx *u, int *key, int n)
-\{ int i = threadIdx.x + blockIdx.x * blockDim.x;
- if(i<n) \{
- int smallestJ = i;
- Tx curP, smallestP, curz=z[i];
- smallestP=Aver(curz,i,i,z);
- for(int j = i+1; j < n; j++) \{
- curP=Aver(curz,i,j,z);
- if(smallestP>curP) \{
- smallestJ = j;
- smallestP = curP;
- \}
- \}
- curP=y[i];
- if(curP > smallestP) t=smallestP;
- else smallestJ=i;
- key[i]=smallestJ;
- u[i]=t;
- \}
-\}
-
-template< typename Tx >
-void MonotonizeData(Tx *y, int n, Tx *u) \{
- thrust::less_equal<int> binary_pred;
- thrust::maximum<Tx> binary_op2;
- thrust::device_vector<Tx> z_d(n+1);
- thrust::device_vector<int> keys_d(n);
- thrust::device_ptr<Tx> y_d(y), u_d(u);
- thrust::fill(u_d, u_d+n, -1e100);
- thrust::fill(keys_d.begin(), keys_d.end(), 0);
-
- thrust::reverse_iterator< typename thrust::device_vector<Tx>::iterator >
- y_reverse_b(y_d+n), y_reverse_end(y_d), z_reverse_b(z_d.end());
+%% %\renewcommand{\baselinestretch}{2}
+
+
+%% \begin{figure}[!hp]
+%% \begin{alltt}
+%% \begin{center}
+%% \begin{minipage}{13cm}\small
+%% template<typename Tx>
+%% __device__ Tx Aver(Tx z,int i,int j, Tx *z) \{return (z-z[j+1])/(j-i+1);\}
+
+%% template<typename Tx>
+%% __global__ void monotonizekernel(Tx *y, Tx *z, Tx *u, int *key, int n)
+%% \{ int i = threadIdx.x + blockIdx.x * blockDim.x;
+%% if(i<n) \{
+%% int smallestJ = i;
+%% Tx curP, smallestP, curz=z[i];
+%% smallestP=Aver(curz,i,i,z);
+%% for(int j = i+1; j < n; j++) \{
+%% curP=Aver(curz,i,j,z);
+%% if(smallestP>curP) \{
+%% smallestJ = j;
+%% smallestP = curP;
+%% \}
+%% \}
+%% curP=y[i];
+%% if(curP > smallestP) t=smallestP;
+%% else smallestJ=i;
+%% key[i]=smallestJ;
+%% u[i]=t;
+%% \}
+%% \}
+
+%% template< typename Tx >
+%% void MonotonizeData(Tx *y, int n, Tx *u) \{
+%% thrust::less_equal<int> binary_pred;
+%% thrust::maximum<Tx> binary_op2;
+%% thrust::device_vector<Tx> z_d(n+1);
+%% thrust::device_vector<int> keys_d(n);
+%% thrust::device_ptr<Tx> y_d(y), u_d(u);
+%% thrust::fill(u_d, u_d+n, -1e100);
+%% thrust::fill(keys_d.begin(), keys_d.end(), 0);
+
+%% thrust::reverse_iterator< typename thrust::device_vector<Tx>::iterator >
+%% y_reverse_b(y_d+n), y_reverse_end(y_d), z_reverse_b(z_d.end());
- thrust::inclusive_scan(y_reverse_b, y_reverse_end, z_reverse_b+1);
+%% thrust::inclusive_scan(y_reverse_b, y_reverse_end, z_reverse_b+1);
- monotonizekernel<<<grid, block>>>(y, thrust::raw_pointer_cast(&z_d[0]),
- u, thrust::raw_pointer_cast(&keys_d[0]), n );
+%% monotonizekernel<<<grid, block>>>(y, thrust::raw_pointer_cast(&z_d[0]),
+%% u, thrust::raw_pointer_cast(&keys_d[0]), n );
- thrust::sort(keys_d.begin(), keys_d.end());
- thrust::inclusive_scan_by_key(keys_d.begin(), keys_d.end(),
- u_d, u_d, binary_pred, binary_op2);
-\}
-\end{minipage}
-\end{center}
-\end{alltt}
-\caption{Fragments of implementation of a parallel version of the MLS algorithm using Thrust library.}
-\label{ch11:algMLS}
-\end{figure}
+%% thrust::sort(keys_d.begin(), keys_d.end());
+%% thrust::inclusive_scan_by_key(keys_d.begin(), keys_d.end(),
+%% u_d, u_d, binary_pred, binary_op2);
+%% \}
+%% \end{minipage}
+%% \end{center}
+%% \end{alltt}
+%% \caption{Fragments of implementation of a parallel version of the MLS algorithm using Thrust library.}
+%% \label{ch11:algMLS}
+%% \end{figure}
+
+\lstinputlisting[label=ch11:algMLS,caption=Fragments of implementation of a parallel version of the MLS algorithm using Thrust library.]{Chapters/chapter11/code4.cu}
\section{Conclusion} \label{ch11:conc}
We presented three GPU-based parallel algorithms for approximating monotone data: monotone quadratic spline, monotone Hermite rational spline and minimum lower sets algorithm for monotonizing noisy data. These tools are valuable in a number of applications that involve large data sets modeled by monotone nonlinear functions.
-The source code of the package monospline is available from \texttt{www.deakin.edu.au/$\sim$ gleb/monospline.html }
+The source code of the package monospline is available from \texttt{www.deakin.edu.au/$\sim$gleb/monospline.html }
\putbib[Chapters/chapter11/biblio11]
