Skip to content
GitLab
Commit dcdf9be0 authored by jml1g18's avatar jml1g18
Browse files

Reorganise directory structure and external dependencies

parent e360ffa4
Expand all Hide whitespace changes
Inline Side-by-side
PIV_2D_wdef.m
View file @ dcdf9be0
......@@ -79,7 +79,8 @@ function piv_result = PIV_2D_wdef(A, B, mask, piv_setup)
'win_x', empty, 'win_y', empty, ...
'ux', empty, 'uy', empty, ...
'nan_mask', empty, 'Q', empty, ...
'peak_mag', empty, 'peak_choice', empty);
'peak_mag', empty, 'peak_choice', empty, ...
'dt', empty);
%% ITERATE OVER PASSES
for pass = 1 : n_passes
......@@ -305,7 +306,8 @@ function piv_result = PIV_2D_wdef(A, B, mask, piv_setup)
'nan_mask', nan_mask, ...
'Q', Q, ...
'peak_mag', peak_mag, ...
'peak_choice', peak_choice);
'peak_choice', peak_choice, ...
'dt', piv_dt);
%% save displacement field from this pass
% into "old" displacement field structure
......
ext/inpaintn/inpaint_nans.m 0 → 100644
View file @ dcdf9be0
function B=inpaint_nans(A,method) % INPAINT_NANS: in-paints over nans in an array % usage: B=INPAINT_NANS(A) % default method % usage: B=INPAINT_NANS(A,method) % specify method used % % Solves approximation to one of several pdes to % interpolate and extrapolate holes in an array % % arguments (input): % A - nxm array with some NaNs to be filled in % % method - (OPTIONAL) scalar numeric flag - specifies % which approach (or physical metaphor to use % for the interpolation.) All methods are capable % of extrapolation, some are better than others. % There are also speed differences, as well as % accuracy differences for smooth surfaces. % % methods {0,1,2} use a simple plate metaphor. % method 3 uses a better plate equation, % but may be much slower and uses % more memory. % method 4 uses a spring metaphor. % method 5 is an 8 neighbor average, with no % rationale behind it compared to the % other methods. I do not recommend % its use. % % method == 0 --> (DEFAULT) see method 1, but % this method does not build as large of a % linear system in the case of only a few % NaNs in a large array. % Extrapolation behavior is linear. % % method == 1 --> simple approach, applies del^2 % over the entire array, then drops those parts % of the array which do not have any contact with % NaNs. Uses a least squares approach, but it % does not modify known values. % In the case of small arrays, this method is % quite fast as it does very little extra work. % Extrapolation behavior is linear. % % method == 2 --> uses del^2, but solving a direct % linear system of equations for nan elements. % This method will be the fastest possible for % large systems since it uses the sparsest % possible system of equations. Not a least % squares approach, so it may be least robust % to noise on the boundaries of any holes. % This method will also be least able to % interpolate accurately for smooth surfaces. % Extrapolation behavior is linear. % % Note: method 2 has problems in 1-d, so this % method is disabled for vector inputs. % % method == 3 --+ See method 0, but uses del^4 for % the interpolating operator. This may result % in more accurate interpolations, at some cost % in speed. % % method == 4 --+ Uses a spring metaphor. Assumes % springs (with a nominal length of zero) % connect each node with every neighbor % (horizontally, vertically and diagonally) % Since each node tries to be like its neighbors, % extrapolation is as a constant function where % this is consistent with the neighboring nodes. % % method == 5 --+ See method 2, but use an average % of the 8 nearest neighbors to any element. % This method is NOT recommended for use. % % % arguments (output): % B - nxm array with NaNs replaced % % % Example: % [x,y] = meshgrid(0:.01:1); % z0 = exp(x+y); % znan = z0; % znan(20:50,40:70) = NaN; % znan(30:90,5:10) = NaN; % znan(70:75,40:90) = NaN; % % z = inpaint_nans(znan); % % % See also: griddata, interp1 % % Author: John D'Errico % e-mail address: woodchips@rochester.rr.com % Release: 2 % Release date: 4/15/06 % I always need to know which elements are NaN, % and what size the array is for any method [n,m]=size(A); A=A(:); nm=n*m; k=isnan(A(:)); % list the nodes which are known, and which will % be interpolated nan_list=find(k); known_list=find(~k); % how many nans overall nan_count=length(nan_list); % convert NaN indices to (r,c) form % nan_list==find(k) are the unrolled (linear) indices % (row,column) form [nr,nc]=ind2sub([n,m],nan_list); % both forms of index in one array: % column 1 == unrolled index % column 2 == row index % column 3 == column index nan_list=[nan_list,nr,nc]; % supply default method if (nargin<2) || isempty(method) method = 0;elseif ~ismember(method,0:5) error 'If supplied, method must be one of: {0,1,2,3,4,5}.'end % for different methodsswitch method case 0 % The same as method == 1, except only work on those % elements which are NaN, or at least touch a NaN. % is it 1-d or 2-d? if (m == 1) || (n == 1) % really a 1-d case work_list = nan_list(:,1); work_list = unique([work_list;work_list - 1;work_list + 1]); work_list(work_list <= 1) = []; work_list(work_list >= nm) = []; nw = numel(work_list); u = (1:nw)'; fda = sparse(repmat(u,1,3),bsxfun(@plus,work_list,-1:1), ... repmat([1 -2 1],nw,1),nw,nm); else % a 2-d case % horizontal and vertical neighbors only talks_to = [-1 0;0 -1;1 0;0 1]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with second partials on row % variable for each element in either list, but only % for those nodes which have a row index > 1 or < n L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % 2nd partials on column index L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 1 % least squares approach with del^2. Build system % for every array element as an unknown, and then % eliminate those which are knowns. % Build sparse matrix approximating del^2 for % every element in A. % is it 1-d or 2-d? if (m == 1) || (n == 1) % a 1-d case u = (1:(nm-2))'; fda = sparse(repmat(u,1,3),bsxfun(@plus,u,0:2), ... repmat([1 -2 1],nm-2,1),nm-2,nm); else % a 2-d case % Compute finite difference for second partials % on row variable first [i,j]=ndgrid(2:(n-1),1:m); ind=i(:)+(j(:)-1)*n; np=(n-2)*m; fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... repmat([1 -2 1],np,1),n*m,n*m); % now second partials on column variable [i,j]=ndgrid(1:n,2:(m-1)); ind=i(:)+(j(:)-1)*n; np=n*(m-2); fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ... repmat([1 -2 1],np,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 2 % Direct solve for del^2 BVP across holes % generate sparse array with second partials on row % variable for each nan element, only for those nodes % which have a row index > 1 or < n % is it 1-d or 2-d? if (m == 1) || (n == 1) % really just a 1-d case error('Method 2 has problems for vector input. Please use another method.') else % a 2-d case L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(nan_list(L,1),1,3), ... repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),n*m,n*m); else fda=spalloc(n*m,n*m,size(nan_list,1)*5); end % 2nd partials on column index L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,3), ... repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),n*m,n*m); end % fix boundary conditions at extreme corners % of the array in case there were nans there if ismember(1,nan_list(:,1)) fda(1,[1 2 n+1])=[-2 1 1]; end if ismember(n,nan_list(:,1)) fda(n,[n, n-1,n+n])=[-2 1 1]; end if ismember(nm-n+1,nan_list(:,1)) fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1]; end if ismember(nm,nan_list(:,1)) fda(nm,[nm,nm-1,nm-n])=[-2 1 1]; end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); % and solve... B=A; k=nan_list(:,1); B(k)=fda(k,k)\rhs(k); end case 3 % The same as method == 0, except uses del^4 as the % interpolating operator. % del^4 template of neighbors talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ... 0 1;0 2;1 -1;1 0;1 1;2 0]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with del^4, but only % for those nodes which have a row & column index % >= 3 or <= n-2 L = find( (all_list(:,2) >= 3) & ... (all_list(:,2) <= (n-2)) & ... (all_list(:,3) >= 3) & ... (all_list(:,3) <= (m-2))); nl=length(L); if nl>0 % do the entire template at once fda=sparse(repmat(all_list(L,1),1,13), ... repmat(all_list(L,1),1,13) + ... repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ... repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % on the boundaries, reduce the order around the edges L = find((((all_list(:,2) == 2) | ... (all_list(:,2) == (n-1))) & ... (all_list(:,3) >= 2) & ... (all_list(:,3) <= (m-1))) | ... (((all_list(:,3) == 2) | ... (all_list(:,3) == (m-1))) & ... (all_list(:,2) >= 2) & ... (all_list(:,2) <= (n-1)))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,5), ... repmat(all_list(L,1),1,5) + ... repmat([-n,-1,0,+1,n],nl,1), ... repmat([1 1 -4 1 1],nl,1),nm,nm); end L = find( ((all_list(:,2) == 1) | ... (all_list(:,2) == n)) & ... (all_list(:,3) >= 2) & ... (all_list(:,3) <= (m-1))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3) + ... repmat([-n,0,n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end L = find( ((all_list(:,3) == 1) | ... (all_list(:,3) == m)) & ... (all_list(:,2) >= 2) & ... (all_list(:,2) <= (n-1))); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3) + ... repmat([-1,0,1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 4 % Spring analogy % interpolating operator. % list of all springs between a node and a horizontal % or vertical neighbor hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1]; hv_springs=[]; for i=1:4 hvs=nan_list+repmat(hv_list(i,:),nan_count,1); k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m); hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]]; end % delete replicate springs hv_springs=unique(sort(hv_springs,2),'rows'); % build sparse matrix of connections, springs % connecting diagonal neighbors are weaker than % the horizontal and vertical springs nhv=size(hv_springs,1); springs=sparse(repmat((1:nhv)',1,2),hv_springs, ... repmat([1 -1],nhv,1),nhv,nm); % eliminate knowns rhs=-springs(:,known_list)*A(known_list); % and solve... B=A; B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs; case 5 % Average of 8 nearest neighbors % generate sparse array to average 8 nearest neighbors % for each nan element, be careful around edges fda=spalloc(n*m,n*m,size(nan_list,1)*9); % -1,-1 L = find((nan_list(:,2) > 1) & (nan_list(:,3) > 1)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % 0,-1 L = find(nan_list(:,3) > 1); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,-1 L = find((nan_list(:,2) < n) & (nan_list(:,3) > 1)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-n+1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % -1,0 L = find(nan_list(:,2) > 1); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,0 L = find(nan_list(:,2) < n); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % -1,+1 L = find((nan_list(:,2) > 1) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n-1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % 0,+1 L = find(nan_list(:,3) < m); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % +1,+1 L = find((nan_list(:,2) < n) & (nan_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(nan_list(L,1),1,2), ... repmat(nan_list(L,1),1,2)+repmat([n+1, 0],nl,1), ... repmat([1 -1],nl,1),n*m,n*m); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); % and solve... B=A; k=nan_list(:,1); B(k)=fda(k,k)\rhs(k); end % all done, make sure that B is the same shape as % A was when we came in. B=reshape(B,n,m); % ==================================================== % end of main function % ==================================================== % ==================================================== % begin subfunctions % ==================================================== function neighbors_list=identify_neighbors(n,m,nan_list,talks_to) % identify_neighbors: identifies all the neighbors of % those nodes in nan_list, not including the nans % themselves % % arguments (input): % n,m - scalar - [n,m]=size(A), where A is the % array to be interpolated % nan_list - array - list of every nan element in A % nan_list(i,1) == linear index of i'th nan element % nan_list(i,2) == row index of i'th nan element % nan_list(i,3) == column index of i'th nan element % talks_to - px2 array - defines which nodes communicate % with each other, i.e., which nodes are neighbors. % % talks_to(i,1) - defines the offset in the row % dimension of a neighbor % talks_to(i,2) - defines the offset in the column % dimension of a neighbor % % For example, talks_to = [-1 0;0 -1;1 0;0 1] % means that each node talks only to its immediate % neighbors horizontally and vertically. % % arguments(output): % neighbors_list - array - list of all neighbors of % all the nodes in nan_list if ~isempty(nan_list) % use the definition of a neighbor in talks_to nan_count=size(nan_list,1); talk_count=size(talks_to,1); nn=zeros(nan_count*talk_count,2); j=[1,nan_count]; for i=1:talk_count nn(j(1):j(2),:)=nan_list(:,2:3) + ... repmat(talks_to(i,:),nan_count,1); j=j+nan_count; end % drop those nodes which fall outside the bounds of the % original array L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m); nn(L,:)=[]; % form the same format 3 column array as nan_list neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn]; % delete replicates in the neighbors list neighbors_list=unique(neighbors_list,'rows'); % and delete those which are also in the list of NaNs. neighbors_list=setdiff(neighbors_list,nan_list,'rows'); else neighbors_list=[]; end
\ No newline at end of file
ext/inpaintn/inpaint_nans_bc.m 0 → 100644
View file @ dcdf9be0
function B=inpaint_nans_bc(A,method,bcclass) % INPAINT_NANS_BC: in-paints over nans in an array, with spherical or toroidal boundary conditions % usage: B=inpaint_nsns_bc(A) % default method % usage: B=inpaint_nsns_bc(A,method) % specify method used % usage: B=inpaint_nsns_bc(A,method,bcclass) % specify class of boundary conditions applied % % Solves approximation to one of several pdes to % interpolate and extrapolate holes in an array. % Depending upon the boundary conditions specified, % the array will effectively be treated as if it lies % on either the surface of a sphere or a toroid. % % arguments (input): % A - nxm array with some NaNs to be filled in % % method - (OPTIONAL) scalar numeric flag - specifies % which approach (or physical metaphor to use % for the interpolation.) All methods are capable % of extrapolation, some are better than others. % There are also speed differences, as well as % accuracy differences for smooth surfaces. % % The methods employed here are a subset of the % methods of the original inpaint_nans. % % methods {0,1} use a simple plate metaphor. % method 4 uses a spring metaphor. % % method == 0 --> (DEFAULT) see method 1, but % this method does not build as large of a % linear system in the case of only a few % NaNs in a large array. % Extrapolation behavior is linear. % % method == 1 --> simple approach, applies del^2 % over the entire array, then drops those parts % of the array which do not have any contact with % NaNs. Uses a least squares approach, but it % does not modify known values. % In the case of small arrays, this method is % quite fast as it does very little extra work. % Extrapolation behavior is linear. % % method == 4 --> Uses a spring metaphor. Assumes % springs (with a nominal length of zero) % connect each node with every neighbor % (horizontally, vertically and diagonally) % Since each node tries to be like its neighbors, % extrapolation is as a constant function where % this is consistent with the neighboring nodes. % % DEFAULT: 0 % % bcclass - (OPTIONAL) character flag, indicating how % the array boundaries will be treated in the % inpainting operation. bcclass may be either % 'sphere' or 'toroid', or any simple contraction % of these words. % % bcclass = 'sphere' --> The first and last rows % of the array will be treated as if they are % at the North and South poles of a sphere. % Adjacent to those rows will be singular % phantom nodes at each pole. % % bcclass = 'toroid' --> The first and last rows % of the array will be treated as if they are % adjacent to ech other. As well, the first and % last columns will be adjacent to each other. % % DEFAULT: 'sphere' % % arguments (output): % B - nxm array with NaNs replaced % % % Example: % [x,y] = meshgrid(0:.01:1); % z0 = exp(x+y); % znan = z0; % znan(20:50,40:70) = NaN; % znan(30:90,5:10) = NaN; % znan(70:75,40:90) = NaN; % % z = inpaint_nans(znan); % % % See also: griddata, interp1 % % Author: John D'Errico % e-mail address: woodchips@rochester.rr.com % Release: 2 % Release date: 4/15/06 % I always need to know which elements are NaN, % and what size the array is for any method [n,m]=size(A); A=A(:); nm=n*m; k=isnan(A(:)); % list those nodes which are known, and which will % be interpolated nan_list=find(k); known_list=find(~k); % how many nans overall nan_count=length(nan_list); % convert NaN indices to (r,c) form % nan_list==find(k) are the unrolled (linear) indices % (row,column) form [nr,nc]=ind2sub([n,m],nan_list); % both forms of index in one array: % column 1 == unrolled index % column 2 == row index % column 3 == column index nan_list=[nan_list,nr,nc]; % supply default method if (nargin<2) || isempty(method) method = 0;elseif ~ismember(method,[0 1 4]) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, method must be one of: {0,1,4}.')end % supply default value for bcclassif (nargin < 3) || isempty(bcclass) bcclass = 'sphere';elseif ~ischar(bcclass) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, bcclass must be ''sphere'' or ''toroid''')else % it was a character string valid = {'sphere' 'toroid'}; % check to see if it is valid [bcclass,errorclass] = validstring(arg,valid); if ~isempty(errorclass) error('INPAINT_NANS_BC:improperargument', ... 'If supplied, bcclass must be ''sphere'' or ''toroid''') endend % choice of methodsswitch method case 0 % The same as method == 1, except only work on those % elements which are NaN, or at least touch a NaN. % horizontal and vertical neighbors only talks_to = [-1 0;0 -1;1 0;0 1]; neighbors_list=identify_neighbors(n,m,nan_list,talks_to); % list of all nodes we have identified all_list=[nan_list;neighbors_list]; % generate sparse array with second partials on row % variable for each element in either list, but only % for those nodes which have a row index > 1 or < n L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); nl=length(L); if nl>0 fda=sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); else fda=spalloc(n*m,n*m,size(all_list,1)*5); end % 2nd partials on column index L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); nl=length(L); if nl>0 fda=fda+sparse(repmat(all_list(L,1),1,3), ... repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ... repmat([1 -2 1],nl,1),nm,nm); end % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list(:,1)),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 1 % least squares approach with del^2. Build system % for every array element as an unknown, and then % eliminate those which are knowns. % Build sparse matrix approximating del^2 for % every element in A. % Compute finite difference for second partials % on row variable first [i,j]=ndgrid(1:n,1:m); ind=i(:)+(j(:)-1)*n; np=n*m; switch bcclass case 'sphere' % we need to have two phantom nodes at the poles np = np + 2; end fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... repmat([1 -2 1],np,1),n*m,n*m); % now second partials on column variable [i,j]=ndgrid(1:n,2:(m-1)); ind=i(:)+(j(:)-1)*n; np=n*(m-2); fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ... repmat([1 -2 1],np,1),nm,nm); % eliminate knowns rhs=-fda(:,known_list)*A(known_list); k=find(any(fda(:,nan_list),2)); % and solve... B=A; B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k); case 4 % Spring analogy % interpolating operator. % list of all springs between a node and a horizontal % or vertical neighbor hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1]; hv_springs=[]; for i=1:4 hvs=nan_list+repmat(hv_list(i,:),nan_count,1); k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m); hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]]; end % delete replicate springs hv_springs=unique(sort(hv_springs,2),'rows'); % build sparse matrix of connections, springs % connecting diagonal neighbors are weaker than % the horizontal and vertical springs nhv=size(hv_springs,1); springs=sparse(repmat((1:nhv)',1,2),hv_springs, ... repmat([1 -1],nhv,1),nhv,nm); % eliminate knowns rhs=-springs(:,known_list)*A(known_list); % and solve... B=A; B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs; end % all done, make sure that B is the same shape as% A was when we came in.B=reshape(B,n,m); end % mainline % ====================================================% end of main function% ====================================================% ====================================================% begin subfunctions% ====================================================function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)% identify_neighbors: identifies all the neighbors of% those nodes in nan_list, not including the nans% themselves%% arguments (input):% n,m - scalar - [n,m]=size(A), where A is the% array to be interpolated% nan_list - array - list of every nan element in A% nan_list(i,1) == linear index of i'th nan element% nan_list(i,2) == row index of i'th nan element% nan_list(i,3) == column index of i'th nan element% talks_to - px2 array - defines which nodes communicate% with each other, i.e., which nodes are neighbors.%% talks_to(i,1) - defines the offset in the row% dimension of a neighbor% talks_to(i,2) - defines the offset in the column% dimension of a neighbor% % For example, talks_to = [-1 0;0 -1;1 0;0 1]% means that each node talks only to its immediate% neighbors horizontally and vertically.% % arguments(output):% neighbors_list - array - list of all neighbors of% all the nodes in nan_list if ~isempty(nan_list) % use the definition of a neighbor in talks_to nan_count=size(nan_list,1); talk_count=size(talks_to,1); nn=zeros(nan_count*talk_count,2); j=[1,nan_count]; for i=1:talk_count nn(j(1):j(2),:)=nan_list(:,2:3) + ... repmat(talks_to(i,:),nan_count,1); j=j+nan_count; end % form the same format 3 column array as nan_list neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn]; % delete replicates in the neighbors list neighbors_list=unique(neighbors_list,'rows'); % and delete those which are also in the list of NaNs. neighbors_list=setdiff(neighbors_list,nan_list,'rows'); else neighbors_list=[];end end % function identify_neighbors function [str,errorclass] = validstring(arg,valid)% validstring: compares a string against a set of valid options% usage: [str,errorclass] = validstring(arg,valid)%% If a direct hit, or any unambiguous shortening is found, that% string is returned. Capitalization is ignored.%% arguments: (input)% arg - character string, to be tested against a list% of valid choices. Capitalization is ignored.%% valid - cellstring array of alternative choices%% Arguments: (output)% str - string - resulting choice resolved from the% list of valid arguments. If no unambiguous% choice can be resolved, then str will be empty.%% errorclass - string - A string argument that explains% the error. It will be one of the following% possibilities:%% '' --> No error. An unambiguous match for arg % was found among the choices. % % 'No match found' --> No match was found among % the choices provided in valid. % % 'Ambiguous argument' --> At least two ambiguous % matches were found among those provided % in valid. % % % Example: % valid = {'off' 'on' 'The sky is falling'} % % % See also: parse_pv_pairs, strmatch, strcmpi % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 3/25/2010 ind = strmatch(lower(arg),lower(valid)); if isempty(ind) % No hit found errorclass = 'No match found'; str = ''; elseif (length(ind) > 1) % Ambiguous arg, hitting more than one of the valid options errorclass = 'Ambiguous argument'; str = ''; return else errorclass = ''; str = valid{ind}; end end % function validstring
\ No newline at end of file
ext/inpaintn/inpaint_nans_demo.m 0 → 100644
View file @ dcdf9be0
%% Surface Fit Artifact Removal
%% Construct the Surface
[x,y] = meshgrid(0:.01:1);
z0 = exp(x+y);
close all
figure
surf(z0)
title 'Original surface'
znan = z0;
znan(20:50,40:70) = NaN;
znan(30:90,5:10) = NaN;
znan(70:75,40:90) = NaN;
figure
surf(znan)
title 'Artifacts (large holes) in surface'
%% In-paint Over NaNs
z = inpaint_nans(znan,3);
figure
surf(z)
title 'Inpainted surface'
figure
surf(z-z0)
title 'Inpainting error surface (Note z-axis scale)'
%% Comapre to GRIDDATA
k = isnan(znan);
zk = griddata(x(~k),y(~k),z(~k),x(k),y(k));
zg = znan;
zg(k) = zk;
figure
surf(zg)
title(['Griddata inpainting (',num2str(sum(isnan(zg(:)))),' NaNs remain)'])
figure
surf(zg-z0)
title 'Griddata error surface'
ext/mmx/build_mmx.m 0 → 100644
View file @ dcdf9be0
function build_mmx(verbose)
% BUILD_MMX - compiles mmx() for different platforms and provides help
% regarding compilation.
%
% BUILD_MMX will try to compile, in this order, 3 different builds of mmx:
% mmx_mkl_single - linked to Intel's single-threaded MKL library (usually fastest)
% mmx_mkl_multi - linked to the multithreaded BLAS/LAPACK libraries that come
% with Matlab.
% mmx_naive - does not link to anything, uses simple C-loops.
%
% The first time BUILD_MMX succeeds, it will compile again to 'mmx', so
% that the mex-file mmx should be the fastest possible build on your
% system.
%
% BUILD_MMX has been tested on Win32, Win64, OSX, Linux 64
%
% %% FOR LINUX OR MAC SYSTEMS:
%
% To properly link to Intel's MKL, user needs to repackage their libraries
% into one single statically linked library. The instructions are as
% follows:
%
%
% Download Intel MKL for Linux here:
% http://software.intel.com/en-us/articles/non-commercial-software-download/
%
% Donwload Intel MKL for Mac here:
% https://registrationcenter.intel.com/RegCenter/AutoGen.aspx?ProductID=1518&AccountID=&EmailID=&ProgramID=&RequestDt=&rm=EVAL&lang=
%
% The Default installation directory for both Linux and Mac will be
% /opt/intel/
% with the MKL libraries in /opt/intel/mkl
%
% %% To build needed static Library
% assuming default installation directory
%
% Run the following commands in Linux/Mac terminal:
%
% sudo -s
% cd /opt/intel/mkl/tools/builder
% cat blas_example_list > blas_lapack_list
% cat lapack_example_list >> blas_lapack_list
%
% For Linux 64 bit:
% make libintel64 interface=ilp64 export=blas_lapack_list name=libsingle_mkl_ilp64 threading=sequential
% For Linux 32 bit:
% make libia32 interface=lp64 export=blas_lapack_list name=libsingle_mkl_32 threading=sequential
%
% For Mac:
% make libuni interface=ilp64 export=blas_lapack_list name=libsingle_mkl_ilp64 threading=sequential
%
% A new libsingle_mkl_ilp64.so, libsingle_mkl_32.so, or
% libsingle_mkl_ilp64.dylib will appear.
% This needs to be copied to Matlab's external libraries directory.
%
% For Mac:
% cp libsingle_mkl_ilp64* MATLAB_ROOT/extern/lib/maci64
%
% For Linux 64 bit:
% cp libsingle_mkl_ilp64* MATLAB_ROOT/extern/lib/glnxa64
% For Linux 32 bit:
% cp libsingle_mkl_32* MATLAB_ROOT/extern/lib/glnx86
%
% Where MATLAB_ROOT is the installation directory of your Matlab.
if nargin == 0
verbose = false;
end
clc
build_names = {'mmx_mkl_single', 'mmx_mkl_multi','mmx_naive'};
built_mmx = false;
arch = computer('arch');
for b = 2
name = build_names{b};
[link, define] = deal({});
[inc_dir, link_dir, Cflags, Lflags] = deal('');
switch arch
case {'win64','win32'}
switch name
case 'mmx_naive'
define = {'WIN_SYSTEM'};
case 'mmx_mkl_multi'
root = matlabroot;
if strcmp(arch,'win32')
inc_dir = [root '\extern\lib\win32\microsoft'];
else
inc_dir = [root '\extern\lib\win64\microsoft'];
end
link = {'libmwblas','libmwlapack'};
define = {'WIN_SYSTEM','USE_BLAS'};
case 'mmx_mkl_single'
root = 'C:\Program Files (x86)\Intel\Composer XE 2011 SP1\mkl';
inc_dir = [root '\include'];
if strcmp(arch,'win32')
link_dir = [root '\lib\ia32'];
link = {'mkl_intel_c','mkl_sequential','mkl_core'};
define = {'WIN_SYSTEM','USE_BLAS','MKL_32'};
else
link_dir = [root '\lib\intel64'];
link = {'mkl_intel_ilp64','mkl_sequential','mkl_core'};
define = {'WIN_SYSTEM','USE_BLAS','MKL_ILP64'};
end
end
case {'glnxa64','glnx86'}
switch name
case 'mmx_naive'
link = {'pthread'};
define = {'UNIX_SYSTEM'};
case 'mmx_mkl_multi'
if strcmp(arch,'glnx86')
inc_dir = [matlabroot '/extern/lib/glnx86'];
else
inc_dir = [matlabroot '/extern/lib/glnxa64'];
end
link = {'mwblas','mwlapack','pthread'};
define = {'UNIX_SYSTEM','USE_BLAS'};
case 'mmx_mkl_single'
root = '/opt/intel/mkl';
inc_dir = [ root '/include'];
if strcmp(arch,'glnx86')
link_dir = [matlabroot '/extern/lib/glnx86'];
link = {'single_mkl_32','pthread'};
define = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_32'};
else
link_dir = [matlabroot '/extern/lib/glnxa64'];
link = {'small_mkl_ilp64','pthread'};
define = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_ILP64'};
end
end
case {'maci64'}
switch name
case 'mmx_naive'
link = {'pthread'};
define = {'UNIX_SYSTEM'};
case 'mmx_mkl_multi'
root = matlabroot;
inc_dir = [root '/extern/lib/maci64'];
link = {'mwblas','mwlapack','pthread'};
define = {'UNIX_SYSTEM','USE_BLAS'};
case 'mmx_mkl_single'
root = '/opt/intel/mkl';
inc_dir = [ root '/include'];
link_dir = [matlabroot '/extern/lib/maci64'];
link = {'single_mkl_ilp64','pthread'};
%link = {'small_mkl_ilp64','pthread'};
define = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_ILP64'};
end
otherwise
error unsupported_architecture
end
if ~isempty(link_dir)
if strcmp(arch,'glnxa64') || strcmp(arch,'maci64')
L_dir = {['LDFLAGS="\$LDFLAGS -L' link_dir ' ' Lflags '"']};
else
L_dir = {['-L' link_dir]};
end
else
L_dir = {};
end
if ~isempty(inc_dir)
if strcmp(arch,'glnxa64') || strcmp(arch,'maci64')
I_dir = {['CXXFLAGS="\$CXXFLAGS -I' inc_dir ' ' Cflags '"']};
else
I_dir = {['-I' inc_dir]};
end
else
I_dir = {};
end
prefix = @(pref,str_array) cellfun(@(x)[pref x],str_array,'UniformOutput',0);
l_link = prefix('-l',link);
D_define = prefix('-D',define);
if verbose
verb = {'-v'};
else
verb = {};
end
try
check_dir(link_dir, link)
check_dir(inc_dir)
clear(name)
command = {verb{:}, I_dir{:}, L_dir{:}, l_link{:}, D_define{:}}; %#ok<*CCAT>
fprintf('==========\nTrying to compile ''%s'', using \n',name);
fprintf('%s, ',command{:})
fprintf('\n')
mex(command{:}, '-output', name, 'mmx.cpp');
fprintf('Compilation of ''%s'' succeeded.\n',name);
if ~built_mmx
fprintf('Compiling again to ''mmx'' target using ''%s'' build.\n',name);
mex(command{:}, '-output','mmx','mmx.cpp');
built_mmx = true;
end
catch err
fprintf('Compilation of ''%s'' failed with error:\n%s\n',name,err.message);
end
end
function check_dir(dir,files)
if ~isempty(dir)
here = cd(dir);
if nargin == 2
for i = 1:size(files)
if isempty(ls(['*' files{i} '.*']))
cd(here);
error('could not find file %s', files{i});
end
end
end
cd(here);
end
ext/mmx/compare_chol_flops.m 0 → 100644
View file @ dcdf9be0
%% === compare timings for CHOL
R = 20; % number of repeats (to increase accuracy)
nn = 1:1:200;
K = length(nn);
time = zeros(K,3,R);
ops = zeros(K,1);
fun = {
@(a)mmx('chol',a,[]),...
@(a)mmx_MMKL('chol',a,[]),...
@(a)mmx_C('chol',a,[])%,...
% @(a)ndfun('chol',a)
};
names = {
'mmx Intel',...
'mmx MKL',...
'mmx Emo'%,...
% 'ndfun'
};
for i = 1:K
n = nn(i);
r1 = n;
c1 = n;
N = round(1e6 / (r1*c1));
A = randn(n,n,N);
A = mmx('s',A,[]);
AA = A;
ops(i) = (n^3/3 + n^2/2 +n/6)*N;
fprintf('cholesky algorithm on %d SPD matrices of size [%dx%d]\n',N,n,n)
for r = 1:R
for f = 1:length(fun)
tic;
C = fun{f}(AA);
time(i,f,r) = toc;
AA = A;
end
end
end
%%
gflops = 1e-6*bsxfun(@times, ops, 1./time);
mFLP = mean(mflops,3);
clf
hold on
cols = get(gca,'ColorOrder');
for i = 1:size(mflops,2)
plot(nn,mFLP(:,i),'color',cols(i,:),'linewidth',3);
end
for i = 1:size(mflops,2)
plot(nn,squeeze(mflops(:,i,:)),'.','color',cols(i,:),'linewidth',2)
end
grid on
set(gca,'xlim',[nn(1) nn(end)])
ylabel('Mflops')
xlabel('dimension')
legend(names{:},'location','northwest')
title('\bf Speed Comparison between mmx and ndfun. (in Mflops, bigger is better)')
drawnow;
ext/mmx/compare_mult_flops.m 0 → 100644
View file @ dcdf9be0
%% === compare timings for MULT
choice = questdlg('Run full comparison (several minutes)?', ...
'Performace Comparison', 'Yes','No','Yes');
if strcmp(choice,'No')
return
end
% number of repeats (to increase measurement accuracy)
R = 10;
% what sizes of (square) matrices do we want measure? (log spacing)
nn = unique(round(exp(linspace(log(1),log(120),100))));
% number of different sizes to try
N = length(nn);
% size of the output C, in Megabytes
Mbytes = 10;
time = zeros(N,3,R);
ops = zeros(N,1);
candidates = {'mmx, single-threaded BLAS', 'mmx_mkl_single', @(a,b)mmx_mkl_single('m',a,b);
'mmx, multi-threaded BLAS', 'mmx_mkl_multi', @(a,b)mmx_mkl_multi('m',a,b);
'mmx, naive loops', 'mmx_naive', @(a,b)mmx_naive('m',a,b);
'mtimesx', 'mtimesx', @(a,b)mtimesx(a,b,'speedomp');
'ndfun', 'ndfun', @(a,b)ndfun('mult',a,b);
'Matlab ''for'' loop', 'matlab_mprod', @(a,b)matlab_mprod(a,b)...
};
funcs = cell(0,0);
for i=1:size(candidates,1)
found_it = ~isempty(which(candidates{i,2}));
if found_it
funcs(end+1,:) = candidates(i,[1 3]); %#ok<SAGROW>
else
fprintf('Could not find %s in your path, ignoring.\n',candidates{i,2});
end
end
nf = length(funcs);
for i = 1:N
n = nn(i);
r1 = n;
c1 = n;
r2 = c1;
c2 = n;
pages = round(1e6*Mbytes / (8*r1*c2));
A = randn(r1,c1,pages);
B = randn(r2,c2,pages);
ops(i) = r1*c2*(2*c1-1)*pages;
fprintf('multiplying %d*%d matrix pairs of dimension [%dx%d]\n',R*nf,pages,n,n)
for r = 1:R
for f = 1:nf %nf:-1:1
tstart = tic;
C = funcs{length(funcs)+1-f,2}(A,B);
time(i,length(funcs)+1-f,r) = toc(tstart);
pause(1/10000);
end
end
end
%% graphics
gflops = 1e-9*bsxfun(@times, ops, 1./time);
gFLP = mean(gflops,3);
clf
hold on
cols = get(gca,'ColorOrder');
for i = 1:size(gflops,2)
plot(nn,gFLP(:,i),'color',cols(i,:),'linewidth',3);
end
for i = 1:size(gflops,2)
plot(nn,squeeze(gflops(:,i,:)),'.','color',cols(i,:))
end
grid on
set(gca,'xlim',[1 max(nn)])
ylabel('\bf Gigaflops')
xlabel('\bf dimension')
legend(funcs{:,1},'location','northwest')
title('\bf Comparison between mmx, ndfun, and mtimesx. (in Gflops, bigger is better)')
%% make PDF
% set(gcf,'color','w')
% export_fig 'comparison' -png
\ No newline at end of file
ext/mmx/html/mmx_web.html 0 → 100644
View file @ dcdf9be0
This diff is collapsed.
ext/mmx/license.txt 0 → 100644
View file @ dcdf9be0
Copyright (c) 2012, Yuval
Copyright (c) 2011, James Tursa
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
ext/mmx/matlab_mprod.m 0 → 100644
View file @ dcdf9be0
function c = matlab_mprod(a,b)
c = zeros(size(a,1),size(b,2),size(a,3));
for i = 1:size(a,3)
c(:,:,i) = a(:,:,i)*b(:,:,i);
end
\ No newline at end of file
ext/mmx/matrix_fun.c 0 → 100644
View file @ dcdf9be0
// ==================================================
// straightforward implementations of matrix multiply
// ==================================================
void multAB(double* C, double* A, double* B, const int rA, const int cA, const int cB) {
int i, j, k;
double *a, *c, tmp;
for( i=0; i<cB; i++ ){
for( k=0; k<cA; k++ ){
tmp = B[i*cA+k];
a = A + k*rA;
c = C + i*rA;
for( j=0; j<rA; j++ ){
c[j] += tmp * a[j];
}
}
}
}
void multAtB(double* C, double* A, double* B, const int rA, const int cA, const int cB) {
int i, j, k;
double *a, *b, *c;
for( i=0; i<cB; i++ ){
for( k=0; k<cA; k++ ){
a = A + k*rA;
b = B + i*rA;
c = C + i*cA + k;
for( j=0; j<rA; j++ ){
(*c) += a[j]*b[j];
}
}
}
}
void multABt(double* C, double* A, double* B, const int rA, const int cA, const int rB) {
int i, j, k;
double *a, *b;
for( j=0; j<cA; j++ ){
a = A + j*rA;
b = B + j*rB;
for( i=0; i<rB; i++ ){
for( k=0; k<rA; k++ ){
C[i*rA + k] += a[k]*b[i];
}
}
}
}
void multAtBt(double* C, double* A, double* B, const int rA, const int cA, const int rB) {
int i, j, k;
double *b, *c, tmp;
for( i=0; i<cA; i++ ){
for( k=0; k<rA; k++ ){
tmp = A[i*rA+k];