text:
function [TrainData, TestData, TrainLabel, TestLabel, numX] = createData_CRA()
    ncFeature = 36;  % features of color
    ntFeature = 51;  % features of text

    % load source domain 1
    A = textread(strcat('data/img/11_color36.txt'));
    B = textread(strcat('data/img/11_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY1(1,i) = 1;
    end
    A = textread(strcat('data/img/15_color36.txt'));
    B = textread(strcat('data/img/15_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY1(1,i+n1) = -1;
    end
    TrainX1 = [Cc;D];
    clear Cc;
    clear D;
    % load source domain 2
    A = textread(strcat('data/img/12_color36.txt'));
    B = textread(strcat('data/img/12_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY2(1,i) = 1;
    end
    A = textread(strcat('data/img/16_color36.txt'));
    B = textread(strcat('data/img/16_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY2(1,i+n1) = -1;
    end
    TrainX2 = [Cc;D];
    clear Cc;
    clear D;
    % load source domain 3
    A = textread(strcat('data/img/13_color36.txt'));
    B = textread(strcat('data/img/13_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY3(1,i) = 1;
    end
    A = textread(strcat('data/img/17_color36.txt'));
    B = textread(strcat('data/img/17_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY3(1,i+n1) = -1;
    end
    TrainX3 = [Cc;D];
    clear Cc;
    clear D;
    % load target domain
    A = textread(strcat('data/img/14_color36.txt'));
    B = textread(strcat('data/img/14_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TestY(1,i) = 1;
    end
    A = textread(strcat('data/img/18_color36.txt'));
    B = textread(strcat('data/img/18_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TestY(1,i+n1) = -1;
    end
    TestX = [Cc;D];

    clear A B Cc D;

    %% normalization
    column = size(TrainX1,2);
    mode_TrainX = sqrt(sum(TrainX1.*TrainX1,1));
    for i = 1 : column
        TrainX1(:,i) = TrainX1(:,i)/mode_TrainX(1,i);
    end
    column = size(TrainX2,2);
    mode_TrainX = sqrt(sum(TrainX2.*TrainX2,1));
    for i = 1 : column
        TrainX2(:,i) = TrainX2(:,i)/mode_TrainX(1,i);
    end
    column = size(TrainX3,2);
    mode_TrainX = sqrt(sum(TrainX3.*TrainX3,1));
    for i = 1 : column
        TrainX3(:,i) = TrainX3(:,i)/mode_TrainX(1,i);
    end
    TrainData{1,1} = TrainX1;
    TrainLabel{1,1} = TrainY1;
    TrainData{1,2} = TrainX2;
    TrainLabel{1,2} = TrainY2;
    TrainData{1,3} = TrainX3;
    TrainLabel{1,3} = TrainY3;

    column = size(TestX,2);
    mode_TestX = sqrt(sum(TestX.*TestX,1));
    for i = 1 : column
        TestX(:,i) = TestX(:,i)/mode_TestX(1,i);
    end
    TestData{1,1} = TestX;
    TestLabel{1,1} = TestY;
    save inputData.mat TrainData TrainLabel TestData TestLabel
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
classdef ChongYaleFaces < handle
% This class is a wrapper over the dataset created by
% Chong You, 2016
% % C. You, D. Robinson, R. Vidal, Scalable Sparse Subspace Clustering by
% Orthogonal Matching Pursuit, CVPR 2016.
% [ExtendedYaleB.mat, https://sites.google.com/site/chongyou1987/] .
% The images are already  cropped.
properties
end
properties(SetAccess=private)
    % the data file from which everything is loaded.
    file_path
    % the raw data loaded from file
    data
    % image height
    image_height
    % image width
    image_width
    % number of subjects in database
    num_subjects
    %  The data set
    Y
    % The labels
    labels
    % The cluster sizes
    cluster_sizes
    % positions for individual clusters
    start_indices
    end_indices
end

properties(Access=private)
    subject_sets_map
end

methods
    function self = ChongYaleFaces()
        root = spx.data_dir;
        self.file_path = fullfile(root, 'faces', 'ExtendedYaleB.mat');
        data = load(self.file_path);
        self.image_height = 48;
        self.image_width = 42;
        self.num_subjects = 38;
        self.labels = data.EYALEB_LABEL;
        self.Y = data.EYALEB_DATA;
        self.cluster_sizes = spx.cluster.cluster_sizes_from_labels(self.labels);
        [self.start_indices, self.end_indices] = spx.cluster.start_end_indices(self.cluster_sizes);
        self.subject_sets_map = cell(1, 10);
    end

    function result = num_total_faces(self)
        result = sum(self.cluster_sizes);
    end

    function result = subject_faces_as_matrices(self, i)
        error('Not implemented.');
        % Returns all faces of a subject as matrices
        result = self.data.I(:, :, i, :);
    end

    function result = subject_faces_as_columns(self, i)
        % Returns all faces of a subject as column vectors
        s= self.start_indices(i);
        e = self.end_indices(i);
        result = self.Y(:, s:e);
    end

    function result = sets_of_n_subjects(self, n)
        if n < 1 || n > 10
            error('Number of subjects must be between 1 and 10');
        end
        if ~isempty(self.subject_sets_map{n})
            result = self.subject_sets_map{n};
            return;
        end
        result1 = nchoosek(1:10, n);
        result2 = nchoosek(11:20, n);
        result3 = nchoosek(21:30, n);
        result = [result1; result2; result3];
        if n <= 8
            result4 = nchoosek(31:38, n);
            result = [result; result4];
        end
        self.subject_sets_map{n} = result;
    end

    function result = num_experiments_of_n_subjects(self, n)
        sets = self.sets_of_n_subjects(n);
        result = size(sets, 1);
    end

    %% experiment_data: Returns data for a particular experiment
    function [Y, cluster_sizes, labels] = experiment_data(self, n, i)
        sets = self.sets_of_n_subjects(n);
        nn  = size(sets, 1);
        if i < 1 || i > nn
            error('Invalid experiment number.');
        end
        indices = sets(i, :);
        Y = [];
        cluster_sizes = self.cluster_sizes(indices);
        for k=1:n
            subject_index = indices(k);
            Y = [Y self.subject_faces_as_columns(subject_index)];
        end
        labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes);
    end

    %% Generates data for specified subjects
    function [Y, cluster_sizes, labels] = data_for_subjects(self, subject_list)
        cluster_sizes = self.cluster_sizes(subject_list);
        Y = [];
        for k=1:numel(subject_list)
            subject_index = subject_list(k);
            Y = [Y self.subject_faces_as_columns(subject_index)];
        end
        labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes);
    end

    function [faces, cluster_sizes, labels] = all_faces(self)
        faces = self.Y;
        cluster_sizes = self.cluster_sizes;
        labels = self.labels;
    end

    function canvas = create_random_canvas(self, rows, cols)
        n = self.num_subjects;
        if nargin < 2
            rows = 4;
        end
        if nargin < 3
            cols = rows;
        end
        k  = rows * cols;
        indices = randperm(n, k);
        faces = self.all_faces;
        % pick first image for each subject
        indices = self.start_indices(indices);
        images = faces(:, indices);
        canvas = spx.graphics.canvas.create_image_grid(images, rows, cols, ...
            self.image_height, self.image_width);
        canvas = uint8(canvas);
    end

end

end

author: indigits

repo: sparse-plex


text:
function [ebsd,options] = loadEBSD_generic(fname,varargin)
% load ebsd data from generic text files
%
% Description
%
% *loadEBSD_generic* loads EBSD data from text or exel files that have a
% column oriented format as
%
%  x1 y1 phi1_1 Phi_1 phi2_1 phase_1
%  x2 y2 phi1_2 Phi_2 phi2_2 phase_2
%  x2 y3 phi1_3 Phi_3 phi2_3 phase_3
%  .      .       .       .
%  .      .       .       .
%  .      .       .       .
%  xM yM phi1_M Phi_M phi2_M phase_m
%
% The assoziation of the columns as Euler angles, phase informationl, etc.
% is specified by the options |ColumnNames| and |Columns|. The files can be
% contain any number of header lines.
%
% Syntax
%   ebsd   = loadEBSD_generic(fname,'ColumnNames',{'x','y','Euler1','Euler2','Euler3','phase'},'cs',cs)
%
% Input
%  fname - file name (text files only)
%  cs - @crystalSymmetry or cell array of @crystalSymmetry
%
% Options
%  ColumnNames       - names of the colums to be imported, mandatory are euler 1, euler 2, euler 3
%  Columns           - postions of the columns to be imported
%  radians           - treat input in radiand
%  delimiter         - delimiter between numbers
%  header            - number of header lines
%  Bunge             - [phi1 Phi phi2] Euler angle in Bunge convention (default)
%  ABG               - [alpha beta gamma] Euler angle in Mathies convention
%  passive           -
%
%
% Example
%
%   fname = fullfile(mtexDataPath,'EBSD','85_829grad_07_09_06.txt');
%   CS = {'notIndexed',...
%          crystalSymmetry('m-3m','mineral','Fe'),...
%          crystalSymmetry('m-3m','mineral','Mg')};
%   SS = specimenSymmetry('triclinic');
%   ebsd = loadEBSD_generic(fname,'CS',CS,'SS',SS, 'ColumnNames', ...
%     {'Index' 'Phase' 'x' 'y' 'Euler1' 'Euler2' 'Euler3' 'MAD' 'BC' 'BS'...
%     'Bands' 'Error' 'ReliabilityIndex'}, 'Bunge')
%
% See also
% ImportEBSDData loadEBSD ebsd_demo

max_num_phases = 40;

try
  % load data
  if ~isnumeric(fname)
    [d,options,header,c] = load_generic(char(fname),varargin{:});
    varargin = options;
  else
    d = fname;
  end

  % no data found
  if size(d,1) < 1 || size(d,2) < 3
    error('Generic interface could not detect any numeric data in %s',fname);
  end

  % no options given -> ask
  if ~check_option(varargin,'ColumnNames')

    options = generic_wizard('data',d(1:end<101,:),'type','EBSD','header',header,'columns',c);
    if isempty(options), ebsd = []; return; end
    varargin = [options,varargin];

  end

  loader = loadHelper(d,varargin{:});

  q      = loader.getRotations();

  % assign phases
  if loader.hasColumn('Phase')

    phase = loader.getColumnData('Phase');

  else

    phase = ones(length(q),1);

  end

  if max(phase)>max_num_phases

    warning('MTEX:toomanyphases', ...
      ['Found more than ' num2str(max_num_phases) '. I''m going to ignore them all.']);
    phase = ones(size(d,1),1);

  end

  % compute unit cells
  if loader.hasColumn({'x' 'y'})

    varargin = [varargin,'unitCell',....
      calcUnitCell(loader.getColumnData({'x' 'y' 'z'}),varargin{:})];

  end

  % return varargin as options
  opt = loader.getOptions('ignoreColumns','Phase');

  % set up EBSD variable
  CSList{1} = 'notIndexed';
  for k = 1:max(1,length(unique(phase))-1)
    CSList{k+1} = crystalSymmetry('432','mineral',['unknown' int2str(k)]);
  end
  CSList = get_option(varargin,'CS',CSList);
  ebsd = EBSD(q,phase,CSList,opt,varargin{:});

catch
  interfaceError(fname)
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function R = rotz(alpha)

    % ROTZ computes a rotation along the z axis (rotation matrix).
    %
    % FORMAT: R = rotz(alpha)
    %
    % INPUT:  - alpha = angle in radians
    %
    % OUTPUT: - R = [3 * 3] rotation matrix along z axis
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

   R = [cos(alpha), -sin(alpha), 0;
        sin(alpha),  cos(alpha), 0;
        0,           0,          1];

end

author: robotology

repo: whole-body-controllers


text:
function xy=mmstream2(x,y,u,v,x0,y0,mark,step)
%MMSTREAM2 Improved 2D Streamlines.
% XY = MMSTREAM2(X,Y,U,V,X0,Y0,Mark,Step) computes streamlines from gradient
% data in matrices U and V.
%
% X and Y can be vectors defining the coordinate axes data where U(i,j) and
% V(i,j) coincide with the coordinates axes points X(j) and Y(i).
% Alternatively, X and Y can be 2D plaid matrices as produced by MESHGRID.
%
% X0 and Y0 are equal length vectors defining coordinates that mark the
% Start, End, or a point On individual streamlines as denoted by the input
% Mark which is 'Start', 'End' or 'On'. If empty or not given Mark='Start'.
%
% Step identifies the normalized step size used. If empty or not given,
% Step = 0.1, i.e., 1/10 of a cell. 0.01 <= Step <= 0.5
%
% XY is a cell array containing streamline data points. XY{k}(:,1) contains
% the x-axis data and XY{k}(:,2) contains the y-axis data for the k-th
% streamline.
%
% Improvements over MATLAB's STREAM2:
% Will find closed streamlines.
% Will selectively go downstream, upstream, or both directions from a point.
% Uses a more accurate integration routine.
% X and Y need not linearly spaced.
%
% STREAMLINE(XY) plots the streamlines on the current axes.

% D.C. Hanselman, University of Maine, Orono, ME 04469
% MasteringMatlab@yahoo.com
% Mastering MATLAB 7
% 2005-06-21
% Revised 2005-08-31, 2007-07-03, 2008-07-29

if nargin<8 || isempty(step)                        % parse input arguments
    step=0.1;
end
if nargin<7 || isempty(mark)
    mark='start';
end
if nargin<6
    error('At Least Six Input Arguments are Required.')
end
if ~isnumeric(step) || numel(step)>1 || step<0.01 || step>0.5
    error('Step Must be a Scalar Between 0.01 and 0.5')
end
if ~ischar(mark) || ~any(lower(mark(1))=='seo')
    error('Mark Must be ''Start'', ''End'', or ''On''')
end
mark=lower(mark(1));
if any(size(u)~=size(v))
    error('U and V Must be the Same Size.')
end
if ndims(u)~=2 || min(size(u))<2
    error('U and V Must be 2D and at Least 2-by-2.')
end
if all(size(x)==size(y)) && min(size(x))>1 % x and y are plaid
    xx=x;                  % save plaid
    yy=y;
    x=x(1,:)';             % create axes vectors
    y=y(:,1);
else                                       % x and y are axes vectors
    x=x(:);                % save axes vectors
    y=y(:);
    [xx,yy]=meshgrid(x,y); % create plaid
end
if any(abs(diff(x))<eps) || any(abs(diff(y))<eps)
    error('X and Y Must Not Have Consecutive Equal Values.')
end
Nx=length(x);
Ny=length(y);
if Nx~=size(u,2)
    error('X Must Have as Many Elements or Columns as U and V Have Columns.')
end
if Ny~=size(u,1)
    error('Y Must Have as Many Elements or Rows as U and V Have Rows.')
end
x0=x0(:);
y0=y0(:);
N0=length(x0);
if N0~=length(y0)
    error('X0 and Y0 Must Contain the Same Number of Elements.')
end

xy=cell(1,N0); % create cells for output

for k=1:N0 % find streamlines given [x0, y0] pairs

    if x0(k)<min(x) || x0(k)>max(x) ||...
            y0(k)<min(y) || y0(k)>max(y)   % point is outside map
        continue
    end
    % figure out what cell [x0,y0] is in
    [idx,idx]=min(abs(x0(k)-xx(:))+abs(y0(k)-yy(:)));
    [i,j]=ind2sub(size(u),idx);
    if j==1
        jlim=[j j+1];
    elseif j==Nx
        jlim=[j-1 j];
    elseif abs(x0(k)-x(j-1))<abs(x0(k)-x(j+1))
        jlim=[j-1 j];
    else
        jlim=[j j+1];
    end
    if i==1
        ilim=[i i+1];
    elseif i==Ny
        ilim=[i-1 i];
    elseif abs(y0(k)-y(i-1))<abs(y0(k)-y(i+1))
        ilim=[i-1 i];
    else
        ilim=[i i+1];
    end
    switch mark
        case 's'
            xy{k}=local_getstream(x,y,u,v,ilim,jlim,x0(k),y0(k),step);
        case 'e'
            xy{k}=local_getstream(x,y,-u,-v,ilim,jlim,x0(k),y0(k),step);
            xy{k}=xy{k}(end:-1:1,:); % make stream flow downhillm
        case 'o'
            [xy{k},ic]=local_getstream(x,y,u,v,ilim,jlim,x0(k),y0(k),step);
            if ~ic % go other direction only if streamline is not closed
                xye=local_getstream(x,y,-u,-v,ilim,jlim,x0(k),y0(k),step);
                xy{k}=[xye(end:-1:2,:);xy{k}]; % make stream flows downhill
            end
    end
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [xy,ic]=local_getstream(x,y,u,v,ilim,jlim,x0,y0,step)
% move from cell to cell gathering a streamline
% xy is a double array
% ic is true if streamline is closed

ic=false;    % logical variable that is true if streamline is closed
xy0=[x0 y0]; % hold starting point so that closed streamlines can be found

xy=zeros(1000,2); % allocate storage
k=1; % pointer to next index to store data

loop2=tic;
while true
    if toc(loop2)>0.1
        disp('Streamline timeout')
        xy(k-1:end,:)=[]; % throw away unused storage
        break
    end
    xlim=x(jlim);
    ylim=y(ilim);
    ulim=u(ilim,jlim);
    vlim=v(ilim,jlim);

    [dxy,ii,jj]=local_cellstream(xlim,ylim,ulim,vlim,x0,y0,xy0,step);
    np=size(dxy,1);
    if k+np-1>size(xy,1);   % allocate more storage if needed
        xy=[xy;zeros(500,2)];%#ok
    end
    xy(k:k+np-1,:)=dxy; % poke data into storage array
    k=k+np;             % update pointer
    jlim=jlim+jj;
    ilim=ilim+ii;
    if jlim(1)<1 || jlim(2)>length(x) || ...
            ilim(1)<1 || ilim(2)>length(y)
        xy(k-1:end,:)=[]; % throw away unused storage
        break
    elseif (ii==0 && jj==0) || np==0
        xy(k:end,:)=[]; % throw away unused storage
        ic=isequal(xy(end,:),xy0);
        break
    end
    x0=dxy(end,1);
    y0=dxy(end,2);
end
%--------------------------------------------------------------------------
function [xy,ii,jj]=local_cellstream(xlim,ylim,ulim,vlim,x0,y0,xy0,step)

% starting at the point [sx, sy] create a streamline until the edge of the
% cell given by xlim and ylim is reached.
%
% xlim = [x1 x2]                ylim = [y1 y2]
%
% ulim = [u(x1,y1) u(x2,y1)     vlim = [v(x1,y1) v(x2,y1)
%         u(x1,y2) u(x2,y2)];           v(x1,y2) v(x2,y2)];

xymin=[xlim(1) ylim(1)]; % [x1 y1]
xymax=[xlim(2) ylim(2)]; % [x2 y2]
xybox=xymax-xymin;       % [x2-x1 y2-y1]
uv=[ulim(:) vlim(:)]';   % [u(x1,y1) u(x1,y2) u(x2,y1) u(x2,y2);
%  v(x1,y1) v(x1,y2) v(x2,y1) v(x2,y2)];
tol=1e-4;
N=round(3/step);
xy=zeros(N,2); % preallocate memory for result
cstep=1.5*step;% step size for closed streamline detection

k=1;
xy(k,:)=[x0 y0]; % first point
ii=0;            % next cell in x
jj=0;            % and y direction

%william mod:
loop1=tic;
while true
    if toc(loop1)>0.1
        xy(k+1:end,:)=[];
        disp('anderer')
        disp('Streamline timeout')
        break
    end
    abk=(xy(k,:)-xymin)./xybox; % normalized current position [alpha beta]
    % compute slopes at current point
    uvk=((1-abk(1))*( (1-abk(2))*uv(:,1) + abk(2)*uv(:,2)) + ...
        abk(1) *( (1-abk(2))*uv(:,3) + abk(2)*uv(:,4)))';      % [uk vk]

    if k>1 && max(abs(uvk))<tol  % at minimum, this stream stops
        xy(k+1:end,:)=[];
        break
    end

    h=min(step*abs(xybox)./(max(abs(uvk),tol)));       % allowable step size

    xy(k+1,:)=xy(k,:) + h*uvk;                    % Forward Euler Prediction

    abi=(xy(k+1,:)-xymin)./xybox;        % normalized position at next point
    uvi=((1-abi(1))*( (1-abi(2))*uv(:,1) + abi(2)*uv(:,2)) + ...
        abi(1) *( (1-abi(2))*uv(:,3) + abi(2)*uv(:,4)))';      % [ui vi]

    xy(k+1,:)=xy(k,:) + h*(uvk + uvi)/2;       % Trapezoidal Rule Correction
    k=k+1;

    if k>2 && norm((xy(k,:)-xy(k-2,:))./xybox,inf)<step/2% stuck inside cell
        xy(k-2:end,:)=[];
        break
    elseif k>2 && norm((xy(k,:)-xy0)./xybox,inf)<=cstep   % closed streamline
        xy(k,:)=xy0; % close the streamline
        xy(k+1:end,:)=[];
        break
    elseif k==N+1                                        % stuck inside cell
        break
    else                                       % check if moved outside cell

        jj=sign(xybox(1))*(-(xy(k,1)<xlim(1)) + (xy(k,1)>xlim(2)));
        ii=sign(xybox(2))*(-(xy(k,2)<ylim(1)) + (xy(k,2)>ylim(2)));
        if jj~=0 || ii~=0                                % point to next cell
            xy(k+1:end,:)=[];
            break
        end
    end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function h = hesschek(net, x, t)
%HESSCHEK Use central differences to confirm correct evaluation of Hessian matrix.
%
%	Description
%
%	HESSCHEK(NET, X, T) takes a network data structure NET, together with
%	input and target data matrices X and T, and compares the evaluation
%	of the Hessian matrix using the function NETHESS and using central
%	differences with the function NETERR.
%
%	The optional return value H is the Hessian computed using NETHESS.
%
%	See also
%	NETHESS, NETERR
%

%	Copyright (c) Ian T Nabney (1996-2001)

w0 = netpak(net);
nwts = length(w0);
h = nethess(w0, net, x, t);

w = w0;
hcent = zeros(nwts, nwts);
h1 =  0.0; h2 =  0.0; h3 =  0.0; h4 = 0.0;
epsilon = 1.0e-4;
fprintf(1, 'Checking Hessian ...\n\n');
for k = 1:nwts;
  for l = 1:nwts;
    if(l == k)
      w(k) = w0(k) + 2.0*epsilon;
      h1 = neterr(w, net, x, t);
      w(k) = w0(k) - 2.0*epsilon;
      h2 = neterr(w, net, x, t);
      w(k) = w0(k);
      h3 = neterr(w, net, x, t);
      hcent(k, k) = (h1 + h2 - 2.0*h3)/(4.0*epsilon^2);
    else
      w(k) = w0(k) + epsilon;
      w(l) = w0(l) + epsilon;
      h1 = neterr(w, net, x, t);
      w(k) = w0(k) - epsilon;
      w(l) = w0(l) - epsilon;
      h2 = neterr(w, net, x, t);
      w(k) = w0(k) + epsilon;
      w(l) = w0(l) - epsilon;
      h3 = neterr(w, net, x, t);
      w(k) = w0(k) - epsilon;
      w(l) = w0(l) + epsilon;
      h4 = neterr(w, net, x, t);
      hcent(k, l) = (h1 + h2 - h3 - h4)/(4.0*epsilon^2);
      w(k) = w0(k);
      w(l) = w0(l);
    end
  end
end

fprintf(1, '   analytical    numerical       delta\n\n');
temp = [h(:), hcent(:), (h(:) - hcent(:))];
fprintf(1, '%12.6f  %12.6f  %12.6f\n', temp');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = showSimilarFrames(varargin)
% SHOWSIMILARFRAMES MATLAB code for showSimilarFrames.fig
%      SHOWSIMILARFRAMES, by itself, creates a new SHOWSIMILARFRAMES or raises the existing
%      singleton*.
%
%      H = SHOWSIMILARFRAMES returns the handle to a new SHOWSIMILARFRAMES or the handle to
%      the existing singleton*.
%
%      SHOWSIMILARFRAMES('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SHOWSIMILARFRAMES.M with the given input arguments.
%
%      SHOWSIMILARFRAMES('Property','Value',...) creates a new SHOWSIMILARFRAMES or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before showSimilarFrames_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to showSimilarFrames_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 15-Mar-2012 16:36:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @showSimilarFrames_OpeningFcn, ...
                   'gui_OutputFcn',  @showSimilarFrames_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

% --- Executes just before showSimilarFrames is made visible.
function showSimilarFrames_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to showSimilarFrames (see VARARGIN)

% Choose default command line output for showSimilarFrames
handles.output = hObject;

handles.maxFrames = 20;
handles.numSimilar = 4;
handles.halfSize = 30*5;
handles.imgX = 1024;
handles.imgY = 1024;
handles.frameNo = handles.maxFrames+1;
handles.isPlaying = false;
handles.fwdLimit = 5;
handles.revLimit = 5;
handles.fps = 3;
handles.maxfps = 10;
handles.tsize = 500;
handles.cache = containers.Map('keyType','char','valueType','any');
handles.display = 'manual';

set(handles.frameSlider,'SliderStep',[1/(2*handles.maxFrames) 3/(2*handles.maxFrames)]);

set(handles.frame_range_edit, 'string', '5');
set(handles.fps_edit, 'string', '5');
set(handles.FrameText,'String','frame:0');

sz = handles.halfSize;
tlsz = handles.maxFrames;

linkaxesGroup = [];

for ndx = 1:handles.numSimilar
  tName = sprintf('top%d',ndx);
  mName = sprintf('middle%d',ndx);
  bName = sprintf('bottom%d',ndx);
  handles.hAxes(1,ndx) = initAxes(handles,handles.(tName),sz,1,ndx);
  handles.hAxes(2,ndx) = initAxes(handles,handles.(mName),sz,2,ndx);
  handles.hAxes(3,ndx) = initAxes(handles,handles.(bName),sz,3,ndx);
  handles.hTimeline(1,ndx) = initTimeline(handles.(['tl_',tName]),tlsz,handles.frameNo);
  handles.hTimeline(2,ndx) = initTimeline(handles.(['tl_',mName]),tlsz,handles.frameNo);
  handles.hTimeline(3,ndx) = initTimeline(handles.(['tl_',bName]),tlsz,handles.frameNo);
  linkaxesGroup = [linkaxesGroup, handles.(tName)];
  linkaxesGroup = [linkaxesGroup, handles.(mName)];
  linkaxesGroup = [linkaxesGroup, handles.(bName)];
end

linkaxes(linkaxesGroup,'xy');
%plot(handles.axes_slider_zero,[0,1],[0,0],'y-','LineWidth',5);
% Update handles structure
guidata(hObject, handles);

% UIWAIT makes showSimilarFrames wait for user response (see UIRESUME)
% uiwait(handles.figure_showSimilarFrames);

% --- Outputs from this function are returned to the command line.
function varargout = showSimilarFrames_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

function SetJLabelData(hObject,obj,HJLabel)
handles = guidata(hObject);
handles.JLDobj = obj;
handles.JLabelHandles = HJLabel;
handles.imgX = HJLabel.guidata.movie_width;
handles.imgY = HJLabel.guidata.movie_height;
guidata(hObject,handles);

%{
function loadmovies(hObject)
% Loads the movies.

handles = guidata(hObject);
pointer = [];
movieNames = handles.JLDobj.expdirs;
for movieNum = 1:length(movieNames)
[pointer(movieNum).readframe,pointer(movieNum).nframes,pointer(movieNum).movie_fid,pointer(movieNum).movieheaderinfo] = ...
  get_readframe_fcn(sprintf('%s/movie.ufmf',movieNames{movieNum}),'interruptible',false);
end
handles.movieNames = movieNames;
handles.pointer = pointer;

guidata(hObject,handles);
%}

function CacheTracksLabeled(hObject,exps)
% Load the tracks.

handles = guidata(hObject);
blk = handles.tsize;
off = handles.maxFrames;
fprintf('Loading ');
if nargin<2
  exps = 1:handles.JLDobj.nexps;
end
for expNum = exps
  fprintf('.');
  trxfile = handles.JLDobj.GetFile('trx',expNum);
  trx = load_tracks(trxfile);

  handles.firstframe{expNum} = [trx.firstframe];
  handles.endframe{expNum} = [trx.endframe];

  % Get labels. Store only those parts of tracks that have been labeled.
  for curFly = 1:length(trx)
    labels = handles.JLDobj.GetLabels(expNum,curFly);
    for curNdx = 1:length(labels.t0s)
      trxOffset = handles.firstframe{expNum}(curFly);
      curblockStart = floor( (labels.t0s(curNdx)-off-trxOffset+1)/blk)+1;
      curblockEnd = ceil( (labels.t1s(curNdx)+off-trxOffset+1)/blk);
      for curBlk = curblockStart:curblockEnd
        idStr = sprintf('%d_%d_%d',expNum,curFly,curBlk);
        if ~isKey(handles.cache,idStr),
          tStart = max((curBlk-1)*blk+1,1);
          tEnd = min(curBlk*blk,trx(curFly).endframe-trxOffset+1);
          tSlice = tStart:tEnd;
          trax = [];
          trax.X = trx(curFly).x(tSlice);
          trax.Y = trx(curFly).y(tSlice);
          trax.theta = trx(curFly).theta(tSlice);
          trax.maj = trx(curFly).a(tSlice);
          trax.min = trx(curFly).b(tSlice);
          handles.cache(idStr) = trax;
        end
      end
    end

  end
end
fprintf(' Done loading tracks\n');

guidata(hObject,handles);

function add_prep_list(hObject,exps)
handles = guidata(hObject);
if ~isempty(handles.JLDobj.allperframefns)
  for i = 1:numel(handles.JLDobj.allperframefns),
    handles.timeline_prop_options{i} = handles.JLDobj.allperframefns{i};
  end
  %handles.d = handles.data.allperframefns(1);
  handles.perframepropis = 1;
  set(handles.timeline_properties,'String',handles.timeline_prop_options,'Value',1);
  %handles.timeline_data_ylims = nan(2,numel(handles.data.allperframefns));
end
guidata(hObject,handles);

function CacheTracks(hObject,expNum,curFly,t0,t1)
handles = guidata(hObject);
blk = handles.tsize;
off = handles.maxFrames;
trxfile = handles.JLDobj.GetFile('trx',expNum);
trx = load_tracks(trxfile);

handles.firstframe{expNum} = [trx.firstframe];
handles.endframe{expNum} = [trx.endframe];
trxOffset = handles.firstframe{expNum}(curFly);
curblockStart = floor( (t0-off-trxOffset+1)/blk)+1;
curblockEnd = ceil( (t1+off-trxOffset+1)/blk);

for curBlk = curblockStart:curblockEnd
  idStr = sprintf('%d_%d_%d',expNum,curFly,curBlk);
  if ~isKey(handles.cache,idStr),
    tSlice = ((curBlk-1)*blk+1):(curBlk*blk);
    tSliceValid = tSlice >= 1 & tSlice<=(trx(curFly).endframe-trxOffset+1);
    trax = [];
    trax.X = nan(1,length(tSlice));
    trax.Y = nan(1,length(tSlice));
    trax.theta = nan(1,length(tSlice));
    trax.maj = nan(1,length(tSlice));
    trax.min = nan(1,length(tSlice));

    validNdx = tSlice(tSliceValid);
    trax.X(tSliceValid) = trx(curFly).x(validNdx);
    trax.Y(tSliceValid) = trx(curFly).y(validNdx);
    trax.theta(tSliceValid) = trx(curFly).theta(validNdx);
    trax.maj(tSliceValid) = trx(curFly).a(validNdx);
    trax.min(tSliceValid) = trx(curFly).b(validNdx);
    handles.cache(idStr) = trax;
  end
end
guidata(hObject,handles);

function setFrames(hObject,frameData)
% Sets the frames to be shown.

handles = guidata(hObject);

handles.frames{1} = frameData.posFrames;
handles.frames{2} = frameData.curFrame;
handles.frames{3} = frameData.negFrames;
guidata(hObject,handles);

readFrames(hObject);
handles = guidata(hObject);
updatePlots(hObject,handles,handles.frameNo);

%zoom to the center
limVal = [handles.halfSize*0.8, handles.halfSize*1.2];
set(handles.middle1,'xLim',limVal,'yLim',limVal);

function readFrames(hObject)
handles = guidata(hObject);
handles.isPlaying = false;
guidata(hObject,handles);
curImg = getRotatedFrame(handles,handles.frames{2});
relTrx = getRelTrx(handles,handles.frames{2});
[labels predictions] = getLabels(handles,handles.frames{2});
ptrx = labels.ptrx;
ntrx = labels.ntrx;

plotPredictionsInAutoMode(handles);
poltLabelsInManualMode(handles);

%for row 2
for t = 1:length(ptrx)
    % tranfer data to image for manual timeline
    timeline{2}(1,t,:) = [0.2 0.2 0.2];
    if ~isnan(ptrx(t)); timeline{2}(1,t,:) = [0.7 0.0 0.0]; end
    if ~isnan(ntrx(t)); timeline{2}(1,t,:) = [0.0 0.0 0.7]; end

    % tranfer data to image for automatic timeline
    if isnan(predictions(t)), timeline{2}(2,t,:) = [0.2 0.2 0.2];
    elseif predictions(t)<1.5, timeline{2}(2,t,:) = [0.7 0.0 0.0];
    else timeline{2}(2,t,:) = [0.0 0.0 0.7];
    end
end

for ndx = 1:handles.numSimilar
    handles.im{2,ndx} = curImg;
    handles.trx{2,ndx} = relTrx;
    set(handles.hAxes(2,ndx).labelPos,'visible', 'on');
    set(handles.hAxes(2,ndx).labelNeg,'visible', 'on');
    set(handles.hAxes(2,ndx).predPos,'visible', 'off');
    set(handles.hAxes(2,ndx).predNeg,'visible', 'off');
    set(handles.hAxes(2,ndx).trax,'XData',relTrx.X,'YData',relTrx.Y);
    set(handles.hTimeline(2,ndx).image, 'CData', timeline{2});
    set(handles.hTimeline(2,ndx).hcurr,'XData',[handles.frameNo,handles.frameNo]);
end

%create property plot
%set(handles.axes_timeline_properties,'YLimMode','auto');
hold(handles.axes_timeline_properties,'on');
handles.htimeline_propdata{2} = plot(handles.axes_timeline_properties,nan,nan,'.-');

%choose the first property on the popup list
s = handles.timeline_prop_options{1};
prop = find(strcmpi(s,handles.JLDobj.allperframefns),1);
handles.perframepropis = prop;
handles.perframeprops = s;

[perframedata,T0,T1] = handles.JLDobj.GetPerFrameData(handles.frames{2}.expNum,handles.frames{2}.flyNum,prop);
TStart = handles.frames{2}.curTime - handles.maxFrames;
IStart = TStart - T0 + 1;
TEnd = handles.frames{2}.curTime + handles.maxFrames;
IEnd = TEnd - T0 + 1;
perframedata = perframedata(IStart:IEnd);
set(handles.htimeline_propdata{2},'XData',-handles.maxFrames:handles.maxFrames,'YData',perframedata,'Color','k', 'LineWidth', 2);

%update text_timeline_prop
handles.cur_rb_row = 2;
handles.cur_rb_col = 1;
perframedata = handles.JLDobj.GetPerFrameData1(handles.frames{2}.expNum,handles.frames{2}.flyNum,prop,handles.frames{2}.curTime);
s = sprintf('%.3f',perframedata);
set(handles.text_timeline_prop,'String',s);

%for row 1 and 3
for row = [1 3]
    for ndx = 1:handles.numSimilar
        relTrx = getRelTrx(handles,handles.frames{row}(ndx));
        handles.im{row,ndx} = getRotatedFrame(handles,handles.frames{row}(ndx));

        handles.trx{row,ndx} = relTrx;
        [labels predictions] = getLabels(handles,handles.frames{row}(ndx));
        ptrx = labels.ptrx;
        ntrx = labels.ntrx;

        % plot for manual timeline
        for t = 1:length(ptrx)
            timeline{row,ndx}(1,t,:) = [0.2 0.2 0.2];
            if ~isnan(ptrx(t)); timeline{row,ndx}(1,t,:) = [0.7 0.0 0.0]; end
            if ~isnan(ntrx(t)); timeline{row,ndx}(1,t,:) = [0.0 0.0 0.7]; end

            % plot for automatic timeline
            if isnan(predictions(t)), timeline{row,ndx}(2,t,:) = [0.2 0.2 0.2];
            elseif predictions(t)<1.5, timeline{row,ndx}(2,t,:) = [0.7 0.0 0.0];
            else timeline{row,ndx}(2,t,:) = [0.0 0.0 0.7];
            end
        end

        set(handles.hAxes(row,ndx).labelPos,'visible', 'on');
        set(handles.hAxes(row,ndx).labelNeg,'visible', 'on');
        set(handles.hAxes(row,ndx).predPos,'visible', 'off');
        set(handles.hAxes(row,ndx).predNeg,'visible', 'off');

        set(handles.hAxes(row,ndx).trax,'XData',relTrx.X,'YData',relTrx.Y);
        set(handles.hTimeline(row,ndx).image, 'CData', timeline{row,ndx});
        fprintf('.');

        %create property plots
        if row == 1
            propColor = 'b';
        else
            propColor = 'r';
        end

        handles.htimeline_propdata{row}(ndx) = plot(handles.axes_timeline_properties,nan,nan,'.-','color', propColor);

        [perframedata,T0,T1] = handles.JLDobj.GetPerFrameData(handles.frames{row}(ndx).expNum,handles.frames{row}(ndx).flyNum,prop);
        TStart = handles.frames{row}(ndx).curTime - handles.maxFrames;
        IStart = TStart - T0 + 1;
        TEnd = handles.frames{row}(ndx).curTime + handles.maxFrames;
        IEnd = TEnd - T0 + 1;
        perframedata = perframedata(IStart:IEnd);
        set(handles.htimeline_propdata{row}(ndx),'XData',-handles.maxFrames:handles.maxFrames,'YData',perframedata);
    end
    fprintf('\n');
end

%plot hcurr for axes_timeline_properties
yLimVal=ylim(handles.axes_timeline_properties);
frameNo = handles.frameNo - handles.maxFrames - 1;
handles.hPropOrig= plot(handles.axes_timeline_properties,[frameNo,frameNo],[yLimVal(1),yLimVal(2)],'g--');
handles.hPropCurr= plot(handles.axes_timeline_properties,[frameNo,frameNo],[yLimVal(1),yLimVal(2)],'k-');
guidata(hObject,handles);

function updatePlots(hObject,handles,frameNo)

  for row = 1:3
    for ndx = 1:handles.numSimilar
      set(handles.hAxes(row,ndx).image,'CData',uint8(handles.im{row,ndx}(:,:,:,frameNo)));
      set(handles.hTimeline(row,ndx).hcurr,'XData', [frameNo,frameNo]);
      set(handles.hPropCurr,'XData',[frameNo- handles.maxFrames - 1,frameNo- handles.maxFrames - 1]);
      updatefly(handles.hAxes(row,ndx).fly,handles.hAxes(row,ndx).flyPred,handles.trx{row,ndx},frameNo, handles.display);
    end
    set(handles.FrameText,'String',sprintf('frame:%d',frameNo-handles.maxFrames-1));
  end

  %update the text_timeline_prop value
  curExperiment = handles.frames{handles.cur_rb_row}(handles.cur_rb_col).expNum;
  curFly = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).flyNum;
  curTime = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).curTime+handles.frameNo-handles.maxFrames-1;
  curProp= handles.perframepropis;
  perframedata = handles.JLDobj.GetPerFrameData1(curExperiment,curFly,curProp,curTime);
  s = sprintf('%.3f',perframedata);
  set(handles.text_timeline_prop,'String',s);

%   guidata(hObject,handles);

function play(hObject)
  handles = guidata(hObject);
  frameNo = handles.frameNo;
  ticLoop = tic; % for loop
  ticUpdate = tic; % for updates;
  while(handles.isPlaying)
    if(frameNo-handles.maxFrames>handles.fwdLimit)
      frameNo = handles.maxFrames+1-handles.revLimit;
    else
      frameNo = frameNo+1;
    end
    handles.frameNo = frameNo;
    guidata(hObject,handles);
    tmp = toc(ticLoop);
    if tmp < 1/handles.fps
      pause(1/handles.fps-tmp);
    end
    ticLoop = tic;
    updatePlots(hObject,handles,frameNo);
    set(handles.frameSlider,'Value', 1-frameNo/(2*handles.maxFrames+1));

    if(toc(ticUpdate)>0.5)
      handles = guidata(hObject);
    end
  end
  guidata(hObject,handles);

function [labels predictions]= getLabels(handles,curFrame)
% Reads labels and predictions.

  curExp = curFrame.expNum;
  curFly = curFrame.flyNum;
  curTime = curFrame.curTime;
  sz = handles.maxFrames;

  tSlice = curTime-sz:curTime+sz;
  ltrx = nan(1,length(tSlice));
  labels.ptrx = nan(1,length(tSlice));
  labels.ntrx = nan(1,length(tSlice));

  curLabel = handles.JLDobj.GetLabels(curExp,curFly);
%   windowNdx = find( (handles.JLDobj.windowdata.exp == curExp) & ...
%       (handles.JLDobj.windowdata.flies == curFly) & ...
%       (handles.JLDobj.windowdata.t == curTime) ,1);

  idxcurr = handles.JLDobj.windowdata.exp == curExp & ...
    all(bsxfun(@eq,handles.JLDobj.windowdata.flies,curFly),2) & ...
    handles.JLDobj.windowdata.t >= (curTime-sz) & ...
    handles.JLDobj.windowdata.t <= (curTime+sz) & ...
    handles.JLDobj.windowdata.isvalidprediction;

  tNdx = handles.JLDobj.windowdata.t(idxcurr);
  predictions = nan(1,length(tSlice));
  predictions(  (tNdx-curTime+sz+1) ) = handles.JLDobj.windowdata.predicted(idxcurr);
  curT0 = curLabel.t0s;
  curT1 = curLabel.t1s;
  curNames = curLabel.names;

  for bnum = 1:length(curNames)
    if strcmpi(curNames{bnum},'None')
      curVal = -1;
    else
      curVal = 1;
    end
    curBoutSlice = curT0(bnum):(curT1(bnum)-1);
    overlap = ismember(tSlice,curBoutSlice);
    ltrx(overlap)=curVal;

  end
  labels.ptrx(ltrx>0) = 1;
  labels.ntrx(ltrx<0) = 1;

function relTrax = getRelTrx(handles,curFrame)
  curExp = curFrame.expNum;
  curFly = curFrame.flyNum;
  curTime = curFrame.curTime;
  sz = handles.maxFrames;

  tSlice = curTime-sz:curTime+sz;
  trax = readCache(handles,curExp,curFly,curTime);

  % Rotate
  curA = trax.theta(handles.maxFrames+1)-pi/2;
  trax.theta = trax.theta-curA;

  rotMat = [cos(curA) sin(curA) ; -sin(curA) cos(curA) ];
  trax.X = trax.X - trax.X(sz+1);
  trax.Y = trax.Y - trax.Y(sz+1);
  R = rotMat*[trax.X; trax.Y];
  trax.X = R(1,:); trax.Y = R(2,:);
  trax.X = trax.X + handles.halfSize+1;
  trax.Y = trax.Y + handles.halfSize+1;
  [trax.labels trax.predictions] = getLabels(handles,curFrame);
  relTrax = trax;

function im = getRotatedFrame(handles,curFrame)
% Reads frames around curFrame and rotates them so that the fly is vertical
  sz = handles.halfSize;
  curExp = curFrame.expNum;
  curFly = curFrame.flyNum;
  curTime = curFrame.curTime;
  curTrx = readCache(handles,curExp,curFly,curTime);
  curX = round(curTrx.X(handles.maxFrames+1));
  curY = round(curTrx.Y(handles.maxFrames+1));
  curA = double(curTrx.theta(handles.maxFrames+1));
  tt = zeros(handles.imgY+4*sz,handles.imgX+4*sz,1);
  im = zeros(2*sz+1,2*sz+1,1,2*handles.maxFrames+1);
  bBoxX = (curX-2*sz:curX+2*sz)+2*sz;
  bBoxY = (curY-2*sz:curY+2*sz)+2*sz;

  pointer = [];
  curMovie = handles.JLDobj.GetFile('movie',curExp);
  [pointer.readframe, pointer.nframes, pointer.movie_fid, pointer.movieheaderinfo] = ...
    get_readframe_fcn(curMovie,'interruptible',false);

  for offset = -handles.maxFrames:handles.maxFrames
    if curTime+offset<handles.firstframe{curExp}(curFly) || ...
        curTime+offset>handles.endframe{curExp}(curFly)
      tt(2*sz+1:end-2*sz,2*sz+1:end-2*sz,1) = ...
        zeros(pointer.movieheaderinfo.nr,pointer.movieheaderinfo.nc);
    else
      s = pointer.readframe(curTime+offset);
      tt(2*sz+1:end-2*sz,2*sz+1:end-2*sz,1) = s(:,:,1);
    end
    timg = tt(bBoxY,bBoxX,:);
    rotI = imrotate(timg,curA*180/pi-90,'bilinear','crop');
    im(:,:,:,offset+handles.maxFrames+1) = rotI(sz+1:end-sz,sz+1:end-sz,:);
  end

  if pointer.movie_fid,
    fclose(pointer.movie_fid);
  end

function trx = readCache(handles,expNum,flyNum,curT)

  blk = handles.tsize;
  sz = handles.maxFrames;
  trxOffset = handles.firstframe{expNum}(flyNum);
  blkNumStart = floor( (curT-sz-trxOffset+1)/blk)+1;
  blkNumEnd = ceil( (curT+sz-trxOffset+1)/blk);
  off = curT-trxOffset+1 - (blkNumStart-1)*blk - sz;
  trx.X = []; trx.Y = []; trx.theta = []; trx.maj = []; trx.min = [];
  for curBlk = blkNumStart:blkNumEnd
    idStr = sprintf('%d_%d_%d',expNum,flyNum,curBlk);
    if ~isKey(handles.cache,idStr)
      CacheTracks(handles.output,expNum,flyNum,curT,curT);
      handles = guidata(handles.output);
    end
    curTrx = handles.cache(idStr);
    trx.X = [trx.X curTrx.X];
    trx.Y = [trx.Y curTrx.Y];
    trx.theta = [trx.theta curTrx.theta];
    trx.maj = [trx.maj curTrx.maj];
    trx.min = [trx.min curTrx.min];
  end

  trx.X = trx.X(off:(off+2*sz));
  trx.Y = trx.Y(off:(off+2*sz));
  trx.theta = trx.theta(off:(off+2*sz));
  trx.maj = trx.maj(off:(off+2*sz));
  trx.min = trx.min(off:(off+2*sz));

function axesH = initAxes(handles,ax,sz,row,col)
  axesH.image = [];
  hold(ax,'on');
%   set(ax,'XLimMode','manual');  xlim(ax,[1 2*sz+1]);
%   set(ax,'YLimMode','manual');  ylim(ax,[1 2*sz+1]);

%   set(ax,'XLimMode','manual');  xlim(ax,[1 sz/5*2+1]);
%   set(ax,'YLimMode','manual');  ylim(ax,[1 sz/5*2+1]);
  axesH.image = imagesc(zeros(2*sz+1),'ButtonDownFcn',{@displayOnJLabel,handles,row,col},'Parent',ax,[0,255]);
  axesH.labelNeg = plot(ax,nan,nan,'Linestyle','-','Color',[0.0 0.0 0.7],'Linewidth',3);
  axesH.labelPos = plot(ax,nan,nan,'Linestyle','-','Color',[0.7 0.0 0.0],'Linewidth',3);
  axesH.predNeg = plot(ax,nan,nan,'Linestyle','-','Color',[0.0 0.0 0.7],'Linewidth',3);
  axesH.predPos = plot(ax,nan,nan,'Linestyle','-','Color',[0.7 0.0 0.0],'Linewidth',3);
  axesH.trax = plot(ax,nan,nan,'Linestyle','-','Marker','.',...
    'Color',[0.1 0.1 0.1],'MarkerSize',4,'Linewidth',0.1);
  axesH.fly = plot(ax,nan,nan,'Linestyle','-','Color',[0.7 0.2 0.2]);
  axesH.flyPred = plot(ax,nan,nan,'Linestyle','-','Color',[0.7 0.2 0.2],...
    'Linewidth',2);
  colormap(ax,'gray');
  axis(ax,'off','image');

function tlH = initTimeline(tl,sz,tcurr)
  tlH.image = [];
  hold(tl,'on');
  %set(tl,'XLimMode','auto');  xlim(tl,[1 2*sz+1]);
  %set(tl,'YLimMode','auto');  ylim(tl,[1 2*sz+1]);
  tlH.image = imagesc(zeros(2,2*sz+1,3),'Parent',tl,[0,255]);
  tlH.horig = plot(tl,[tcurr,tcurr],[0.5, 2.5],'g--','HitTest','off','linewidth',2);
  tlH.hcurr = plot(tl,[tcurr,tcurr],[0.5, 2.5],'y-','HitTest','off','linewidth',2);
  colormap(tl,'gray');
  axis(tl,'on','image');
  set(tl,'XTick',[]);
  set(tl,'YTick',[]);

function updatefly(h1,h2,trx,t,displayChoice)
% Coped from JCtrax.

% draw an isosceles triangle with center (x,y)
% with rotation theta
% with height maj*4
% with base min*4

x = trx.X(t);
y = trx.Y(t);
theta = trx.theta(t);
maj = trx.maj(t);
min = trx.min(t);

% isosceles triangle not yet rotated or centered
pts = [-maj*2,-min*2,
       -maj*2,min*2,
       maj*2,0];

% rotate
costheta = cos(theta);
sintheta = sin(theta);
R = [costheta,sintheta;-sintheta,costheta];
pts = pts*R;

% translate
pts(:,1) = pts(:,1) + x;
pts(:,2) = pts(:,2) + y;

switch displayChoice
    case 'automatic'
        if isnan(trx.predictions(t)), colr = [0.2 0.2 0.2];
        elseif trx.predictions(t)<1.5, colr = [0.7 0.0 0.0];
        else colr = [0.0 0.0 0.7];
        end
    case 'manual'
        colr = [0.2 0.2 0.2];
        if ~isnan(trx.labels.ptrx(t)); colr = [0.7 0.0 0.0]; end
        if ~isnan(trx.labels.ntrx(t)); colr = [0.0 0.0 0.7]; end
    otherwise
        colr = [0.2 0.2 0.2];
end

% plot for labels.
set(h1,'xdata',pts([1 2],1),'ydata',pts([1 2],2));
set(h1,'Color',colr);
set(h2,'xdata',pts([2 3 1], 1),'ydata',pts([2 3 1],2));
set(h2,'Color',colr);

% --- Executes on slider movement.
function frameSlider_Callback(hObject, eventdata, handles)
% hObject    handle to frameSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
v = 1-get(hObject,'Value');
handles.frameNo = round(v*(2*handles.maxFrames))+1;
guidata(hObject,handles);
updatePlots(hObject,handles,handles.frameNo);

% --- Executes during object creation, after setting all properties.
function frameSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to frameSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in playButton.
function playButton_Callback(hObject, eventdata, handles)
% hObject    handle to playButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of playButton
if(get(hObject,'Value'));
  handles.isPlaying = true;
  guidata(hObject,handles);
  play(hObject);
else
  handles.isPlaying = false;
  guidata(hObject,handles);
end

% --- Executes when user attempts to close figure_showSimilarFrames.
function figure_showSimilarFrames_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure_showSimilarFrames (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure

handles.JLDobj.frameFig = [];
delete(hObject);

function fps_edit_Callback(hObject, eventdata, handles)
% hObject    handle to fps_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of fps_edit as text
%        str2double(get(hObject,'String')) returns contents of fps_edit as a double
v = str2num(get(hObject,'String'));
if (isempty(v))||(v >handles.maxfps|| v < 0)
    error('frame range requires 1 argument that is a number between 0 and maxfps ');
end

handles.fps = round(v);
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function fps_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to fps_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function frame_range_edit_Callback(hObject, eventdata, handles)
% hObject    handle to frame_range_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of frame_range_edit as text
%        str2double(get(hObject,'String')) returns contents of frame_range_edit as a double

frameRanges = str2double(get(hObject,'String'));
if (isempty(frameRanges))||(frameRanges >handles.maxFrames||frameRanges < 0)
    warnStr = sprintf('Frame range requires 1 argument that is a number between 0 and %d, please reinput a value.\n',handles.maxFrames);
    warndlg(warnStr,'frame range is out of range');
else

    handles.revLimit = round(frameRanges);
    handles.fwdLimit = handles.revLimit;

    guidata(hObject,handles);
end

% --- Executes during object creation, after setting all properties.
function frame_range_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to frame_range_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in timeline_properties.
function timeline_properties_Callback(hObject, eventdata, handles)
% hObject    handle to timeline_properties (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns timeline_properties contents as cell array
%        contents{get(hObject,'Value')} returns selected item from
%        timeline_properties

v = get(hObject,'Value');
s = handles.timeline_prop_options{v};

prop = find(strcmpi(s,handles.JLDobj.allperframefns),1);
handles.perframepropis = prop;
handles.perframeprops = s;

[perframedata,T0,T1] = handles.JLDobj.GetPerFrameData(handles.frames{2}.expNum,handles.frames{2}.flyNum,prop);
TStart = handles.frames{2}.curTime - handles.maxFrames;
IStart = TStart - T0 + 1;
TEnd = handles.frames{2}.curTime + handles.maxFrames;
IEnd = TEnd - T0 + 1;
perframedata = perframedata(IStart:IEnd);
set(handles.htimeline_propdata{2},'XData',-handles.maxFrames:handles.maxFrames,'YData',perframedata);

%reset Ymax and Ymin value for current time and original time
set(handles.hPropOrig,'YData',[0,0]);
set(handles.hPropCurr,'XData',[0,0],'YData',[0,0]);

%for row 1 and 3
for row = [1 3]
    for ndx = 1:handles.numSimilar
        [perframedata,T0,T1] = handles.JLDobj.GetPerFrameData(handles.frames{row}(ndx).expNum,handles.frames{row}(ndx).flyNum,prop);
        TStart = handles.frames{row}(ndx).curTime - handles.maxFrames;
        IStart = TStart - T0 + 1;
        TEnd = handles.frames{row}(ndx).curTime + handles.maxFrames;
        IEnd = TEnd - T0 + 1;
        perframedata = perframedata(IStart:IEnd);
        set(handles.htimeline_propdata{row}(ndx),'XData',-handles.maxFrames:handles.maxFrames,'YData',perframedata);
    end
    fprintf('\n');
end

%plot hcurr for axes_timeline_properties
yLimVal=ylim(handles.axes_timeline_properties);
frameNo = handles.frameNo - handles.maxFrames - 1;
set(handles.hPropOrig,'YData',[yLimVal(1),yLimVal(2)]);
set(handles.hPropCurr,'XData',[frameNo,frameNo],'YData',[yLimVal(1),yLimVal(2)]);

%update the text_timeline_prop value
curExperiment = handles.frames{handles.cur_rb_row}(handles.cur_rb_col).expNum;
curFly = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).flyNum;
curTime = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).curTime+handles.frameNo-handles.maxFrames-1;
curProp= handles.perframepropis;
perframedata = handles.JLDobj.GetPerFrameData1(curExperiment,curFly,curProp,curTime);
s = sprintf('%.3f',perframedata);
set(handles.text_timeline_prop,'String',s);

guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function timeline_properties_CreateFcn(hObject, eventdata, handles)
% hObject    handle to timeline_properties (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes when selected object is changed in uipanel6.
function uipanel6_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in uipanel6
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)
switch get(eventdata.NewValue,'Tag') % Get Tag of selected object.
    case 'rb_automatic'
        for row = 1:3
            for ndx = 1:handles.numSimilar

                % Code for when rb_top1 is selected.
                handles.display = 'automatic';
                set(handles.hAxes(row,ndx).predPos,'visible', 'on');
                set(handles.hAxes(row,ndx).predNeg,'visible', 'on');
                set(handles.hAxes(row,ndx).labelPos,'visible', 'off');
                set(handles.hAxes(row,ndx).labelNeg,'visible', 'off');

            end
        end

    case 'rb_manual'
        for row = 1:3
            for ndx = 1:handles.numSimilar
                handles.display = 'manual';
                set(handles.hAxes(row,ndx).labelPos,'visible', 'on');
                set(handles.hAxes(row,ndx).labelNeg,'visible', 'on');
                set(handles.hAxes(row,ndx).predPos,'visible', 'off');
                set(handles.hAxes(row,ndx).predNeg,'visible', 'off');
            end
        end
    otherwise
        % Code for when there is no match.
end
updatePlots(hObject,handles,handles.frameNo);
guidata(hObject,handles);

function displayOnJLabel(obj,eventdata,handles,row,col)

switch get(gcf,'selectiontype')
    case 'normal'
        % LEFT CLICK
    case 'extend'
        % SHIFT-CLICK LEFT (or L/R simultaneous)
    case 'alt'
        % CTRL-CLICK LEFT (or RIGHT-CLICK)
    case 'open'
        % DOUBLE-CLICK

        handles=guidata(obj);
        JLabelHandles = handles.JLabelHandles;
        figure(JLabelHandles.figure_JLabel);
        if row==2
            JLabelHandles.guidata.ts(1) =999; % in order to update plot in SetCurrentFrame
            handles.curExp = handles.frames{2}.expNum;
            [JLabelHandles,~] = JLabel('SetCurrentMovie',JLabelHandles,handles.curExp);
            handles.curFly = handles.frames{2}.flyNum;
            JLabelHandles = JLabel('SetCurrentFlies',JLabelHandles,handles.curFly);
            handles.curFrame = handles.frames{2}.curTime;
            JLabelHandles = JLabel('SetCurrentFrame',JLabelHandles,1,handles.curFrame,obj);
            guidata(obj, handles);
        else
            handles.curExp = handles.frames{row}(col).expNum;
            [JLabelHandles,~] = JLabel('SetCurrentMovie',JLabelHandles,handles.curExp);
            handles.curFly = handles.frames{row}(col).flyNum;
            JLabelHandles = JLabel('SetCurrentFlies',JLabelHandles,handles.curFly);
            handles.curFrame = handles.frames{row}(col).curTime;
            JLabelHandles = JLabel('SetCurrentFrame',JLabelHandles,1,handles.curFrame,obj);
            guidata(obj, handles);
        end

end

% --- Executes when selected object is changed in chooseFrame.
function chooseFrame_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in chooseFrame
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)

old_rb_row = handles.cur_rb_row;
old_rb_col = handles.cur_rb_col;

switch get(eventdata.NewValue,'Tag') % Get Tag of selected object.
    case 'rb_top1'
        % Code for when rb_top1 is selected.
        handles.cur_rb_row = 1;
        handles.cur_rb_col = 1;
    case 'rb_top2'
        % Code for when radiobutton2 is selected.
        handles.cur_rb_row = 1;
        handles.cur_rb_col = 2;
    case 'rb_top3'
        handles.cur_rb_row = 1;
        handles.cur_rb_col = 3;

    case 'rb_top4'
        handles.cur_rb_row = 1;
        handles.cur_rb_col = 4;

    case 'rb_middle1'
        handles.cur_rb_row = 2;
        handles.cur_rb_col = 1;
    case 'rb_bottom1'
        handles.cur_rb_row = 3;
        handles.cur_rb_col = 1;
    case 'rb_bottom2'
        handles.cur_rb_row = 3;
        handles.cur_rb_col = 2;

    case 'rb_bottom3'
        handles.cur_rb_row = 3;
        handles.cur_rb_col = 3;
    case 'rb_bottom4'
        handles.cur_rb_row = 3;
        handles.cur_rb_col = 4;

    otherwise
        % Code for when there is no match.
end

%change the color of the old lines to blue
set(handles.htimeline_propdata{old_rb_row}(old_rb_col),'LineWidth',0.5);

%change the color of the new line to red
set(handles.htimeline_propdata{handles.cur_rb_row}(handles.cur_rb_col),'LineWidth',2);

%update the text_timeline_prop value
curExperiment = handles.frames{handles.cur_rb_row}(handles.cur_rb_col).expNum;
curFly = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).flyNum;
curTime = handles.frames{ handles.cur_rb_row}(handles.cur_rb_col).curTime+handles.frameNo-handles.maxFrames-1;
curProp= handles.perframepropis;
perframedata = handles.JLDobj.GetPerFrameData1(curExperiment,curFly,curProp,curTime);
s = sprintf('%.3f',perframedata);
set(handles.text_timeline_prop,'String',s);

guidata(hObject,handles);

function plotPredictionsInAutoMode(handles)
%for row 2

relTrx = getRelTrx(handles,handles.frames{2});
[labels predictions] = getLabels(handles,handles.frames{2});
ptrx = nan(size(predictions));
ntrx = nan(size(predictions));

for i=1:length(ptrx)
    if predictions(i)<1.5
        ptrx(i) = 1;
    else
        ntrx(i) = 1;
    end
end

%To fill in the gap between the lablePos and labelNeg
diffPred = diff(predictions);
changeIdex = find(diffPred);

for i = 1:length(changeIdex)
    if diffPred(changeIdex(i))<0 %negative -> positive
        ntrx(changeIdex(i)+1) = 1;
    else %positive -> negative
        ptrx(changeIdex(i)+1) = 1;
    end
end

for ndx = 1:handles.numSimilar
    handles.trx{2,ndx} = relTrx;
    set(handles.hAxes(2,ndx).predPos,'XData',relTrx.X.*ptrx,'Ydata',relTrx.Y.*ptrx);
    set(handles.hAxes(2,ndx).predNeg,'XData',relTrx.X.*ntrx,'Ydata',relTrx.Y.*ntrx);
end

%for row 1 and 3
for row = [1 3]
    for ndx = 1:handles.numSimilar
        relTrx = getRelTrx(handles,handles.frames{row}(ndx));
        handles.trx{row,ndx} = relTrx;
        [labels predictions] = getLabels(handles,handles.frames{row}(ndx));

        ptrx = nan(size(predictions));
        ntrx = nan(size(predictions));

        for i=1:length(ptrx)
            if predictions(i)<1.5
                ptrx(i) = 1;
            else
                ntrx(i) = 1;
            end
        end

        diffPred = diff(predictions);
        changeIdex = find(diffPred);

        for i = 1:length(changeIdex)
            if diffPred(changeIdex(i))<0 %negative -> positive
                ntrx(changeIdex(i)+1) = 1;
            else %positive -> negative
                ptrx(changeIdex(i)+1) = 1;
            end
        end

        set(handles.hAxes(row,ndx).predPos,'XData',relTrx.X.*ptrx,'YData',relTrx.Y.*ptrx);
        set(handles.hAxes(row,ndx).predNeg,'XData',relTrx.X.*ntrx,'YData',relTrx.Y.*ntrx);
    end
end

function poltLabelsInManualMode(handles)

%for row 2
relTrx = getRelTrx(handles,handles.frames{2});
[labels predictions] = getLabels(handles,handles.frames{2});
ptrx = labels.ptrx;
ntrx = labels.ntrx;

for ndx = 1:handles.numSimilar
    handles.trx{2,ndx} = relTrx;
    set(handles.hAxes(2,ndx).labelPos,'XData',relTrx.X.*ptrx,'Ydata',relTrx.Y.*ptrx);
    set(handles.hAxes(2,ndx).labelNeg,'XData',relTrx.X.*ntrx,'Ydata',relTrx.Y.*ntrx);
end

%for row 1 and 3
for row = [1 3]
    for ndx = 1:handles.numSimilar
        relTrx = getRelTrx(handles,handles.frames{row}(ndx));
        handles.trx{row,ndx} = relTrx;
        [labels predictions] = getLabels(handles,handles.frames{row}(ndx));
        ptrx = labels.ptrx;
        ntrx = labels.ntrx;

        set(handles.hAxes(row,ndx).labelPos,'XData',relTrx.X.*ptrx,'YData',relTrx.Y.*ptrx);
        set(handles.hAxes(row,ndx).labelNeg,'XData',relTrx.X.*ntrx,'YData',relTrx.Y.*ntrx);
    end
end

% --- Executes on key press with focus on figure_showSimilarFrames and none of its controls.
function figure_showSimilarFrames_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to figure_showSimilarFrames (see GCBO)
% eventdata  structure with the following fields (see FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
switch eventdata.Key,

  case 'leftarrow',
      go_previous_frame(hObject, handles);

  case 'rightarrow',
      go_next_frame(hObject, handles);

end

function go_next_frame(hObject, handles)
% hObject    handle to menu_go_next_frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% TODO: make this work with multiple preview axes

frameNo = handles.frameNo;
if(frameNo-handles.maxFrames>handles.fwdLimit)
    frameNo = handles.maxFrames+1-handles.revLimit;
else
    frameNo = frameNo+1;
end

handles.frameNo = frameNo;

updatePlots(hObject,handles,frameNo);
set(handles.frameSlider,'Value', 1-frameNo/(2*handles.maxFrames+1));

guidata(hObject,handles);

% --------------------------------------------------------------------
function go_previous_frame(hObject, handles)
% hObject    handle to menu_go_previous_frame (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

frameNo = handles.frameNo;
if(handles.maxFrames - frameNo>= handles.revLimit - 1)
    frameNo = handles.maxFrames + handles.fwdLimit + 1;
else
    frameNo = frameNo-1;
end

handles.frameNo = frameNo;

updatePlots(hObject,handles,frameNo);
set(handles.frameSlider,'Value', 1-frameNo/(2*handles.maxFrames+1));

guidata(hObject,handles);

% --- Executes on mouse press over figure background.
function figure_showSimilarFrames_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to figure_showSimilarFrames (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function zoomOut_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to zoomOut (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
oldLimVal = get(handles.middle1,'xLim');
mag = (oldLimVal(2) - oldLimVal(1))/(handles.halfSize*0.4);

if mag<=4
    mag = mag +1;
else
    mag=5;
end
limVal = [handles.halfSize*(1-mag/5), handles.halfSize*(1+mag/5)];
set(handles.middle1,'xLim',limVal,'yLim',limVal);

% --------------------------------------------------------------------
function zoomIn_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to zoomIn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
oldLimVal = get(handles.middle1,'xLim');
mag = (oldLimVal(2) - oldLimVal(1))/(handles.halfSize*0.4);

if mag>=2
    mag = mag-1;
else
    mag=1;
end
limVal = [handles.halfSize*(1-mag/5), handles.halfSize*(1+mag/5)];
set(handles.middle1,'xLim',limVal,'yLim',limVal);

author: kristinbranson

repo: JAABA


text:
% Laser sensor's parameters
function lidar = SetLidarParameters()

lidar.angle_min = -2.351831;
lidar.angle_max =  2.351831;
lidar.angle_increment = 0.004363;
lidar.npoints   = 1079;
lidar.range_min = 0.023;
lidar.range_max = 60;
lidar.scan_time = 0.025;
lidar.time_increment  = 1.736112e-05;
lidar.angles = (lidar.angle_min : lidar.angle_increment : lidar.angle_max)';

author: meyiao

repo: LaserSLAM

save_path: github-repos/MATLAB/meyiao-LaserSLAM


text:
function mapVals = rmOverlapMap(model, mask, X, Y);
% Produce a single pRF / stimulus overlap mapVals given some low-level inputs.
%
%  mapVals = rmOverlapMap(model, mask, X, Y);
%
% model should be a pRF / RM Model struct (or a path to an RM file).
%
% mask should be a single, 2D, binary stimulus mask, sampled at the
% retinotopic points represented in X and Y.
%
% X and Y can either be the range of x and y points at which the mask was
% sampled, or matrices (e.g. produced my MESHGRID) matching the size of
% mask. Direction conventions: +x is right, -x is left; +y is up, -y is
% down.
%
% Returns a vector of mapVals values representing the proportion of each
% voxel's pRF which is contained within the mask.
%
% Currently implemented only for Volume/Gray views models.
%
% SEE ALSO: rmComputeOverlapMaps.
%
% ras, 11/2007.
if nargin < 4, error('Need all 4 input arguments.');	end
if iscell(model),	model = model{1};		end

if ~isequal( size(mask), size(X) ) | ~isequal( size(mask), size(X) )
	% assume x and y ranges were input
	xRange = unique(X);
	yRange = unique(Y);
	[X Y] = meshgrid(xRange, yRange);
end

% get indices in the mask
I = find(mask);

%% main loop, create pRF for each voxel, convolve w/ mask
% Doing this one voxel at a time takes a LONG time (hours).
% But trying to make all the pRFs at once may cause an
% out-of-memory error. So, we split the difference: do it
% in batches of size n (chosen to work on my ~2005 Windows box;
% adjust as needed).
n = 50;  % compute n pRFs at once
nVoxels = length(model.x0);
verbose = prefsVerboseCheck;
if verbose
	hwait = mrvWaitbar(0, 'Computing pRF Overlap Map...');
end

% go
for v = 1:n:nVoxels
	rng = [0:n-1] + v;  % take this bunch of voxels first
	rng = rng( rng < nVoxels );  % out-of-range check

	% grab params for this batch
	x0 = model.x0(rng);
	y0 = model.y0(rng);
	sigma = model.sigma.major(rng);

	% make pRFs for this batch
	pRFs = rfGaussian2D(X, Y, sigma, sigma, 0, x0, y0);

	% loop across pRFs, convolving each with the mask
	for ii = 1:length(rng)
		mapVals(rng(ii)) = sum(pRFs(I,ii)) ./ sum(pRFs(:,ii));
	end

	if verbose, mrvWaitbar(v/nVoxels, hwait); end
end

if verbose,	close(hwait); end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [grandavg] = ft_freqgrandaverage(cfg, varargin)

% FT_FREQGRANDAVERAGE computes the average powerspectrum or time-frequency spectrum
% over multiple subjects
%
% Use as
%   [grandavg] = ft_freqgrandaverage(cfg, freq1, freq2, freq3...)
%
% The input data freq1..N are obtained from either FT_FREQANALYSIS with
% keeptrials=no or from FT_FREQDESCRIPTIVES. The configuration structure
% can contain
%   cfg.keepindividual = 'yes' or 'no' (default = 'no')
%   cfg.foilim         = [fmin fmax] or 'all', to specify a subset of frequencies (default = 'all')
%   cfg.toilim         = [tmin tmax] or 'all', to specify a subset of latencies (default = 'all')
%   cfg.channel        = Nx1 cell-array with selection of channels (default = 'all'),
%                        see FT_CHANNELSELECTION for details
%   cfg.parameter      = string or cell-array of strings indicating which
%                        parameter(s) to average. default is set to
%                        'powspctrm', if it is present in the data.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure. For this particular function, the input should be
% specified as a cell-array.
%
% See also FT_TIMELOCKGRANDAVERAGE, FT_FREQANALYSIS, FT_FREQDESCRIPTIVES,
% FT_FREQBASELINE

% FIXME averaging coherence is not possible if inputs contain different amounts of data (i.e. chan/freq/time)

% Copyright (C) 2005-2006, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
for i=1:length(varargin)
  varargin{i} = ft_checkdata(varargin{i}, 'datatype', 'freq', 'feedback', 'no');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.keepindividual = ft_getopt(cfg, 'keepindividual', 'no');
cfg.channel        = ft_getopt(cfg, 'channel',    'all');
cfg.foilim         = ft_getopt(cfg, 'foilim',     'all');
cfg.toilim         = ft_getopt(cfg, 'toilim',     'all');
cfg.parameter      = ft_getopt(cfg, 'parameter',  []);

if isempty(cfg.parameter) && isfield(varargin{1}, 'powspctrm')
  cfg.parameter = 'powspctrm';
elseif isempty(cfg.parameter)
  ft_error('you should specify a valid parameter to average');
end

if ischar(cfg.parameter)
  cfg.parameter = {cfg.parameter};
end

Nsubj    = length(varargin);
dimord   = varargin{1}.dimord;
hasfreq  = ~isempty(strfind(varargin{i}.dimord, 'freq')); % this should always be true
hastime  = ~isempty(strfind(varargin{i}.dimord, 'time'));
hasrpt   = ~isempty(strfind(varargin{i}.dimord, 'rpt'));
hastap   = ~isempty(strfind(varargin{i}.dimord, 'tap'));

% check whether the input data is suitable
if hasrpt
  ft_error('the input data of each subject should be an average, use FT_FREQDESCRIPTIVES first');
end
if hastap
  ft_error('multiple tapers in the input are not supported');
end

if ischar(cfg.foilim) && strcmp(cfg.foilim, 'all')
  fbeg = -inf;
  fend =  inf;
else
  fbeg = cfg.foilim(1);
  fend = cfg.foilim(2);
end

if ischar(cfg.toilim) && strcmp(cfg.toilim, 'all')
  tbeg = -inf;
  tend =  inf;
else
  tbeg = cfg.toilim(1);
  tend = cfg.toilim(2);
end

% select the data in all inputs
for k=1:numel(cfg.parameter)

  % determine which channels, frequencies and latencies are available for all inputs
  for i=1:Nsubj
    cfg.channel = ft_channelselection(cfg.channel, varargin{i}.label);
    if hasfreq
      fbeg = max(fbeg, varargin{i}.freq(1  ));
      fend = min(fend, varargin{i}.freq(end));
    end
    if hastime
      tbeg = max(tbeg, varargin{i}.time(1  ));
      tend = min(tend, varargin{i}.time(end));
    end
  end
  cfg.foilim = [fbeg fend];
  cfg.toilim = [tbeg tend];

  % pick the selections
  for i=1:Nsubj
    if ~isfield(varargin{i}, cfg.parameter{k})
      ft_error('the field %s is not present in data structure %d', cfg.parameter{k}, i);
    end
    [dum, chansel] = match_str(cfg.channel, varargin{i}.label);
    varargin{i}.label = varargin{i}.label(chansel);

    if hasfreq
      freqsel = nearest(varargin{i}.freq, fbeg):nearest(varargin{i}.freq, fend);
      varargin{i}.freq = varargin{i}.freq(freqsel);
    end
    if hastime
      timesel = nearest(varargin{i}.time, tbeg):nearest(varargin{i}.time, tend);
      varargin{i}.time = varargin{i}.time(timesel);
    end
    % select the overlapping samples in the power spectrum
    switch dimord
      case 'chan_freq'
        varargin{i}.(cfg.parameter{k}) = varargin{i}.(cfg.parameter{k})(chansel,freqsel);
      case 'chan_freq_time'
        varargin{i}.(cfg.parameter{k}) = varargin{i}.(cfg.parameter{k})(chansel,freqsel,timesel);
      case {'rpt_chan_freq' 'rpttap_chan_freq' 'subj_chan_freq'}
        varargin{i}.(cfg.parameter{k}) = varargin{i}.(cfg.parameter{k})(:,chansel,freqsel);
      case {'rpt_chan_freq_time' 'rpttap_chan_freq_time' 'subj_chan_freq_time'}
        varargin{i}.(cfg.parameter{k}) = varargin{i}.(cfg.parameter{k})(:,chansel,freqsel,timesel);
      otherwise
        ft_error('unsupported dimord');
    end
  end % for i = subject
end % for k = parameter

% determine the size of the data to be averaged
dim = cell(1,numel(cfg.parameter));
for k=1:numel(cfg.parameter)
  dim{k} = size(varargin{1}.(cfg.parameter{k}));
end

% give some feedback on the screen
if strcmp(cfg.keepindividual, 'no')
  for k=1:numel(cfg.parameter)
    ft_info('computing average %s over %d subjects\n', cfg.parameter{k}, Nsubj);
  end
else
  for k=1:numel(cfg.parameter)
    ft_info('not computing average, but keeping individual %s for %d subjects\n', cfg.parameter{k}, Nsubj);
  end
end

% allocate memory to hold the data and collect it
for k=1:numel(cfg.parameter)
  if strcmp(cfg.keepindividual, 'no')
    tmp = zeros(dim{k});
    for s=1:Nsubj
      tmp = tmp + varargin{s}.(cfg.parameter{k})./Nsubj; % do a weighted running sum
    end
  elseif strcmp(cfg.keepindividual, 'yes')
    tmp = zeros([Nsubj dim{k}]);
    for s=1:Nsubj
      tmp(s,:,:,:,:) = varargin{s}.(cfg.parameter{k});
    end
  end
  grandavg.(cfg.parameter{k}) = tmp;
end

% collect the output data
grandavg.label = varargin{1}.label;
grandavg.freq  = varargin{1}.freq;
if isfield(varargin{1}, 'time')
  % remember the time axis
  grandavg.time = varargin{1}.time;
end
if isfield(varargin{1}, 'labelcmb')
  grandavg.labelcmb = varargin{1}.labelcmb;
end
if isfield(varargin{1}, 'grad')
  ft_warning('discarding gradiometer information because it cannot be averaged');
end
if isfield(varargin{1}, 'elec')
  ft_warning('discarding electrode information because it cannot be averaged');
end
if strcmp(cfg.keepindividual, 'yes')
  grandavg.dimord = ['subj_',varargin{1}.dimord];
elseif strcmp(cfg.keepindividual, 'no')
  grandavg.dimord = varargin{1}.dimord;
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   varargin
ft_postamble provenance grandavg
ft_postamble history    grandavg
ft_postamble savevar    grandavg

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function chrpak_test057 ( )

%*****************************************************************************80
%
%% TEST057 tests OCT_TO_I4 and I4_TO_OCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2007
%
%  Author:
%
%    John Burkardt
%
  test_num = 3;

  i4_test = [ 21, -32, 1776 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST057\n' );
  fprintf ( 1, '  OCT_TO_I4, octal->integer.\n' );
  fprintf ( 1, '  I4_TO_OCT, integer->octal\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        I4       OCT        I4\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    i4 = i4_test(test);
    oct = i4_to_oct ( i4 );
    j4 = oct_to_i4 ( oct );

    fprintf ( 1, '  %8d  %8s  %8d\n', i4, oct, j4 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function wseg = pb2wseg(pb, maxsp)

nsp = Inf;
c = 1;
pb = max(pb, [], 3);
while nsp > maxsp
    if c > 1
        wseg = watershed(medfilt2(pb, [c c]));
    else
        wseg = watershed(pb);
    end
    nsp = max(wseg(:));
    c = c + 2;
end
wseg = uint16(wseg);

author: Cloud-CV

repo: object-proposals


text:
function [C,Atp,Btp]=etprod(Cidx,A,Aidx,B,Bidx)
% tprod wrapper to make calls more similar to Einstein Summation Convention
%
% [C,Atp,Btp]=etprod(Cidx,A,Aidx,B,Bidx);
% Wrapper function for tprod to map between Einstein summation
% convetion (ESC) and tprod's numeric calling convention e.g.
% 1) Matrix Matrix product:
%   ESC: C_ij = A_ik B_kj <=> C = etprod('ij',A,'ik',B,'kj');
% 2) Vector outer product
%   ESC: C_ij = A_i B_j   <=> C = etprod('ij',A,'i',B,'j');   % A,B col vec
%                             C = etprod('ij',A,' i',B,' j'); % A,B row vec
% N.B. use spaces ' ' to indicate empty/ignored *singlenton* dimensions
% 3) Matrix vector product
%   ESC: C_i  = A_ik B_k  <=> C = etprod('i',A,'ik',B,'k');
% 4) Spatial Filtering
%   ESC: FX_fte = A_cte B_cf <=> C = etprod('fte',A,'cte',B,'cf')
% OR:
%   C = etprod({'feat','time','epoch'},A,{'ch','time','epoch'},B,{'ch','feat'})
%
% Inputs:
%  Cidx -- the list of dimension labels for the output
%  A    -- [n-d] array of the A values
%  Aidx -- [ndims(A) x 1] (array, string, or cell array of strings)
%          list of dimension labels for A array
%  B    -- [m-d] array of the B values
%  Bidx -- [ndims(B) x 1] (array, string or cell array of strings)
%          list of dimension labels for B array
% Outputs:
%  C    -- [p-d] array of output values.  Dimension labels are as in Cidx
%  Atp  -- [ndims(A) x 1] A's dimspec as used in the core tprod call
%  Btp  -- [ndims(B) x 1] B's dimspec as used in the core tprod call
%
% See Also: tprod, tprod_testcases
%
% Copyright 2006-     by Jason D.R. Farquhar (jdrf@zepler.org)
% Permission is granted for anyone to copy, use, or modify this
% software and accompanying documents for any uncommercial
% purposes, provided this copyright notice is retained, and note is
% made of any changes that have been made. This software and
% documents are distributed without any warranty, express or
% implied
if ( iscell(Aidx)~=iscell(Bidx) || iscell(Cidx)~=iscell(Aidx) )
 error('Aidx,Bidx and Cidx cannot be of different types, all cells or arrays');
end

Atp = zeros(size(Aidx));
Btp = zeros(size(Bidx));
% Map inner product dimensions, to unique *negative* index
for i=1:numel(Aidx)
   if ( iscell(Aidx) ) Bmatch = strcmp(Aidx{i}, Bidx);
   else                Bmatch =       (Aidx(i)==Bidx);
   end
   if ( any(Bmatch) ) Btp(Bmatch)=-i; Atp(i)=-i; end;
end
% Spaces/empty values in the input become 0's, i.e. ignored dimensions
if ( iscell(Aidx) )
   Btp(strcmp(' ',Bidx))=0;Btp(strcmp('',Bidx))=0;
   Atp(strcmp(' ',Aidx))=0;Atp(strcmp('',Aidx))=0;
else
   Btp(' '==Bidx)=0;
   Atp(' '==Aidx)=0;
end

% Map to output position numbers, to correct *positive* index
for i=1:numel(Cidx);
   if ( iscell(Aidx) )
      Atp(strcmp(Cidx{i}, Aidx))=i;
      Btp(strcmp(Cidx{i}, Bidx))=i;
   else
      Atp(Cidx(i)==Aidx)=i;
      Btp(Cidx(i)==Bidx)=i;
   end
end

% now do the tprod call.
global LOG; if ( isempty(LOG) ) LOG=0; end; % N.B. set LOG to valid fd to log
if ( LOG>0 )
   fprintf(LOG,'tprod(%s,[%s], %s,[%s])\n',mxPrint(A),sprintf('%d ',Atp),mxPrint(B),sprintf('%d ',Btp));
end
C=tprod(A,Atp,B,Btp,'n');
return;

function [str]=mxPrint(mx)
sz=size(mx);
if ( isa(mx,'double') ) str='d'; else  str='s'; end;
if ( isreal(mx)) str=[str 'r']; else str=[str 'c']; end;
str=[str ' [' sprintf('%dx',sz(1:end-1)) sprintf('%d',sz(end)) '] '];
return;

%----------------------------------------------------------------------------
function testCase();
A=randn(10,10); B=randn(10,10);
C2 = tprod(A,[1 -2],B,[-2 2]);
C  = etprod('ij',A,'ik',B,'kj'); mad(C2,C)
C  = etprod({'i' 'j'},A,{'i' 'k'},B,{'k' 'j'});
A=randn(100,100);B=randn(100,100,4);
C3 = tprod(A,[-1 -2],B,[-1 -2 3]);
C3 = tprod(B,[-1 -2 1],A,[-1 -2]);
C3 = tprod(A,[-1 -2],B,[-1 -2 1]);
C  = etprod('3',A,'12',B,'123');
C  = etprod([3],A,[1 2],B,[1 2 3]);
C  = etprod({'3'},A,{'1' '2'},B,{'1' '2' '3'})


text:
function [ptheta] = tapas_sem_prosa_ptheta()
%% Returns the standard priors of the model.
%
% Input
%
% Output
% ptheta -- Structure containing the priors. The prior distribution is assumed
%           to be log Gaussian, so that the prior are the means and covariance
%           matrix. It is assumed that the covariance is diagonal so only the
%           eigenvalues are returned. ptheta.jm is a projection matrix. It can
%           be replaced with a rank deficient matrix in order to project the
%           samples to a lower dimensional space.

%
% aponteeduardo@gmail.com
% copyright (C) 2015
%

DIM_THETA = tapas_sem_prosa_ndims();
[ptheta] = tapas_sem_prosa_gaussian_priors();

% Projection matrix
ptheta.jm = eye(DIM_THETA);

% Likelihood function and priors

ptheta.name = 'prosa';
ptheta.llh = [];
ptheta.lpp = @tapas_sem_prosa_lpp;
ptheta.prepare = @tapas_sem_prepare_gaussian_ptheta;
ptheta.sample_priors = @tapas_sem_sample_gaussian_uniform_priors;
ptheta.ndims = tapas_sem_prosa_ndims();
ptheta.npars = 2; % It has two sets of parameters.

end

author: translationalneuromodeling

repo: tapas


text:
function filters = make_filters(radii, gtheta)

d = 2;

filters = cell(numel(radii), numel(gtheta));
for r = 1:numel(radii),
    for t = 1:numel(gtheta),

        ra = radii(r);
        rb = ra / 4;
        theta = gtheta(t);

        ra = max(1.5, ra);
        rb = max(1.5, rb);
        ira2 = 1 / ra^2;
        irb2 = 1 / rb^2;
        wr = floor(max(ra, rb));
        wd = 2*wr+1;
        sint = sin(theta);
        cost = cos(theta);

        % 1. compute linear filters for coefficients
        % (a) compute inverse of least-squares problem matrix
        filt = zeros(wd,wd,d+1);
        xx = zeros(2*d+1,1);
        for u = -wr:wr,
            for v = -wr:wr,
                ai = -u*sint + v*cost; % distance along major axis
                bi = u*cost + v*sint; % distance along minor axis
                if ai*ai*ira2 + bi*bi*irb2 > 1, continue; end % outside support
                xx = xx + cumprod([1;ai+zeros(2*d,1)]);
            end
        end
        A = zeros(d+1,d+1);
        for i = 1:d+1,
            A(:,i) = xx(i:i+d);
        end

        % (b) solve least-squares problem for delta function at each pixel
        for u = -wr:wr,
            for v = -wr:wr,
                ai = -u*sint + v*cost; % distance along major axis
                bi = u*cost + v*sint; % distance along minor axis
                if (ai*ai*ira2 + bi*bi*irb2) > 1, continue; end % outside support
                yy = cumprod([1;ai+zeros(d,1)]);
                filt(v+wr+1,u+wr+1,:) = A\yy;
            end
        end

        filters{r,t}=filt;
    end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function slopeplot(f,xs,X,Xplot,n)
%SLOPEPLOT    Plot of one-dimensional function with slope expansion
%
%  slopeplot(f,xs,X,Xplot,n)
%
%Simple functions, which can be written in one string, can be entered
%  directly. The unknown must be 'x'. E.g.,
%
%  slopeplot('sqrt(abs(x))',2,infsup(-1,1));            or
%  slopeplot('sin(x)',0,infsup(-1,1));                  or
%  slopeplot('x.*x',1,infsup(-2,3));
%
%For a one-dimensional function f to be called by u=f(x) consisting of
%  operators overloaded by slope operators such that for a column vector
%  v, f(v) computes the vector f(v(i)), i=1:length(v).
%  Then the function f is expanded within X with respect to xs and
%  plotted together with the slopes in the interval Xplot).
%
%Parameter Xplot is optional; default is  hull(xs,X) * midrad(1,0.2);
%Parameter n is optional (number of grid points for function plot); default is 100
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if prod(size(xs))~=1
    error('slopeplot only for one-dimensional functions')
  end

  if nargin==3 | isempty(Xplot)
    Xplot = hull(xs,X);
    Xplot = Xplot + midrad(0,0.2*diam(Xplot));
  end

  if nargin<=4
    n = 100;
  end

  % slope evaluation
  x = slopeinit(xs,X);
  eval([ 'u = ' f ';' ])

  % plot function
  v = linspace(Xplot.inf,Xplot.sup,n);
  x = v;
  eval([ 'y = ' f ';' ])
  plot(x,y,'k')
  hold on

  % plot expansion range X
  A = axis;
  dy = A(4)-A(3);
  axis([ A(1:2) A(3)-.1*dy A(4)+.1*dy ])
  H = (A(4)-A(3))/200;         % height of box for X
  infX = inf(X);
  supX = sup(X);
  XX = [ infX supX supX infX ];
  patch(XX,A(3)+[-H -H H H],'b')

  % plot expansion point xs
  if isa(xs,'intval')
    infxs = inf(xs);
    supxs = sup(xs);
  else
    infxs = xs;
    supxs = xs;
  end
  xxs = [ infxs supxs supxs infxs ];
  patch(xxs,A(3)+2*[-H -H H H],'r')
  plot([infX infxs supxs],A(3)*ones(1,3),'b')

  % lines infX,supX,xs to function
  x = infX; eval(['y = ' f ';']); plot([infX infX],[A(3) y],'b:')
  x = supX; eval(['y = ' f ';']); plot([supX supX],[A(3) y],'b:')
  x = infxs; eval(['y = ' f ';']); plot([infxs infxs],[A(3) y],'r:')
  x = supxs; eval(['y = ' f ';']); plot([supxs supxs],[A(3) y],'r:')
  x = infxs; eval(['y = ' f ';']); plot(infxs,y,'r.')
  x = supxs; eval(['y = ' f ';']); plot(supxs,y,'r.')

  % slopes
  plot( v , inf(u.c) + inf(u.s)*(v-infxs) , 'b-' );
  plot( v , inf(u.c) + sup(u.s)*(v-infxs) , 'b-' );
  plot( v , sup(u.c) + inf(u.s)*(v-supxs) , 'b-' );
  plot( v , sup(u.c) + sup(u.s)*(v-supxs) , 'b-' );

  % add text X and xs
  if infxs+.5*(supxs-infxs)<infX+.5*(supX-infX)
    set(gcf,'DefaultTextColor','blue')
    text(supX,A(3)-8*H,'\bf X')
    set(gcf,'DefaultTextColor','red')
    text(infxs,A(3)-8*H,'\bf xs')
  else
    set(gcf,'DefaultTextColor','blue')
    text(infX,A(3)-8*H,'\bf X')
    set(gcf,'DefaultTextColor','red')
    text(supxs,A(3)-8*H,'\bf xs')
  end

  % add title
  set(gcf,'DefaultTextColor','default')
  strX = infsup(X);
  strX(isspace(strX)) = '';
  if isa(xs,'intval')
    strxs = infsup(xs);
  else
    strxs = sprintf('%g',xs);
  end
  strxs(isspace(strxs)) = '';
  title(['slope expansion of ' '' f '' ' in X=' strX ...
         ' with respect to xs=' strxs ])

  hold off

  setround(rndold)

repo: OpenMAS


text:
function ch = rocholhFactorise(v);

% ROCHOLHFACTORISE Rank one Cholesky factorise.

% ROCHOL

% Return structure which represents factorisation of I+vv^T

if size(v, 2) ~= 1
  error('v should be a row vector');
end
t = 0;
ch.n = size(v, 1);
ch.v = v;
ch.s = zeros(size(v));
ch.u = zeros(size(v));
t = ch.v(1)*ch.v(1);
ch.s(1) = inf;
ch.u(1) = 1/ch.v(1);
for i = 2:ch.n
  tnew = t + ch.v(i)*ch.v(i);
  ch.s(i) = sqrt(tnew/t);
  ch.u(i) = ch.v(i)/tnew;
  t = tnew;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%ESTIMATERIGIDTRANSFORM  Computes an optimal affine transformation between two 2D point sets
%
%     M = cv.estimateRigidTransform(src, dst)
%     M = cv.estimateRigidTransform(src, dst, 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ First input 2D point set stored in a cell array of 2-element
%   vectors `{[x,y], ...}`, or first image (8-bit numeric array, 1- or
%   3-channels).
% * __dst__ Second input 2D point set of the same size and the same type as
%   `src`, or second image (8-bit numeric array, 1- or 3-channels).
%
% ## Output
% * __M__ output 2x3 affine transformation `[A|b]` matrix (see below).
%
% ## Options
% * __FullAffine__ If true, the function finds an optimal affine transformation
%   with no additional resrictions (6 degrees of freedom). Otherwise, the
%   class of transformations to choose from is limited to combinations of
%   translation, rotation, and uniform scaling (4 degrees of freedom).
%   default false
%
% The function finds an optimal affine transform `[A|b]` (a 2x3 floating-point
% matrix) that approximates best the affine transformation between:
%
% * Two point sets
% * Two raster images. In this case, the function first finds some features in
%   the `src` image and finds the corresponding features in `dst` image. After
%   that, the problem is reduced to the first case.
%
% In case of point sets, the problem is formulated as follows: you need to
% find a 2x2 matrix `A` and 2x1 vector `b` so that:
%
%     [A*|b*] = argmin_{[A|b]} sum_{i}(|| dst{i} - A*src{i}' - b ||^2)
%
% where `src{i}` and `dst{i}` are the i-th points in `src` and `dst`,
% respectively. `[A|b]` can be either arbitrary (when `FullAffine=true`) or
% have a form of:
%
%     [ a11, a12, b1;
%      -a12, a11, b2 ]
%
% when `FullAffine=false`.
%
% See also: cv.estimateAffine2D, cv.estimateAffinePartial2D,
%  cv.getAffineTransform, cv.getPerspectiveTransform, cv.findHomography,
%  imregtform
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sca=RLSsca(X,loc,p)
%Syntax: sca=RLSsca(X,loc,p)
%___________________________
%
% Calculates the Reweighted Least Squares (RLS) scale parameter
% of the columns of a matrix X.
%
% sca is the RLS estimated vector of scales.
% X is the matrix with the data sets.
% loc is the location vector. Its default value is the LMS location.
% p is the number of parameters. Its default value is 1.
%
% Reference:
% Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
%
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
%
% University e-mail: al10@ukc.ac.uk (until December 2002)
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Sep 3, 2001.

if nargin<1 | isempty(X)==1
   error('Not enough input arguments.');
else
   % X must be 2-dimensional
   if ndims(X)>2
      error('Invalid data set.');
   end
   % If X is a row vector make it a column vector
   if size(X,1)==1
      X=X';
   end

   % n is the length of the data set
   n=size(X,1);
end

% For a single data point there is no need to proceed
if n==1
   sca=0;
else

   % Check loc
   if nargin<2 | isempty(loc)==1
      % If you don't give loc, it is supposed the LMS location
      loc=LMSloc(X);
   else
      % loc must be a scalar/vector
      if min(size(loc))>1
         error('loc must be a scalar/vector.');
      end
      % loc must have the length of the columns of X
      if length(loc)~=size(X,2);
         error('loc must have the length of the columns of X.');
      end
   end

   % Check p
   if nargin<3 | isempty(p)==1
      % If p is omitted give it the value of 1
      p=1;
   else
      % p must be a scalar
      if sum(size(p))>2
         error('p must be a scalar.');
      end
      % p must be an integrer greater than 1
      if round(p)-p~=0 | p<1
         error('p must be an integrer greater than or equal to 1.');
      end
   end

   % Calculate the scale parameter for every column
   for i=1:size(X,2)

      % The differences from the location
      r=X(:,i)-loc(i);

      % Estimate the preliminary scale parameter
      s=LMSsca(r,0,p);
      if s==0
         w=1:n;
      else
         w=find(abs(r)/s<=2.5);
      end
      sca(i)=sqrt(sum(r(w).^2)/(length(w)-p));
   end
end

author: Sable

repo: mcbench-benchmarks


text:
% File testLUSOL.m
%
% Script for testing various LUSOL factorizations
% on the merged_S.mat matrix from 5 Dec 2006
% (S = 62177 x 75644 916437).
%
% 25 Jan 2008: First experiments on this particular S.
%              All used options.FactorTol = 2.
%              S is badly scaled, but several options
%              find the rank to be 61833 (rank deficiency 344).

load merged_S
[m,n] = size(S)
nnzS  = nnz(S)

options = lusolSet;
options.FactorTol = 2.0;
options.Pivoting = 'TPP';

%-----------------------------------------------------------------------------
% Try cheapest method:  L*U = S'
% It returns rank 61833.
disp(' ')
disp('Factor S(transpose)')
ST = S';
tic
[L1,U1,p1,q1,options] = lusolFactor(ST,options);
toc

% m       75644 >n       62177  Elems   916437  Amax   8.0E+05  Density   0.02
% Singular(m>n)  rank    61833  n-rank     344  nsing      344
% Merit    95.3  lenL   306685  L+U    1051242  Cmpressns    0  Incres   14.71
% Utri      556  lenU   744557  Ltol  2.00E+00  Umax   8.0E+05  Ugrwth 1.0E+00
% Ltri     3633  dense1      0  Lmax  2.00E+00
% bump    71455  dense2      0  DUmax  2.5E+03  DUmin  5.0E-01  condU  5.0E+03

rank1 = options.Rank;
rows1 = q1(1:rank1);
S1    = S(rows1,:);  % These should be independent rows of S.

%-----------------------------------------------------------------------------
% Try Rook Pivoting.
% S is too badly scaled for this to be efficient.
% Find column and row scales first.
disp(' ')
disp('Scale S now')
iprint  = 1;
scltol  = 0.9;
tic
[cscale,rscale] = gmscal(S,iprint,scltol);
disp(' ')
toc

% Apply scale factors to S.
C = spdiags(cscale,0,n,n);   Cinv = spdiags(1./cscale,0,n,n);
R = spdiags(rscale,0,m,m);   Rinv = spdiags(1./rscale,0,m,m);
SS = Rinv*S*Cinv;  % Scaled S
%-----------------------------------------------------------------------------

% Factor scaled SS with rook pivoting.
% It returns rank 61833.
disp(' ')
disp('Factor scaled S now with rook pivoting')
options.Pivoting = 'TRP';
tic
[L2,U2,p2,q2,options] = lusolFactor(SS,options);
toc

% m       62177 <n       75644  Elems   916437  Amax   1.0E+00  Density   0.02
% Singular(m<n)  rank    61833  n-rank   13811  nsing    13811
% MerRP    39.5  lenL   210899  L+U    1617359  Cmpressns    5  Incres   76.48
% Utri     2966  lenU  1406460  Ltol  2.00E+00  Umax   2.0E+00  Ugrwth 2.0E+00
% Ltri      207  dense1      0  Lmax  2.00E+00  Akmax  0.0E+00  Agrwth 0.0E+00
% bump    59004  dense2      0  DUmax  2.0E+00  DUmin  2.3E-05  condU  8.5E+04

rank2 = options.Rank;
rows2 = p2(1:rank2);
S2    = S(rows2,:);  % These should be independent rows of S.

%-----------------------------------------------------------------------------
% Factor scaled SS with rook pivoting.
% It returns rank 61833 also.
disp(' ')
disp('Factor scaled S(transpose) now with rook pivoting')
tic
[L3,U3,p3,q3,options] = lusolFactor(SS',options);
toc

% m       75644 >n       62177  Elems   916437  Amax   1.0E+00  Density   0.02
% Singular(m>n)  rank    61833  n-rank     344  nsing      344
% Merit    12.7  lenL   571936  L+U    1015445  Cmpressns    0  Incres   10.80
% Utri      556  lenU   443509  Ltol  2.00E+00  Umax   3.2E+01  Ugrwth 3.2E+01
% Ltri     2959  dense1      0  Lmax  2.00E+00
% bump    72129  dense2      0  DUmax  2.0E+01  DUmin  5.4E-05  condU  3.8E+05

rank3 = options.Rank;
rows3 = q3(1:rank3);
S3    = S(rows3,:);  % These should be independent rows of S.

disp(' ')
disp('rows1, rows2, rows3 are 3 sets of independent rows of S')
disp('   S1,    S2,    S3 are those submatrices of S')
disp(' ')

author: opencobra

repo: cobratoolbox


text:
function x = cvx_r2c( x, dim )

%
% Quick exit for complex quantities
%

if ~isreal( x ),
    error( 'Matrix must be real.' );
end

%
% Determine expansion dimension
%

sx = size( x );
if nargin < 2,
    dim = [ find( sx > 1 ), 1 ];
    dim = dim( 1 );
elseif ~isnumeric( dim ) || dim <= 0 || dim ~= floor( dim ),
    error( 'Second argument must be a dimension.' );
end
nd = length( sx );
if nd < dim || rem( sx( dim ), 2 ) ~= 0,
    error( 'The size of the array along the key dimension must be even.' );
end

%
% Extract the real and imaginary halves
%

[ ndxs{ 1 : nd } ] = deal( ':' );
ndxs{ dim } = 1 : 2 : sx( dim );
xr = cvx_subsref( x, ndxs{:} );
ndxs{ dim } = 2 : 2 : sx( dim );
xi = cvx_subsref( x, ndxs{:} );

%
% Combine
%

x = xr + 1i * xi;

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST