1 /*****************************************************************************
2 * *
3 * Small (Matlab/Octave) Toolbox for Kriging *
4 * *
5 * Copyright Notice *
6 * *
7 * Copyright (C) 2015 CentraleSupelec *
8 * Copyright (C) 2013 SUPELEC *
9 * *
10 * Author: Julien Bect <julien.bect@centralesupelec.fr> *
11 * *
12 * Copying Permission Statement *
13 * *
14 * This file is part of *
15 * *
16 * STK: a Small (Matlab/Octave) Toolbox for Kriging *
17 * (http://sourceforge.net/projects/kriging) *
18 * *
19 * STK is free software: you can redistribute it and/or modify it under *
20 * the terms of the GNU General Public License as published by the Free *
21 * Software Foundation, either version 3 of the License, or (at your *
22 * option) any later version. *
23 * *
24 * STK is distributed in the hope that it will be useful, but WITHOUT *
25 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
26 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public *
27 * License for more details. *
28 * *
29 * You should have received a copy of the GNU General Public License *
30 * along with STK. If not, see <http://www.gnu.org/licenses/>. *
31 * *
32 ****************************************************************************/
33
34 #include "stk_mex.h"
35
36
gpquadform_pairwise(double * x,double * y,double * rx,double * ry,double * h,size_t m,size_t dim)37 static void gpquadform_pairwise
38 (
39 double* x, double* y, double* rx, double* ry,
40 double* h, size_t m, size_t dim
41 )
42 {
43 size_t i, kx, ky;
44 double diff, lambda, rxi, ryi;
45
46 for (i = 0; i < m; i++)
47 {
48 /* compute distance between x[i,:] and y[i,:] */
49 lambda = 0.0;
50 for (kx = i, ky = i; kx < dim * m; kx += m, ky += m)
51 {
52 diff = x[kx] - y[ky];
53 rxi = rx[kx];
54 ryi = ry[ky];
55 lambda += (diff * diff) / (rxi * rxi + ryi * ryi);
56 }
57
58 /* store the result in h */
59 h[i] = lambda;
60 }
61 }
62
63
compute_qxy_matrixy(const mxArray * x,const mxArray * y,const mxArray * rx,const mxArray * ry)64 mxArray* compute_qxy_matrixy
65 (
66 const mxArray* x,
67 const mxArray* y,
68 const mxArray* rx,
69 const mxArray* ry
70 )
71 {
72 size_t d, m;
73 mxArray *h;
74
75 if((!stk_is_realmatrix(x)) || (!stk_is_realmatrix(y)) ||
76 (!stk_is_realmatrix(rx)) || (!stk_is_realmatrix(ry)))
77 mexErrMsgTxt("Input arguments should be real-valued double-precision array.");
78
79 /* Check that the all input arguments have the same number of columns */
80 d = mxGetN(x);
81 if ((mxGetN(y) != d) || (mxGetN(rx) != d) || (mxGetN(ry) != d))
82 mexErrMsgTxt("All input arguments should have the same number of columns.");
83
84 /* Read the number of rows of x and y */
85 if (mxGetM(y) != (m = mxGetM(x)))
86 mexErrMsgTxt("x and y should have the same number of rows.");
87
88 /* Check that rx and ry have the appropriate number of rows */
89 if (mxGetM(rx) != m)
90 mexErrMsgTxt("x and rx should have the same number of rows.");
91 if (mxGetM(ry) != m)
92 mexErrMsgTxt("y and ry should have the same number of rows.");
93
94 /* Create a matrix for the return argument */
95 h = mxCreateDoubleMatrix(m, 1, mxREAL);
96
97 /* Do the actual computations in a subroutine */
98 gpquadform_pairwise(mxGetPr(x), mxGetPr(y), mxGetPr(rx), mxGetPr(ry), mxGetPr(h), m, d);
99
100 return h;
101 }
102
103
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * prhs[])104 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
105 {
106 if (nlhs > 1) /* Check number of output arguments */
107 mexErrMsgTxt("Too many output arguments.");
108
109 if (nrhs != 4) /* Check number of input arguments */
110 mexErrMsgTxt("Incorrect number of input arguments.");
111
112 plhs[0] = compute_qxy_matrixy(prhs[0], prhs[1], prhs[2], prhs[3]);
113 }
114