function d2img = draw_2Dskeleton_MSRA(tline,frameId,folderId,modelId,coord_pixel)

    db_path = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/MSRA/cvpr15_MSRAHandGestureDB/';
    jointNum =  21;
    cubicSz = 200;
    imgWidth = 320;
    imgHeight = 240;
    folder_list = ["1","2","3","4","5","6","7","8","9","I","IP","L","MP","RP","T","TIP","Y"];

    coord_pixel = squeeze(coord_pixel);
    rgb_img = ones(imgHeight,imgWidth,3);
    refDepths = zeros(1,jointNum);
    line_width = 4;

    bin_name = strcat(db_path,'P',num2str(modelId),'/',folder_list(folderId),'/',sprintf('%06d',frameId-1),'_depth.bin');
    splitted = strsplit(tline);

    for jid = 1:jointNum
        refDepths(1,jid) = -str2num(splitted{(jid-1)*3+3});
    end
    refDepth = (min(refDepths(:)) + max(refDepths(:)))/2;

    fp_bin = fopen(bin_name,'r');
    img_info = fread(fp_bin,6,'int');
    left = img_info(3);
    top = img_info(4);
    right = img_info(5)-1;
    bottom = img_info(6)-1;

    img = fread(fp_bin,[right-left+1 bottom-top+1],'float');
    img = permute(img,[2,1]);
    img(img==0) = refDepth+cubicSz/2;
    fclose(fp_bin);

    img = img - refDepth;
    img = img/(cubicSz/2);
    img(img>1) = 1;
    img(img<-1) = -1;

    img = img + 1;
    img = img/2;

    rgb_img(:,:,1) = 255/255;
    rgb_img(:,:,2) = 240/255;
    rgb_img(:,:,3) = 204/255;

    rgb_img(top:bottom,left:right,1) = img*255/255;
    rgb_img(top:bottom,left:right,2) = img*240/255;
    rgb_img(top:bottom,left:right,3) = img*204/255;

    f = figure;
    set(f, 'visible', 'off');
    imshow(rgb_img);
    hold on;

    plot([coord_pixel(1,1),coord_pixel(1,18)],[coord_pixel(2,1),coord_pixel(2,18)],'Color',[255/255,153/255,153/255],'LineWidth',line_width) %wrist to thumb mcp
    plot([coord_pixel(1,18),coord_pixel(1,19)],[coord_pixel(2,18),coord_pixel(2,19)],'Color',[255/255,102/255,102/255],'LineWidth',line_width) %thumb mcp to thumb pip
    plot([coord_pixel(1,19),coord_pixel(1,20)],[coord_pixel(2,19),coord_pixel(2,20)],'Color',[255/255,51/255,51/255],'LineWidth',line_width) %thumb pip to thumb dip
    plot([coord_pixel(1,20),coord_pixel(1,21)],[coord_pixel(2,20),coord_pixel(2,21)],'Color',[255/255,0/255,0/255],'LineWidth',line_width) %thumb dip to thumb tip

    plot([coord_pixel(1,1),coord_pixel(1,2)],[coord_pixel(2,1),coord_pixel(2,2)],'Color',[153/255,255/255,153/255],'LineWidth',line_width) %wrist to index mcp
    plot([coord_pixel(1,2),coord_pixel(1,3)],[coord_pixel(2,2),coord_pixel(2,3)],'Color',[102/255,255/255,102/255],'LineWidth',line_width) %index mcp to index pip
    plot([coord_pixel(1,3),coord_pixel(1,4)],[coord_pixel(2,3),coord_pixel(2,4)],'Color',[51/255,255/255,51/255],'LineWidth',line_width) %index pip to index dip
    plot([coord_pixel(1,4),coord_pixel(1,5)],[coord_pixel(2,4),coord_pixel(2,5)],'Color',[0/255,255/255,0/255],'LineWidth',line_width) %index dip to index tip

    plot([coord_pixel(1,1),coord_pixel(1,6)],[coord_pixel(2,1),coord_pixel(2,6)],'Color',[255/255,204/255,153/255],'LineWidth',line_width) %wrist to middle mcp
    plot([coord_pixel(1,6),coord_pixel(1,7)],[coord_pixel(2,6),coord_pixel(2,7)],'Color',[255/255,178/255,102/255],'LineWidth',line_width) %middle mcp to middle pip
    plot([coord_pixel(1,7),coord_pixel(1,8)],[coord_pixel(2,7),coord_pixel(2,8)],'Color',[255/255,153/255,51/255],'LineWidth',line_width) %middle pip to middle dip
    plot([coord_pixel(1,8),coord_pixel(1,9)],[coord_pixel(2,8),coord_pixel(2,9)],'Color',[255/255,128/255,0/255],'LineWidth',line_width) %middle dip to middle tip

    plot([coord_pixel(1,1),coord_pixel(1,10)],[coord_pixel(2,1),coord_pixel(2,10)],'Color',[153/255,204/255,255/255],'LineWidth',line_width) %wrist to ring mcp
    plot([coord_pixel(1,10),coord_pixel(1,11)],[coord_pixel(2,10),coord_pixel(2,11)],'Color',[102/255,178/255,255/255],'LineWidth',line_width) %ring mcp to ring pip
    plot([coord_pixel(1,11),coord_pixel(1,12)],[coord_pixel(2,11),coord_pixel(2,12)],'Color',[51/255,153/255,255/255],'LineWidth',line_width) %ring pip to ring dip
    plot([coord_pixel(1,12),coord_pixel(1,13)],[coord_pixel(2,12),coord_pixel(2,13)],'Color',[0/255,128/255,255/255],'LineWidth',line_width) %ring dip to ring tip

    plot([coord_pixel(1,1),coord_pixel(1,14)],[coord_pixel(2,1),coord_pixel(2,14)],'Color',[255/255,153/255,255/255],'LineWidth',line_width) %wrist to little mcp
    plot([coord_pixel(1,14),coord_pixel(1,15)],[coord_pixel(2,14),coord_pixel(2,15)],'Color',[255/255,102/255,255/255],'LineWidth',line_width) %little mcp to little pip
    plot([coord_pixel(1,15),coord_pixel(1,16)],[coord_pixel(2,15),coord_pixel(2,16)],'Color',[255/255,51/255,255/255],'LineWidth',line_width) %little pip to little dip
    plot([coord_pixel(1,16),coord_pixel(1,17)],[coord_pixel(2,16),coord_pixel(2,17)],'Color',[255/255,0/255,255/255],'LineWidth',line_width) %little dip to little tip

    colorList = [
    230/255 230/255 0/255;

    153/255 255/255 153/255;
    102/255 255/255 102/255;
    51/255 255/255 51/255;
    0/255 255/255 0/255;

    255/255 204/255 153/255;
    255/255 178/255 102/255;
    255/255 153/255 51/255;
    255/255 128/255 0/255;

    153/255 204/255 255/255;
    102/255 178/255 255/255;
    51/255 153/255 255/255;
    0/255 128/255 255/255;

    255/255 153/255 255/255;
    255/255 102/255 255/255;
    255/255 51/255 255/255;
    255/255 0/255 255/255;

    255/255 153/255 153/255;
    255/255 102/255 102/255;
    255/255 51/255 51/255;
    255/255 0/255 0/255;
    ];

    scatter(coord_pixel(1,:),coord_pixel(2,:),100,colorList,'filled');

    set(gca,'Units','normalized','Position',[0 0 1 1]);  %# Modify axes size
    set(gcf,'Units','pixels','Position',[200 200 2*imgWidth 2*imgHeight]);  %# Modify figure size

    frame = getframe(gcf);
    framedata = frame.cdata;

    xmin = min(coord_pixel(1,:));
    xmax = max(coord_pixel(1,:));
    ymin = min(coord_pixel(2,:));
    ymax = max(coord_pixel(2,:));

    len = max(xmax-xmin+1,ymax-ymin+1) + 70;
    xcenter = (xmin + xmax)/2;
    ycenter = (ymin + ymax)/2;

    xmin = max(round(xcenter - len/2),1);
    xmax = min(round(xmin + len),imgWidth);
    ymin = max(round(ycenter - len/2),1);
    ymax = min(round(ymin + len),imgHeight);

    framedata = framedata(2*ymin:2*ymax,2*xmin:2*xmax,:);

    hold off;
    close(f);

    d2img = framedata;

end


text:
% This routine do:
% 1. Discretize schrodinger equation over the specified set of points:
%     v(x)=2m*V(x)/hbar^2; ee=2m*E/hbar^2.
%     [    d^2       2m        ]           2m*E
%     [ - -----  + ------ V(x) ] phi(x) = ------  phi(x)
%     [   d x^2    hbar^2      ]          hbar^2
%
% 2. Find eigen values and eigen vectors.
% 3. Plot the found eigen values and eigen vectors.
%
% Usage:
% [ee,ev] = qm1d_fast(NPTS,NSTM,a,b,f_pot_handle)
% NPTS - number of points for discretization of schrodinger equation
% NSTM - number of eigen values and eigen vector to find
% a - the start point of te interval for x
% b - the end point of the interval for y
% f_pot_handle - handle to function which defines the potential
%
% Examples:
%
% 1. Harmonic oscillator - symetric at x=L/2
% NPTS=1000;
% NSTM=5;
% L=10d0;
% f_pot = @(x) {4d0*(x-L/2).^2};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% 2. Harmonic oscillator - symetric at x=0
% NPTS=1000;
% NSTM=5;
% L=10d0;
% f_pot = @(x) {4d0*x.^2};
% qm1d_fast(NPTS,NSTM,-L/2,L/2,f_pot);
%
% 3. Barrier - defined with heaviside function
% NPTS=1000;
% NSTM=5;
% L=100d0;
% bb=3d0;
% aa=1d0;
% f_pot = @(x) {heaviside(x-aa)-heaviside(x-bb)};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% If you don't have heaviside defined in your Matlab version define this function as:
% function [y]=heaviside(x)
% y=zeros(1,length(x))
% ipos = find(x>0);
% y(ipos)=1;
%
% 4.Barrier
% NPTS=1000;
% NSTM=5;
% L=100d0;
% bb=3d0;
% aa=1d0;
% pot_dat = [zeros(1,int16(aa/L*NPTS)), ones(1,int16((bb-aa)/L*NPTS)),  zeros(1,int16((L-bb)/L*NPTS))];
% f_pot = @(x) {pot_dat(int16(x/L*(NPTS-1)))};
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
% 5. Harmonic oscillator - with potential defined in file
% Define the potential in file:
% NPTS=1000;
% L=10d0;
% dx=L/(NPTS-1);
% pot_dat=[[0:NPTS-1]*dx; 4d0*([0:NPTS-1]-((NPTS-1)/2)).^2*dx^2]';
% save('pot.dat','-ascii','pot_dat');
% clear
%
% pot_dat=load('pot.dat');
% f_pot = @(x) {spline(pot_dat(:,1),pot_dat(:,2),x)};
% %It is not obligatory the NPTS to be the same as the number of points for the potential
% NPTS=10000;
% NSTM=5;
% L=10d0;
% qm1d_fast(NPTS,NSTM,0,L,f_pot);
%
%
% See also: http://iffwww.iff.kfa-juelich.de/~ekoch/DFT/qm1d.html
function [ee,ev] = qm1d_fast(NPTS,NSTM,a,b,f_pot_handle)
%pot_dat=load(pot_filename);
j=1:NPTS; % indexes for main diagonal
L=b-a;
h=L/(NPTS-1); % space step
x=j*h+a;
V=cell2mat(f_pot_handle(x));

main_diag=2/h^2+V(j);
sub_diag=-1/h^2*ones(1,NPTS-1);

[ee,ev] = trideigs(main_diag,sub_diag,'I',1,NSTM);

% average spacing of energy levels (for adjusting scale of ev)
de0=(ee(NSTM)-ee(1))/(NSTM-1);

% plot potential
plot(x,V(j),'r'); hold on;

% plot eign vectors
de=0.15*sqrt(NPTS)*de0;
for n=1:NSTM
    plot(x,ee(n)+de*ev(:,n)); hold on;
    plot(x,ones(length(j),1)*ee(n)); hold on;
end

% set up plotting options
xlim([0 NPTS]*h+a); ylim([min(V) ee(NSTM)+de0]);

author: Sable

repo: mcbench-benchmarks


text:
function ls_symbfct(P,opts)
% SYMBFCT     - Symbolic factorization for sparse Cholesky factorization.
%            symbfct(P) does a minimum degree ordering and a symbolic
%            factorization for a symmetric square sparse matrix P, and
%            passes the results as global variables.
%            It calls two Fortran MEX files.

% Yin Zhang, January, 1995
% Department of Mathematics and Statistics
% University of Maryland  Baltimore County
% Modified J.Currie AUT May 2013

% PATCH by Jos F. Sturm, May 1998.
% 1. Use global variable "OUTFID" for sending messages. (NOT USED)
% 2. Don't call symfct.mex* if P is diagonal, since that would
%    cause a segmentation violation.

global probData

verb = opts.verb;

% ----- check input matrix P -----
m = size(P,1);
if(size(P,2) ~= m || ~issparse(P) || isempty(P))
   error('Matrix must be square, sparse and nonempty.');
end;
Pdiag = diag(diag(P)); P = P - Pdiag;

% ----- min-degree ordering -----
t0 = tic;
if(verb), fprintf('  Minimum-Degree Ordering ... '); end
[probData.PERM, probData.INVP] = ls_ordmmd(P);
if(verb), fprintf('Done. [%1.3fs]\n',toc(t0)); end

% ----- symbalic factorization -----
if(nnz(P) > 0)
    t0 = tic;
    if(verb), fprintf('  Symbolic Factorization ... '); end
    [probData.XLNZ,probData.NNZL,probData.XSUPER,probData.XLINDX,probData.LINDX,...
     probData.SNODE,probData.SPLIT,probData.TMPSIZ] = ls_symfct(P,probData.PERM,probData.INVP,opts.cache_size);
    if(verb), fprintf('Done. [%1.3fs]\n',toc(t0)); end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% vectorTensorProductArray
% Below is a demonstration of the features of the |vectorTensorProductArray| function

%%
clear; close all; clc;

%% Syntax
% |[AB]=vectorTensorProductArray(a,b);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

% display montage
function MU_dispMontage(Temp,Event,Disp_Matrix,handles)

        % Montage slice selection checking
        if str2num(get(handles.Sli_from,'String'))>str2num(get(handles.Sli_to,'String'))
                errordlg('No slice is chosen for creating montage !');
                return;
        end
        if str2num(get(handles.Sli_from,'String'))<1
                set(handles.Sli_from,'String','1')
                return;
        end
        if str2num(get(handles.Sli_to,'String'))>handles.V.Layer
                set(handles.Sli_to,'String',num2str(handles.V.Layer))
                return;
        end

        % Montage layout initialization
        [Row,Column,Layer]=size(Disp_Matrix);
        ind=0;
        b_flag=0;
        for i=1:str2num(get(handles.Mont_row,'String'))
                for j=1:str2num(get(handles.Mont_col,'String'))
                        DMatrix((i-1)*Row+1:i*Row,(j-1)*Column+1:j*Column)=Disp_Matrix(:,:,str2num(get(handles.Sli_from,'String'))+ind);
                        ind=ind+1;
                        if ind>str2num(get(handles.Sli_to,'String'))-str2num(get(handles.Sli_from,'String'))
                                b_flag=1;
                                break;
                        end
                end
                if b_flag==1
                        break;
                end
        end

        %-------------------------Windows Contrast Change
        Color_map=handles.V.Color_map;
        Contrast_Low=handles.V.C_lower;
        Contrast_High=handles.V.C_upper;
        Contrast_Change=0;
        point=[0 0];
        point2=[0 0];
        Contrast_Interval=(Contrast_High-Contrast_Low)/100;
        %-------------------------End

        delete(handles.Cre_Mont);
        figure ('KeyReleaseFcn',@WindowKeyRelease,'KeyPressFcn',@WindowKeyPress, 'WindowScrollWheelFcn',@MouseScrollWheel,...
                'WindowButtonUpFcn',@MouseUp,'WindowButtonDownFcn', @MouseClick,'Name',['Display ' handles.V.Current_matrix ' Montage']);
        imagesc(DMatrix,[Contrast_Low Contrast_High]);
        colormap(Color_map);
        colorbar;

        function MouseScrollWheel (Temp, Event)
                if Contrast_Change==-1
                    Contrast_Low=Contrast_Low + Event.VerticalScrollCount*Contrast_Interval;
                end
                if Contrast_Change==1
                    Contrast_High=Contrast_High + Event.VerticalScrollCount*Contrast_Interval;
                end

                imagesc (DMatrix,[Contrast_Low Contrast_High]);
                colormap(Color_map);
                colorbar;

        end

        function MouseClick(Temp,Event)
                tpoint=get(gca,'currentpoint');
                point=[round(tpoint(1)),round(tpoint(3))];

        end

        function MouseUp(Temp,Event)
                tpoint2=get(gca,'currentpoint');
                point2=[round(tpoint2(1)),round(tpoint2(3))];

                if sum(abs(point2-point))~=0
                    rectangle('Position',[point(1),point(2),point2(1)-point(1),point2(2)-point(2)],'EdgeColor',[0 1 0]);
                    DispS=DMatrix(point(2):point2(2),point(1):point2(1));
                    DispS=DispS(:);
                    tmean=mean(double(DispS));
                    tstd=std(double(DispS));
                    xlabel(['ROI mean:' num2str(tmean) '  ROI std:' num2str(tstd) '  RSD(%):' num2str(abs(tstd./tmean)*100)],'FontSize',18);
                end
        end

        function WindowKeyPress(Temp,Event)
                %------------------------------Windows Contrast Change
                if Event.Character=='l'
                    Contrast_Change=-1;
                elseif Event.Character=='u'
                    Contrast_Change=1;
                end
                %------------------------------End
        end

        function WindowKeyRelease(Temp,Event)
                Contrast_Change=0;
        end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function test_solve_kron
warning('This function is deprecated')

N1 = 150;
x1 = 1:N1;
K1 = gp_cov_pp(log(10), x1, x1);
I1 = speye(N1);

N2 = 100;
x2 = 1:N2;
K2 = gp_cov_pp(log(10), x2, x2);
I2 = speye(N2);

% Full matrix approach
%K = kron(K1, K2);
%K_noise = kron(K1+I1, K2+I2);
%K_posterior = K - K * (K_noise \ K);

%
% Solve: kron(K1+I1,K2+I2) \ kron(K1,K2)
%

%K = kron(K2,K1);
%X_full = K * (kron(K2+I2,K1+I1) \ K);

% Kronecker approach
tic
C1 = K1 * ((K1 + I1) \ K1);
C2 = K2 * ((K2 + I2) \ K2);
X_kron = diag(C1) * diag(C2)';
%X_kron = dot(K1,C1,1)' * dot(K2,C2,1);
toc

% Conjugate gradient approach
%
% Computing the variances is veeery expensive...
tic
afun = @(x) reshape(kronprod(K1+I1,K2+I2,reshape(x,N1,N2)), N1*N2, 1);
X_cg = zeros(N1,N2);
for n1=1:N1
  for n2=1:N2
    fprintf('Progress: %d %%\n', floor(100*n1/N1));
    x = pcg(afun, kron(K2(:,n2),K1(:,n1)), [], 2);
    X_cg(n1,n2) = K1(n1,:) * reshape(x,N1,N2) * K2(:,n2);
  end
end
toc

full([X_kron(:), X_cg(:)])
%full([X_kron(:), X_cg(:), diag(X_full)])

function Y = kronprod(A, B, X)
% kron(B', A)
Y = A*X*B;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [x, c, funVal, ValueL]=tree_LogisticR(A, y, z, opts)
%
%%
% Function tree_LogisticR
%      Logistic Loss with the
%           tree structured group Lasso Regularization
%
%% Problem
%
%  min  f(x,c) = - sum_i weight_i * log (p_i) + z * sum_j w_j ||x_{G_j}||
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  G_j's are nodes with tree structure
%
%  The tree overlapping group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x( opts.ind(1,j):opts.ind(2,j) ) denotes x_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and n.
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x( opts.G ( opts.ind(1,j):opts.ind(2,j) ) ) denotes x_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind(1:2,:) are within 1 and length(opts.G).
%
% The following example shows how G and ind works:
%
% G={ {1, 2}, {4, 5}, {3, 6}, {7, 8},
%     {1, 2, 3, 6}, {4, 5, 7, 8},
%     {1, 2, 3, 4, 5, 6, 7, 8} }.
%
% ind={ [1, 2, 100]', [3, 4, 100]', [5, 6, 100]', [7, 8, 100]',
%       [9, 12, 100]', [13, 16, 100]', [17, 24, 100]' },
%
% where each node has a weight of 100.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        Regularization parameter (z >=0)
%  opts-      optional inputs (default value: opts=[])
%             !!For tr_LogisticR, we require that opts.ind is specified.!!
%
%% Output parameters:
%  x-         The obtained weight of size n x 1
%  c-         The obtained intercept (scalar)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on October 3, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z<0)
    error('\n z should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to the function 'sll_opts'
%                 for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly.
%     This implicit normalization is suggested for the sparse matrix,
%                                    but not for the dense matrix.
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind
if (~isfield(opts,'ind'))
    error('\n In tree_LeastR, the field .ind should be specified');
else
    ind=opts.ind;

    if (size(ind,1)~=3)
        error('\n Check opts.ind');
    end
end

% The parameter 'weight' contains the weight for each training sample.
% See the definition of the problem above.
% The summation of the weights for all the samples equals to 1.
if (isfield(opts,'sWeight'))
    sWeight=opts.sWeight;

    if ( length(sWeight)~=2 || sWeight(1) <=0 || sWeight(2) <= 0)
        error('\n Check opts.sWeight, which contains two positive values');
    end

    % we process the weight, so that the summation of the weights for all
    % the samples is 1.

    p_flag=(y==1);                  % the indices of the postive samples
    m1=sum(p_flag) * sWeight(1);    % the total weight for the positive samples
    m2=sum(~p_flag) * sWeight(2);   % the total weight for the positive samples

    weight(p_flag,1)=sWeight(1)/(m1+m2);
    weight(~p_flag,1)=sWeight(2)/(m1+m2);
else
    weight=ones(m,1)/m;             % if not specified, we apply equal weight
end

GFlag=0;
% if GFlag=1, we shall apply general_altra
if (isfield(opts,'G'))
    GFlag=1;

    G=opts.G;
    if (max(G) >n || max(G) <1)
        error('\n The input G is incorrect. It should be within %d and %d',1,n);
    end
end

%% Starting point initialization

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(weight(p_flag));         % the total weight for the positive samples
m2=1-m1;                        % the total weight for the positive samples

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
%     if (z<0 || z>1)
%         error('\n opts.rFlag=1, and z should be in [0,1]');
%     end

    computedFlag=0;
    if (isfield(opts,'lambdaMax'))
        if (opts.lambdaMax~=-1)
            lambda=z*opts.lambdaMax;
            computedFlag=1;
        end
    end

    if (~computedFlag)

        % we compute ATb for computing lambda_max, when the input z is a ratio
        b(p_flag,1)=m2;  b(~p_flag,1)=-m1;
        b=b.*weight;

        % compute AT b
        if (opts.nFlag==0)
            ATb =A'*b;
        elseif (opts.nFlag==1)
            ATb= A'*b - sum(b) * mu';  ATb=ATb./nu;
        else
            invNu=b./nu;               ATb=A'*invNu-sum(invNu)*mu';
        end

        % compute lambda_max
        if (GFlag==0)
            lambda_max=findLambdaMax(ATb, n, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax(ATb, n, G, ind, size(ind,2));
        end

        % As .rFlag=1, we set lambda as a ratio of lambda_max
        lambda=z*lambda_max;

    end
end

% The following is for computing lambdaMax
% we use opts.lambdaMax=-1 to show that we need the computation.
%
% One can use this for setting up opts.lambdaMax
if (isfield(opts,'lambdaMax'))

    if (opts.lambdaMax==-1)

        % we compute ATb for computing lambda_max, when the input z is a ratio
        b(p_flag,1)=m2;  b(~p_flag,1)=-m1;
        b=b.*weight;

        % compute AT b
        if (opts.nFlag==0)
            ATb =A'*b;
        elseif (opts.nFlag==1)
            ATb= A'*b - sum(b) * mu';  ATb=ATb./nu;
        else
            invNu=b./nu;               ATb=A'*invNu-sum(invNu)*mu';
        end

        % compute lambda_max
        if (GFlag==0)
            lambda_max=findLambdaMax(ATb, n, ind, size(ind,2));
        else
            lambda_max=general_findLambdaMax(ATb, n, G, ind, size(ind,2));
        end

        x=lambda_max;
        funVal=lambda_max;
        ValueL=lambda_max;

        return;
    end
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1); c=log(m1/m2);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,1);
    end

    if isfield(opts,'c0')
        c=opts.c0;
    else
        c=log(m1/m2);
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);    Ax=Ax./nu;
end

%% The main program

if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    %% The Armijo Goldstein line search schemes + accelearted gradient descent

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb );

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v=s-g/L;  c= sc- gc/L;

            % tree overlapping group Lasso projection
            ind_work(1:2,:)=ind(1:2,:);
            ind_work(3,:)=ind(3,:) * (lambda / L);

            if (GFlag==0)
                x=altra(v, n, ind_work, size(ind_work,2));
            else
                x=general_altra(v, n, G, ind_work, size(ind_work,2));
            end

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );

            r_sum= (v'*v + (c-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x;

        % compute the regularization part
        if (GFlag==0)
            tree_norm=treeNorm(x, n, ind, size(ind,2));
        else
            tree_norm=general_treeNorm(x, n, G, ind, size(ind,2));
        end

        % function value = loss + regularizatioin
        funVal(iterStep)=fun_x + lambda * tree_norm;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
end

author: thomaskuestner


text:
function L = watershedsegment(bw)

D = bwdist(~bw);
D = -D;
D(~bw) = -Inf;
L = watershed(D);

author: kristinbranson

repo: JAABA


text:
% This file standard contains degradation data for the growth of orange trees as used in the article
% An approximation to the log-likelihood function in the Non-Linear Mixed Effects model
% by Pinheiro & Bates (1995,2001)

% Every row contains the degradation or growth observations for one unit/product/object/system
% Replace the data with your own observations

 Y=[30 58 87 115 120 142 145; 33 69 111 156 172 203 203; 30 51 75 108 115 139 140;
          32 62 112 167 179 209 214; 30 49 81 125 142 174 177];

 % Do not change!!!!
 matrixsize_Y=size(Y);
 M=matrixsize_Y(1);
 N=matrixsize_Y(1)*matrixsize_Y(2);


text:
function d = ipdm(data1,varargin)
% ipdm: Inter-Point Distance Matrix
% usage: d = ipdm(data1)
% usage: d = ipdm(data1,data2)
% usage: d = ipdm(data1,prop,value)
% usage: d = ipdm(data1,data2,prop,value)
%
% Arguments: (input)
%  data1 - array of data points, each point is one row. p dimensional
%          data will be represented by matrix with p columns.
%          If only data1 is provided, then the distance matrix
%          is computed between all pairs of rows of data1.
%
%          If your data is one dimensional, it MUST form a column
%          vector. A row vector of length n will be interpreted as
%          an n-dimensional data set.
%
%  data2 - second array, supplied only if distances are to be computed
%          between two sets of points.
%
%
% Class support: data1 and data2 are assumed to be either
% single or double precision. I have not tested this code to
% verify its success on integer data of any class.
%
%
% Additional parameters are expected to be property/value pairs.
% Property/value pairs are pairs of arguments, the first of which
% (properties) must always be a character string. These strings
% may be shortened as long as the shortening is unambiguous.
% Capitalization is ignored. Valid properties for ipdm are:
%
%  'Metric', 'Subset', 'Limit', 'Result'
%
%  'Metric' - numeric flag - defines the distance metric used
%          metric = 2 --> (DEFAULT) Euclidean distance = 2-norm
%                         The standard distance metric.
%
%          metric = 1 --> 1-norm = sum of absolute differences
%                         Also sometimes known as the "city block
%                         metric", since this is the sum of the
%                         differences in each dimension.
%
%          metric = inf --> infinity-norm = maximum difference
%                         over all dimensions. The name refers
%                         to the limit of the p-norm, as p
%                         approaches infinity.
%
%          metric = 0 --> minimum difference over all dimensions.
%                         This is not really a useful norm in
%                         practice.
%
%          Note: while other distance metrics exist, IMHO, these
%          seemed to be the common ones.
%
%
%  'Result' - A string variable that denotes the style of returned
%          result. Valid result types are 'Array', 'Structure'.
%          Capitalization is ignored, and the string may be
%          shortened if you wish.
%
%          result = 'Array' --> (DEFAULT) A matrix of all
%                         interpoint distances will be generated.
%                         This array may be large. If this option
%                         is specified along with a minimum or
%                         maximum value, then those elements above
%                         or below the limiting values will be
%                         set as -inf or +inf, as appropriate.
%
%                         When any of 'LargestFew', 'SmallestFew',
%                         or 'NearestNeighbor' are set, then the
%                         resulting array will be a sparse matrix
%                         if 'array' is specified as the result.
%
%          result = 'Structure' --> A list of all computed distances,
%                         defined as a structure. This structure
%                         will have fields named 'rowindex',
%                         'columnindex', and 'distance'.
%
%                         This option will be useful when a subset
%                         criterion for the distances has been
%                         specified, since then the distance matrix
%                         may be very sparsely populated. Distances
%                         for pairs outside of the criterion will
%                         not be returned.
%
%
%  'Subset' - Character string, any of:
%
%          'All', 'Maximum', 'Minimum', 'LargestFew', 'SmallestFew',
%          'NearestNeighbor', 'FarthestNeighbor', or empty
%
%          Like properties, capitalization is ignored here, and
%          any unambiguous shortening of the word is acceptable.
%
%          DEFAULT = 'All'
%
%          Some interpoint distance matrices can be huge. Often
%          these matrices are too large to be fully retained in
%          memory, yet only the pair of points with the largest
%          or smallest distance may be needed. When only some
%          subset of the complete set of distances is of interest,
%          these options allow you to specify which distances will
%          be returned.
%
%          If 'result' is defined to be an array, then a sparse
%          matrix will be returned for the 'LargestFew', 'SmallestFew',
%          'NearestNeighbor', and 'FarthestNeighbor' subset classes.
%          'Minimum' and 'Maximum' will yield full matrices by
%          default. If a structure is specified, then only those
%          elements which have been identified will be returned.
%
%          Where a subset is specified, its limiting value is
%          specified by the 'Limit' property. Call that value k.
%
%
%          'All' -->     (DEFAULT) Return all interpoint distances
%
%          'Minimum' --> Only look for those distances above
%                        the cutoff k. All other distances will
%                        be returned as -inf.
%
%          'Maximum' --> Only look for those distances below
%                        the cutoff k. All other distances will
%                        be returned as +inf.
%
%          'SmallestFew' --> Only return the subset of the k
%                        smallest distances. Where only one data
%                        set is provided, only the upper triangle
%                        of the inter-point distance matrix will
%                        be generated since that matrix is symmetric.
%
%          'LargestFew' --> Only return the subset of the k
%                        largest distances. Where only one data
%                        set is provided, only the upper triangle
%                        of the inter-point distance matrix will
%                        be generated since that matrix is symmetric.
%
%          'NearestNeighbor' --> Only return the single nearest
%                        neighbor in data2 to each point in data1.
%                        No limiting value is required for this
%                        option. If multiple points have the same
%                        nearest distance, then return the first
%                        such point found. With only one input set,
%                        a point will not be its own nearest
%                        neighbor.
%
%                        Note that exact replicates in a single set
%                        will cause problems, since a sparse matrix
%                        is returned by default. Since they will have
%                        a zero distance, they will not show up in
%                        the sparse matrix. A structure return will
%                        show those points as having a zero distance
%                        though.
%
%          'FarthestNeighbor' --> Only return the single farthest
%                        neighbor to each point. No limiting value
%                        is required for this option. If multiple
%                        points have the same farthest distance,
%                        then return the first such point found.
%
%
%  'Limit' - scalar numeric value or []. Used only when some
%           Subset is specified.
%
%          DEFAULT = []
%
%
%  'ChunkSize' - allows a user with lower RAM limits
%          to force the code to only grab smaller chunks of RAM
%          at a time (where possible). This parameter is specified
%          in bytes of RAM. The default is 32 megabytes, or 2^22
%          elements in any piece of the distance matrix. Only some
%          options will break the problem into chunks, thus as long
%          as a full matrix is expected to be returned, there seems
%          no reason to break the problem up into pieces.
%
%          DEFAULT = 2^25
%
%
% Arguments: (output)
%  d     - array of interpoint distances, or a struct wth the
%          fields {'rowindex', 'columnindex', 'distance'}.
%
%          d(i,j) represents the distance between point i
%          (from data1) and point j (from data2).
%
%          If only one (n1 x p) array is supplied, then d will
%          be an array of size == [n1,n1].
%
%          If two arrays (of sizes n1 x p and n2 x p) then d
%          will be an array of size == [n1,n2].
%
%
% Efficiency considerations:
%  Where possible, this code will use bsxfun to compute its
%  distances.
%
%
% Example:
%  Compute the interpoint distances between all pairs of points
%  in a list of 5 points, in 2 dimensions and using Euclidean
%  distance as the distance metric.
%
%  A = randn(5,2);
%  d = ipdm(A,'metric',2)
%  d =
%            0       2.3295       3.2263       2.0263       2.8244
%       2.3295            0       1.1485      0.31798       1.0086
%       3.2263       1.1485            0       1.4318       1.8479
%       2.0263      0.31798       1.4318            0       1.0716
%       2.8244       1.0086       1.8479       1.0716            0
%
% (see the demo file for many other examples)
%
% See also: pdist
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 2/26/08

% Default property values
params.Metric = 2;
params.Result = 'array';
params.Subset = 'all';
params.Limit = [];
params.ChunkSize = 2^25;

% untangle the arguments
if nargin<1
  % if called with no arguments, then the user probably
  % needs help. Give it to them.
  help ipdm
  return
end

% were two sets of data provided?
pvpairs = {};
if nargin==1
  % only 1 set of data provided
  dataflag = 1;
  data2 = [];
else
  if ischar(varargin{1})
    dataflag = 1;
    data2 = [];
    pvpairs = varargin;
  else
    dataflag = 2;
    data2 = varargin{1};
    if nargin>2
      pvpairs = varargin(2:end);
    end
  end
end

% get data sizes for later
[n1,dim] = size(data1);
if dataflag == 2
  n2 = size(data2,1);
end

% Test the class of the input variables
if ~(isa(data1,'double') || isa(data1,'single')) || ...
    ((dataflag == 2) && ~(isa(data2,'double') || isa(data2,'single')))
  error('data points must be either single or double precision variables.')
end

% do we need to process any property/value pairs?
if nargin>2
  params = parse_pv_pairs(params,pvpairs);

  % check for problems in the properties

  % was a legal Subset provided?
  if ~isempty(params.Subset) && ~ischar(params.Subset)
    error('If provided, ''Subset'' must be character')
  elseif isempty(params.Subset)
    params.Subset = 'all';
  end
  valid = {'all','maximum','minimum','largestfew','smallestfew', ...
    'nearestneighbor','farthestneighbor'};
  ind = find(strncmpi(params.Subset,valid,length(params.Subset)));
  if (length(ind)==1)
    params.Subset = valid{ind};
  else
    error(['Invalid Subset: ',params.Subset])
  end

  % was a limit provided?
  if ~ismember(params.Subset,{'all','nearestneighbor','farthestneighbor'}) && ...
      isempty(params.Limit)
    error('No limit provided, but a Subset that requires a limit value was specified')
  end
  % check the limit values for validity
  if length(params.Limit)>1
    error('Limit must be scalar or empty')
  end

  switch params.Subset
    case {'largestfew', 'smallestfew'}
      % must be at least 1, and an integer
      if (params.Limit<1) || (round(params.Limit)~=params.Limit)
        error('Limit must be a positive integer for LargestFew or NearestFew')
      end
  end

  % was a legal Result provided?
  if isempty(params.Result)
    params.result = 'Array';
  elseif ~ischar(params.Result)
    error('If provided, ''Result'' must be character or empty')
  end
  valid = {'array','structure'};
  ind = find(strncmpi(params.Result,valid,length(params.Result)));
  if (length(ind)==1)
    params.Result = valid{ind};
  else
    error(['Invalid Result: ',params.Subset])
  end

  % check for the metric
  if isempty(params.Metric)
    params.Metric = 2;
  elseif (length(params.Metric)~=1) || ~ismember(params.Metric,[0 1 2 inf])
    error('If supplied, ''Metric'' must be a scalar, and one of [0 1 2 inf]')
  end
end % if nargin>2

% If Metric was given as 2, but the dimension is only 1, then it will
% be slightly faster (and equivalent) to use the 1-norm Metric.
if (dim == 1) && (params.Metric == 2)
  params.Metric = 1;
end

% Can we use bsxfun to compute the interpoint distances?
% Older Matlab releases will not have bsxfun, but if it is
% around, it will ne both faster and less of a memory hog.
params.usebsxfun = (5==exist('bsxfun','builtin'));

% check for dimension mismatch if 2 sets
if (dataflag==2) && (size(data2,2)~=dim)
  error('If 2 point sets provided, then both must have the same number of columns')
end

% Total number of distances to compute, in case I must do it in batches
if dataflag==1
  n2 = n1;
end
ntotal = n1*n2;

% FINALLY!!! Compute inter-point distances
switch params.Subset
  case 'all'
    % The complete set of interpoint distances. There is no need
    % to break this into chunks, since we must return all distances.
    % If that is too much to compute in memory, then it will fail
    % anyway when we try to store the result. bsxfun will at least
    % do the computation efficiently.

    % One set or two?
    if dataflag == 1
      d = distcomp(data1,data1,params);
    else
        d = distcomp(data1,data2,params);
    end

    % Must we return it as a struct?
    if params.Result(1) == 's'
      [rind,cind] = ndgrid(1:size(d,1),1:size(d,2));
      ds.rowindex = rind(:);
      ds.columnindex = cind(:);
      ds.distance = d(:);
      d = ds;
    end

  case {'minimum' 'maximum'}
    % There is no reason to break this into pieces if the result
    % sill be filled in the end with +/- inf. Only break it up
    % if the final result is a struct.
    if ((ntotal*8)<=params.ChunkSize) || (params.Result(1) == 'a')
      % its small enough to do it all at once

      % One set or two?
      if dataflag == 1
        d = distcomp(data1,data1,params);
      else
          d = distcomp(data1,data2,params);
      end

      % Must we return it as a struct?
      if params.Result(1) == 'a'
        % its an array, fill the unwanted distances with +/- inf
        if params.Subset(2) == 'i'
          % minimum
          d(d<=params.Limit) = -inf;
        else
          % maximum
          d(d>=params.Limit) = +inf;
        end
      else
        % a struct will be returned
        if params.Subset(2) == 'i'
          % minimum
          [dist.rowindex,dist.columnindex] = find(d>=params.Limit);
        else
          % maximum
          [dist.rowindex,dist.columnindex] = find(d<=params.Limit);
        end
        dist.distance = d(dist.rowindex + n1*(dist.columnindex-1));
        d = dist;
      end

    else
      % we need to break this into chunks. This branch
      % will always return a struct.

      % this is the number of rows of data1 that we will
      % process at a time.
      bs = floor(params.ChunkSize/(8*n2));
      bs = min(n1,max(1,bs));

      % Accumulate the result into a cell array. Do it this
      % way because we don't know in advance how many elements
      % that we will find satisfying the minimum or maximum
      % limit specified.
      accum = cell(0,1);

      % now loop over the chunks
      batch = 1:bs;
      while ~isempty(batch)

        % One set or two?
        if dataflag == 1
          dist = distcomp(data1(batch,:),data1,params);
        else
          dist = distcomp(data1(batch,:),data2,params);
        end

        % big or small as requested
        if ('i'==params.Subset(2))
          % minimum value specified
          [I,J,V] = find(dist>=params.Limit);
        else
          % maximum limit
          [I,J] = find(dist<=params.Limit);
          I = I(:);
          J = J(:);
          V = dist(I + (J-1)*length(batch));
          I = I + (batch(1)-1);
        end

        % and stuff them into the cell structure
        if ~isempty(V)
          accum{end+1,1} = [I,J,V(:)]; %#ok
        end

        % increment the batch
        batch = batch + bs;
        if batch(end)>n1
          batch(batch>n1) = [];
        end

      end

      % convert the cells into one flat array
      accum = cell2mat(accum);

      if isempty(accum)
        d.rowindex = [];
        d.columnindex = [];
        d.distance = [];
      else
        % we found something

        % sort on the second column, to put them in a reasonable order
        accum = sortrows(accum,[2 1]);

        d.rowindex = accum(:,1);
        d.columnindex = accum(:,2);
        d.distance = accum(:,3);
      end

    end

  case {'smallestfew' 'largestfew'}
    % find the k smallest/largest distances. k is
    % given by params.Limit

    % if only 1 set, params.Limit must be less than n*(n-1)/2
    if dataflag == 1
      params.Limit = min(params.Limit,n1*(n1-1)/2);
    end

    % is this a large problem?
    if ((ntotal*8) <= params.ChunkSize)
      % small potatoes

      % One set or two?
      if dataflag == 1
        dist = distcomp(data1,data1,params);
        % if only one data set, set the diagonal and
        % below that to +/- inf so we don't find it.
        temp = find(tril(ones(n1,n1),0));
        if params.Subset(1) == 's'
          dist(temp) = inf;
        else
          dist(temp) = -inf;
        end
      else
        dist = distcomp(data1,data2,params);
      end

      % sort the distances to find those we need
      if ('s'==params.Subset(1))
        % smallestfew
        [val,tags] = sort(dist(:),'ascend');
      else
        % largestfew
        [val,tags] = sort(dist(:),'descend');
      end
      val = val(1:params.Limit);
      tags = tags(1:params.Limit);

      % recover the row and column index from the linear
      % index returned by sort in tags.
      [d.rowindex,d.columnindex] = ind2sub([n1,size(dist,2)],tags);

      % create the matrix as a sparse one or a struct?
      if params.Result(1)=='a'
        % its an array, so make the array sparse.
        d = sparse(d.rowindex,d.columnindex,val,n1,size(dist,2));
      else
        % a structure
        d.distance = val;
      end

    else
      % chunks

      % this is the number of rows of data1 that we will
      % process at a time.
      bs = floor(params.ChunkSize/(8*n2));
      bs = min(n1,max(1,bs));

      % We need to find the extreme cases. There are two possible
      % algorithms, depending on how many total elements we will
      % search for.
      % 1. Only a very few total elements.
      % 2. A relatively large number of total elements, forming
      %    a significant fraction of the total set.
      %
      % Case #1 would suggest to retain params.Limit numberr of
      % elements from each batch, then at the end, sort them all
      % to find the best few. Case #2 will result in too many
      % elements to retain, so we must distinguish between these
      % alternatives.
      if (8*params.Limit*n1/bs) <= params.ChunkSize
        % params.Limit is small enough to fall into case #1.

        % Accumulate the result into a cell array. Do it this
        % way because we don't know in advance how many elements
        % that we will find satisfying the minimum or maximum
        % limit specified.
        accum = cell(0,1);

        % now loop over the chunks
        batch = (1:bs)';
        while ~isempty(batch)
          % One set or two?
          if dataflag == 1
            dist = distcomp(data1(batch,:),data1,params);
            k = find(tril(ones(length(batch),n2),batch(1)-1));
            if ('s'==params.Subset(1))
              dist(k) = inf;
            else
              dist(k) = -inf;
            end
          else
            dist = distcomp(data1(batch,:),data2,params);
          end

          % big or small as requested, keeping only the best
          % params.Limit number of elements
          if ('s'==params.Subset(1))
            % minimum value specified
            [tags,tags] = sort(dist(:),1,'ascend');
            tags = tags(1:bs);
            [I,J] = ndgrid(batch,1:n2);
            ijv = [I(tags),J(tags),dist(tags)];
          else
            % maximum limit
            [tags,tags] = sort(dist(:),1,'descend');
            tags = tags(1:bs);
            [I,J] = ndgrid(batch,1:n2);
            ijv = [I(tags),J(tags),dist(tags)];
          end
          % and stuff them into the cell structure
          accum{end+1,1} = ijv; %#ok

          % increment the batch
          batch = batch + bs;
          if batch(end)>n1
            batch(batch>n1) = [];
          end
        end

        % convert the cells into one flat array
        accum = cell2mat(accum);

        % keep only the params.Limit best of those singled out
        accum = sortrows(accum,3);
        if ('s'==params.Subset(1))
          % minimum value specified
          accum = accum(1:params.Limit,:);
        else
          % minimum value specified
          accum = accum(end + 1 - (1:params.Limit),:);
        end
        d.rowindex = accum(:,1);
        d.columnindex = accum(:,2);
        d.distance = accum(:,3);

        % create the matrix as a sparse one or a struct?
        if params.Result(1)=='a'
          % its an array, so make the array sparse.
          d = sparse(d.rowindex,d.columnindex,d.distance,n1,size(dist,2));
        end

      else
        % params.Limit forces us into the domain of case #2.
        % Here we cannot retain params.Limit elements from each chunk.
        % so we will grab each chunk and append it to the best elements
        % found so far, then filter out the best after each chunk is
        % done. This may be slower than we want, but its the only way.
        ijv = zeros(0,3);

        % loop over the chunks
        batch = (1:bs)';
        while ~isempty(batch)
          % One set or two?
          if dataflag == 1
            dist = distcomp(data1(batch,:),data1,params);
            k = find(tril(ones(length(batch),n2),batch(1)-1));
            if ('s'==params.Subset(1))
              dist(k) = inf;
            else
              dist(k) = -inf;
            end
          else
            dist = distcomp(data1(batch,:),data2,params);
          end

          [I,J] = ndgrid(batch,1:n2);
          ijv = [ijv;[I(:),J(:),dist(:)]]; %#ok

          % big or small as requested, keeping only the best
          % params.Limit number of elements
          if size(ijv,1) > params.Limit
            if ('s'==params.Subset(1))
              % minimum value specified
              [tags,tags] = sort(ijv(:,3),1,'ascend');
            else
              [tags,tags] = sort(ijv(:,3),1,'ascend');
            end
            ijv = ijv(tags(1:params.Limit),:);
          end

          % increment the batch
          batch = batch + bs;
          if batch(end)>n1
            batch(batch>n1) = [];
          end
        end

        % They are fully trimmed down. stuff a structure
        d.rowindex = ijv(:,1);
        d.columnindex = ijv(:,2);
        d.distance = ijv(:,3);

        % create the matrix as a sparse one or a struct?
        if params.Result(1)=='a'
          % its an array, so make the array sparse.
          d = sparse(d.rowindex,d.columnindex,d.distance,n1,size(dist,2));
        end

      end

    end

  case {'nearestneighbor' 'farthestneighbor'}
    % find the closest/farthest neighbor for every point

    % is this a large problem? Or a 1-d problem?
    if dim == 1
      % its a 1-d nearest/farthest neighbor problem. we can
      % special case these easily enough, and all the distance
      % metric options are the same in 1-d.

      % first split it into the farthest versus nearest cases.
      if params.Subset(1) == 'f'
        % farthest away

        % One set or two?
        if dataflag == 1
          [d2min,minind] = min(data1);
          [d2max,maxind] = max(data1);
        else
          [d2min,minind] = min(data2);
          [d2max,maxind] = max(data2);
        end

        d.rowindex = (1:n1)';
        d.columnindex = repmat(maxind,n1,1);
        d.distance = repmat(d2max,n1,1);

        % which endpoint was further away?
        k = abs((data1 - d2min)) >= abs((data1 - d2max));
        if any(k)
          d.columnindex(k) = minind;
          d.distance(k) = d2min;
        end

      else
        % nearest. this is mainly a sort and some fussing around.
        d.rowindex = (1:n1)';
        d.columnindex = ones(n1,1);
        d.distance = zeros(n1,1);

        % One set or two?
        if dataflag == 1
          % if only one data point, then we are done
          if n1 == 2
            % if exactly two data points, its trivial
            d.columnindex = [2 1];
            d.distance = repmat(abs(diff(data1)),2,1);
          elseif n1>2
            % at least three points. do a sort.
            [sorted_data,tags] = sort(data1);

            % handle the first and last points separately
            d.columnindex(tags(1)) = tags(2);
            d.distance(tags(1)) = sorted_data(2) - sorted_data(1);
            d.columnindex(tags(end)) = tags(end-1);
            d.distance(tags(end)) = sorted_data(end) - sorted_data(end-1);

            ind = (2:(n1-1))';

            d1 = sorted_data(ind) - sorted_data(ind-1);
            d2 = sorted_data(ind+1) - sorted_data(ind);

            k = d1 < d2;
            d.distance(tags(ind(k))) = d1(k);
            d.columnindex(tags(ind(k))) = tags(ind(k)-1);
            k = ~k;
            d.distance(tags(ind(k))) = d2(k);
            d.columnindex(tags(ind(k))) = tags(ind(k)+1);
          end % if n1 == 2
        else
          % Two sets of data. still really a sort and some fuss.
          if n2 == 1
            % there is only one point in data2
            d.distance = abs(data1 - data2);
            % d.columnindex is already set correctly
          else
            % At least two points in data2
            % We need to sort all the data points together, but also
            % know which points from each set went where. ind12 and
            % bool12 will help keep track.
            ind12 = [1:n1,1:n2]';
            bool12 = [zeros(n1,1);ones(n2,1)];
            [sorted_data,tags] = sort([data1;data2]);

            ind12 = ind12(tags);
            bool12 = bool12(tags);

            % where did each point end up after the sort?
            loc1 = find(~bool12);
            loc2 = find(bool12);

            % for each point in data1, what is the (sorted) data2
            % element which appears most nearly to the left of it?
            cs = cumsum(bool12);
            leftelement = cs(loc1);

            % any points which fell below the minimum element in data2
            % will have a zero for the index of the element on their
            % left. fix this.
            leftelement = max(1,leftelement);

            % likewise, any point greater than the max in data2 will
            % have an n2 in left element. this too will be a problem
            % later, so fix it.
            leftelement = min(n2-1,leftelement);

            % distance to the left hand element
            dleft = abs(sorted_data(loc1) - sorted_data(loc2(leftelement)));
            dright = abs(sorted_data(loc1) - sorted_data(loc2(leftelement+1)));

            % find the points which are closer to the left element in data2
            k = (dleft < dright);
            d.distance(ind12(loc1(k))) = dleft(k);
            d.columnindex(ind12(loc1(k))) = ind12(loc2(leftelement(k)));
            k = ~k;
            d.distance(ind12(loc1(k))) = dright(k);
            d.columnindex(ind12(loc1(k))) = ind12(loc2(leftelement(k)+1));

          end % if n2 == 1
        end % if dataflag == 1
      end % if params.Subset(1) == 'f'

      % create the matrix as a sparse one or a struct?
      if params.Result(1)=='a'
        % its an array, so make the array sparse.
        d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
      end

    elseif (ntotal>1000) && (((params.Metric == 0) && (params.Subset(1) == 'n')) || ...
        ((params.Metric == inf) && (params.Subset(1) == 'f')))
      % nearest/farthest neighbour in n>1 dimensions, but for an
      % infinity norm metric. Reduce this to a sequence of
      % 1-d problems, each of which will be faster in general.
      % do this only if the problem is moderately large, since
      % we must overcome the extra overhead of the recursive
      % calls to ipdm.

      % do the first dimension
      if dataflag == 1
        d = ipdm(data1(:,1),data1(:,1),'subset',params.Subset,'metric',params.Metric,'result','struct');
      else
        d = ipdm(data1(:,1),data2(:,1),'subset',params.Subset,'metric',params.Metric,'result','struct');
      end

      % its slightly different for nearest versus farthest here
      % now, loop over dimensions
      for i = 2:dim
        if dataflag == 1
          di = ipdm(data1(:,i),data1(:,i),'subset',params.Subset,'metric',params.Metric,'result','struct');
        else
          di = ipdm(data1(:,i),data2(:,i),'subset',params.Subset,'metric',params.Metric,'result','struct');
        end

        % did any of the distances change?
        if params.Metric == 0
          % the 0 norm, with nearest neighbour, so take the
          % smallest distance in any dimension.
          k = d.distance > di.distance;
        else
          % inf norm. so take the largest distance across dimensions
          k = d.distance < di.distance;
        end

        if any(k)
          d.distance(k) = di.distance(k);
          d.columnindex(k) = di.columnindex(k);
        end
      end

      % create the matrix as a sparse one or a struct?
      if params.Result(1)=='a'
        % its an array, so make the array sparse.
        d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
      end

    elseif ((ntotal*8) <= params.ChunkSize)
      % None of the other special cases apply, so do it using brute
      % force for the small potatoes problem.

      % One set or two?
      if dataflag == 1
        dist = distcomp(data1,data1,params);
      else
        dist = distcomp(data1,data2,params);
      end

      % if only one data set and if a nearest neighbor
      % problem, set the diagonal to +inf so we don't find it.
      if (dataflag==1) && (n1>1) && ('n'==params.Subset(1))
        diagind = (1:n1) + (0:n1:(n1^2-1));
        dist(diagind) = +inf;
      end

      if ('n'==params.Subset(1))
        % nearest
        [val,j] = min(dist,[],2);
      else
        % farthest
        [val,j] = max(dist,[],2);
      end

      % create the matrix as a sparse one or a struct?
      if params.Result(1)=='a'
        % its an array, so make the array sparse.
        d = sparse((1:n1)',j,val,n1,size(dist,2));
      else
        % a structure
        d.rowindex = (1:n1)';
        d.columnindex = j;
        d.distance = val;
      end

    else

      % break it into chunks
      bs = floor(params.ChunkSize/(8*n2));
      bs = min(n1,max(1,bs));

      % pre-allocate the result
      d.rowindex = (1:n1)';
      d.columnindex = zeros(n1,1);
      d.distance = zeros(n1,1);

      % now loop over the chunks
      batch = 1:bs;
      while ~isempty(batch)

        % One set or two?
        if dataflag == 1
          dist = distcomp(data1(batch,:),data1,params);
        else
          dist = distcomp(data1(batch,:),data2,params);
        end

        % if only one data set and if a nearest neighbor
        % problem, set the diagonal to +inf so we don't find it.
        if (dataflag==1) && (n1>1) && ('n'==params.Subset(1))
          diagind = 1:length(batch);
          diagind = diagind + (diagind-2+batch(1))*length(batch);
          dist(diagind) = +inf;
        end

        % big or small as requested
        if ('n'==params.Subset(1))
          % nearest
          [val,j] = min(dist,[],2);
        else
          % farthest
          [val,j] = max(dist,[],2);
        end

        % and stuff them into the result structure
        d.columnindex(batch) = j;
        d.distance(batch) = val;

        % increment the batch
        batch = batch + bs;
        if batch(end)>n1
          batch(batch>n1) = [];
        end

      end

      % did we need to return a struct or an array?
      if params.Result(1) == 'a'
        % an array. make it a sparse one
        d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
      end

    end % if dim == 1

end  % switch params.Subset

% End of mainline

% ======================================================
% begin subfunctions
% ======================================================
function d = distcomp(set1,set2,params)
% Subfunction to compute all distances between two sets of points
dim = size(set1,2);
% can we take advantage of bsxfun?
% Note: in theory, there is no need to loop over the dimensions. We
% could Just let bsxfun do ALL the work, then wrap a sum around the
% outside. In practice, this tends to create large intermediate
% arrays, especially in higher numbers of dimensions. Its also when
% we might gain here by use of a vectorized code. This will only be
% a serious gain when the number of points is relatively small and
% the dimension is large.
if params.usebsxfun
  % its a recent enough version of matlab that we can
  % use bsxfun at all.
  n1 = size(set1,1);
  n2 = size(set2,1);
  if (dim>1) && ((n1*n2*dim)<=params.ChunkSize)
    % its a small enough problem that we might gain by full
    % use of bsxfun
    switch params.Metric
      case 2
        d = sum(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim])).^2,3);
      case 1
        d = sum(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),3);
      case inf
        d = max(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),[],3);
      case 0
        d = min(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),[],3);
    end
  else
    % too big, so that the ChunkSize will have been exceeded, or just 1-d
    if params.Metric == 2
      d = bsxfun(@minus,set1(:,1),set2(:,1)').^2;
    else
      d = abs(bsxfun(@minus,set1(:,1),set2(:,1)'));
    end
    for i=2:dim
      switch params.Metric
        case 2
          d = d + bsxfun(@minus,set1(:,i),set2(:,i)').^2;
        case 1
          d = d + abs(bsxfun(@minus,set1(:,i),set2(:,i)'));
        case inf
          d = max(d,abs(bsxfun(@minus,set1(:,i),set2(:,i)')));
        case 0
          d = min(d,abs(bsxfun(@minus,set1(:,i),set2(:,i)')));
      end
    end
  end
else
  % Cannot use bsxfun. Sigh. Do things the hard (and slower) way.
  n1 = size(set1,1);
  n2 = size(set2,1);
  if params.Metric == 2
    % Note: While some people might use a different Euclidean
    % norm computation based on expanding the square of the
    % difference of two numbers, that computation is inherantly
    % inaccurate when implemented in floating point arithmetic.
    % While it might be faster, I won't use it here. Sorry.
    d = (repmat(set1(:,1),1,n2) - repmat(set2(:,1)',n1,1)).^2;
  else
    d = abs(repmat(set1(:,1),1,n2) - repmat(set2(:,1)',n1,1));
  end
  for i=2:dim
    switch params.Metric
      case 2
        d = d + (repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)).^2;
      case 1
        d = d + abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1));
      case inf
        d = max(d,abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)));
      case 0
        d = min(d,abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)));
    end
  end
end
% if 2 norm, then we must sqrt at the end
if params.Metric==2
  d = sqrt(d);
end

% ==============================================================
%    end main ipdm
%    begin included function - parse_pv_pairs
% ==============================================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params.
  % Use setfield for comptability issues with older releases.
  params = setfield(params,p_i,v_i); %#ok

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function spline_test205 ( )

%*****************************************************************************80
%
%% TEST205 tests SPLINE_LINEAR_INT, SPLINE_LINEAR_INTSET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST205\n' );
  fprintf ( 1, '  SPLINE_LINEAR_INTSET is given some interval endpoints,\n' );
  fprintf ( 1, '  and a value associated with each interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  It determines a linear spline, with breakpoints\n' );
  fprintf ( 1, '  at the centers of each interval, whose integral\n' );
  fprintf ( 1, '  over each interval is equal to the given value.\n' );

  int_x(1:n+1) = [ 0.0, 1.0, 4.0, 5.0, 10.0 ]';
  int_v(1:n) = [ 10.0, 2.0, 8.0, 27.5 ]';

  r8vec_print ( n+1, int_x, '  The interval end points:' );
  r8vec_print ( n, int_v, '  The desired interval integral values:' );

  [ data_x, data_y ] = spline_linear_intset ( n, int_x, int_v );

  r8vec_print ( n, data_x, '  The spline break points:' );
  r8vec_print ( n, data_y, '  The spline data values: ' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  As a check, call SPLINE_LINEAR_INT to compute\n' );
  fprintf ( 1, '  the integral of the spline over each interval,\n' );
  fprintf ( 1, '  and compare to the desired value.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     A       B    Desired      Computed\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    a = int_x(i);
    b = int_x(i+1);
    value = spline_linear_int ( n, data_x, data_y, a, b );
    fprintf ( 1, '%6f  %6f  %12f  %12f\n', a, b, int_v(i), value );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p_ = ViewCurveSlope(X,p)

[J,K]=size(X);

% constrain probabilities to sum to one...
Aeq = ones(1,J);
beq=1;

% ...constrain the expectation...
V= X(:,14)-X(:,13);
v=.0005;

Aeq=[Aeq
    V'];
beq=[beq
    v];

A=[];
b=[];

% ...compute posterior probabilities
p_ = EntropyProg(p,A,b,Aeq ,beq);


text:
MATLAB_PATH = '/usr/local/matlab/';

global Start_path;
Start_path = '/home.dokt/martid1/export/svoboda/Matlab';
setpaths;

screen = get(0, 'ScreenSize');
xwidth = screen(4)/2;
ywidth = screen(4)/2;
set(0, 'defaultfigureposition', [screen(3)-xwidth*1.1 5 xwidth*1.1-4 ywidth*0.9]);
clear screen xwidth ywidth

set(0,'FormatSpacing',         'compact',...
      'DefaultFigureColormap' , (0:1/63:1)'*[1 1 1],...
      'DefaultFigureMenu', 'none' );
close

flops(0);

author: strawlab

repo: MultiCamSelfCal


text:
clear;clc;
util;
VAEdata = load([input_folder, filesep, data_filename]);% load the generated data

tic
dataset = {};
for i = 1:length(VAEdata.boxes)
    data.leafreps = VAEdata.boxes{i}';
    data.relposreps = VAEdata.rplist{i}';
    data.kids = VAEdata.kids{i};
    [ ndata ] = convertPytorchData( data );
    ndata = adjustReconOffset( ndata );%adjust the one-hot vectors
    dataset{length(dataset)+1} = ndata;
end

rdataset = {};
for i = 1:length(dataset)
    fprintf([num2str(i),'\n']);
    p_data = dataset{i};

    % adjust the relpos between floor and wallnodes
    rplist = p_data.relposreps{length(p_data.relposreps)};
    rplist(1,:) = [-1,0.5,0.5,0.5];
    rplist(24:end,1:4) = 0;
    rplist(25,1) = 1;
    rplist(28,2:end) = 1;
    p_data.relposreps{length(p_data.relposreps)} = rplist;
	try
        [ rdata ] = reconOffset2Scene( p_data, 7, 6 , i, labelset_filename, objperlabel_filename, objfolder_path);
        [ rdata ] = optimize3dScene( rdata ); % generate 3D scenes from 2D

        % plot and save the topviews
        plotRoom2(rdata.obblist, rdata.labellist, 0);
        set(gcf,'position',[100 100 900 900])
        saveas(gcf, [output_img_folder, filesep,num2str(i),'.jpg']);
        close gcf;
catch
end
    % save the object meshes in the scenes
    % scene_folder = [output_obj_folder, filesep, 'scene-', num2str(i)];
    % mkdir(scene_folder);
    % recon3dscene( rdata, scene_folder, 7, objfolder_path);
    rdataset{i} = rdata;
end
toc

repo: GRAINS


text:
function lambda = ortega_eigenvalues ( n, u, v, d )

%*****************************************************************************80
%
%% ORTEGA_EIGENVALUES returns the eigenvalues of the ORTEGA matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Ortega,
%    Generation of Test Matrices by Similarity Transformations,
%    Communications of the ACM,
%    Volume 7, 1964, pages 377-378.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    2 <= N.
%
%    Input, real U(N,1), V(N,1), vectors which define the matrix.
%    U'V must not equal -1.0.  If, in fact, U'V = 0, and U, V and D are
%    integers, then the matrix, inverse, eigenvalues, and eigenvectors
%    will be integers.
%
%    Input, real D(N,1), the desired eigenvalues.
%
%    Output, real LAMBDA(N,1), the determinant.
%
  lambda = zeros ( n, 1 );

  lambda(1:n,1) = d(:);

  return
end


text:
function z = mrdivide( x, y )

%   Disciplined convex/geomtric programming information for MRDIVIDE:
%      The MRDIVIDE operation X/Y is quite often employed with Y as a
%      scalar. In that case, it is equivalent to the RDIVIDE operation
%      X./Y, and must obey the same rules as outlined in the help for
%      CVX/RDIVIDE.
%
%      When Y is a matrix, the MRDIVIDE operation X/Y is equivalent to
%      X*inv(Y) for both DCP and DGP purposes. The inv() operation is
%      not supported for non-constant expressions, so Y must be both
%      constant and nonsingular. The resulting matrix multiplication
%      must obey the same rules as outlined in the help for CVX/MTIMES.

sz = size( y );
if all( sz == 1 ),
    z = rdivide( x, y, '/' );
elseif length( sz ) > 2,
    cvx_throw( 'Inputs must be 2-D, or at least one input must be scalar.' );
elseif sz( 1 ) ~= sz( 2 ) && length( sz ) == 2,
    cvx_throw( 'Non-square matrix divisors are not supported in CVX.' );
elseif ~cvx_isconstant( y ),
    cvx_throw( 'Matrix divisors must be constant.' );
else
    z = mtimes( x, cvx_constant( y ), '/' );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% SYNTAX:
% yAvgStd = hmrS_SessAvgStd2_oldversion_nirs(yAvgRuns, ySum2Runs, mlActRuns, nTrialsRuns)
%
% UI NAME:
% Run_Average_Standard_Deviation
%
% DESCRIPTION:
% Calculate avearge HRF standard deviation of all runs for one subject.
%
% INPUTS:
% yAvgRuns:
% ySum2Runs:
% mlActRuns:
% nTrialsRuns:
%
% OUTPUTS:
% yAvgStdOut: the standard deviation across runs
%
% USAGE OPTIONS:
% Run_Average_Standard_Deviation_on_Concentration_Data:  dcAvgStd  = hmrS_SessAvgStd2_oldversion_nirs(dcAvgRuns, dcSum2Runs, mlActRuns, nTrialsRuns)
% Run_Average_Standard_Deviation_on_Delta_OD_Data:       dodAvgStd = hmrS_SessAvgStd2_oldversion_nirs(dodAvgRuns, dodSum2Runs, mlActRuns, nTrialsRuns)
%

function yAvgStdOut = hmrS_SessAvgStd2_oldversion_nirs(yAvgRuns, ySum2Runs, mlActRuns, nTrialsRuns)

yAvgStdOut = DataClass().empty();

nDataBlks = length(yAvgRuns{1});
nTrials_tot = cell(nDataBlks,1);

for iBlk = 1:nDataBlks

    grp1 = [];

    for iRun = 1:length(yAvgRuns)

        yAvgStdOut(iBlk) = DataClass();

        yAvg    = yAvgRuns{iRun}(iBlk).GetDataTimeSeries('reshape');
        if isempty(yAvg)
            continue;
        end

        yAvgStd = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), size(yAvg,4));
        ySum2   = ySum2Runs{iRun}(iBlk).GetDataTimeSeries('reshape');
        tHRF    = yAvgRuns{iRun}(iBlk).GetTime();
        nTrials = nTrialsRuns{iRun}{iBlk};
        if isempty(mlActRuns{iRun})
            mlActRuns{iRun} = cell(length(nDataBlks),1);
        end

        %
        datatype  = yAvgRuns{iRun}(iBlk).GetDataTypeLabel();
        if strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))
            ml    = yAvgRuns{iRun}(iBlk).GetMeasListSrcDetPairs('reshape');
        elseif strcmp(datatype{1}, 'HRF dOD')
            ml    = yAvgRuns{iRun}(iBlk).GetMeasList('reshape');
        end
        if isempty(mlActRuns{iRun}{iBlk})
            mlActRuns{iRun}{iBlk} = ones(size(ml,1),1);
        end
        mlAct = mlActRuns{iRun}{iBlk}(1:size(ml,1));

        nCond = size(nTrials,2);
        yAvgStdOut(iBlk).SetTime(tHRF);

        if strcmp(datatype{1}, 'HRF dOD')

            if isempty(grp1)
                grp1 = zeros(size(yAvg,1), size(yAvg,2), nCond);
                grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), nCond);
                nTrials_tot{iBlk} = zeros(size(yAvg,2), nCond);
            end

            lstChInc = find(mlAct==1);
            for iC = 1:nCond
                nT = nTrials(iC);
                if nT>0
                    if iRun==1
                        grp1(:,lstChInc,iC) = yAvg(:,lstChInc,iC) * nT;
                        grp1Sum2(:,lstChInc,iC) = ySum2(:,lstChInc,iC);
                        nTrials_tot{iBlk}(lstChInc,iC) = nT;
                    else
                        for iCh=1:length(lstChInc) %size(yAvg,2)
                            % Make sure 3rd arg to interp1 is column vector to guarauntee interp1 output is column vector
                            % which matches grp1 dimensions when adding the two.
                            grp1(:,lstChInc(iCh),iC) = grp1(:,lstChInc(iCh),iC) + interp1(tHRF,yAvg(:,lstChInc(iCh),iC),tHRF(:)) * nT;
                            grp1Sum2(:,lstChInc(iCh),iC) = grp1Sum2(:,lstChInc(iCh),iC) + interp1(tHRF,ySum2(:,lstChInc(iCh),iC),tHRF(:));
                            nTrials_tot{iBlk}(lstChInc(iCh),iC) = nTrials_tot{iBlk}(lstChInc(iCh),iC) + nT;
                        end
                    end
                end
            end

            if ~isempty(grp1)
                for iC = 1:size(grp1,3)
                    for iCh = 1:size(grp1,2)
                        yAvg(:,iCh,iC) = grp1(:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC);

                        % We want to distinguish between no trials and 1 trial:
                        % If there are no trials, we have no HRF data and no std which
                        % the first case will calculate as opposed to one trial (2nd case)
                        % where we have all zeros.
                        if(nTrials_tot{iBlk}(iCh,iC)~=1)
                            yAvgStd(:,iCh,iC) = ( (1/(nTrials_tot{iBlk}(iCh,iC)-1))*grp1Sum2(:,iCh,iC) - (nTrials_tot{iBlk}(iCh,iC)/(nTrials_tot{iBlk}(iCh,iC)-1))*(grp1(:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC)).^2).^0.5 ;
                        else
                            yAvgStd(:,iCh,iC) = zeros(size(grp1Sum2(:,iCh,iC)));
                        end

                        %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                        %%%% Now we can set channel descriptors for avg and standard deviation
                        if iRun == length(yAvgRuns)
                            yAvgStdOut(iBlk).AddChannelDod(ml(iCh,1), ml(iCh,2), ml(iCh,4), iC);
                        end
                    end

                    %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                    %%%% Now we can set channel descriptors for avg and standard deviation
                    if iRun == length(yAvgRuns)
                        yAvgStdOut(iBlk).AppendDataTimeSeries(yAvgStd(:,:,iC));
                    end
                end
            end

        elseif strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))

            if isempty(grp1)
                grp1 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
                grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
                nTrials_tot{iBlk} = zeros(size(yAvg,3), nCond);
            end

            lstChInc = find(mlAct==1);
            for iC = 1:1:nCond
                nT = nTrials(iC);
                if nT>0
                    if iRun==1
                        grp1(:,:,lstChInc,iC) = yAvg(:,:,lstChInc,iC) * nT;
                        grp1Sum2(:,:,lstChInc,iC) = ySum2(:,:,lstChInc,iC);
                        nTrials_tot{iBlk}(lstChInc,iC) = nT;
                    else
                        for iCh=1:length(lstChInc) %size(yAvg,3)
                            for iHb=1:size(yAvg,2)
                                % Make sure 3rd arg to interp1 is column vector to guarauntee interp1 output is column vector
                                % which matches grp1 dimensions when adding the two.
                                grp1(:,iHb,lstChInc(iCh),iC) = grp1(:,iHb,lstChInc(iCh),iC) + interp1(tHRF,yAvg(:,iHb,lstChInc(iCh),iC),tHRF(:)) * nT;
                                grp1Sum2(:,iHb,lstChInc(iCh),iC) = grp1Sum2(:,iHb,lstChInc(iCh),iC) + interp1(tHRF,ySum2(:,iHb,lstChInc(iCh),iC),tHRF(:));
                            end
                            nTrials_tot{iBlk}(lstChInc(iCh),iC) = nTrials_tot{iBlk}(lstChInc(iCh),iC) + nT;
                        end
                    end
                end
            end

            if ~isempty(grp1)
                for iC = 1:size(grp1,4)
                    for iCh = 1:size(grp1,3)
                        yAvg(:,:,iCh,iC) = grp1(:,:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC);

                        % We want to distinguish between no trials and 1 trial:
                        % If there are no trials, we have no HRF data and no std which
                        % the first case will calculate as opposed to one trial (2nd case)
                        % where we have all zeros.
                        if(nTrials_tot{iBlk}(iCh,iC)~=1)
                            yAvgStd(:,:,iCh,iC) = ( (1/(nTrials_tot{iBlk}(iCh,iC)-1))*grp1Sum2(:,:,iCh,iC) - (nTrials_tot{iBlk}(iCh,iC)/(nTrials_tot{iBlk}(iCh,iC)-1))*(grp1(:,:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC)).^2).^0.5 ;
                        else
                            yAvgStd(:,:,iCh,iC) = zeros(size(grp1Sum2(:,:,iCh,iC)));
                        end

                        %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                        %%%% Now we can set channel descriptors for avg and standard deviation
                        if iRun == length(yAvgRuns)
                            yAvgStdOut(iBlk).AddChannelHbO(ml(iCh,1), ml(iCh,2), iC);
                            yAvgStdOut(iBlk).AddChannelHbR(ml(iCh,1), ml(iCh,2), iC);
                            yAvgStdOut(iBlk).AddChannelHbT(ml(iCh,1), ml(iCh,2), iC);
                        end
                    end

                    %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                    %%%% Now we can set channel descriptors for avg and standard deviation
                    if iRun == length(yAvgRuns)
                        yAvgStdOut(iBlk).AppendDataTimeSeries(yAvgStd(:,:,:,iC));
                    end
                end
            end
        end
    end
end

author: BUNPC


text:
function datatype = fmf_get_datatype( data_format )

if strcmp(data_format,'MONO8'),
  datatype = 'uint8';
else,
  if strcmp(data_format,'YUV422'),
    datatype = 'uint16';
  else,
    error( 'Only MONO8 supported for now' );
  end
end


text:
function varargout = DataSelectorGUI(varargin)
%DataSelectorGUI MATLAB code for DataSelectorGUI.fig
%
%   See also OlyDat.DataSelector.

% Edit the above text to modify the response to help DataSelectorGUI

% Last Modified by GUIDE v2.5 31-May-2013 23:29:24

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @DataSelectorGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @DataSelectorGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% dsGui = DataSelectorGUI(DataSelectorObj)
% DataSelectorObj needs to have its fDataPuller field set.
function DataSelectorGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
handles.qb = varargin{1};

pbWidth = 55;
dx = 5;
x = 10;
y = 8;
handles.pbEdit = uicontrol('tag','pbEdit','callback',@pbEditCbk,'keyPressFcn',@zlclPbKeyPressFcn,...
    'Units','pixels','position', [x y pbWidth 20],'string','Replace',...
    'FontName','default','FontUnits','pixels','FontSize',13,'parent',handles.pnlOuter);
handles.pbAppend = uicontrol('tag','pbAppend','callback',@pbAppendCbk,'keyPressFcn',@zlclPbKeyPressFcn,...
    'Units','pixels','position',[x+pbWidth+dx y pbWidth 20],'string','Add',...
    'FontName','default','FontUnits','pixels','FontSize',13,'parent',handles.pnlOuter);
handles.pbInsert = uicontrol('tag','pbInsert','callback',@pbInsertCbk,'keyPressFcn',@zlclPbKeyPressFcn,...
    'Units','pixels','position',[x+2*(pbWidth+dx) y pbWidth 20],'string','Insert',...
    'FontName','default','FontUnits','pixels','FontSize',13,'parent',handles.pnlOuter);
handles.pbDelete = uicontrol('tag','pbDelete','callback',@pbDeleteCbk,'keyPressFcn',@zlclPbKeyPressFcn,...
    'Units','pixels','position',[x+3*(pbWidth+dx) y pbWidth 20],'string','Delete',...
    'FontName','default','FontUnits','pixels','FontSize',13,'parent',handles.pnlOuter);
handles.pbDeleteAll = uicontrol('tag','pbDeleteAll','callback',@pbDeleteAllCbk,'keyPressFcn',@zlclPbKeyPressFcn,...
    'Units','pixels','position',[x+4*(pbWidth+dx) y pbWidth 20],'string','Reset',...
    'FontName','default','FontUnits','pixels','FontSize',13,'parent',handles.pnlOuter);

guidata(hObject, handles);
set(hObject,'Name',sprintf('OlyDat DataSelector: %s',handles.qb.AssayName));

if handles.qb.fDataPullerHasOptions
    set(handles.pumPullerOptions,'String',handles.qb.fDataPuller.options,'Visible','on');
else
    set(handles.pumPullerOptions,'Visible','off');
end
if handles.qb.fDataPuller.tfPullFSDefault
    set(handles.cbPullFileSys,'Value',1);
else
    set(handles.cbPullFileSys,'Value',0);
end

% tab order
uistack(handles.pnlQuery,'top');
uistack(handles.pumField,'top');
uistack(handles.pumCompare,'top');
uistack(handles.etValue,'top');
uistack(handles.lbValidValues,'top');
uistack(handles.pbEdit,'top');
uistack(handles.pbAppend,'top');
uistack(handles.pbInsert,'top');
uistack(handles.pbDelete,'top');
uistack(handles.pbDeleteAll,'top');
uistack(handles.pbLoad,'top');
uistack(handles.pbSave,'top');

zlclUpdateFileSysControls(handles);

centerfig(hObject);

function varargout = DataSelectorGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

function pbSave_Callback(hObject,~,handles) %#ok<*INUSL,*DEFNU>
[fname pname] = uiputfile('*.qry','Save to file',handles.qb.getCurrentQuerySet);
if isequal(fname,0) || isequal(pname,0)
    return;
end
fullfname = fullfile(pname,fname);
handles.qb.setCurrentQuerySet(fullfname);
handles.qb.fileWrite();

function pbLoad_Callback(hObject,~,handles)
[fname pname] = uigetfile('*.qry','Select query file',handles.qb.getCurrentQuerySet);
if isequal(fname,0) || isequal(pname,0)
    return;
end
fullfname = fullfile(pname,fname);
handles.qb.fileRead(fullfname);
handles.qb.setCurrentQuerySet(fullfname);

function pumField_Callback(hObject, eventdata, handles)
handles.qb.newFieldSelection;

function pumCompare_Callback(hObject, eventdata, handles)
handles.qb.newCompareSelection;

function pbEditCbk(src,~,handles) %#ok<*INUSD>
if nargin < 3
    handles = guidata(src);
end
handles.qb.replaceQuery;

function pbAppendCbk(src,~,handles)
if nargin < 3
    handles = guidata(src);
end
handles.qb.appendQuery;

function pbInsertCbk(src,~,handles)
if nargin < 3
    handles = guidata(src);
end
handles.qb.insertQuery;

function pbDeleteCbk(src,~,handles)
if nargin < 3
    handles = guidata(src);
end
handles.qb.deleteQuery;

function pbDeleteAllCbk(src,~,handles)
if nargin < 3
    handles = guidata(src);
end
handles.qb.deleteAllQueries;
set(handles.pumField,'Value',1);
set(handles.pumCompare,'Value',1);
set(handles.etValue,'String','');
set(handles.lbValidValues,'Value',1);
handles.qb.newFieldSelection();

function pbRunQuery_Callback(hObject,~,handles)
handles.qb.runQuery();

function etValue_KeyPressFcn(hObject, eventdata, handles)
switch eventdata.Key
    case 'return'
        drawnow;
        handles.qb.appendQuery;
end

function lbValidValues_Callback(hObject, eventdata, handles)
switch get(gcf,'SelectionType')
    case 'open' % captures double-click and return
        handles.qb.appendQuery;
end

% Generic key-press function for pushbuttons that translates 'return' into
% 'press'
function zlclPbKeyPressFcn(hObject,eventdata)
handles = guidata(hObject);
c = get(handles.figure1,'CurrentCharacter');
if strcmp(c,sprintf('\r'))
    cbk = get(hObject,'Callback');
    cbk(hObject,[]);
end

function pbRunQuery_KeyPressFcn(hObject, eventdata, handles)
zlclPbKeyPressFcn(hObject,eventdata);

function pbLoad_KeyPressFcn(hObject, eventdata, handles)
zlclPbKeyPressFcn(hObject,eventdata);

function pbSave_KeyPressFcn(hObject, eventdata, handles)
zlclPbKeyPressFcn(hObject,eventdata);

function zlclUpdateFileSysControls(handles)
tf = logical(get(handles.cbPullFileSys,'Value'));
if tf
    enab = 'on';
else
    enab = 'off';
end
set(handles.etDataDir,'enable',enab);
set(handles.pbBrowseFS,'enable',enab);

function cbPullFileSys_Callback(hObject, eventdata, handles)
zlclUpdateFileSysControls(handles);

function pbBrowseFS_Callback(hObject, eventdata, handles)
d = uigetdir([],'Select data directory');
if isequal(d,0)
    % user canceled
elseif exist(d,'dir')~=7
    warning('DataSelectorGUI:dirNotFound','Directory ''%s'' not found.',d);
else
    set(handles.etDataDir,'String',d);
end