1    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2%%
3%% REGRESSION MODULE FOR RESPONSE SURFACE METHODS
4%%
5%% Implemented from the article:
6%% "A CONSERVATIVE MESH-FREE APPROACH
7%% FOR FLUID-STRUCTURE INTERFACE PROBLEMS"
8%% by Quaranta, Masarati and Mantegazza.
9%%
10%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11
12%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13%%
14%%  Building the regression coefficients I
15%%  and their derivatives
16%%  given the number of interpolant points
17%%
18%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20
21
22function [I,varargout] = f_regression_derivatives_kdtree(opoints, inpoints, tree, basis, n_points, phi_ord, varargin)
23
24%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25%%
26%% Inputs:
27%%  - opoints:  the desired output point of the regression.
28%%  - inpoints: the vectors of known point positions arranged as rows.
29%%  - tree: pointer to kdtree structure built using kdtree function if
30%%                  empty a new tree is built
31%%  - basis: the type of basis desired:
32%%          1 --> first-order polynomial
33%%          2 --> full second-order polynomial
34%%          3 --> full third-order polynomial
35%%          4 --> first-order polynomial + quadratic terms
36%%          5 --> sine and cosine
37%%  - n_points: number of points to be used in the local support for the interpolation.
38%%  - phi_ord: the order of weight desired. For example, in 3 dimensions:
39%%          0 --> (1 - r)^2+
40%%          2 --> (1 - r)^4+ (4r + 1)
41%%          4 --> (1 - r)^6+ (35/5 r^2 + 18/3 r + 1)
42%%          6 --> (1 - r)^8+ (32 r^3 + 25 r^2 + 8 r + 1)
43%%  - Weight: the weight of the trim points (Z)  optional
44%%
45%% Outputs:
46%%  - I: the regression matrix
47%%  - I_i: cell of sparse matrices for the derivatives of the regression coefficients on the i-th direction.
48%%  - Index: the known points that are taken for the regression.
49%%
50%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52
53% Total number of known points:
54n_pt = size(inpoints,1) ;
55%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56% total number of query points
57n_qpt = size(opoints,1) ;
58% space dimension
59dim = size(inpoints,2);
60if (dim ~= size(opoints,2))
61    error('Incopatible dimension between input points and output points');
62end
63
64if nargout > 1
65    comp_derivatives = 1;
66else
67    comp_derivatives = 0;
68end
69if (nargin > 6)
70    useW = 1;
71    Weight = varargin{1};
72else
73    useW = 0;
74end
75
76%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77
78%% FINDING THE n_points NEAREST POINTS:
79%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80Index = zeros(n_qpt, n_points+1);
81dist  = zeros(n_qpt, n_points+1);
82if isempty(tree)
83    tree = kd_buildtree(inpoints,0);
84end
85if (length(tree) > n_points)
86    for i = 1 : n_qpt
87        idx  = (kd_knn(tree,opoints(i,:),n_points+1,0))';
88        Index(i,:) = idx;
89        for j = 1 : n_points + 1
90            dist(i,j) = sqrt(sum((opoints(i,:) - inpoints(Index(i,j),:)).^2));
91        end
92    end
93elseif (length(tree) > n_points - 1)
94    for i = 1 : n_qpt
95        idx  = (kd_knn(tree,opoints(i,:),n_points,0))';
96        Index(i,1:length(idx)) = idx;
97        Index(i,length(idx)+1:end) = idx(end);
98        for j = 1 : n_points + 1
99            dist(i,j) = sqrt(sum((opoints(i,:) - inpoints(Index(i,j),:)).^2));
100        end
101    end
102else
103    error('MLS Interpolation: There are not enought points reduce the number of points required');
104end
105
106
107%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108% structure for the sparse I matrix
109r = zeros(n_qpt*n_points,1);
110c = zeros(n_qpt*n_points,1);
111h = zeros(n_qpt*n_points,1);
112if comp_derivatives
113    rd = cell(dim,1);
114    cd = cell(dim,1);
115    hd = cell(dim,1);
116    for i = 1 : dim
117        rd{i} = zeros(n_qpt*n_points,1);
118        cd{i} = zeros(n_qpt*n_points,1);
119        hd{i} = zeros(n_qpt*n_points,1);
120    end
121end
122coef = 0;
123
124p_dim = 0;
125for i = 0 : basis
126    p_dim = p_dim + nchoosek(dim + i - 1, i);
127end
128P = zeros(size(Index,2) - 1, p_dim);
129
130for stp = 1 : n_qpt
131
132%% BUILDING MATRIX P:
133%%%%%%%%%%%%%%%%%%%%%
134
135    for i = 1 : size(Index,2) - 1
136
137        P(i,:) = f_polyn_basis(inpoints(Index(stp,i),:),basis, opoints(stp,:)) ;
138
139    end
140
141%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142%P = P';
143%% BUILDING THE MATRIX FOR WEIGHTS:
144%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145    if  (abs(dist(stp,end) - dist(stp,end-1)) > 100*eps)
146        r_max = 0.5 * (dist(stp,end) + dist(stp,end-1)) ;
147    else
148        r_max = 1.1 * dist(stp,end);
149    end
150
151    if useW
152        phi = Weight(i)*f_weight_function_givendist(dist(stp,1:end-1),dim,r_max,phi_ord) ;
153    else
154        phi = f_weight_function_givendist(dist(stp,1:end-1),dim,r_max,phi_ord) ;
155    end
156
157    % scaling proposto da Levin
158    %Scale = diag(sqrt( phi ./ (phi + 1) ));
159    %Scale = eye(n_points);
160    %P = Scale * P ;
161    %Phi =   diag(phi + 1) ;
162    Phi =   diag(phi) ;
163%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164
165%% BUILDING MATRICES A and b:
166%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
168    A = P' * Phi * P ;
169    b = P' * Phi ;
170
171%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
172
173%% FINDING THE REGRESSION VALUE FOR THE CHOSEN POINT:  f(point) = p'(point) * a(point)
174%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
175
176%% Building p:
177    p_point = zeros(1, size(P, 2));
178    p_point(1) = 1;
179
180    tol = max(size(A)) * norm(A) * eps*10^0 ;
181    rank_A=rank(A,tol) ;
182    size_A=size(A,1) ;
183    %svd(A)
184    if rank_A < size_A
185       a = pinv(A,tol) * b ;  % it's the interpolation coefficients vector
186    else
187       a = A \ b ;            % it's the interpolation coefficients vector
188    end
189    h(coef+1:coef+n_points) = p_point * a ;
190    r(coef+1:coef+n_points) = stp*ones(n_points,1);
191    c(coef+1:coef+n_points) = Index(stp,1:n_points)';
192    %if (sum(h(coef+1:coef+n_points)) < 0.999) || (sum(h(coef+1:coef+n_points)) > 1.001)
193    %    disp('INTERPOLATION Warning:');
194    %    disp(['Interpolation point ', num2str(stp)]);
195    %    disp('the problem may be ill conditioned please try to change the interpolation parameters');
196    %end
197
198
199%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
200%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
201%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202    if comp_derivatives
203%% BUILDING REGRESSION DERIVATIVES: %%
204%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
205        p_i = zeros(dim, size(P, 2));
206        phi_i = zeros(size(phi));
207        Phi_i = zeros(size(Phi,1), size(Phi,2), dim);
208        b_i = zeros(size(b,1), size(b,2), dim);
209        gamma_i = zeros(size(a,1),size(a,2), dim);
210
211
212%% Loop for each direction:
213        for i = 1 : dim
214
215%% BUILDING DERIVATIBES FOR THE BASE p_i:
216%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
217            p_i(i,:) = f_polyn_basis_derivatives(opoints(stp,:),i,basis, opoints(stp,:)) ;
218
219
220%% BUILDING DERIVATIVES FOR THE WEIGHTS phi_i:
221%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222            adist = mean(dist(stp,1:n_points));
223            for j = 1 : size(Index,2) - 1
224                phi_i(j) = f_weight_function_derivatives(opoints(stp,:),inpoints(Index(stp,j),:),i,dist(stp,j),dim,r_max,adist,phi_ord) ;
225            end
226
227% figure;
228% plot(vars(Index),phi,'*')
229% grid on
230% figure;
231% plot(vars(Index),phi_i,'*')
232% grid on
233
234            Phi_i(:,:,i) = diag(phi_i);
235
236%% BUILDING A_i and B_i:
237%%%%%%%%%%%%%%%%%%%%%%%%
238            b_i(:,:,i) = P' * Phi_i(:,:,i);
239
240%% BUILDING DERIVATIVES FOR GAMMA gamma_i:
241%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
242            if rank_A < size_A
243                gamma_i(:,:,i) = pinv(A,tol) * b_i(:,:,i) ;
244            else
245                gamma_i(:,:,i) = A \ b_i(:,:,i) ;
246            end
247
248%% BUILDING DERIVATIVES FOR THE REGRESSION COEFFICIENTS I_i:
249%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
250            rd{i}(coef+1:coef+n_points) = stp*ones(n_points,1);
251            cd{i}(coef+1:coef+n_points) = Index(stp,1:n_points)';
252            hd{i}(coef+1:coef+n_points) = p_i(i,:) * a + p_point * ( gamma_i(:,:,i) - gamma_i(:,:,i)* P * a) ;
253
254        end
255    end
256    coef = coef + n_points;
257end
258I = sparse(r, c, h, n_qpt, n_pt);
259if comp_derivatives
260    I_i = cell(dim,1);
261    for i = 1 : dim
262        I_i{i} = sparse(rd{i}, cd{i}, hd{i}, n_qpt, n_pt);
263    end
264    varargout{1} = I_i;
265end
266