Initial commit

readme

@@ -0,0 +1,7 @@
+Contains all my matlab toolboxes. 
+
+Note that these toolboxes are not anymore developped, and you should go online to 
+	www.numerical-tours.com
+for updated codes.
+
+Copyright (c) 2013 Gabriel Peyre

startup.m

@@ -0,0 +1,38 @@
+% set the paths for all the tooboxes
+%
+%   Copyright (c) 2008 Gabriel Peyre
+
+% return;
+
+% base directory for the toolboxes
+rep = 'Desktop/dev/matlab/matlabcompanion/toolbox/';
+rep = '/Volumes/gpeyre 1/matlab/toolbox/';
+
+
+% return;
+
+if not(rep(1)=='/');
+
+% retrieve the user name
+s = pwd();
+f = strfind(s, '/');
+s = s(1:f(3));
+rep = [s rep];
+
+end
+
+rep_old = pwd;
+cd(rep);
+
+A = dir('toolbox_*');
+for i=1:length(A)
+    name = A(i).name;
+    if not(strcmp(name(end-2:end), 'zip'))
+        path(path, [rep name]);
+    end
+end
+
+% set default font size larger
+set_thicklines(2);
+
+cd(rep_old);

toolbox_alpert/build_alpert_matrix_1d.m

@@ -0,0 +1,114 @@
+function V = build_alpert_matrix(x,k)
+
+% build_alpert_matrix - build 1D Alpert basis
+%   adapted for the sampling provided.
+%
+%   V = build_alpert_matrix(pos,k);
+%
+%   'pos' is a 1D vector, the location of the points.
+%   'k' is the number of vanishing moments (1=>Haar, 2=>linear basis ...).
+%   'V' is an orthogonal matrix, each column is a basis vector.
+%
+%   Copyright (c) 2004 Gabriel Peyré
+
+n = length(x);
+x = x(:);
+
+if n<k
+    k=n;    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% first sort the sampling, and record the order !!
+[x,I] = sort(x);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% find a regroupement
+options.ptmax = k;
+[part,B,G] = dichotomic_grouping(x',options);
+
+P = length(G);          % number of groups
+
+% we have got nbr.packets = 2^(J-1)
+J = log2(P)+1;
+
+si = [0, cumsum(G)]+1;    % si(i) is the index of the 1st point of ith group
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% initialization of moments matrices
+for i = 1:P
+    % build moment matrix
+    xs = x( si(i):(si(i+1)-1) );        % selected points
+
+    kk = G(i);                  % nbr of point in this bin, ie. 2k in regular case
+    % for bookkeeping we need 2*k polynmials
+    [xJ,xI] = meshgrid(0:2*k-1,0:kk-1);
+    Mi{i} = xs(xI+1).^xJ;    
+    % for orthogonalization the matrix must be square
+    [xJ,xI] = meshgrid(0:kk-1, 0:kk-1);
+    M = xs(xI+1).^xJ;    
+    [Ui{i},R] = qr(M);    Ui{i} = transpose(Ui{i});
+end
+
+% initialization of first matrix
+U = zeros(n,n);
+for i = 1:P
+    selj = si(i):(si(i+1)-1);                   % selected column
+    % to keep : upper part is of size n-P*k
+    offs = si(i)-(i-1)*k;   % offset on row
+    long = G(i)-k;          % length on row
+    U( offs + (0:long-1), selj )       = Ui{i}((k+1):end, :);   % we keep the G(i)-k last
+    % to retransform : lower part is of size P*k
+    offs = n - P*k+1 + (i-1)*k;
+    U( offs+(0:k-1), selj )   = Ui{i}(1:k, :);   % the first k doesn't have vanishing moments, keep them   
+end
+
+for j=2:J   % for each scale
+
+    % at this scale, we have nj = P/2^(j-1) groups
+    nj = P/2^(j-1) ; % n/(2^j*k);
+    mj = nj*2*k;     % total length of the blocks
+
+    % update each sub matrix
+    for i = 1:nj
+
+        M = zeros(2*k, 2*k);
+        M(1:k,:)        = Ui{2*i-1}(1:k,:)  *   Mi{2*i-1};         % Ui^U is just k first row
+        M((k+1):2*k,:)  = Ui{2*i}(1:k,:)    *   Mi{2*i};
+        MMi{i} = M;
+        [UUi{i},R] = qr(MMi{i});    UUi{i} = transpose(UUi{i});
+
+    end
+    Mi = MMi;
+    Ui = UUi;
+
+    % lower part of the multiplicative matrix
+    UU = zeros(mj,mj);
+    for i = 1:nj
+        UU( k*(i-1)+(1:k), 2*k*(i-1)+(1:2*k) )          = mlow( Ui{i} );
+        UU( mj/2+k*(i-1)+(1:k), 2*k*(i-1)+(1:2*k) )     = mup( Ui{i} );    
+    end
+
+    % multiplicative matrix
+    Uj = eye(n);
+    Uj( (end-mj+1):end, (end-mj+1):end ) = UU;
+
+    % update vectors
+    U = Uj*U;
+end
+
+V = U';
+
+% reorder the index of the vectors !
+II = reverse_permutation(I);
+V = V(II,:);
+
+
+% extract uper part of the matrix
+function MU = mup(M)
+MU = M(1:end/2, :);
+
+% extract lower part of the matrix
+function ML = mlow(M)
+ML = M((end/2+1):end, :);

toolbox_alpert/build_alpert_matrix_2d.m

@@ -0,0 +1,156 @@
+function V = build_alpert_matrix_2d(pos,k, options)
+
+% build_alpert_matrix_2d - build 2D Alpert basis adapted to the sampling provided.
+%
+%   V = build_alpert_matrix_2d(pos,k);
+%
+%   The basis is hierarchical with respect to X direction ...
+%
+%   'pos' is a 2D vector, pos(:,i) is the ith point.
+%   'k' is the number of vanishing moments (1=>Haar, 2=>linear basis ...).
+%   'V' is an orthogonal matrix, each column is a basis vector.
+%
+%   Copyright (c) 2004 Gabriel Peyré
+
+if nargin<3
+    options.null = 0;
+end
+
+if isfield(options, 'part_type')
+    part_type = options.part_type;
+else
+    part_type = '1axis';
+end
+
+if isfield(options, 'part')
+    part = options.part;
+else
+    part = [];
+end
+
+use_scaling = 1;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% number of points
+n = size(pos,2);
+if n==0
+    V = [];
+    return;
+end
+
+if n<k^2
+    % special case, not enough data
+    k = floor(sqrt(n));
+end
+
+k2 = k^2;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% find a regroupement
+if isempty(part)
+    clear options;
+    options.ptmax = k2;
+    options.part_type = part_type;
+    [part,B,G] = dichotomic_grouping(pos(1,:),options);
+end
+
+P = length(G);          % number of groups
+
+% we have got nbr.packets = 2^(J-1)
+J = log2(P)+1;
+
+si = [0, cumsum(G)]+1;    % si(i) is the index of the 1st point of ith group
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% initialization of moments matrices
+for i = 1:P
+
+    seli = part{i};    
+    % build moment matrix
+    xs = pos( 1,seli );        % selected points, x component
+    ys = pos( 2,seli );        % selected points, y component
+
+    kk = length(seli);                   % nbr of point in this bin, ie. 2k^2 in regular case
+    M = zeros( kk, 2*k2 );
+    for ii=1:kk
+        j = 0;
+        % first part    : polynoms X^j1*Y^j2 for 0<=j1,j2<k
+        % second part   : polynoms X^j1*Y^j2 for k<=j1<2k and 0<=j2<k
+        for j1=1:2*k
+            for j2=1:k
+                j = j+1;
+                M(ii,j) = xs(ii)^(j1-1) * ys(ii)^(j2-1);        
+            end
+        end
+    end
+
+    Mi{i} = M;    
+    % orthogonalize
+    [Ui{i},R] = qr(Mi{i});    Ui{i} = transpose(Ui{i});
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% initialization of first matrix
+U = zeros(n,n);
+for i = 1:P
+    selj = part{i};                   % selected column
+    % to keep : upper part is of size n-P*k^2
+    offs = si(i)-(i-1)*k2;   % offset on row
+    long = G(i)-k2;          % length on row
+    U( offs + (0:long-1), selj )       = Ui{i}((k2+1):end, :);   % we keep the G(i)-k^2 last
+    % to retransform : lower part is of size P*k2
+    offs = n - P*k2+1 + (i-1)*k2;
+    if offs>0
+        U( offs+(0:k2-1), selj )   = Ui{i}(1:k2, :);   % the first k^2 doesn't have vanishing moments, keep them   
+    else
+        warning('Pbm empty bin.');
+        % bug ...
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+for j=2:J   % for each scale
+
+    % at this scale, we have nj = P/2^(j-1) groups
+    nj = P/2^(j-1) ; % n/(2^j*k);
+    mj = nj*2*k2;     % total length of the blocks
+
+    % update each sub matrix
+    for i = 1:nj
+
+        M = zeros(2*k2, 2*k2);
+        M(1:k2,:)           = Ui{2*i-1}(1:k2,:)  *   Mi{2*i-1};         % Ui^U is just k first row
+        M((k2+1):2*k2,:)    = Ui{2*i}(1:k2,:)    *   Mi{2*i};
+        MMi{i} = M;
+        [UUi{i},R] = qr(MMi{i});    UUi{i} = transpose(UUi{i});
+
+    end
+    Mi = MMi;
+    Ui = UUi;
+
+    % lower part of the multiplicative matrix
+    UU = zeros(mj,mj);
+    for i = 1:nj
+        UU( k2*(i-1)+(1:k2), 2*k2*(i-1)+(1:2*k2) )          = mlow( Ui{i} );
+        UU( mj/2+k2*(i-1)+(1:k2), 2*k2*(i-1)+(1:2*k2) )     = mup( Ui{i} );    
+    end
+
+    % multiplicative matrix
+    Uj = eye(n);
+    Uj( (end-mj+1):end, (end-mj+1):end ) = UU;
+
+    % update vectors
+    U = Uj*U;
+end
+
+V = U';
+
+
+% extract uper part of the matrix
+function MU = mup(M)
+MU = M(1:end/2, :);
+
+% extract lower part of the matrix
+function ML = mlow(M)
+ML = M((end/2+1):end, :);

toolbox_alpert/compile_mex.m

@@ -0,0 +1,2 @@
+% compile mex file
+mex mex/perform_moment_transform.cpp mex/perform_moment_transform_mex.cpp