text:
function[d,x,y,t,z,w]=polysmooth_presort(varargin)
%POLYSMOOTH_PRESORT  Sort arguments to POLYSMOOTH in case of missing data.
%
%   POLYSMOOTH_PRESORT is a low-level function called by POLYSMOOTH.
%
%   [DS,XS,YS,TS,ZS,WS]=POLYSMOOTH_PRESORT(DS,XS,YS,TS,ZS,WS,P,TAU,STR),
%   where all the input fields are in the form output by TWODSORT or
%   SPHERESORT, performs several prelimary processing steps.
%
%   First, if the temporal bandwidth TAU is nonempty, any datapoints
%   outside of temporal window ABS(TS./TAU)>1 are set to NaNs.
%
%   After this, the data field, ZS, is checked to see if it contains any
%   NaNs. If it does, then all fields are sorted by distance DS along
%   third dimension, thus moving any interior NaN value to the bottom.
%
%   In the case that STR='pop', indicating a fixed population algorithm,
%   then finally all input fields are additionally truncated to length P
%   along their third dimension.
%
%   See also POLYSMOOTH.
%
%   Usage: [ds,xs,ys,zs,ws]=polysmooth_presort(ds,xs,ys,zs,ws,P,tau,str);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2018 J.M. Lilly --- type 'help jlab_license' for details

t=[];
z=[];
w=[];
B=[];
tau=[];
varstr='bandwidth';

d=varargin{1};
x=varargin{2};
y=varargin{3};
if nargin>=4
    t=varargin{4};
end
if nargin>=5
    z=varargin{5};
end
if nargin>=6
    w=varargin{6};
end
if nargin>=7
    B=varargin{7};
end
if nargin>=8
    tau=varargin{8};
end
if nargin>=9
    varstr=varargin{9};
end
%figure,jpcolor(squeeze(z))

if ~isempty(tau)
    if length(tau(:))~=1
        tau=vrep(tau,size(z,3),3);
    end
    z(abs(t./tau)>1)=nan;
    %d(abs(t./tau)>1)=nan;
end

%figure,jpcolor(squeeze(z))
%d(~isfinite(z))=nan;
%x(~isfinite(z))=nan;
%y(~isfinite(z))=nan;

if ~isempty(d)
    if ~strcmpi(varstr(1:3),'pop')
        if ~isempty(B)
            if length(B(:))~=1
                B=vrep(B,size(z,3),3);
            end
            z(abs(d./B)>1)=nan;
        end
    end
    %figure,jpcolor(squeeze(z))

    if ~isempty(z)
        z(~isfinite(d))=nan;  %make sure pattern of nans matches
        d(~isfinite(z))=nan;  %make sure pattern of nans matches
    end

    %bool will be true if current point is finite, but previous is not
    bool1=isnan(d);
    bool2=isnan(vshift(d,-1,3));bool2(:,:,1)=false;
    bool=(sum(~bool1&bool2,3)>0)&~(sum(bool1,3)==size(bool1,3));

    if anyany(bool)
        disp('POLYSMOOTH_PRESORT detecting interior NaNs; sorting to compress these.')
        [d,kk]=sort(d,3);
        ii=vrep([1:size(d,1)]',[size(d,2) size(d,3)],[2 3]);
        jj=vrep([1:size(d,2)],[size(d,1) size(d,3)],[1 3]);
        index=sub2ind(size(d),ii,jj,kk);
        %Can't use vindex here because i want to preserve the size
        x=x(index);
        y=y(index);
        if ~isempty(z),z=z(index);end
        if ~isempty(t),t=t(index);end
        if ~isempty(w),w=w(index);end
        %vindex(x,y,t,z,w,index,3);
        disp('Sorting complete.')
    end
    %1,vsize(d,x,y,t,z,w)

    if strcmpi(varstr(1:3),'pop')
        if ~isempty(B)
            vindex(d,x,y,t,z,w,1:min(maxmax(B),size(d,3)),3);
        end
    end
    %2,vsize(d,x,y,t,z,w)
    %figure,jpcolor(squeeze(d))
    numgood=squeeze(sum(sum(isfinite(d),2),1));
    %figure,plot(numgood)
    index=find(numgood~=0,1,'last');
    %index
    vindex(d,x,y,t,z,w,1:index,3);
    %3,vsize(d,x,y,t,z,w)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
clear all
close all

d1=3;
d2=5;
N=100;

disp('Model order selection');

%m=min([d1,d2])
m=3;

% Generate true factor matrices
W1=10*randn(d1,m);
W2=10*randn(d2,m);

% Observation noise covariance
sig=1;
E1=sig*randn(d1,N);
E2=sig*randn(d2,N);

if m==0
    X1=E1;
    X2=E2;
else
    Z=randn(m,N);
    X1=W1*Z+E1;
    X2=W2*Z+E2;
end

for i=1:4,
    CVA = spm_cva_prob (X1,X2,i-1);
    L(i)=CVA.L;
    bic(i)=CVA.bic;
    aic(i)=CVA.aic;
end
figure
plot([0:3],L);
hold on
plot([0:3],bic,'r');
plot([0:3],aic,'g');
xlabel('Number of Canonical Vectors');
legend('LogLike','BIC','AIC');

author: spm


text:
function matrix = genMatEllPPIFE3D(pde,mesh,fem,femI,femIF,PPtype,sig)
%% Generate global matrices and load vector of PPIFEM for 3D ellipic eq
%     -div(A grad u)  = f,    x\in \Omega
%      where A is a piecewise constant on Omega^+ and Omega^-.
% INPUTS:
% pde --- given data function from equation, e.g.
%         pde.A --- diffusion coefficient
%         pde.f --- right hand side function
%         pde.gD --- Dirichlet boundary value function
%         pde.one --- constant function 1.
% mesh --- mesh structure.
% fem --- global degree of freedom of FEM
% femI --- quadrature info on interface cells
% femIF --- quadrature info on interface faces, required in PPIFE.
% PPtype --- partial penalty type possible value
%            'N' : Nonsymmetric PPIFE (e = 1)
%            'S' : Symmetric PPIFE (e = -1)
%            'I' : Incomplete PPIFE (e = 0)
% OUTPUTS:
% matrix.S --- stiffness matrix (w/o boundary condition)
% matrix.E --- consistence matrix (w/o boundary condition)
% matrix.P --- penalty matrix (w/o boundary condition)
% matrix.A --- final FEM matrix (after boundary condition)
% matrix.rhsF --- load vector (w/o boundary condition)
% matrix.rhsE --- consistence vector (=0 if no interface face on boundary)
% matrix.rhsP --- penalty vector (=0 if no interface face on boundary)
% matrix.f --- final RHS matrix (after boundary condition)

% Last Modified: 08/07/2020 by Xu Zhang

%% 1. Stiffness Matrix
S = globMatrixIFE3DStiff(pde.A,mesh,femI,fem,fem);
% This matrix generator is a little faster than doing it separately.
% Sx = globMatrixIFE3D(pde.A,mesh,femI,fem,[1,0,0],fem,[1,0,0]);
% Sy = globMatrixIFE3D(pde.A,mesh,femI,fem,[0,1,0],fem,[0,1,0]);
% Sz = globMatrixIFE3D(pde.A,mesh,femI,fem,[0,0,1],fem,[0,0,1]);
% S = Sx+Sy+Sz;

%% 2. Face Matrices
d1 = 0; j1 = 1; d2 = 1; j2 = 0;
E = globMatrixIFE3DFace(pde.one, pde.A, fem, fem, femIF, d1,j1, d2,j2);

d1 = 0; j1 = 1; d2 = 0; j2 = 1;
P = globMatrixIFE3DFace(pde.one, pde.one,fem,fem, femIF, d1,j1, d2,j2);

%% 3. Generate the Right Hand Side Vector
rhsF = globRHSIFE3D(pde.f, mesh, fem, femI, [0,0,0]);
rhsE = globRHSIFE3DFace(pde.A, pde.exactu, fem, femIF, 1);
rhsP = globRHSIFE3DFace(pde.one, pde.exactu, fem, femIF, 0);

%% 4. Dirichlet Boundary Conditions
if strcmp(PPtype,'N')
    e = 1;
elseif strcmp(PPtype,'S')
    e = -1;
elseif strcmp(PPtype,'I')
    e = 0;
end
h = mesh.p(2,1) - mesh.p(1,1);
Atotal = S - E + e*E' + (sig/h)*P;
rhs = rhsF + e*rhsE + (sig/h)*rhsP;
tu = feval(pde.gD,fem.p(:,1),fem.p(:,2),fem.p(:,3));
ub = tu;
ub(fem.mapper) = 0;
rhsB = Atotal*ub;
A = Atotal(fem.mapper,fem.mapper);
f = rhs(fem.mapper) - rhsB(fem.mapper);

%% 5. Outputs
matrix = struct('A',A,'f',f,'S',S,'E',E,'P',P,'rhsF',rhsF,'rhsE',rhsE,...
    'rhsP',rhsP);

repo: ifem


text:
function value = s_begin ( s1, s2 )

%*****************************************************************************80
%
%% S_BEGIN is TRUE if one string matches the beginning of the other.
%
%  Discussion:
%
%    The strings are compared, ignoring blanks and capitalization.
%
%  Example:
%
%     S1              S2      S_BEGIN
%
%    'Bob'          'BOB'     TRUE
%    '  B  o b '    ' bo b'   TRUE
%    'Bob'          'Bobby'   TRUE
%    'Bobo'         'Bobb'    FALSE
%    ' '            'Bob'     FALSE    (Do not allow a blank to match
%                                       anything but another blank string.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character S1(*), S2(*), the strings to be compared.
%
%    Output, logical S_BEGIN, is TRUE if the strings match up to
%    the end of the shorter string, ignoring case.
%
  len1 = s_len_trim ( s1 );
  len2 = s_len_trim ( s2 );
%
%  If either string is blank, then both must be blank to match.
%  Otherwise, a blank string matches anything, which is not
%  what most people want.
%
  if ( len1 == 0 || len2 == 0 )

    if ( len1 == 0 && len2 == 0 )
      value = 1;
    else
      value = 0;
    end

    return

  end

  i1 = 0;
  i2 = 0;
%
%  Find the next nonblank in S1.
%
  while ( 1 )

    while ( 1 )

      i1 = i1 + 1;

      if ( len1 < i1 )
        value = 1;
        return
      end

      if ( s1(i1) ~= ' ' )
        break
      end

    end
%
%  Find the next nonblank in S2.
%
    while ( 1 )

      i2 = i2 + 1;

      if ( len2 < i2 )
        value = 1;
        return
      end

      if ( s2(i2) ~= ' ' )
        break
      end

    end
%
%  If the characters match, get the next pair.
%
    if ( ~ch_eqi ( s1(i1), s2(i2) ) )
      break
    end

  end

  value = 0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ values, times, vertices ] = mne_label_time_courses(labelfile,stcfile)
%
% function [ values, times ] = mne_label_time_courses(labelfile,stcfile)
%
% Extract the time courses corresponding to a label file from an stc file
%
% labelfile - The name of the label file
% stcfile   - The name of the stc file (must be on the same subject and
% hemisphere as the stc file
%
% values    - The time courses
% times     - The time points
% vertices  - The vertices corresponding to the time points
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%
me='MNE:mne_label_time_courses';

if nargin ~= 2
    error(me,'Incorrect number of arguments');
end
try
    stc = mne_read_stc_file(stcfile);
catch
    error(me,'%s',mne_omit_first_line(lasterr));
end

try
    lab = mne_read_label_file(labelfile);
catch
    error(me,'%s',mne_omit_first_line(lasterr));
end

[vertices,ia,ib] = intersect(double(stc.vertices),double(lab.vertices));
if length(vertices) == 0
    error(me,'No vertices match the label in the stc file');
end

values = stc.data(ia,:);
times  = zeros(1,size(stc.data,2));
for k = 0:length(times)-1
    times(k+1) = stc.tmin + k*stc.tstep;
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s = resampwor(p,m,n)
%RESAMWOR Random resampling without replacement
%
%   Description:
%   S = RESAMWOR(P) returns a new set of indices according to the
%   probabilities P without replacemnt. P is array of probabilities,
%   which are not necessarily normalized, though they must be
%   non-negative, and not all zero. The size of S is the size of P.
%
%   S = RESAMWOR(P,M,N) returns M by N matrix.
%
%   S = RESAMWOR(P,M) returns M by M matrix.
%
%   See also RESAMPSIM, RESAMPRES, RESAMPSTR, RESAMPDET
%
% Copyright (c) 2003-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin<2
    [m,n] = size(p);
elseif nargin==2
    n = m;
end
if m*n>numel(p)
  error('In resampling without replacment M*N has to be smaller than numel(P)')
end
r=rand([m,n]);
s=zeros([m,n]);
for i=1:m*n
  pc=cumsum(p(:));
  pc=pc./pc(end);
  s(i)=binsgeq(pc,r(i));
  p(s(i))=0;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% StackExchange Signal Processing Q63449
% https://dsp.stackexchange.com/questions/63449
% Deconvolution of an Image Acquired by a Square Uniform Detector
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     25/01/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

%% Simulation Parameters

imageFileName   = 'Lenna256.png';
kernelRadius    = 1;
convShape       = CONVOLUTION_SHAPE_VALID;

paramLambda = 0.005;

maxSize = 256;

%% Generate Data

mA = im2double(imread(imageFileName));
numRows = size(mA, 1);
numCols = size(mA, 2);

imgSize = min(numRows, numCols);
maxSize = min(imgSize, maxSize);
mA = mA(1:imgSize, 1:imgSize, :);
mA = imresize(mA, [maxSize, maxSize]);

mK = ones((2 * kernelRadius) + 1);
% mK = rand((2 * kernelRadius) + 1);
mK = mK / sum(mK(:));

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        convShapeString = 'full';
    case(CONVOLUTION_SHAPE_SAME)
        convShapeString = 'same';
    case(CONVOLUTION_SHAPE_VALID)
        convShapeString = 'valid';
end

hFigure     = figure();
hAxes       = axes();
hImageObj   = imshow(mA);
set(get(hAxes, 'Title'), 'String', {['Input Image - Lenna']}, ...
    'FontSize', fontSizeTitle);

mB = conv2(mA, mK, convShapeString);

hFigure     = figure();
hAxes       = axes();
hImageObj   = imshow(mB);
set(get(hAxes, 'Title'), 'String', {['Sensor Image - Lenna']}, ...
    'FontSize', fontSizeTitle);

%% Solution by Linear Algebra

mKK = CreateConvMtx2D(mK, maxSize, maxSize, convShape);
% Basically:
% vB = mKK * mA(:);
vA = mKK \ mB(:);
% vA = pinv(full(mKK)) * mB(:);

vA = ((mKK.' * mKK) + (paramLambda * speye(maxSize * maxSize))) \ (mKK.' * mB(:));

mAA = reshape(vA, maxSize, maxSize); %<! Restored Image

hFigure     = figure();
hAxes       = axes();
hImageObj   = imshow(mAA);
set(get(hAxes, 'Title'), 'String', {['Estimated Image - Lenna']}, ...
    'FontSize', fontSizeTitle);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [traj] = initTraj(NbSteps)
traj.Rot = zeros(3,3,NbSteps);
traj.phi = zeros(1,NbSteps);
traj.theta = zeros(1,NbSteps);
traj.psi = zeros(1,NbSteps);
traj.v = zeros(3,NbSteps);
traj.x = zeros(3,NbSteps);
traj.omega_b = zeros(3,NbSteps);
traj.a_b = zeros(3,NbSteps);
end

author: mbrossar


text:
% Init Script
% Remarks:
%   1.  If there is a predefined variable called 'subStreamNumberDefault'
%       it will be used to set the SubStream of the Random Number
%       Generator.
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.005     02/08/2019  Royi Avital
%   *   Added Octave Comptability.
% - 1.0.004     05/07/2018  Royi Avital
%   *   Added case for 'subStreamNumberDefault = 0'.
% - 1.0.003     29/06/2018  Royi Avital
%   *   Clearing all variables but 'subStreamNumberDefault'.
%   *   Setting 'subStreamNumber' only if 'subStreamNumberDefault' not
%       defined.
% - 1.0.002     19/04/2018  Royi Avital
%   *   Added 'markerSizeMedium' (Matches 'markerSizeNormal').
% - 1.0.001     24/11/2016  Royi Avital
%   *   Added figure positions options.
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%
%
%% General Parameters

close('all');
clearvars('-except', 'subStreamNumberDefault');
clc();

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

% See https://stackoverflow.com/questions/2246579
isOctave = (exist('OCTAVE_VERSION', 'builtin') ~= 0);

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 0400, 0300];
figPosMedium    = [100, 100, 0800, 0600];
figPosLarge     = [100, 100, 0960, 0720];
figPosXLarge    = [100, 100, 1100, 0825];
figPosX2Large   = [100, 100, 1200, 0900];
figPosX3Large   = [100, 100, 1400, 1050];
figPosDefault   = [100, 100, 0560, 0420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeMedium    = 8;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

% https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html
% https://www.mathworks.com/matlabcentral/answers/160332
mColorOrder = get(groot, 'DefaultAxesColorOrder');

if(isOctave == FALSE)
    randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');

    if(exist('subStreamNumberDefault', 'var') && (subStreamNumberDefault ~= 0))
        subStreamNumber = subStreamNumberDefault;
    else
        subStreamNumber = round(sum(clock()));
    end
    set(randomNumberStream, 'Substream', subStreamNumber);
    RandStream.setGlobalStream(randomNumberStream);
end


text:
function max_pr = trace_fit_extreme(C,fr,t_int,fac)

if nargin < 4 || isempty(fac); fac = 1; end
if nargin < 3 || isempty(t_int); t_int = 0.25; end
if nargin < 2 || isempty(fr); fr = 30; end

Np = round(t_int*fr);
[K,T] = size(C);
bas = zeros(K,1);
sn = zeros(K,1);
for i = 1:K
    [~,density,xmesh] = kde(C(i,:));
    [~,ind] = max(density);
    bas(i) = xmesh(ind);
    sn(i) = std(C(i,C(i,:)<bas(i)))/sqrt(1-2/pi);
end

mu = norminv(1-1/T);
sig = norminv(1-exp(-1)/T) - mu;

z = bsxfun(@times, 1./(fac*sn(:)), bsxfun(@minus, C, bas));  % normalized z scores
z_pr = exp(-exp(-(z-mu)/sig));

filt_z = filter(ones(1,Np),1,log(z_pr),[],2)/sqrt(Np);
max_pr = exp(max(filt_z,[],2));

author: flatironinstitute

repo: CaImAn-MATLAB


text:
% drawLineSegments draws Lines using given set of points.
% Input     - MapSize : The size of EdgeMap which can be expressed as
%                       [ROW COL]
%             LineSegments : The set of points which will be drawn.
%                           EdgeContour can be used as input as well.
function drawLineSegments(MapSize, LineSegments)

%% Initialize
NumberOfSegments = length(LineSegments);

% Default setting for plot.
hold on;
minx = 1; maxx = MapSize(2);
miny = 1; maxy = MapSize(1);
axis([minx maxx miny maxy]);
axis off;
set(gca, 'Ydir', 'reverse')

% Use different colors to distinguish each edge.
ColorSet = [
0.00 0.00 1.00 % Data 1 - blue
0.00 1.00 0.00 % Data 2 - green
1.00 0.00 0.00 % Data 3 - red
0.00 1.00 1.00 % Data 4 - cyan
1.00 0.00 1.00 % Data 5 - magenta
0.75 0.75 0.00 % Data 6 - RGB
0.25 0.25 0.25 % Data 7
0.75 0.25 0.25 % Data 8
0.95 0.95 0.00 % Data 9
0.25 0.25 0.75 % Data 10
0.75 0.75 0.75 % Data 11
0.00 0.50 0.00 % Data 12
0.76 0.57 0.17 % Data 13
0.54 0.63 0.22 % Data 14
0.34 0.57 0.92 % Data 15
1.00 0.10 0.60 % Data 16
0.88 0.75 0.73 % Data 17
0.10 0.49 0.47 % Data 18
0.66 0.34 0.65 % Data 19
0.99 0.41 0.23 % Data 20
];

%% Draw each Segement.
for i=1:NumberOfSegments
    Segment = LineSegments{i};
    x = Segment(:,2);
    y = Segment(:,1);
    plot(x, y, 'LineWidth', 2, 'Color', ColorSet(mod(i,20)+1,:));
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [f,Ph]=haralick_n(qs,nL,q)
% Haralick textures measurements

Ph = coocurrance_alldir_mod(q);
% last row an column corresponds to outside ROI!
Ph=Ph(1:end-1,1:end-1);
nL=nL-1;

%compute features
R=sum(Ph(:));
Ph=Ph/R;
% Energy (1)
f(1)=sum(sum(Ph.^2));
% Contrast (2)
f(2)=0.0;
for n=0:nL-1
   temp=0;
   for i=1:nL
      for j=1:nL
         if (abs(i-j) == n)
            temp=temp+Ph(i,j);
         end
      end
   end
   f(2)=f(2)+n^2*temp;
end
% Correlation
%using symmetry Ph'=Ph!
Px=sum(Ph);
Py=sum(Ph');
vec=[1:nL];
ux=sum(Px .*vec);
uy=sum(Py .*vec);

varx=sum(Px .* vec.^2)-ux^2;
sigx=sqrt(varx);
vary=sum(Py .* vec.^2)-uy^2;
sigy=sqrt(vary);
u=vec*Ph(i,j)*vec';
f(3)=(u-ux*uy)/(sigx*sigy);

%Entropy (3)
f(4)=-sum(sum(Ph.*log(Ph+realmin))); % log????

% variance
f(5)=0;
for i=1:nL
   for j=1:nL
      f(5)=f(5)+(i-u)^2*Ph(i,j);
   end
end
% P(x+y)
f(6)=0;  % sum of entropy...
for k=2:2*nL
   temp=0.0;
   for i=1:nL
      for j=1:nL
         if ((i+j) == k)
            temp=temp+Ph(i,j);
         end
      end
   end
   Pxpy(k)=temp;
   f(6)=f(6)-temp*log(temp+realmin);
end

% inverse different moment..
f(7)=0;
f(8)=0;  %Homogeneity (4)
for i=1:nL
   for j=1:nL
      temp1=1/(1+(i-j)^2)*Ph(i,j);
      temp2=1/(1+abs(i-j))*Ph(i,j);
      f(7)=f(7)+temp1;
      f(8)=f(8)+temp2;
   end
end

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function kern = rbfwhiteKernExpandParam(kern, params)

% RBFWHITEKERNEXPANDPARAM Create kernel structure from RBF-WHITE kernel's
% parameters.
% FORMAT
% DESC returns a RBF-WHITE kernel structure filled with the parameters in
% the given vector. This is used as a helper function to enable parameters
% to be optimised in, for example, the NETLAB optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : rbfwhiteKernParamInit, rbfwhiteKernExtractParam, kernExpandParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

kern.inverseWidth = params(1);
kern.variance = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [Model2] = train()
    [X] = LoadImages();
    [Y] = LoadLabels();
    [M,P] = GetMP(X,Y);

    field1 = 'M';
    field2 = 'P';
    Model2 = struct(field1,M,field2,P);
    save('Model2.mat','Model2');
end


text:
function varargout = log(varargin)
%LOG (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/LOG CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
         % Try to detect logsumexp construction etc
        varargout{1} = check_for_special_cases(varargin{:});
        % Nope, then just define this logarithm
        if isempty(varargout{1})
            varargout{1} = InstantiateElementWise(mfilename,varargin{:});
        end

    case 'char'

        X = varargin{3};
        F = (X >= 1e-8);

        operator = struct('convexity','concave','monotonicity','increasing','definiteness','none','model','callback');
        operator.convexhull = @convexhull;
        operator.bounds = @bounds;
        operator.domain = [0 inf];
        operator.derivative = @(x)(1./(abs(x)+eps));
        operator.inverse = @(x)(exp(x));

        varargout{1} = F;
        varargout{2} = operator;
        varargout{3} = X;

    otherwise
        error('SDPVAR/LOG called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
if xL <= 0
    % The variable is not bounded enough yet
    L = -inf;
else
    L = log(xL);
end
if xU < 0
    % This is an infeasible problem
    L = inf;
    U = -inf;
else
    U = log(xU);
end

function [Ax, Ay, b, K] = convexhull(xL,xU)
K = [];
if xL <= 0
    fL = inf;
else
    fL = log(xL);
end
fU = log(xU);
dfL = 1/(xL);
dfU = 1/(xU);
%xM = (xU - xL)/(fU-fL);
xM = (xL + xU)/2;
fM = log(xM);
dfM = 1/xM;

[Ax,Ay,b] = convexhullConcave(xL,xM,xU,fL,fM,fU,dfL,dfM,dfU);
remove = isinf(b) | isinf(Ax) | isnan(b);
if any(remove)
    remove = find(remove);
    Ax(remove)=[];
    b(remove)=[];
    Ay(remove)=[];
end

function f = check_for_special_cases(x)
f = [];
% Check for log(1+x)
base = getbase(x);
if all(base(:,1)==1)
    f = slog(x-1);
    return;
end
% Check if user is constructing log(sum(exp(x)))
if base(1)~=0
    return
end
if ~all(base(2:end)==1)
    return
end
modelst = yalmip('extstruct',getvariables(x));
if isempty(modelst)
    return;
end
if length(modelst)==1
    models{1} = modelst;
else
    models = modelst;
end
% LOG(DET(X))
if length(models)==1
    if strcmp(models{1}.fcn,'det_internal')
        n = length(models{1}.arg{1});
        try
            f = logdet(reshape(models{1}.arg{1},sqrt(n),sqrt(n)));
        catch
        end
        return
    end
end
% LOG(EXP(x1)+...+EXP(xn))
for i = 1:length(models)
    if ~strcmp(models{i}.fcn,'exp')
        return
    end
end
p = [];
for i = 1:length(models)
    p = [p;models{i}.arg{1}];
end
f = logsumexp(p);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
x0=2;
[x,y]=fminsearch(@fun5,x0)

author: Eurus-Holmes


text:
function [kld, llh1, llh2] = GLG_kld( theta1, theta2, sz )

% Kullback-Leibler divergence between GLG models
% The KLD is computed by simulating observations from model 1 and
% evaluating the composite log-likelihood of these observations under
% model 1 and model 2.
%
% Syntax:
%   [kld, llh1, llh2] = GLG_kld( theta1, theta2, sz )
%
%
% Input:
%   theta1 : Parameters of model 1.
%
%   theta2 : Parameters of model 2.
%
%   sz     : Size of image simulated from model 1
%
%
% Output:
%   kld    : The Kullback-Leibler divergence from model 1 to model 2
%
%   llhX   : The likelihood of the simulation from model 1 under model X
%
%
% See also: GLG_SIMULATION

if ~isequal( size(theta1), size(theta2) )
    error('Models must be of the same size');
end

% Simulate from model1
tree = GLG_simulation( theta1, sz );

% Evaluate likelihood for each direction and level
llh1 = GLG_llh( theta1, tree );
llh2 = GLG_llh( theta2, tree );

% Sum composite likelihood and directions into KLD
kld = sum( llh1(:) - llh2(:) );

end

author: Sable

repo: mcbench-benchmarks


text:
function prettifiedJSON = prettyJson(rawJSON)
% prettyJson
%   Read in raw JSON text that is generated by encodejson function and
%   usually in one line, then convert to pretty format with proper
%   indentations and line breaks
%
% Input:
%   rawJSON             input char array in JSON format
%
% Output:
%   prettifiedJSON      prettified char array in JSON format
%
% NOTE: This function is to optimize displaying effect of JSON string through
% recognizing general delimited characters, but may change the file content
% and cause problem in several uncommon cases. Please avoid the use of this
% function on Matlab structures that contain elements with any of the three
% character combinations: `":[` and `",` and `"],`.
%
% Usage: prettifiedJSON = prettyJson(rawJSON)
%

% handle input
if isstring(rawJSON)
    % if input is supplied as string
    rawJSON = char(rawJSON);
elseif ~ischar(rawJSON)
    % if input is not supplied as char array
    error('The input data type is Not correct!');
end
jsonStr = rawJSON;

% process input string to pretty JSON format
jsonStr = regexprep(jsonStr,'^{','{\n\t','lineanchors');  % left  curly  bracket
jsonStr = regexprep(jsonStr,'}$','\n}\n','lineanchors');  % right curly  bracket
jsonStr = strrep(jsonStr,'":[','":[\n\t\t');              % left  square bracket
jsonStr = strrep(jsonStr,'",','",\n\t\t');                % internal elements in array
jsonStr = strrep(jsonStr,'],','],\n\t');                  % last element in array
jsonStr = compose(jsonStr);                               % translate escape-characters

prettifiedJSON = jsonStr{1};                              % convert from cell to char

end

author: SysBioChalmers

repo: Human-GEM


text:
%PLOT_CIRCLE Draw a circle
%
% plot_circleC, R, OPTIONS) draws a circle on the current plot with
% centre C=[X,Y] and radius R.  If C=[X,Y,Z] the circle is drawn in the
% XY-plane at height Z.
%
% If C (2xN) then N circles are drawn.  If R (1x1) then all
% circles have the same radius or else R (1xN) to specify the radius of
% each circle.
%
% H = plot_circle(...) as above but return handles. For multiple
% circles H is a vector of handles, one per circle.
%
% Options::
% 'edgecolor'   the color of the circle's edge, Matlab color spec
% 'fillcolor'   the color of the circle's interior, Matlab color spec
% 'alpha'       transparency of the filled circle: 0=transparent, 1=solid
% 'alter',H     alter existing circles with handle H
%
% - For an unfilled circle:
%   - any standard MATLAB LineStyle such as 'r' or 'b---'.
%   - any MATLAB LineProperty options can be given such as 'LineWidth', 2.
% - For a filled circle any MATLAB PatchProperty options can be given.
%
% Example::
%
%          H = plot_circle([3 4]', 2, 'r');  % draw red circle
%          plot_circle([3 4]', 3, 'alter', H); % change the circle radius
%          plot_circle([3 4]', 3, 'alter', H, 'LineColor', 'k'); % change the color
%
% Notes::
% - The 'alter' option can be used to create a smooth animation.
% - The circle(s) is added to the current plot irrespective of hold status.
%
% See also PLOT_ELLIPSE, PLOT_BOX, PLOT_POLY.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math
function handles = plot_circle(centre, rad, varargin)

    opt.fillcolor = [];
    opt.alpha = 1;
    opt.edgecolor = 'k';
    opt.alter = [];

    [opt,arglist] = tb_optparse(opt, varargin);

    if ~isempty(opt.alter) & ~ishandle(opt.alter)
        error('SMTB:plot_circle:badarg', 'argument to alter must be a valid graphic object handle');
    end

    holdon = ishold;
    hold on

	n = 50;
	th = [0:n]'/n*2*pi;

    if length(rad) == 1
        rad = rad*ones(numcols(centre),1);
    end
    if length(centre) == 2 || length(centre) == 3
        centre = centre(:);
    end

    for i=1:numcols(centre)
        x = rad(i)*cos(th) + centre(1,i);
        y = rad(i)*sin(th) + centre(2,i);
        if numrows(centre) > 2
            % plot 3D data
            z = ones(size(x))*centre(3,i);
            if isempty(opt.alter)
                h(i) = plot3(x, y, z, varargin{:});
            else
                set(opt.alter(i), 'xdata', x, 'ydata', y, 'zdata', z, arglist{:});
            end
        else
            % plot 2D data
            if isempty(opt.fillcolor)
                if isempty(opt.alter)
                    h(i) = plot(x, y, arglist{:});
                else
                    set(opt.alter(i), 'xdata', x, 'ydata', y, arglist{:});
                end
            else
                if isempty(opt.alter)
                    h(i) = patch(x, y, 0*y, 'FaceColor', opt.fillcolor, ...
                        'FaceAlpha', opt.alpha, 'EdgeColor', opt.edgecolor, arglist{:});
                else
                    set(opt.alter(i), 'xdata', x, 'ydata', y, arglist{:});
                end

            end
        end
    end

    if holdon == 0
        hold off
    end

    if nargout > 0
        handles = h;
    end

author: petercorke

repo: spatialmath-matlab


text:
function partition_count_values_test ( )

%*****************************************************************************80
%
%% PARTITION_COUNT_VALUES_TEST demonstrates the use of PARTITION_COUNT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PARTITION_COUNT_VALUES_TEST:\n' );
  fprintf ( 1, '  PARTITION_COUNT_VALUES returns values of \n' );
  fprintf ( 1, '  the integer partition count function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N         P(N)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, fn ] = partition_count_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %10d\n', n, fn );

  end

  return
end


text:
function dq = DQconj(dv,WHICHCONJ)

% DQCONJ   Dual quaternion conjugate
%
%    DQ = DQCONJ(DV) returns the dual quaternion conjugate of the dual
%      quaternion DV. DV (resp. DQ) is a 8-vector representing a dual
%      quaternion or an array 8*N (column i represents dual quaternion i)
%      where N is the number of dual quaternions. The default type of
%      conjugate is 'point' (see below). DQ has the same size as DV.
%
%    DQ = DQCONJ(DV,WHICHCONJ) returns the type of dual quaternion conjugate
%       corresponding to WHICHCONJ (DV = Q0+eps*Q1):
%         - 'point': mixed conjugate: dV* = Q0*-eps*Q1* (used for point
%         transformations).
%         - 'line': quaternion conjugate: dV* = Q0*+eps*Q1* (used for line
%         transformations).
%         - 'other': dual number conjugate: dV* = Q0-eps*Q1
%
% See also QCONJ, DQMULT

if nargin < 2
    WHICHCONJ = 'point'; % default is 'mixed' conjugate' (point transformation)
end

sdv = size(dv);
if sdv == [1 8], dv = dv'; sdv = size(dv); end

% wrong size
if sdv(1) ~= 8
    error('DualQuaternion:DQconj:wrongsize',...
        '%d rows in array dv. It should be 8.',sdv(1));
end

n = sdv(2);
dq = zeros(8,n);
switch WHICHCONJ
    case 'point' % % case of point transformation
        dq(1:4,:) = Qconj(dv(1:4,:));
        dq(5:8,:) = -Qconj(dv(5:8,:));
    case 'line' % case of line transformation
        dq(1:4,:) = Qconj(dv(1:4,:));
        dq(5:8,:) = Qconj(dv(5:8,:));
    case 'other' % dual number conjugate
        dq(1:4,:) = dv(1:4,:);
        dq(5:8,:) = -dv(5:8,:);
    otherwise
        error('DualQuaternion:DQconj:wrongsize',...
            'WHICHCONJ must be ''point'', ''line'' or ''other'' ');
end


text:
function theta=Theta_ADMM(Y,W,H,X,L,alfa,beta,theta,rho,maxiter)

%%
%  Solve the JPLAY's subproblem: theta, using ADMM

%% Initializing Setting
epsilon = 1e-6;
iter=0;

G=zeros(size(theta));
Q=zeros(size(theta*X));
P=zeros(size(theta*X));
% M=zeros(size(theta*X));
lamda1=zeros(size(theta*X));
lamda2=zeros(size(theta));
lamda3=zeros(size(theta*X));
lamda4=zeros(size(theta*X));
% lamda5=zeros(size(theta*X));

stop = false;
mu=1e-3;
% rho=2;%1.4:Houston2018 2:Houston2013
mu_bar=1e+6;

GL=(X*L*X');%Graph Laplacian

  while ~stop && iter < maxiter+1

      iter=iter+1;
      %solve theta
      theta=(mu*(H*X')+lamda1*X'+mu*G+lamda2+mu*(Q*X')+lamda3*X'+mu*(P*X')+lamda4*X')/(3*mu*(X*X')+beta*GL+mu*eye(size(X*X')));
      %solve H
      H=(alfa*(W'*W)+(G*G')+mu*eye(size(W'*W)))\(alfa*(W'*Y)+(G*X)+mu*(theta*X)-lamda1);
      %solve G
      G=((H*H')+mu*eye(size(H*H')))\(mu*theta-lamda2+(H*X'));
      %solve Q
      Q=max(theta*X-(lamda3/mu),0);
      %solve P
      MidV=theta*X-lamda4/mu;
      for i=1:size(MidV,2)
          if norm(MidV(:,i))<=1
             P(:,i)=MidV(:,i);
          else
             P(:,i)=MidV(:,i)/norm(MidV(:,i));
          end
      end
%       M=max(abs(theta*(X-XW)-lamda5/mu)-(beta*0.01/mu),0).*sign(theta*(X-XW)-lamda5/mu);
     %update Lagrange multipliers
     lamda1=lamda1+mu*(H-theta*X);
     lamda2=lamda2+mu*(G-theta);
     lamda3=lamda3+mu*(Q-theta*X);
     lamda4=lamda4+mu*(P-theta*X);
%      lamda5=lamda5+mu*(M-theta*(X-XW));
     %update penalty parameter
     mu=min(mu*rho,mu_bar);
     %computer errors
     r_H=norm(H-theta*X,'fro');
     r_G=norm(G-theta,'fro');
     r_Q=norm(Q-theta*X,'fro');
     r_P=norm(P-theta*X,'fro');
%      r_M=norm(M-theta*(X-XW),'fro');
     %check the convergence conditions
     if r_H<epsilon&&r_G<epsilon&&r_Q<epsilon&&r_P<epsilon%&&r_M<epsilon
         stop = true;
         break;
     end
  end
end

author: BehnoodRasti


text:
function [q] = MeshQuality(EToV,VX,VY)
%
% Purpose: Assess the triangle mesh quality using a mesh quality indicators.
%          q=2r/R, where r is the inradius and R is the circumradius of a
%          triangle.
%
%          Function is useful for detecting degenerate triangles.
%
%          EToV : Element-To-Vertice table
%          VX   : x-table for vertices
%          VY   : y-table for vertices
%
% By Allan P. Engsig-Karup, apek@imm.dtu.dk.

% Compute side lengths of triangles
lx = VX(EToV(:,[1 2 3])) - VX(EToV(:,[3 1 2]));
ly = VY(EToV(:,[1 2 3])) - VY(EToV(:,[3 1 2]));
l = sqrt(lx.^2+ly.^2);
a = l(:,1); b = l(:,2); c = l(:,3);
q = (b+c-a).*(c+a-b).*(a+b-c)./(a.*b.*c);

N = 100;
x = linspace(0,1,N+1);
dx = x(2)-x(1);
subdivision = x(1:N)+dx/2;
count = zeros(1,N);
K = size(EToV,1);
for i = 1 : length(subdivision)
    count(i) = length(find(q>subdivision(i)-dx/2 & q<=subdivision(i)+dx/2))/K*100;
end
figure
bar(subdivision,count,1)
minlimit = find(q<0.7);
if isempty(minlimit)
    axis([0.7 1 -1 max(count)*1.05])
else
    axis([min(q) 1 -1 max(count)*1.05])
    hold on
    plot([0.7 0.7],[0 100],'r--')
end
xlabel('Mesh quality')
ylabel('Percentage of elements')
colormap(cool)


text:
function FB = FbMake( dim, flag, show )
% Various 1D/2D/3D filterbanks (hardcoded).
%
% USAGE
%  FB = FbMake( dim, flag, [show] )
%
% INPUTS
%  dim     - dimension
%  flag    - controls type of filterbank to create
%          - if d==1
%            1: gabor filter bank for spatiotemporal stuff
%          - if d==2
%            1: filter bank from Serge Belongie
%            2: 1st/2nd order DooG filters.  Similar to Gabor filterbank.
%            3: similar to Laptev&Lindberg ICPR04
%            4: decent seperable steerable? filterbank
%            5: berkeley filterbank for textons papers
%            6: symmetric DOOG filters
%          - if d==3
%            1: decent seperable steerable filterbank
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%
% EXAMPLE
%  FB = FbMake( 2, 1, 1 );
%
% See also FBAPPLY2D
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(show) ); show=0; end

% create FB
switch dim
  case 1
    FB = FbMake1D( flag );
  case 2
    FB = FbMake2D( flag );
  case 3
    FB = FbMake3d( flag );
  otherwise
    error( 'dim must be 1 2 or 3');
end

% display
FbVisualize( FB, show );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake1D( flag )
switch flag
  case 1  %%% gabor filter bank for spatiotemporal stuff
    omegas = 1 ./ [3 4 5 7.5 11];
    sigmas =      [3 4 5 7.5 11];
    FB = FbMakegabor1D( 15, sigmas, omegas );

  otherwise
    error('none created.');
end

function FB = FbMakegabor1D( r, sigmas, omegas )
for i=1:length(omegas)
  [feven,fodd]=filterGabor1d(r,sigmas(i),omegas(i));
  if( i==1 ); FB=repmat(feven,[2*length(omegas) 1]); end
  FB(i*2-1,:)=feven; FB(i*2,:)=fodd;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake2D( flag )

switch flag
  case 1  %%% filter bank from Berkeley / Serge Belongie
    r=15;
    FB = FbMakegabor( r, 6, 3, 3, sqrt(2) );
    FB2 = FbMakeDOG( r, .6, 2.8, 4);
    FB = cat(3, FB, FB2);
    %FB = FB(:,:,1:2:36); %include only even symmetric filters
    %FB = FB(:,:,2:2:36); %include only odd symmetric filters

  case 2 %%% 1st/2nd order DooG filters.  Similar to Gabor filterbank.
    FB = FbMakeDooG( 15, 6, 3, 5, .5) ;

  case 3 %%% similar to Laptev&Lindberg ICPR04
    % Wierd filterbank of Gaussian derivatives at various scales
    % Higher order filters probably not useful.
    r = 9;   dims=[2*r+1 2*r+1];
    sigs = [.5 1 1.5 3]; % sigs = [1,1.5,2];

    derivs = [];
    %derivs = [ derivs; 0 0 ]; % 0th order
    %derivs = [ derivs; 1 0; 0 1 ]; % first order
    %derivs = [ derivs; 2 0; 0 2; 1 1 ]; % 2nd order
    %derivs = [ derivs; 3 0; 0 3; 1 2; 2 1 ]; % 3rd order
    %derivs = [ derivs; 4 0; 0 4; 1 3; 3 1; 2 2 ]; % 4th order
    derivs = [ derivs; 0 1; 0 2; 0 3; 0 4; 0 5]; % 0n order
    derivs = [ derivs; 1 0; 2 0; 3 0; 4 0; 5 0]; % n0 order
    cnt=1;  nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, [sigs(s) sigs(s)], derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 length(sigs)*nderivs]); end
        FB(:,:,cnt) = dG; cnt=cnt+1;
        %dG = filterDoog( dims, [sigs(s)*3 sigs(s)], derivs(i,:), 0 );
        %FB(:,:,cnt) = dG; cnt=cnt+1;
        %dG = filterDoog( dims, [sigs(s) sigs(s)*3], derivs(i,:), 0 );
        %FB(:,:,cnt) = dG; cnt=cnt+1;
      end
    end

  case 4 % decent seperable steerable? filterbank
    r = 9;   dims=[2*r+1 2*r+1];
    sigs = [.5 1.5 3];
    derivs = [1 0; 0 1; 2 0; 0 2];
    cnt=1;  nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, [sigs(s) sigs(s)], derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 length(sigs)*nderivs]); end
        FB(:,:,cnt) = dG; cnt=cnt+1;
      end
    end
    FB2 = FbMakeDOG( r, .6, 2.8, 4);
    FB = cat(3, FB, FB2);

  case 5  %%% berkeley filterbank for textons papers
    FB = FbMakegabor( 7, 6, 1, 2, 2 );

  case 6  %%% symmetric DOOG filters
    FB = FbMakeDooGSym( 4, 2, [.5 1] );

  otherwise
    error('none created.');
end

function FB = FbMakegabor( r, nOrient, nScales, lambda, sigma )
% multi-scale even/odd gabor filters. Adapted from code by Serge Belongie.
cnt=1;
for m=1:nScales
  for n=1:nOrient
    [F1,F2]=filterGabor2d(r,sigma^m,lambda,180*(n-1)/nOrient);
    if(m==1 && n==1); FB=repmat(F1,[1 1 nScales*nOrient*2]); end
    FB(:,:,cnt)=F1;  cnt=cnt+1;   FB(:,:,cnt)=F2;  cnt=cnt+1;
  end
end

function FB = FbMakeDooGSym( r, nOrient, sigs )
% Adds symmetric DooG filters.  These are similar to gabor filters.
cnt=1; dims=[2*r+1 2*r+1];
for s=1:length(sigs)
  Fodd = -filterDoog( dims, [sigs(s) sigs(s)], [1 0], 0 );
  Feven = filterDoog( dims, [sigs(s) sigs(s)], [2 0], 0 );
  if(s==1); FB=repmat(Fodd,[1 1 length(sigs)*nOrient*2]); end
  for n=1:nOrient
    theta = 180*(n-1)/nOrient;
    FB(:,:,cnt) = imrotate( Feven, theta, 'bil', 'crop' );  cnt=cnt+1;
    FB(:,:,cnt) = imrotate( Fodd,  theta, 'bil', 'crop' );  cnt=cnt+1;
  end
end

function FB = FbMakeDooG( r, nOrient, nScales, lambda, sigma )
% 1st/2nd order DooG filters.  Similar to Gabor filterbank.
% Defaults: nOrient=6, nScales=3, lambda=5, sigma=.5,
cnt=1; dims=[2*r+1 2*r+1];
for m=1:nScales
  sigma = sigma * m^.7;
  Fodd = -filterDoog( dims, [sigma lambda*sigma^.6], [1,0], 0 );
  Feven = filterDoog( dims, [sigma lambda*sigma^.6], [2,0], 0 );
  if(m==1); FB=repmat(Fodd,[1 1 nScales*nOrient*2]); end
  for n=1:nOrient
    theta = 180*(n-1)/nOrient;
    FB(:,:,cnt) = imrotate( Feven, theta, 'bil', 'crop' );  cnt=cnt+1;
    FB(:,:,cnt) = imrotate( Fodd,  theta, 'bil', 'crop' );  cnt=cnt+1;
  end
end

function FB = FbMakeDOG( r, sigmaStr, sigmaEnd, n )
% adds a serires of difference of Gaussian filters.
sigs = sigmaStr:(sigmaEnd-sigmaStr)/(n-1):sigmaEnd;
for s=1:length(sigs)
  FB(:,:,s) = filterDog2d(r,sigs(s),2); %#ok<AGROW>
  if( s==1 ); FB=repmat(FB,[1 1 length(sigs)]); end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function FB = FbMake3d( flag )

switch flag
  case 1 % decent seperable steerable filterbank
    r = 25; dims=[2*r+1 2*r+1 2*r+1];
    sigs = [.5 1.5 3];
    derivs = [0 0 1; 0 1 0; 1 0 0; 0 0 2; 0 2 0; 2 0 0];
    cnt=1; nderivs = size(derivs,1);
    for s=1:length(sigs)
      for i=1:nderivs
        dG = filterDoog( dims, repmat(sigs(s),[1 3]), derivs(i,:), 0 );
        if(s==1 && i==1); FB=repmat(dG,[1 1 1 nderivs*length(sigs)]); end
        FB(:,:,:,cnt) = dG; cnt=cnt+1;
      end
    end

  otherwise
    error('none created.');
end

author: zouchuhang

repo: LayoutNet


text:
% CIRCLES3D  Create a data set containing 2 circles in 3 dimensions.
%
%   DATA = CIRCLES3D(N)
%
%	Creates a data set containing N points in 3 dimensions.
%
% If N is a vector of sizes, exactly N(I) objects are generated
% for class I, I = 1,2.Default: N = [50 50].
%
% See also DATASETS, PRDATASETS

% Copyright: E. Pekalska, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: circles3d.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function data = circles3d(N)
			if nargin< 1, N = [50 50]; end
	N = genclass(N,ones(1,2)/2);

	n2a = N(1);
	n2b = N(2);
	ha = 0:(2*pi/n2a):2*pi*(n2a/(n2a+1)); ha = ha';
	hb = 0:(2*pi/n2b):2*pi*(n2b/(n2b+1)); hb = hb';

	a = [ sin(ha) cos(ha) zeros(n2a,1) ];
	b = [ sin(hb) cos(hb) ones(n2b,1)  ];

	data = prdataset([a;b],genlab(N));
	data = setname(data,'3D Circles');

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
 classdef matRad_MinMaxDVH < DoseConstraints.matRad_DoseConstraint
    % matRad_MinMaxDVH Implements a MinMaxDVH constraint
    %   See matRad_DoseConstraint for interface description
    %
    % References
    %   -
    %
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Copyright 2020 the matRad development team.
    %
    % This file is part of the matRad project. It is subject to the license
    % terms in the LICENSE file found in the top-level directory of this
    % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
    % of the matRad project, including this file, may be copied, modified,
    % propagated, or distributed except according to the terms contained in the
    % LICENSE file.
    %
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'DVH constraint';
        parameterNames = {'d^{ref}', 'V^{min}', 'V^{max}'};
        %parameterIsDose = logical([1 0 0]);
        parameterTypes = {'dose','numeric','numeric'};
    end

    properties
        voxelScalingRatio = 1;
        referenceScalingVal = 0.01;
        parameters = {30,0,100};
    end

    methods
        function obj = matRad_MinMaxDVH(dRef,vMin,vMax)

            %If we have a struct in first argument
            if nargin == 1 && isstruct(dRef)
                inputStruct = dRef;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseConstraints.matRad_DoseConstraint(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct
                if nargin == 3 && isscalar(vMax)
                    obj.parameters{3} = vMax;
                end

                if nargin >= 1 && isscalar(dRef)
                    obj.parameters{1} = dRef;
                end

                if nargin >= 2 && isscalar(vMin)
                    obj.parameters{2} = vMin;
                end
            end
        end

        %Overloads the struct function to add constraint specific
        %parameters
        function s = struct(obj)
            s = struct@DoseConstraints.matRad_DoseConstraint(obj);
            s.voxelScalingRatio = 1;
            s.referenceScalingVal = 0.01;
        end

        function cu = upperBounds(obj,n)
            cu = obj.parameters{3} / 100;
        end
        function cl = lowerBounds(obj,n)
            cl = obj.parameters{2} / 100;
        end
        %% Calculates the Constraint Function value
        function cDose = computeDoseConstraintFunction(obj,dose)

            %Fast DVH point calculation
            cDose = sum(dose >= obj.parameters{1})/numel(dose);

            %cDose = 100 * cDose; %In Percent

            % alternative constraint calculation 3/4 %
            % % get reference Volume
            % refVol = cst{j,6}(k).volume/100;
            %
            % % calc deviation
            % deviation = d_i - d_ref;
            %
            % % calc d_ref2: V(d_ref2) = refVol
            % d_ref2 = matRad_calcInversDVH(refVol,d_i);
            %
            % % apply lower and upper dose limits
            % if isequal(cst{j,6}(k).type, 'max DVH constraint')
            %    deviation(d_i < d_ref | d_i > d_ref2) = 0;
            % elseif isequal(cst{j,6}(k).type, 'min DVH constraint')
            %    deviation(d_i > d_ref | d_i < d_ref2) = 0;
            % end
            %
            % %c = sum(deviation);                              % linear deviation
            % %c = deviation'*deviation;                        % square devioation
            % c = (1/size(cst{j,4},1))*(deviation'*deviation); % square deviation with normalization
            % %c = (deviation).^2'*(deviation).^2;               % squared square devioation
            % alternative constraint calculation 3/4 %
        end

        %% Calculates the Constraint jacobian
        function cDoseJacob  = computeDoseConstraintJacobian(obj,dose)
            %logistic approximation

            %Do we really need to sort two times?
            dose_sort = sort(dose);

            % calculate scaling
            NoVoxels     = max(obj.voxelScalingRatio*numel(dose),10);
            absDiffsort  = sort(abs(obj.parameters{1} - dose_sort));

            deltaDoseMax = absDiffsort(min(ceil(NoVoxels/2),numel(dose)));

            % calclulate DVHC scaling
            DVHCScaling = min((log(1/obj.referenceScalingVal-1))/(2*deltaDoseMax),250);

            d_diff = dose - obj.parameters{1};

            cDoseJacob = (2/numel(dose))*DVHCScaling*exp(2*DVHCScaling*d_diff)./(exp(2*DVHCScaling*d_diff)+1).^2;

            % alternative constraint calculation 4/4 %
            % % get reference Volume
            % refVol = cst{j,6}(k).volume/100;
            %
            % % calc deviation
            % deviation = d_i - d_ref;
            %
            % % calc d_ref2: V(d_ref2) = refVol
            % d_ref2 = matRad_calcInversDVH(refVol,d_i);
            %
            % % apply lower and upper dose limits
            % if isequal(cst{j,6}(k).type, 'max DVH constraint')
            %      deviation(d_i < d_ref | d_i > d_ref2) = 0;
            % elseif isequal(cst{j,6}(k).type, 'min DVH constraint')
            %      deviation(d_i > d_ref | d_i < d_ref2) = 0;
            % end
            %
            % %jacobVec = ones(size(cst{j,4}));             % linear deviation
            % %jacobVec = 2*deviation;                      % square deviation
            % jacobVec = (1/size(cst{j,4},1))*2*deviation; % square deviation with normalization
            % %jacobVec = 4*(deviation).^3;                  % squared square devioation
            % alternative constraint calculation 4/4 %
        end
    end

end

repo: matRad


text:
function [gJacobMat] = gJacob(p_pc_c, calibParams)
%GJACOBMAT Returns the 4x3 Jacobian of the observation model with respect
%to the feature position in the camera frame

%Redefine local shorthand notation
x = p_pc_c(1);
y = p_pc_c(2);
z = p_pc_c(3);

c_u = calibParams.c_u;
c_v = calibParams.c_v;
f_u = calibParams.f_u;
f_v = calibParams.f_v;
b = calibParams.b;

%Evaluate the Jacobian
gJacobMat = [f_u/z, 0, -(1/z^2)*(f_u*x )
             0, f_v/z, -(1/z^2)*(f_v*y )
             f_u/z, 0, -(1/z^2)*(f_u*(x-b) )
             0, f_v/z, -(1/z^2)*(f_v*y )];
end

repo: MSCKF


text:
%ReplicationInitiationDuration
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 9/12/2010
classdef ReplicationInitiationDuration
    methods (Static)
        %inputFilePath is a path to .mat files containing replication initiation
        %durations generated by
        %edu.stanford.covert.cell.sim.process.ReplicationInitiation_Test.meanDur
        %ation
        function run(inputFilePath, outputFileName)
            import edu.stanford.covert.cell.sim.analysis.ReplicationInitiationDuration;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %load data
            [durations, constants] = ReplicationInitiationDuration.loadDurations(inputFilePath);

            %statistics
            [vals, ~, means, stds, mins, maxs, meds, p10, p25, p75, p90] = ...
                ReplicationInitiationDuration.calculateStatistics(durations, constants);
            colLabels = {'Mean DnaA Copy Number', 'Site Cooperativity', 'State Cooperativity', 'No. DnaA Boxes', ...
                'Mean', 'Std. Dev', 'Min', '10%', '25%', 'Median', '75%', '90%', 'Max'};
            if nargin == 1
                PrintUtil.printToStdIO(num2cell([vals means stds mins p10 p25 meds p75 p90 maxs]), colLabels);
            else
                PrintUtil.printToFile(num2cell([vals means stds mins p10 p25 meds p75 p90 maxs]), ...
                    colLabels, [outputFileName '.xls'], 'Statistics')
            end

            %plots
            if nargin == 2, [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle(); end

            if nargin == 1, [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle(); end
            cla(axesHandle);
            tfs = ismember(constants(:, [1 3 4]), constants(1, [1 3 4]), 'rows');
            ReplicationInitiationDuration.plotScatterPlot(axesHandle, durations(tfs), constants(tfs, :));
            if nargin == 2, saveas(figHandle, [outputFileName '-ScatterPlot.pdf']); end

            if nargin == 1, [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle(); end
            cla(axesHandle);
            tfs = ismember(constants, constants(1, :), 'rows');
            ReplicationInitiationDuration.plotHistogram(axesHandle, durations(tfs), constants(tfs, :));
            if nargin == 2, saveas(figHandle, [outputFileName '-Histogram.pdf']); end

            if nargin == 1, [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle(); end
            cla(axesHandle);
            ReplicationInitiationDuration.plotDistributions(axesHandle, durations, constants);
            if nargin == 2, saveas(figHandle, [outputFileName '-Distributions.pdf']); end

            if nargin == 1, [axesHandle, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle(); end
            cla(axesHandle);
            ReplicationInitiationDuration.plotMeanVsStd(axesHandle, durations, constants);
            if nargin == 2, saveas(figHandle, [outputFileName '-MeanVsStd.pdf']); end

            if nargin == 2, close(figHandle); end
        end
    end

    methods (Static)
        function [durations, constants] = loadDurations(filePath)
            fileNames = dir([filePath filesep '*.mat']);
            durations = zeros(0, 1);
            constants = zeros(0, 4);
            path = fileparts([filePath filesep]);
            for i = 1:size(fileNames, 1)
                [~, parameters] = regexp(fileNames(i).name, ...
                    'ReplicationInitiationDuration-(?<meanDnaACopyNumber>\d+)-(?<siteCooperativity>\d+\.\d+)-(?<stateCooperativity>\d+\.\d+)-(?<nDnaABoxes>\d+)\.mat', 'match', 'names');

                data = load([path filesep fileNames(i).name]);
                durations = [durations; data.durations]; %#ok<AGROW>
                constants = [constants; repmat([...
                    str2double(parameters.meanDnaACopyNumber) ...
                    str2double(parameters.siteCooperativity) ...
                    str2double(parameters.stateCooperativity) ...
                    str2double(parameters.nDnaABoxes) ...
                    ], numel(data.durations), 1)]; %#ok<AGROW>
            end

            durations = durations / 3600;
        end

        function [vals, cnts, means, stds, mins, maxs, meds, p10, p25, p75, p90] = ...
                calculateStatistics(durations, constants)
            durationBins = 0:2:20;

            vals = unique(constants, 'rows');
            cnts = zeros(numel(durationBins), size(vals, 1));

            means = zeros(size(vals, 1), 1);
            stds = zeros(size(vals, 1), 1);
            mins = zeros(size(vals, 1), 1);
            maxs = zeros(size(vals, 1), 1);
            meds = zeros(size(vals, 1), 1);
            p25 = zeros(size(vals, 1), 1);
            p75 = zeros(size(vals, 1), 1);
            p90 = zeros(size(vals, 1), 1);
            p10 = zeros(size(vals, 1), 1);

            durations = min(durations, 3*23825/3600);

            for i = 1:size(vals, 1)
                tmpDurations = durations(ismember(constants, vals(i, :), 'rows'), 1);
                cnts(:,i) = histc(tmpDurations, durationBins);
                means(i) = mean(tmpDurations);
                stds(i) = std(tmpDurations);
                mins(i) = min(tmpDurations);
                maxs(i) = max(tmpDurations);
                meds(i) = median(tmpDurations);
                p25(i) = quantile(tmpDurations, 0.25);
                p75(i) = quantile(tmpDurations, 0.75);
                p90(i) = quantile(tmpDurations, 0.90);
                p10(i) = quantile(tmpDurations, 0.10);
            end
        end
    end

    %plots
    methods (Static)
        function plotScatterPlot(axesHandle, durations, constants)
            if isempty(axesHandle)
                axesHandle = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            end

            plot(axesHandle, constants(:, 2), durations, '.');
            xlabel(axesHandle, 'Site Cooperativity Constant', 'fontsize', 14);
            ylabel(axesHandle, 'Replication Initiation Time (h)', 'fontsize', 14);
            title(axesHandle, sprintf('DnaA %d, State %f, DnaA Boxes %d', ...
                constants(1), constants(3), constants(4)), 'fontsize', 14);
        end

        function plotHistogram(axesHandle, durations, constants)
            if isempty(axesHandle)
                axesHandle = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            end

            hist(axesHandle, durations);
            ylabel(axesHandle, 'Frequency', 'fontsize', 14);
            xlabel(axesHandle, 'Replication Initiation Time (h)', 'fontsize', 14);
            title(axesHandle, sprintf('DnaA %d, State %f, DnaA Boxes %d', ...
                constants(1), constants(3), constants(4)), 'fontsize', 14);
        end

        function plotDistributions(axesHandle, durations, constants)
            if isempty(axesHandle)
                axesHandle = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            end

            [vals, cnts] = edu.stanford.covert.cell.sim.analysis.ReplicationInitiationDuration.calculateStatistics(...
                durations, constants);

            colors = [
                1.00 0 0;
                0.75 0 0;
                0.50 0 0;
                0.25 0 0;
                0 0 0;
                0 0 0.25;
                0 0 0.50;
                0 0 0.75;
                0 0 1];

            h = plot(axesHandle, (0:2:20)', cnts);
            for i = 1:size(vals, 1)
                set(h(i), 'color', colors(mod(i-1, size(colors, 1))+1, :));
            end
            legend(h, cellfun(@(val1, val2, val3, val4) sprintf('%d %.2f %.2f %d', val1, val2, val3, val4), ...
                num2cell(vals(:, 1)), ...
                num2cell(vals(:, 2)), ...
                num2cell(vals(:, 3)), ...
                num2cell(vals(:, 4)), ...
                'UniformOutput', false), 'Location', 'BestOutside');
            ylabel(axesHandle, 'Frequency', 'fontsize', 14);
            xlabel(axesHandle, 'Replication Initiation Time (h)', 'fontsize', 14);
            ylim(axesHandle, [0 max(cnts(:)) + eps]);
        end

        function plotMeanVsStd(axesHandle, durations, constants)
            if isempty(axesHandle)
                axesHandle = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            end

            [vals, ~, means, stds] = edu.stanford.covert.cell.sim.analysis.ReplicationInitiationDuration.calculateStatistics(...
                durations, constants);

            vals3 = unique(vals(:, 3));
            markerSizes = repmat([10 18 8 12 24]', ceil(numel(vals3)/5), 1);
            markerSizes = markerSizes(1:numel(vals3));
            vals2 = unique(vals(:, 2));
            markerStyles = repmat(('.+*oxsdv^<')', ceil(numel(vals2)/10), 1);
            markerStyles = markerStyles(1:numel(vals2));
            vals1 = unique(vals(:, 1));
            colors = repmat([
                1 0 0
                0 0 0
                0 0 1
                0.75 0 0
                0.50 0 0
                0.25 0 0
                0 0 0.25
                0 0 0.50
                0 0 0.75
                ], ceil(numel(vals1)/9), 1);
            colors = colors(1:numel(vals1), :);

            hold(axesHandle, 'on');
            for i = 1:size(vals, 1)
                h(i) = plot(axesHandle, means(i), stds(i), '.', ...
                    'MarkerSize', markerSizes(ismember(vals3, vals(i, 3))), ...
                    'Marker', markerStyles(ismember(vals2, vals(i, 2))), ...
                    'Color', colors(ismember(vals1, vals(i, 1)), :));
            end
            legend(h, cellfun(@(val1, val2, val3, val4) sprintf('%d %.2f %.2f %d', val1, val2, val3, val4), ...
                num2cell(vals(:, 1)), ...
                num2cell(vals(:, 2)), ...
                num2cell(vals(:, 3)), ...
                num2cell(vals(:, 4)), ...
                'UniformOutput', false), 'Location', 'BestOutside');
            xlim(axesHandle, [min(means) max(means)]+[-1 1]*0.1*range(means));
            ylim(axesHandle, [min(stds) max(stds)]+[-1 1]*0.1*range(stds));
            ylabel(axesHandle, 'Duration standard deviation (h)', 'fontsize', 14);
            xlabel(axesHandle, 'Mean duration (h)', 'fontsize', 14);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function hypercube_grid_test03 ( )

%*****************************************************************************80
%
%% HYPERCUBE_GRID_TEST03 tests HYPERCUBE_GRID on a three dimensional example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 August 2014
%
%  Author:
%
%    John Burkardt
%
  m = 3;

  a = [ -1.0, -1.0, -1.0 ];
  b = [ +1.0, +1.0, +1.0 ];
  c = [ 1, 1, 1 ];
  ns = [ 3, 3, 3 ];

  n = prod ( ns(1:m) );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HYPERCUBE_GRID_TEST03\n' );
  fprintf ( 1, '  Create a grid using HYPERCUBE_GRID.\n' );
  fprintf ( 1, '  Use the same parameters in every dimension.\n' );
  fprintf ( 1, '  Spatial dimension M = %d\n', m );
  fprintf ( 1, '  Number of grid points N = %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    NS     C      A         B\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '  %4d  %4d  %4d  %8.4f  %8.4f\n', i, ns(i), c(i), a(i), b(i) );
  end

  x = hypercube_grid ( m, n, ns, a, b, c );
  r8mat_transpose_print ( m, n, x, '  Grid points:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function isUniformized = isUniformStr(structNum,planC)
%function isUniformized = isUniformStr(structnum,planC)
%
%This function checks whether the strucutre structNum has been uniformized
%
%APA, 04/23/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC','var')
    global planC
end
indexS = planC{end};
[scanSet,relStrNum] = getStructureAssociatedScan(structNum, planC);

if relStrNum<=52
    cellNum = 1;
    bitsM = planC{indexS.structureArray}(scanSet).bitsArray;
    bitGetM = bitget(bitsM, structNum);
else
    cellNum = ceil((relStrNum-52)/8)+1;
    bitsM = planC{indexS.structureArrayMore}(scanSet).bitsArray{cellnum};
    bitGetM = bitget(bitsM, relStrNum-52-(cellNum-2)*8);
end

if any(bitGetM)
    isUniformized = 1;
else
    isUniformized = 0;
end

return;


text:
function FeatureDetection(wkdir, dataset)
% Detect and save DoG keypoints

disp('Detecting keypoints...');

dataset_dir = [wkdir 'Dataset/' dataset '/'];

feature_dir = [wkdir 'Features/' dataset '/'];
if exist(feature_dir, 'dir') == 0
    mkdir(feature_dir);
end

pairs_gts = dlmread([dataset_dir 'pairs_with_gt.txt']);
pairs_which_dataset = importdata([dataset_dir 'pairs_which_dataset.txt']);

pairs = pairs_gts(:,1:2);
l_pairs = pairs(:,1);
r_pairs = pairs(:,2);

num_pairs = size(pairs,1);
for idx = 1 : num_pairs
    l = l_pairs(idx);
    r = r_pairs(idx);

    I1 = imread([dataset_dir pairs_which_dataset{idx} 'Images/' sprintf('%.8d.jpg', l)]);
    I2 = imread([dataset_dir pairs_which_dataset{idx} 'Images/' sprintf('%.8d.jpg', r)]);

    if size(I1,3) == 3
        I1gray = rgb2gray(I1);
    else
        I1gray = I1;
    end

    if size(I2,3) == 3
        I2gray = rgb2gray(I2);
    else
        I2gray = I2;
    end

    path_l = [feature_dir sprintf('%.4d_l', idx)];
    path_r = [feature_dir sprintf('%.4d_r', idx)];

    if exist([path_l '.keypoints'], 'file') == 2 && exist([path_r '.keypoints'], 'file') == 2
        continue;
    end

    keypoints_l = vl_sift(single(I1gray))';
    write_keypoints([path_l '.keypoints'], keypoints_l);

    keypoints_r = vl_sift(single(I2gray))';
    write_keypoints([path_r '.keypoints'], keypoints_r);
end

disp('Finished.');
end

author: JiawangBian

repo: FM-Bench


text:
function an=random_unit_vector(varargin)
% random_unit_vector
% random_unit_vector(n)
% random_unit_vector(m,n)
% random_unit_vector([m n])
% random_unit_vector('double')
% random_unit_vector(n,'double')
% random_unit_vector(m,n,'double')
% random_unit_vector([m n],'double')
% random_unit_vector('single')
% random_unit_vector(n,'single')
% random_unit_vector(m,n,'single')
% random_unit_vector([m n],'single')

% m - dimentionarity
% n - number of unit vectors

md=3; % default m

isg=false; % if single
switch nargin
    case 0
        % if no inputs
        n=1;
        m=md;
    case 1
        %
        i1=varargin{1};
        if ischar(i1)
            if strcmpi(i1,'single')
                isg=true;
            end
        else
            if length(i1)==1
                m=md;
                n=i1;
            else
                m=i1(1);
                n=i1(2);
            end
        end
    case 2
        i1=varargin{1};
        i2=varargin{2};
        if ischar(i2)
            if length(i1)==1
                m=md;
                n=i1;
            else
                m=i1(1);
                n=i1(2);
            end
            if strcmpi(i2,'single')
                isg=true;
            end
        else
            m=i1;
            n=i2;
        end
    case 3
        m=varargin{1};
        n=varargin{2};
        if strcmpi(varargin{3},'single')
            isg=true;
        end

end

% simple case of 1d
if m==1
    if isg
        an=single(2*(randn(1,n)>0)-1);
    else
        an=2*(randn(1,n)>0)-1;
    end
    return;
end

if isg
    v=randn(m,n,'single');
else
    v=randn(m,n);
end

% normalize:
if isg
    an=zeros(m,n,'single');
else
    an=zeros(m,n);
end
for nc=1:n
    while 1
        v2=v(:,nc)'*v(:,nc);
        if v2>1e-10 % too small values must be excluded
            % because it will have discretization errors
            an(:,nc)=v(:,nc)/sqrt(v2);
            break;
        else
            % otherwise repeat random generation
            if isg
                v(:,nc)=randn(m,1,'single');
            else
                v(:,nc)=randn(m,1);
            end
        end
    end
end


text:
function ExecFlag=DorfExec(handles)

global VObj;
global VMag;
global VVar;
global VSeq;
global VCtl;
global VSig;
global VCoi;

%preserve VObj VMag & VCoi
VTmpObj=VObj;
VTmpMag=VMag;
VTmpCoi=VCoi;

handles.Simuh=guidata(handles.Simuh.SimuPanel_figure);
DoDisableButton([],[],handles.Simuh);
%% Do spin execution
try
    % Read tab parameters
    fieldname=fieldnames(handles.Attrh0);
    for i=1:length(fieldname)/2
        switch get(handles.Attrh0.(fieldname{i*2}),'Style')
            case 'edit'
                eval(['RF.' fieldname{i*2} '=[' get(handles.Attrh0.(fieldname{i*2}),'String') '];']);
            case 'popupmenu'
                TAttr=get(handles.Attrh0.(fieldname{i*2}),'String');
                eval(['RF.' fieldname{i*2} '=''' TAttr{get(handles.Attrh0.(fieldname{i*2}),'Value')}  ''';']);
        end
    end

    % Prescan config
    DoPreScan(handles.Simuh);

    % Create SpinWatcher VOrf
    if strcmp(RF.Freq_Flag,'on') % spatial-spectral analysis
        if strcmp(RF.Spat_Flag,'on') % 2D spatial analysis
            error('Freq_Flag is incompatible with Spat_Flag! Turn Spat_Flag to off.');
        end
        set(handles.SliProfile_uipanel,'Title','Spin Response : Spatial-Spectral');
        Freq=linspace(RF.FreqDownLimit,RF.FreqUpLimit,RF.FreqRes)/VCtl.B0;

        VOrf=VObj;
        VOrf.Gyro=RF.Gyro;
        VOrf.ChemShift=Freq;
        VOrf.SpinNum=1;
        VOrf.TypeNum=length(Freq);
        VOrf.XDim=1;
        VOrf.YDim=1;
        VOrf.ZDim=RF.ZSpin;
        VOrf.XDimRes=1e-3;
        VOrf.YDimRes=1e-3;
        VOrf.ZDimRes=RF.ZSpinGap;
        VOrf=rmfield(VOrf, 'Rho');
        VOrf=rmfield(VOrf, 'T1');
        VOrf=rmfield(VOrf, 'T2');
        VOrf=rmfield(VOrf, 'T2Star');
        VOrf=rmfield(VOrf, 'Type');
        for i=1:VOrf.TypeNum
            VOrf.Rho(:,:,1:RF.ZSpin,i)=RF.Rho(1);
            VOrf.T1(:,:,1:RF.ZSpin,i)=RF.T1(1);
            VOrf.T2(:,:,1:RF.ZSpin,i)=RF.T2(1);
            VOrf.T2Star(:,:,1:RF.ZSpin,i)=0;
        end
        VOrf=rmfield(VOrf, 'Mx');
        VOrf=rmfield(VOrf, 'My');
        VOrf=rmfield(VOrf, 'Mz');
        VOrf.Mx= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.My= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.Mz= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        for i=1:VOrf.TypeNum
            VOrf.Mz(:,:,:,:,i)= ones(1,1,RF.ZSpin,VOrf.SpinNum,1).*(RF.Rho(1)./VOrf.SpinNum);
        end
    elseif strcmp(RF.Spat_Flag,'on') % 2D spatial analysis
        set(handles.SliProfile_uipanel,'Title','Spin Response : 2D Spatial');
        VOrf=VObj;
        VOrf.Gyro=RF.Gyro;
        VOrf.ChemShift=RF.ChemShift(1);
        VOrf.SpinNum=1;
        VOrf.TypeNum=1;
        VOrf.XDim=RF.XSpin;
        VOrf.YDim=RF.YSpin;
        VOrf.ZDim=1;
        VOrf.XDimRes=RF.XSpinGap;
        VOrf.YDimRes=RF.YSpinGap;
        VOrf.ZDimRes=1e-3;
        VOrf=rmfield(VOrf, 'Rho');
        VOrf=rmfield(VOrf, 'T1');
        VOrf=rmfield(VOrf, 'T2');
        VOrf=rmfield(VOrf, 'T2Star');
        VOrf=rmfield(VOrf, 'Type');

        VOrf.Rho(1:RF.YSpin,1:RF.XSpin,1,1)=RF.Rho(1);
        VOrf.T1(1:RF.YSpin,1:RF.XSpin,1,1)=RF.T1(1);
        VOrf.T2(1:RF.YSpin,1:RF.XSpin,1,1)=RF.T2(1);
        VOrf.T2Star(1:RF.YSpin,1:RF.XSpin,1,1)=0;

        VOrf=rmfield(VOrf, 'Mx');
        VOrf=rmfield(VOrf, 'My');
        VOrf=rmfield(VOrf, 'Mz');
        VOrf.Mx= zeros(VOrf.YDim,VOrf.XDim,1,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.My= zeros(VOrf.YDim,VOrf.XDim,1,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.Mz= zeros(VOrf.YDim,VOrf.XDim,1,VOrf.SpinNum,VOrf.TypeNum);
        for i=1:VOrf.TypeNum
            VOrf.Mz(:,:,:,:,i)= ones(VOrf.YDim,VOrf.XDim,1,VOrf.SpinNum,1).*(RF.Rho(i)./VOrf.SpinNum);
        end
    else  % normal 1D analysis
        set(handles.SliProfile_uipanel,'Title','Spin Response : 1D Spatial');
        VOrf=VObj;
        VOrf.Gyro=RF.Gyro;
        VOrf.ChemShift=RF.ChemShift;
        VOrf.SpinNum=1;
        VOrf.TypeNum=RF.TypeNum;
        VOrf.XDim=1;
        VOrf.YDim=1;
        VOrf.ZDim=RF.ZSpin;
        VOrf.XDimRes=1e-3;
        VOrf.YDimRes=1e-3;
        VOrf.ZDimRes=RF.ZSpinGap;
        VOrf=rmfield(VOrf, 'Rho');
        VOrf=rmfield(VOrf, 'T1');
        VOrf=rmfield(VOrf, 'T2');
        VOrf=rmfield(VOrf, 'T2Star');
        VOrf=rmfield(VOrf, 'Type');
        for i=1:VOrf.TypeNum
            VOrf.Rho(:,:,1:RF.ZSpin,i)=RF.Rho(i);
            VOrf.T1(:,:,1:RF.ZSpin,i)=RF.T1(i);
            VOrf.T2(:,:,1:RF.ZSpin,i)=RF.T2(i);
            VOrf.T2Star(:,:,1:RF.ZSpin,i)=0;
        end
        VOrf=rmfield(VOrf, 'Mx');
        VOrf=rmfield(VOrf, 'My');
        VOrf=rmfield(VOrf, 'Mz');
        VOrf.Mx= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.My= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        VOrf.Mz= zeros(1,1,VOrf.ZDim,VOrf.SpinNum,VOrf.TypeNum);
        for i=1:VOrf.TypeNum
            VOrf.Mz(:,:,:,:,i)= ones(1,1,RF.ZSpin,VOrf.SpinNum,1).*(RF.Rho(i)./VOrf.SpinNum);
        end
    end

    if strcmp(RF.Spat_Flag,'on') % 2D analysis
        % Gradient Grid
        VMrf=VMag;
        VMrf=rmfield(VMrf, 'Gxgrid');
        VMrf=rmfield(VMrf, 'Gygrid');
        VMrf=rmfield(VMrf, 'Gzgrid');
        VMrf=rmfield(VMrf, 'dB0');
        VMrf=rmfield(VMrf, 'dWRnd');
        VMrf=rmfield(VMrf, 'FRange');

        [VMrf.Gxgrid,VMrf.Gygrid,VMrf.Gzgrid]=meshgrid((-RF.XCenter+1)*RF.XSpinGap:RF.XSpinGap:(RF.XSpin-RF.XCenter)*RF.XSpinGap, ...
                                                       (-RF.YCenter+1)*RF.YSpinGap:RF.YSpinGap:(RF.YSpin-RF.YCenter)*RF.YSpinGap,1);
        VMrf.dB0=ones([RF.YSpin RF.XSpin 1])*RF.dB0;
        VMrf.dWRnd=zeros(size(VOrf.Mz));
        VMrf.FRange=ones(size(VOrf.Rho));

        % Coil
        VCrf=VCoi;
        VCrf.TxCoilmg  = ones([RF.YSpin RF.XSpin 1]);  % Unit B1+ for rf analysis
        VCrf.TxCoilpe  = zeros([RF.YSpin RF.XSpin 1]);
        VCrf.TxCoilNum = 1;

    else % spatial-spectral analysis or 1D analysis
        % Gradient Grid
        VMrf=VMag;
        VMrf=rmfield(VMrf, 'Gxgrid');
        VMrf=rmfield(VMrf, 'Gygrid');
        VMrf=rmfield(VMrf, 'Gzgrid');
        VMrf=rmfield(VMrf, 'dB0');
        VMrf=rmfield(VMrf, 'dWRnd');
        VMrf=rmfield(VMrf, 'FRange');

        [VMrf.Gxgrid,VMrf.Gygrid,VMrf.Gzgrid]=meshgrid(0,0,(-RF.ZCenter+1)*RF.ZSpinGap:RF.ZSpinGap:(RF.ZSpin-RF.ZCenter)*RF.ZSpinGap); % assume constant unit gradient in GzSS
        VMrf.dB0=ones([1 1 RF.ZSpin])*RF.dB0;
        VMrf.dWRnd=zeros(size(VOrf.Mz));
        VMrf.FRange=ones(size(VOrf.Rho));

        % Coil
        VCrf=VCoi;
        VCrf.TxCoilmg  = ones([1 1 RF.ZSpin]);  % Unit B1+ for rf analysis
        VCrf.TxCoilpe  = zeros([1 1 RF.ZSpin]);
        VCrf.TxCoilNum = 1;

    end
    %% Spin execution
    VObj=VOrf;
    VMag=VMrf;
    VCoi=VCrf;

    % Generate Pulse line
    if isfield(handles,'rfNode')
        for parai=1:length(handles.rfNode.Attributes)
            if ~isempty(handles.rfNode.Attributes(parai).Value)
                switch handles.rfNode.Attributes(parai).Value(1)
                    case '$'
                        eval(['AttributeOpt={' handles.rfNode.Attributes(parai).Value(3:end) '};']);
                        eval(['p.' handles.rfNode.Attributes(parai).Name '=AttributeOpt{' handles.rfNode.Attributes(parai).Value(2) '};']);
                    otherwise
                        if ~strcmp(handles.rfNode.Attributes(parai).Name,'Notes')
                            eval(['p.' handles.rfNode.Attributes(parai).Name '=' handles.rfNode.Attributes(parai).Value ';']);
                        else
                            eval(['p.' handles.rfNode.Attributes(parai).Name '=''' handles.rfNode.Attributes(parai).Value ''';']);
                        end
                end

            else
                eval(['p.' handles.rfNode.Attributes(parai).Name '=[];']);
            end
        end
        VVar.TRCount=1; % may used for rfUser
        eval(['[rfAmp,rfPhase,rfFreq,rfCoil,rfTime]=' handles.rfNode.Name '(p);']);
    else
        errordlg('Please choose a rf wave type (e.g. rfSinc)!');
        ExecFlag=0;
        %recover VObj VMag VCoi
        VObj=VTmpObj;
        VMag=VTmpMag;
        VCoi=VTmpCoi;
        return
    end

    if strcmp(RF.Spat_Flag,'off')
        if isfield(handles,'GzNode')
            for parai=1:length(handles.GzNode.Attributes)
                if ~isempty(handles.GzNode.Attributes(parai).Value)
                    switch handles.GzNode.Attributes(parai).Value(1)
                        case '$'
                            eval(['AttributeOpt={' handles.GzNode.Attributes(parai).Value(3:end) '};']);
                            eval(['p.' handles.GzNode.Attributes(parai).Name '=AttributeOpt{' handles.GzNode.Attributes(parai).Value(2) '};']);
                        otherwise
                            if ~strcmp(handles.GzNode.Attributes(parai).Name,'Notes')
                                eval(['p.' handles.GzNode.Attributes(parai).Name '=' handles.GzNode.Attributes(parai).Value ';']);
                            else
                                eval(['p.' handles.GzNode.Attributes(parai).Name '=''' handles.GzNode.Attributes(parai).Value ''';']);
                            end
                    end

                else
                    eval(['p.' handles.GzNode.Attributes(parai).Name '=[];']);
                end
            end
            VVar.TRCount=1; % may used for GzUser
            eval(['[GzAmp,GzTime]=' handles.GzNode.Name '(p);']);
        else
            GzAmp=[0 RF.ConstantGrad RF.ConstantGrad 0];
            GzTime=[0 VCtl.MinUpdRate max(rfTime)-VCtl.MinUpdRate max(rfTime)];
            warndlg('Please choose a GzSS gradient type (e.g. GzSelective), otherwise constant gradient (i.e. ConstantGrad) will be used!');
        end

    else
        if isfield(handles,'GyNode')
            for parai=1:length(handles.GyNode.Attributes)
                if ~isempty(handles.GyNode.Attributes(parai).Value)
                    switch handles.GyNode.Attributes(parai).Value(1)
                        case '$'
                            eval(['AttributeOpt={' handles.GyNode.Attributes(parai).Value(3:end) '};']);
                            eval(['p.' handles.GyNode.Attributes(parai).Name '=AttributeOpt{' handles.GyNode.Attributes(parai).Value(2) '};']);
                        otherwise
                            if ~strcmp(handles.GyNode.Attributes(parai).Name,'Notes')
                                eval(['p.' handles.GyNode.Attributes(parai).Name '=' handles.GyNode.Attributes(parai).Value ';']);
                            else
                                eval(['p.' handles.GyNode.Attributes(parai).Name '=''' handles.GyNode.Attributes(parai).Value ''';']);
                            end
                    end

                else
                    eval(['p.' handles.GyNode.Attributes(parai).Name '=[];']);
                end
            end
            VVar.TRCount=1; % may used for GyUser
            eval(['[GyAmp,GyTime]=' handles.GyNode.Name '(p);']);
        else
            GyAmp=[0 RF.ConstantGrad RF.ConstantGrad 0];
            GyTime=[0 VCtl.MinUpdRate max(rfTime)-VCtl.MinUpdRate max(rfTime)];
            warndlg('Please choose a GyPE gradient type, otherwise constant gradient (i.e. ConstantGrad) will be used!');
        end

        if isfield(handles,'GxNode')
            for parai=1:length(handles.GxNode.Attributes)
                if ~isempty(handles.GxNode.Attributes(parai).Value)
                    switch handles.GxNode.Attributes(parai).Value(1)
                        case '$'
                            eval(['AttributeOpt={' handles.GxNode.Attributes(parai).Value(3:end) '};']);
                            eval(['p.' handles.GxNode.Attributes(parai).Name '=AttributeOpt{' handles.GxNode.Attributes(parai).Value(2) '};']);
                        otherwise
                            if ~strcmp(handles.GxNode.Attributes(parai).Name,'Notes')
                                eval(['p.' handles.GxNode.Attributes(parai).Name '=' handles.GxNode.Attributes(parai).Value ';']);
                            else
                                eval(['p.' handles.GxNode.Attributes(parai).Name '=''' handles.GxNode.Attributes(parai).Value ''';']);
                            end
                    end

                else
                    eval(['p.' handles.GxNode.Attributes(parai).Name '=[];']);
                end
            end
            VVar.TRCount=1; % may used for GxUser
            eval(['[GxAmp,GxTime]=' handles.GxNode.Name '(p);']);
        else
            GxAmp=[0 RF.ConstantGrad RF.ConstantGrad 0];
            GxTime=[0 VCtl.MinUpdRate max(rfTime)-VCtl.MinUpdRate max(rfTime)];
            warndlg('Please choose a GxR gradient type, otherwise constant gradient (i.e. ConstantGrad) will be used!');
        end
    end

    if strcmp(RF.Spat_Flag,'on') % 2D analysis
        SEt=[0 max([rfTime, GxTime, GyTime])]; % SEt time landmark must start from 0
        SEflag=repmat([0 0 0 0 0 0]',[1 2]);
        rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTime))]);
        Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTime))]);
        Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTime))]);

        ts=[SEt rfTime GyTime GxTime];
        flags=[SEflag rfflag Gyflag Gxflag];
        [ts,ind]=sort(ts);
        uts=unique(ts);
        flags=flags(:,ind);

        rfAmp(abs(rfAmp)<eps)=0;
        VSeq.rfAmpLine=rfAmp;
        VSeq.rfPhaseLine=rfPhase;
        VSeq.rfFreqLine=rfFreq;
        VSeq.rfCoilLine=ones(size(rfCoil)); % default single-Tx
        VSeq.GzAmpLine=0;
        VSeq.GyAmpLine=GyAmp;
        VSeq.GxAmpLine=GxAmp;
        VSeq.ADCLine=0;
        VSeq.ExtLine=0;
        VSeq.utsLine=uts;
        VSeq.tsLine=ts;
        VSeq.flagsLine=flags;

        % Initialize Mx, My, Mz, Muts in VSig
        VSig.Mx= zeros(VObj.YDim,VObj.XDim,1,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.My= zeros(VObj.YDim,VObj.XDim,1,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.Mz= zeros(VObj.YDim,VObj.XDim,1,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.Muts=zeros(max(size(VSeq.utsLine)),1);

    else % spatial-spectral analysis or 1D analysis

        SEt=[0 max(max(rfTime),max(GzTime))]; % SEt time landmark must start from 0
        SEflag=repmat([0 0 0 0 0 0]',[1 2]);
        rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTime))]);
        Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTime))]);

        ts=[SEt rfTime GzTime];
        flags=[SEflag rfflag Gzflag];
        [ts,ind]=sort(ts);
        uts=unique(ts);
        flags=flags(:,ind);

        rfAmp(abs(rfAmp)<eps)=0;
        VSeq.rfAmpLine=rfAmp;
        VSeq.rfPhaseLine=rfPhase;
        VSeq.rfFreqLine=rfFreq;
        VSeq.rfCoilLine=ones(size(rfCoil)); % default single-Tx
        VSeq.GzAmpLine=GzAmp;
        VSeq.GyAmpLine=0;
        VSeq.GxAmpLine=0;
        VSeq.ADCLine=0;
        VSeq.ExtLine=0;
        VSeq.utsLine=uts;
        VSeq.tsLine=ts;
        VSeq.flagsLine=flags;

        % Initialize Mx, My, Mz, Muts in VSig
        VSig.Mx= zeros(1,1,VObj.ZDim,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.My= zeros(1,1,VObj.ZDim,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.Mz= zeros(1,1,VObj.ZDim,1,VObj.TypeNum, max(size(VSeq.utsLine)));
        VSig.Muts=zeros(max(size(VSeq.utsLine)),1);

    end

    VSig.Mx(:,:,:,:,:,1)=VObj.Mx;
    VSig.My(:,:,:,:,:,1)=VObj.My;
    VSig.Mz(:,:,:,:,:,1)=VObj.Mz;
    VSig.Muts(1)=0;

    % Simulation Process
    try
        VCtl.CS=double(VObj.ChemShift*VCtl.B0);
        VCtl.RunMode=int32(1); % Spin scan
        VCtl.MaxThreadNum=int32(handles.Simuh.CPUInfo.NumThreads);
        DoDataTypeConv(handles.Simuh);
        DoScanAtCPU; % global (VSeq,VObj,VCtl,VMag,VCoi,VVar,VSig) are needed
    catch me
        error_msg{1,1}='ERROR!!! rf execution process aborted.';
        error_msg{2,1}=me.message;
        errordlg(error_msg);
        ExecFlag=0;
        %recover VObj VMag VCoi
        VObj=VTmpObj;
        VMag=VTmpMag;
        VCoi=VTmpCoi;
        return;
    end

    %modify sequence line for display purpose
    if strcmp(RF.Spat_Flag,'on') % 2D analysis
        handles.rfAmp=[0 rfAmp 0];
        handles.rfPhase=[0 rfPhase 0];
        handles.rfFreq=[0 rfFreq 0];
        handles.rfTime=[SEt(1) rfTime SEt(2)];
        handles.GyAmp=[0 GyAmp 0];
        handles.GyTime=[SEt(1) GyTime SEt(2)];
        handles.GxAmp=[0 GxAmp 0];
        handles.GxTime=[SEt(1) GxTime SEt(2)];
    else % spatial-spectral analysis or 1D analysis
        handles.rfAmp=[0 rfAmp 0];
        handles.rfPhase=[0 rfPhase 0];
        handles.rfFreq=[0 rfFreq 0];
        handles.rfTime=[SEt(1) rfTime SEt(2)];
        handles.GzAmp=[0 GzAmp 0];
        handles.GzTime=[SEt(1) GzTime SEt(2)];
    end

    if strcmp(RF.Freq_Flag,'on')
        handles.Mx=permute(VSig.Mx,[3, 5, 6, 1, 2, 4]).*double(VObj.SpinNum);
        handles.My=permute(VSig.My,[3, 5, 6, 1, 2, 4]).*double(VObj.SpinNum);
        handles.Mz=permute(VSig.Mz,[3, 5, 6, 1, 2, 4]).*double(VObj.SpinNum);
        handles.Mx=handles.Mx(:,:,:,1,1,1);
        handles.My=handles.My(:,:,:,1,1,1);
        handles.Mz=handles.Mz(:,:,:,1,1,1);
        handles.Muts=VSig.Muts;
        handles.Freq=Freq;
    elseif strcmp(RF.Spat_Flag,'on')
        handles.Mx=permute(VSig.Mx,[1, 2, 6, 3, 4, 5]).*double(VObj.SpinNum);
        handles.My=permute(VSig.My,[1, 2, 6, 3, 4, 5]).*double(VObj.SpinNum);
        handles.Mz=permute(VSig.Mz,[1, 2, 6, 3, 4, 5]).*double(VObj.SpinNum);
        handles.Mx=handles.Mx(:,:,:,1,1,1);
        handles.My=handles.My(:,:,:,1,1,1);
        handles.Mz=handles.Mz(:,:,:,1,1,1);
        handles.Muts=VSig.Muts;
    else
        handles.Mx=VSig.Mx.*double(VObj.SpinNum);
        handles.My=VSig.My.*double(VObj.SpinNum);
        handles.Mz=VSig.Mz.*double(VObj.SpinNum);
        handles.Muts=VSig.Muts;
    end

    handles.Freq_Flag=RF.Freq_Flag;
    handles.Spat_Flag=RF.Spat_Flag;
    handles.Gxgrid=VMag.Gxgrid;
    handles.Gygrid=VMag.Gygrid;
    handles.Gzgrid=VMag.Gzgrid;
catch me
    error_msg{1,1}='ERROR!!! rf execution process aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    ExecFlag=0;
    %recover VObj VMag VCoi
    VObj=VTmpObj;
    VMag=VTmpMag;
    VCoi=VTmpCoi;
    return;
end

%recover VObj VMag VCoi
VObj=VTmpObj;
VMag=VTmpMag;
VCoi=VTmpCoi;

guidata(handles.rfDesignPanel_figure, handles);
ExecFlag=1;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%% FUNCTION Least_TGL
% L21 Joint Feature Learning with Least Squares Loss.
%
%% OBJECTIVE
% argmin_W { sum_i^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%            + opts.rho_L2 * \|W\|_2^2 + rho1 * \|W\|_{2,1} }
%
%% INPUT
% X: {n * d} * t - input matrix
% Y: {n * 1} * t - output matrix
% rho1: L2,1-norm group Lasso parameter.
% optional:
%   opts.rho_L2: L2-norm parameter (default = 0).
%
%% OUTPUT
% W: model: d * t
% funcVal: function value vector.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Evgeniou, A. and Pontil, M. Multi-task feature learning, NIPS 2007.
%   [2] Liu, J. and Ye, J. Efficient L1/Lq Norm Regularization, Technical
%       Report, 2010.
%
%% RELATED FUNCTIONS
%  Least_L21, init_opts

%% Code starts here
function [W, funcVal] = Least_L21(X, Y, rho1, opts)

if nargin <3
    error('\n Inputs: X, Y, rho1, should be specified!\n');
end
X = multi_transpose(X);

if nargin <4
    opts = [];
end

% initialize options.
opts=init_opts(opts);

if isfield(opts, 'rho_L2')
    rho_L2 = opts.rho_L2;
else
    rho_L2 = 0;
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

% initialize a starting point
if isfield(opts,'W0')
    W0=opts.W0;
    if (nnz(size(W0)-[dimension, task_num]))
        error('\n Check the input .W0');
    end
elseif opts.init==2
    W0 = zeros(dimension, task_num);
elseif opts.init == 0
    XY = cell(task_num, 1);
    W0_prep = [];
    for t_idx = 1: task_num
        XY{t_idx} = X{t_idx}*Y{t_idx};
        W0_prep = cat(2, W0_prep, XY{t_idx});
    end
    W0 = W0_prep;
end

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

Wz= W0;
Wz_old = W0;

t = 1;
t_old = 0;

iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;

    % compute function value and gradients of the search point
    gWs  = gradVal_eval(Ws);
    Fs   = funVal_eval (Ws);

    while true
        Wzp = FGLasso_projection(Ws - gWs/gamma, rho1 / gamma);
        Fzp = funVal_eval  (Wzp);

        delta_Wzp = Wzp - Ws;
        r_sum = norm(delta_Wzp, 'fro')^2;
        %         Fzp_gamma = Fs + trace(delta_Wzp' * gWs)...
        %             + gamma/2 * norm(delta_Wzp, 'fro')^2;
        Fzp_gamma = Fs + sum(sum(delta_Wzp.* gWs))...
            + gamma/2 * norm(delta_Wzp, 'fro')^2;

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

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Wz = Wzp;

    funcVal = cat(1, funcVal, Fzp + nonsmooth_eval(Wz, rho1));

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

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;

% private functions

    function [X] = FGLasso_projection (D, lambda )
        % l2.1 norm projection.
        X = repmat(max(0, 1 - lambda./sqrt(sum(D.^2,2))),1,size(D,2)).*D;
    end

% smooth part gradient.
    function [grad_W] = gradVal_eval(W)
        if opts.pFlag
            grad_W = zeros(zeros(W));
            parfor i = 1:task_num
                grad_W (i, :) = X{i}*(X{i}' * W(:,i)-Y{i});
            end
        else
            grad_W = [];
            for i = 1:task_num
                grad_W = cat(2, grad_W, X{i}*(X{i}' * W(:,i)-Y{i}) );
            end
        end
        grad_W = grad_W+ rho_L2 * 2 * W;
    end

% smooth part function value.
    function [funcVal] = funVal_eval (W)
        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        else
            for i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        end
        funcVal = funcVal + rho_L2 * norm(W,'fro')^2;
    end

    function [non_smooth_value] = nonsmooth_eval(W, rho_1)
        non_smooth_value = 0;
        if opts.pFlag
            parfor i = 1 : size(W, 1)
                w = W(i, :);
                non_smooth_value = non_smooth_value ...
                    + rho_1 * norm(w, 2);
            end
        else
            for i = 1 : size(W, 1)
                w = W(i, :);
                non_smooth_value = non_smooth_value ...
                    + rho_1 * norm(w, 2);
            end
        end
    end
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function c8_sqrt_test ( )

%*****************************************************************************80
%
%% C8_SQRT_TEST tests C8_SQRT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
  test_num = 10;
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_SQRT_TEST\n' );
  fprintf ( 1, '  C8_SQRT computes the principal square root of a C8.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '            C1=random            C2=C8_SQRT(C1)         C3=C2*C2\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ c1, seed ] = c8_uniform_01 ( seed );
    c2 = c8_sqrt ( c1 );
    c3 = c1 * c1;

    fprintf ( 1, '  %10f  %10f    %10f  %10f    %10f  %10f\n', ...
      real ( c1 ), imag ( c1 ), ...
      real ( c2 ), imag ( c2 ), ...
      real ( c3 ), imag ( c3 ) );

  end

  return
end


text:
function k = invcmpndKernCompute(kern, x, x2)

% INVCMPNDKERNCOMPUTE Compute the INVERSE-PRECISION-CMPND kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the kernel that is the inverse of
% the sum of the precisions of 1 or more kernels, given
% inputs associated with rows and columns.
% That is, K(x,x2) = inv(inv(K1(x,x2)) + inv(K2(x,x2)) + ... )
% Notice that the individual kernels Ki are each allowed to take a different
% subset of the inputs x, x2 (defined in the field kern.comp{i}.index).
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel parameters for the kernel that is the inverse of
% the sum of the precisions of 1 or more kernels
% given a design matrix of inputs.
% That is, K(x,x2) = inv(inv(K1(x1,x2)) + inv(K2(x1,x2)) + ... )
% Notice that the individual kernels Ki are each allowed to take a different
% subset of the inputs x, x2 (defined in the field kern.comp{i}.index).
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : invcmpndKernParamInit, kernCompute, kernCreate, invcmpndKernDiagCompute
%
% COPYRIGHT : Andreas C. Damianou, 2012

% KERN

% !!! TODO: There are several precomputations and tricks to speed this up...!

if nargin > 2
  i = 1;
  if ~isempty(kern.comp{i}.index)
    % only part of the data is involved in the kernel.
    k = kernCompute(kern.comp{i}, ...
                         x(:, kern.comp{i}.index), ...
                         x2(:, kern.comp{i}.index));
  else
    % all the data is involved with the kernel.
    k = kernCompute(kern.comp{i}, x, x2);
  end
  k = pdinv(k); %%%
  for i = 2:length(kern.comp)
    if ~isempty(kern.comp{i}.index)
      % only part of the data is involved in the kernel.
      k_cur = kernCompute(kern.comp{i}, ...
                           x(:, kern.comp{i}.index), ...
                           x2(:, kern.comp{i}.index));
    else
      % all the data is involved with the kernel.
      k_cur = kernCompute(kern.comp{i}, x, x2);
    end
    k = k + pdinv(k_cur);
  end
else
  i = 1;
  if ~isempty(kern.comp{i}.index)
    % only part of the data is involved with the kernel.
    k  = kernCompute(kern.comp{i}, x(:, kern.comp{i}.index));
  else
    % all the data is involved with the kernel.
    k  = kernCompute(kern.comp{i}, x);
  end
  k = pdinv(k); %%%
  for i = 2:length(kern.comp)
    if ~isempty(kern.comp{i}.index)
      % only part of the data is involved with the kernel.
      k_cur = kernCompute(kern.comp{i}, x(:, kern.comp{i}.index));
    else
      % all the data is involved with the kernel.
      k_cur = kernCompute(kern.comp{i}, x);
    end
    k = k + pdinv(k_cur);
  end
end

k = pdinv(k);


text:
function [Gamma,Gammasum,Xi,LL] = hsinference_batched(T,hmm,residuals,XX)

N = length(T);
K = hmm.K;

setstateoptions;
rangeK = 1:K;
ndim = size(residuals,2);

% Cache shared results for use in obslike
for k = rangeK
    hmm.cache.train{k} = train;
    if k == 1 && strcmp(train.covtype,'uniquediag') || strcmp(train.covtype,'shareddiag')
        ldetWishB = 0;
        PsiWish_alphasum = 0;
        for n = 1:ndim
            if ~regressed(n), continue; end
            ldetWishB = ldetWishB+0.5*log(hmm.Omega.Gam_rate(n));
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.Omega.Gam_shape);
        end
        C = hmm.Omega.Gam_shape ./ hmm.Omega.Gam_rate;
    elseif k == 1 && strcmp(train.covtype,'uniquefull') || strcmp(train.covtype,'sharedfull')
        ldetWishB = 0.5*logdet(hmm.Omega.Gam_rate(regressed,regressed));
        PsiWish_alphasum = 0;
        for n = 1:sum(regressed)
            PsiWish_alphasum = PsiWish_alphasum+psi(hmm.Omega.Gam_shape/2+0.5-n/2);
        end
        PsiWish_alphasum=PsiWish_alphasum*0.5;
        C = hmm.Omega.Gam_shape * hmm.Omega.Gam_irate;
        for kk = rangeK
            [hmm.cache.WishTrace{kk},hmm.cache.codevals] = computeWishTrace(hmm,regressed,XX,C,kk);
        end
    elseif strcmp(train.covtype,'diag')
        ldetWishB=0;
        PsiWish_alphasum = 0;
        for n = 1:ndim
            if ~regressed(n), continue; end
            ldetWishB = ldetWishB+0.5*log(hmm.state(k).Omega.Gam_rate(n));
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.state(k).Omega.Gam_shape);
        end
        C = hmm.state(k).Omega.Gam_shape ./ hmm.state(k).Omega.Gam_rate;
        if episodic, iC = hmm.state(k).Omega.Gam_rate / hmm.state(k).Omega.Gam_shape; end
    elseif strcmp(train.covtype,'full')
        ldetWishB = 0.5*logdet(hmm.state(k).Omega.Gam_rate(regressed,regressed));
        PsiWish_alphasum = 0;
        for n = 1:sum(regressed)
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.state(k).Omega.Gam_shape/2+0.5-n/2);
        end
        C = hmm.state(k).Omega.Gam_shape * hmm.state(k).Omega.Gam_irate;
        if episodic, iC = hmm.state(k).Omega.Gam_rate / hmm.state(k).Omega.Gam_shape; end
        [hmm.cache.WishTrace{k},hmm.cache.codevals] = computeWishTrace(hmm,regressed,XX,C,k);
    end
    % Set up cache
    if ~isfield(train,'distribution') || strcmp(train.distribution,'Gaussian')
        hmm.cache.ldetWishB{k} = ldetWishB;
        hmm.cache.PsiWish_alphasum{k} = PsiWish_alphasum;
        hmm.cache.C{k} = C;
    end
end

Gamma = cell(N,1);
LL = cell(N,1);
Xi = cell(N,1);

for j = 1:N

    t0 = sum(T(1:j-1)) - order*(j-1);
    if order > 0
        R = [zeros(order,size(residuals,2)); residuals(t0+1:t0+T(j)-order,:)];
    else
        R = residuals(t0+1:t0+T(j),:);
    end

    Tj = size(R,1);

    L = obslike([],hmm,residuals,XX,hmm.cache);
    L(L<realmin) = realmin;
    L(L>realmax) = realmax;
    L = embeddata_batched(L,Tj,hmm.train.embeddedlags_batched);
    L = permute(mean(L,2),[1,3,2]);

    [Gamma{j},Xij,LL{j}] = fb_Gamma_inference_sub(L,hmm.P,hmm.Pi,size(L,1),order);

    Xi{j} = reshape(Xij,size(Xij,1),K,K);

end

Gamma = cell2mat(Gamma);
Xi = cell2mat(Xi);
LL = cell2mat(LL);
Gammasum = sum(Gamma);

end

author: OHBA-analysis

repo: HMM-MAR