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