text:
function test_bug3337

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

%%
data = [];
data.label = {'1', '2'};
for i=1:3
  data.trial{i} = randn(2, 1000);
  data.time{i} = (1:1000)/1000;
  data.sampleinfo(i,:) = [(i-1)*1000+1 i*1000];
end

data.trial{1}(:,30:40)  = nan;
data.trial{1}(:,80:90)  = nan; % 2nd segment within one trial
data.trial{2}(1,100)    = nan; % one sample
data.trial{3}(1,:)      = nan; % whole trial

cfg = [];
[cfg, artifact] = ft_artifact_nan(cfg, data);

assert(size(artifact,1)==4)
assert(size(artifact,2)==2)

correct = [
    30          40
    80          90
  1100        1100
  2001        3000
  ];

assert(isequal(artifact, correct));

%%

cfg.artfctdef.reject = 'nan';
datanan = ft_rejectartifact(cfg, data);
assert(isequaln(data.trial{1}, datanan.trial{1}));

%%

cfg.artfctdef.minaccepttim = 0;
cfg.artfctdef.reject = 'partial';
datapartial = ft_rejectartifact(cfg, data);
assert(numel(datapartial.trial)==5);


text:
function [Population,z,znad] = EnvironmentalSelection(Population,W,N,z,znad)
% The environmental selection of theta-DEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,N);
    St = find(FrontNo<=MaxFNo);

    %% Normalization
    [PopObj,z,znad] = Normalization(Population(St).objs,z,znad);

    %% theta-non-dominated sorting
    tFrontNo = tNDSort(PopObj,W);

    %% Selection
    MaxFNo    = find(cumsum(hist(tFrontNo,1:max(tFrontNo)))>=N,1);
    LastFront = find(tFrontNo==MaxFNo);
    LastFront = LastFront(randperm(length(LastFront)));
    tFrontNo(LastFront(1:sum(tFrontNo<=MaxFNo)-N)) = inf;
    Next      = St(tFrontNo<=MaxFNo);
    % Population for next generation
    Population = Population(Next);
end

function tFrontNo = tNDSort(PopObj,W)
% Do theta-non-dominated sorting

    N  = size(PopObj,1);
    NW = size(W,1);

    %% Calculate the d1 and d2 values for each solution to each weight
    normP  = sqrt(sum(PopObj.^2,2));
    Cosine = 1 - pdist2(PopObj,W,'cosine');
    d1     = repmat(normP,1,size(W,1)).*Cosine;
    d2     = repmat(normP,1,size(W,1)).*sqrt(1-Cosine.^2);

    %% Clustering
    [~,class] = min(d2,[],2);

    %% Sort
    theta = zeros(1,NW) + 5;
    theta(sum(W>1e-4,2)==1) = 1e6;
    tFrontNo = zeros(1,N);
    for i = 1 : NW
        C = find(class==i);
        [~,rank] = sort(d1(C,i)+theta(i)*d2(C,i));
        tFrontNo(C(rank)) = 1 : length(C);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [objout] = convm2java(objin)

if iscell(objin),
    objin = cell2mat(objin);
end

if isnumeric(objin)
   objout = m2j(objin);
elseif isstruct(objin)
    fields = fieldnames(objin);
    vals = struct2cell(objin);
    [x, y, k] = size(objin);
    objout = matlab.jmstruct.JMStruct;
    objout.setFields(fields);
    objout.setDimX(x);
    objout.setDimY(y);
    for I=1:numel(vals)
        tmp = convm2java(vals(I));
        vals(I) = tmp;
    end
    objout.setValues(vals);
else
    objout = cellstr(objin);
end;
return;

function [jarray] = m2j(marray)
    jarray = matlab.jmarray.JMArray;
    [x, y, k] = size(marray);
    jarray.setDimX(x);
    jarray.setDimY(y);
    jarray.setRealpart(real(marray));
    if isreal(marray),
        jarray.setReal(true);
    else
        jarray.setReal(false);
        jarray.setImagpart(imag(marray));
    end
return;

author: Sable

repo: mcbench-benchmarks


text:
%DEMO_SVI_CLASSIFIC  Classification problem demonstration for SVI GP
%
%  Description
%    The demonstration program is based on synthetic two class data
%    used by B.D. Ripley (Pattern Recognition and Neural Networks,
%    1996}. The data consists of 2-dimensional vectors that are
%    divided into two classes, labeled 0 or 1. Each class has a
%    bimodal distribution generated from equal mixtures of Gaussian
%    distributions with identical covariance matrices. A Bayesian
%    approach is used to find the decision line and predict the
%    classes of new data points. The result can be compared to the
%    ones from the DEMO_CLASSIFIC.
%
%    The probability of y being one is assumed to be
%
%      p(y=1|f) = normcdf(f)
%
%    The latent values f are given a zero mean Gaussian process
%    prior. This implies that at the observed input locations
%    latent values have prior
%
%      f ~ N(0, K),
%
%    where K is the covariance matrix, whose elements are given as
%    K_ij = k(x_i, x_j | th). The function k(x_i, x_j | th) is
%    covariance function and th its parameters.
%
%    Here we demonstarte use of stochastic variational inference
%    methods to find the posterior of the latent values and parameters.
%    With these we can make predictions on the class probability of
%    future observations. See Hensman et. al. (2013) for the
%    detailed treatment.
%
%  See also
%    DEMO_SVI_REGRESSION, DEMO_CLASSIFIC
%
%
%  References:
%    Hensman, J., Fusi, N. and Lawrence, N. D. (2013). Gaussian
%    processes for big data. arXiv preprint arXiv:1309.6835.

% Copyright (c) 2014 Tuomas Sivula

% 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.

% This demonstration is based on the dataset used in the book Pattern
% Recognition and Neural Networks by B.D. Ripley (1996), Cambridge
% University Press.

% Training data
S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];
[n, nin] = size(x);

% Test data
xt1=repmat(linspace(min(x(:,1)),max(x(:,1)),20)',1,20);
xt2=repmat(linspace(min(x(:,2)),max(x(:,2)),20)',1,20)';
xt=[xt1(:) xt2(:)];

% Create likelihood function
lik = lik_probit();
%lik = lik_logit();

% Create covariance functions
gpcf = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 10);
% Set the prior for the parameters of covariance functions
pl = prior_t();
pm = prior_sqrtt('s2',0.5);
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl,'magnSigma2_prior', pm);

% Create the GP structure (type is by default FULL)
fprintf('SVI GP classification model with probit likelihood\n')
gp = gp_set('lik', lik, 'cf', gpcf, ...
  'latent_method', 'SVI', 'jitterSigma2', 1e-6);

% Select 20 inducing inputs by clustering 10 from both training class
% inputs
fprintf(['Select 20 inducing inputs by clustering 10 points ', ...
  'from both training classes ...'])
Sw = warning('off','stats:kmeans:EmptyCluster');
[~,X_u1] = kmeans(x(y==1,:), 10,'Start','uniform',...
    'EmptyAction','singleton');
[~,X_u2] = kmeans(x(y==-1,:), 10,'Start','uniform',...
    'EmptyAction','singleton');
warning(Sw);
X_u = [X_u1 ; X_u2];
fprintf(' done\n')

% Optimise
maxi = 1000; % The maximum number of iteration rounds
gp = svigp(gp,x,y,'X_u',X_u,'maxiter',maxi,'mu2',1e-7);
% Make predictions
[Eft, Varft, lpyt, Eyt, Varyt] = ...
    gpsvi_pred(gp, x, y, xt, 'yt', ones(size(xt,1),1) );

% Visualise predictive probability p(ystar = 1) with grayscale
figure, hold on;
n_pred=size(xt,1);
h1=pcolor(reshape(xt(:,1),20,20),reshape(xt(:,2),20,20),reshape(exp(lpyt),20,20));
set(h1, 'edgealpha', 0), set(h1, 'facecolor', 'interp')
colormap(repmat(linspace(1,0,64)', 1, 3).*repmat(ones(1,3), 64,1))
axis([-inf inf -inf inf]), %axis off
plot(x(y==-1,1),x(y==-1,2),'o', 'markersize', 8, 'linewidth', 2);
plot(x(y==1,1),x(y==1,2),'rx', 'markersize', 8, 'linewidth', 2);
plot(gp.X_u(:,1),gp.X_u(:,2),'gs', 'markersize', 8, 'linewidth', 1);
set(gcf, 'color', 'w'), title('predictive probability, training cases and inducing inputs with SVIGP', 'fontsize', 14)

author: gpstuff-dev

repo: gpstuff

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


text:
function conf = rfcn_config_ohem(varargin)
% conf = rfcn_config(varargin)
% --------------------------------------------------------
% D&T implementation
% Modified from MATLAB  R-FCN (https://github.com/daijifeng001/R-FCN/)
% and Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2017, Christoph Feichtenhofer
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------
%

     ip = inputParser;

    %% training
    % whether use gpu
    ip.addParamValue('use_gpu',         gpuDeviceCount > 0, ...
                                                        @islogical);
    % Image scales -- the short edge of input image
    ip.addParamValue('scales',          600,            @ismatrix);
    % Max pixel size of a scaled input image
    ip.addParamValue('max_size',        1000,           @isscalar);
    % Images per batch
    ip.addParamValue('ims_per_batch',   2,              @isscalar);
    % Minibatch size, set as -1 if using all the rois
    ip.addParamValue('batch_size',      -1,            @isscalar);
    % Fraction of minibatch that is foreground labeled (class > 0),
    % which is disabled when batch_size = -1
    ip.addParamValue('fg_fraction',     -1,           @isscalar);
    % Overlap threshold for a ROI to be considered foreground (if >= fg_thresh)
    ip.addParamValue('fg_thresh',       0.5,            @isscalar);
    % Overlap threshold for a ROI to be considered background (class = 0 if
    % overlap in [bg_thresh_lo, bg_thresh_hi))
    ip.addParamValue('bg_thresh_hi',    0.5,            @isscalar);
    ip.addParamValue('bg_thresh_lo',    0.0,            @isscalar);
    % mean image, in RGB order
    ip.addParamValue('image_means',     128,            @ismatrix);
    % Use horizontally-flipped images during training?
    ip.addParamValue('use_flipped',     true,           @islogical);
    % Vaild training sample (IoU > bbox_thresh) for bounding box regresion
    ip.addParamValue('bbox_thresh',     0.5,            @isscalar);
    % Whether to perform class agnostic bbox regression
    ip.addParamValue('bbox_class_agnostic', true,  @islogical);

    % random seed
    ip.addParamValue('rng_seed',        6,              @isscalar);

    %% testing
    ip.addParamValue('test_scales',     600,            @isscalar);
    ip.addParamValue('test_max_size',   1000,           @isscalar);
    ip.addParamValue('test_nms',        0.3,            @isscalar);
    ip.addParamValue('test_binary',     false,          @islogical);

    ip.parse(varargin{:});
    conf = ip.Results;

    % if image_means is a file, load it
    if ischar(conf.image_means)
        s = load(conf.image_means);
        s_fieldnames = fieldnames(s);
        assert(length(s_fieldnames) <= 2);
        conf.image_means = s.(s_fieldnames{1});
        if length(s_fieldnames) == 2
          conf.image_std = s.(s_fieldnames{2});
        end
    end

end

author: feichtenhofer

repo: Detect-Track


text:
function [ x, tan, status ] = p00_limit ( problem, option, nvar, x1, tan1, ...
  x2, tan2, lim )

%*****************************************************************************80
%
%% P00_LIMIT seeks a limit point.
%
%  Discussion:
%
%    For a given index 1 <= LIM <= NVAR, a limit point X is a point which
%    satisfies F(X) = 0 and TAN(X)(LIM) = 0, that is, X is a point on the
%    solution curve, and the LIM-th component of the tangent vector at X
%    is zero.
%
%    This function may be called if a limit point has been bracketed,
%    that is, if X1 and X2 are points on the curve with the property that
%    there is a change in sign in the LIM-th component of the tangent
%    vector between X1 and X2.
%
%    The function carries out an iteration seeking a point X between
%    X1 and X2 for which the LIM-th tangent component is zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X1(NVAR), TAN1(NVAR), a point on the curve,
%    and its tangent vector.
%
%    Input, real X2(NVAR), TAN2(NVAR), a second point on the curve,
%    and its tangent vector.
%
%    Input, integer LIM, the index of the entry of TAN which
%    we are seeking to zero.
%
%    Output, real X(NVAR), TAN(NVAR), the computed limit point
%    and its tangent vector.
%
%    Output, integer STATUS.
%    nonnegative, the limit point was computed in STATUS steps.
%    negative, the limit point could not be computed.
%
  VERBOSE = 0;
%
%  Use a fixed parameter index, but do NOT use LIM.
%
  x = x2 - x1;
  x(lim) = 0.0;
  par_index = r8vec_amax_index ( nvar, x );
%
%  Start the zero finding process.
%
  a = 0.0;
  b = 1.0;
  tol = sqrt ( eps );
  arg = 0.0;
  status_zero = 0;
  value = 0.0;

  status = 0;

  while ( 1 )

    [ arg, status_zero ] = zero_rc ( a, b, tol, value, status_zero );

    if ( status_zero < 0 )
      status = -1;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P00_LIMIT - Fatal error!\n' );
      fprintf ( 1, '  ZERO_RC returned an error flag.\n' );
      break
    end

    if ( arg == 0.0 )

      x = x1;
      tan = tan1;

    elseif ( arg == 1.0 )

      x = x2;
      tan = tan2;

    else

      x = ( 1.0 - arg ) * x1 ...
          +       arg   * x2;

      [ x, status_newton ] = p00_newton ( problem, option, nvar, x, par_index );

      if ( status_newton < 0 )
        status = -2;
        fprintf ( 1, '\n' );
        fprintf ( 1, 'P00_LIMIT - Fatal error!\n' );
        fprintf ( 1, '  ZERO_RC returned an error flag.\n' );
        break
      end

      tan = p00_tan ( problem, option, nvar, x );

    end

    value = tan(lim);

    if ( VERBOSE )
      fprintf ( 1, '  %8d  %14e  %14e\n', status_zero, arg, value );
    end

    status = status + 1;

    if ( status_zero == 0 )
      break
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function [y,dy] = splineTransformation2D(w,x,varargin)
%
% computes y = dy*w and returns dy = kron(I2,Q2,Q1), where
% Q{i}(:,1) = spline(x(:,1));
%
% if no argumanets are given, the parameters for the identity map are returned.
%
% required inputs are: p (number of spline coefficients), omega and m
%
% see also transformations/contents.m, trafo.m
%==============================================================================

function [y,dy] = splineTransformation2D(w,x,varargin)

% the persitent variable stores the matrix
% Q(x) = kron( I_2 , Q2, Q1 );
% p is number of spline coefficients, m is size of grid, omega is domain
persistent Q p m omega

% p = []; m = []; omega = [];
for k=1:2:length(varargin), % overwrite default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

if nargin==0
    help(mfilename)
    runMinimalExample;
    return;
else
  y = mfilename('fullfile');
  dy = zeros(2*prod(p),1);             % parameterization of identity
  if ischar(w),  Q  = []; w = []; end; % reset Q
  if isempty(w), return;          end;
end;

if isempty(w) || (size(Q,1) ~= numel(x)) || (size(Q,2) ~= numel(w)),
  % it is assumed that x is a cell centered grid, extract xi1 and xi2
  dim = size(omega,2)/2
  n   = numel(x)/dim
  if n == prod(m)
    q = m;
  elseif n == prod(m+1)
    q = m+1
  else
    error('can not handle this grid')
  end;

  x  = reshape(x,[q,2]);
  Q1 = getQ1d(omega(1:2),q(1),p(1),x(:,1,1));
  Q2 = getQ1d(omega(3:4),q(2),p(2),x(1,:,2));
  Q  = kron(speye(2),kron(sparse(Q2),sparse(Q1)));
  if nargout == 0, return; end;
end;

y = x(:) + Q*w;
dy = Q;

%------------------------------------------------------------------------------

function Q = getQ1d(omega,m,p,xi)
Q  = zeros(m,p); xi = reshape(xi,[],1);
for j=1:p,
  cj=zeros(p,1); cj(j) = 1;
  Q(:,j) = splineInter(cj,omega,xi);
end;

%------------------------------------------------------------------------------

function runMinimalExample
fprintf('%s: minimal example\n',mfilename)

omega = [0,10,0,8]; m = [8,9]; p = [5,6];
w = zeros([p,2]);  w(3,3,1) = 0.05; w(3,4,2) = -0.1;
x = getCellCenteredGrid(omega,m);
y = feval(mfilename,w(:),x,'omega',omega,'m',m,'p',p,'Q',[]);
FAIRfigure(1); clf;
plotGrid(x,omega,m,'color','r'); axis image; hold on;
plotGrid(y,omega,m,'color','b'); axis image; hold off;

fctn = @(w) feval(mfilename,w(:),x,'omega',omega,'m',m,'p',p,'Q',[]);
w = w + randn(size(w));
checkDerivative(fctn,w(:),'fig',2);

%==============================================================================


text:
%%***************************************************
%% min sum_k bk*yk
%% s.t. sum yk*Hk  <= 0
%%      y1 = 1
%% Hk = -hankel(ek) if   1 <= k <= n
%%    = -hankel(0,e(k-n+1)) if n+1 <= k <=2*n-1
%%
%% [blk,At,C,b] = sdphankel(n);
%%***************************************************

   function [blk,At,C,b] = sdphankel(n);

   randn('seed',0);
   tmp = randn(n,n);
   tmp = tmp+tmp';
   X{1} = tmp + norm(tmp,'fro')*speye(n,n);
%%
   for k = 1:n
      ek = zeros(n,1); ek(k) = -1;
      AA{k} = sparse(hankel(ek));
   end
   zz = zeros(n,1);
   for k = n+1:2*n-1
      ek = zeros(n,1); ek(k-n+1) = -1;
      AA{k} = sparse(hankel(zz,ek));
   end
   blk{1,1} = 's'; blk{1,2} = n;
   At = svec(blk,AA,1);
   C{1} = spconvert([n n 0]);
   b = AXfun(blk,At,[],X);
%%
   blk{2,1} = 'u'; blk{2,2} = 1;
   ee = zeros(1,2*n-1); ee(1) = 1;
   At{2,1} = ee;
   C{2,1} = 1;
%%***************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% The COBRAToolbox: testFEA.m
%
% Purpose:
%     - testFEA tests the Flux Enrichemnt Analysis
%     function
%
% Author:
%     - Marouen BEN GUEBILA 09/02/2017

global CBTDIR

%Test requirements
requiredToolboxes = {'bioinformatics_toolbox','statistics_toolbox'};
prepareTest('toolboxes',requiredToolboxes);

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testFEA'));
cd(fileDir);

% load the model and reference data
load('testDataFEA.mat');
model = getDistributedModel('ecoli_core_model.mat');

% run FEA
resultCellFtest = FEA(model, 1:10, 'subSystems');

% assert equality of test results and reference data
assert(isequal(resultCellFtest, resultCellF));

% check when the groups argument is not a string
try
    resultCellFtest = FEA(1:10, model, 0);
catch ME
    assert(length(ME.message) > 0)
end

% check less than 3 input arguments
try
    resultCellFtest = FEA(model, 1:10);
catch ME
    assert(length(ME.message) > 0)
end

% check when the rxnSet is not a vector
try
    resultCellFtest = FEA(model, [1:10; 1:10], 'subSystems');
catch ME
    assert(length(ME.message) > 0)
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [Aver,p,ip] = hor_to_ver(Ahor)
%hor_to_ver   reorder matrix from horizontal to vertical
%   [Aver,p,ip] = hor_to_ver(Ahor)
%   input
%          Ahor    matrix from horizontal ordering of square grid
%   output
%          Aver    matrix from vertical ordering of square grid
%          p       horizontal to vertical permutation, xv = x(ip)
%          ip      vertical to horizontal permutation, y = yv(p)
%
%   IFISS function: HCE; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage

%To multiply a vector x by Aver:
%   xv = x(ip), yv = Aver*xv, y = yv(p);

N = length(Ahor);
n = sqrt(N);
p  = zeros(N,1);
ip = zeros(N,1);

for j=1:n,
   for i=1:n,
      hor = (j-1)*n+i;
      ver = (i-1)*n+j;
      ip(ver) = hor;
      p(hor) = ver;
   end
end

Aver = Ahor(ip,ip);


text:
function [orient,cano_rot]=cosmo_fmri_orientation(ds)
% get orientation of a dataset
%
% [orient,cano_rot]=cosmo_fmri_orientation(ds)
%
% Inputs:
%    ds             fmri dataset struct;
%
% Output:
%    orient         three letter string indicating the orientation of this
%                   dataset
%    cano_rot       canonical rotation matrix (relative to LPI)
%
% Example:
%     ds=cosmo_synthetic_dataset();
%     [orient,cano_rot]=cosmo_fmri_orientation(ds);
%     disp(orient)
%     %|| LPI
%     disp(cano_rot)
%     %|| 1 0 0 0
%     %|| 0 1 0 0
%     %|| 0 0 1 0
%     %|| 0 0 0 1
%
% Notes:
%   - there are 3!*3^2 valid orientations, these are:
%         'SAR'  'SAL'  'SPR'  'SPL'  'IAR'  'IAL'  'IPR'  'IPL'
%         'SRA'  'SLA'  'SRP'  'SLP'  'IRA'  'ILA'  'IRP'  'ILP'
%         'ASR'  'ASL'  'PSR'  'PSL'  'AIR'  'AIL'  'PIR'  'PIL'
%         'ARS'  'ALS'  'PRS'  'PLS'  'ARI'  'ALI'  'PRI'  'PLI'
%         'RAS'  'LAS'  'RPS'  'LPS'  'RAI'  'LAI'  'RPI'  'LPI'
%         'RSA'  'LSA'  'RSP'  'LSP'  'RIA'  'LIA'  'RIP'  'LIP'
%     For example, 'LPI' (used in Talairach/MNI) means that
%       * the first dimension goes from left to right
%       * the second dimension goes from posterior to anterior
%       * the third dimension goes from inferior to superior
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    mx=get_affine_matrix(ds);

    rot_orig=mx(1:3,1:3);

    % normalize the rows
    rot_row_norms=sqrt(sum(rot_orig.^2,2));
    rot=bsxfun(@rdivide,rot_orig,rot_row_norms);

    cano_rot=get_canonical_matrix(rot);

    orient=get_orientation(cano_rot);

function cano_rot=get_canonical_matrix(rot)
    % find the three major axes of the rotation matrix
    visited=false(3);
    cano_rot=zeros(4);
    cano_rot(4,4)=1;

    for dim=1:3
        max_v=0;
        for row=1:3
            for col=1:3
                v=abs(rot(row,col));
                if ~visited(row,col) && v>max_v
                    max_row=row;
                    max_col=col;
                    max_v=v;
                end
            end
        end
        if max_v==0
            error('Illegal matrix');
        end
        visited(max_row,:)=true;
        visited(:,max_col)=true;

        cano_rot(max_row,max_col)=sign(rot(max_row,max_col));
    end

    assert(all(visited(:)));
    assert(all(sum(cano_rot~=0,1)==1));
    assert(all(sum(cano_rot~=0,2)==1));

function orient=get_orientation(cano_rot)
    labs=['LR';'PA';'IS'];

    orient='   ';
    for k=1:3
        col=find(cano_rot(1:3,k)~=0);
        assert(numel(col)==1);
        v=cano_rot(col,k);
        orient(k)=labs(col,1+(v<0));
    end

function mx=get_affine_matrix(ds)
    if cosmo_isfield(ds,'a.vol.mat')
        mx=ds.a.vol.mat;
    elseif isnumeric(ds) && size(ds,1)>=3 && size(ds,2)>=3
        mx=ds;
    else
        error('cannot find affine transformation matrix');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [count, failures] = estimate_failures(trajectory, sequence)
% estimate_failures Computes number of failures score
%
% Scans the trajectory for instances of tracker failure and returns the number of failures.
%
% Input:
% - trajectory (cell): A trajectory as a cell array of regions.
% - sequence (cell): A valid sequence descriptor.
%
% Output:
% - count (integer): Number of failures
% - failures (integer vector): Indices of frames where the tracker failed
%

failures = find(cellfun(@(x) numel(x) == 1 && x == 2, trajectory, 'UniformOutput', true));

count = numel(failures);

author: votchallenge

repo: toolkit-legacy


text:
function compatibleStatus = isCompatible(solverName, printLevel, specificSolverVersion)
% determine the compatibility status of a solver based on the file compatMatrix.rst
%
% USAGE:
%    compatibleStatus = isCompatible(solverName, printLevel, specificSolverVersion)
%
% INPUT:
%    solverName:            Name of the solver
%    printLevel:            verbose level (default: 0)
%    specificSolverVersion: string with specific solver version (example: '12.7.1' or '6.5.1')
%
% OUTPUT:
%    compatibleStatus:      compatibility status
%
%                             * 0: not compatible with the COBRA Toolbox (tested)
%                             * 1: compatible with the COBRA Toolbox (tested)
%                             * 2: unverified compatibility with the COBRA Toolbox (not tested)
%
% .. Author: - Laurent Heirendt, Ronan Fleming
%

    global CBTDIR

    % define a default printLevel
    if nargin < 2
        printLevel = 1;
    end

    % determine specific version
    if nargin < 3
        specificSolverVersion = '';
    end

    % define the same interface for tomlab and quadMinos
    if strcmp(solverName, 'tomlab_snopt')
        solverName = 'tomlab_cplex';
    end

    if strcmp(solverName, 'quadMinos')
        solverName = 'dqqMinos';
    end

    if strcmp(solverName, 'lindo_old')
        solverName = 'lindo_legacy';
    end

    % default compatibility status (untested)
    compatibleStatus = 2;

    % check if the solver and the matlab version are compatible
    compatMatrixFile = [CBTDIR filesep 'docs' filesep 'source' filesep 'installation' filesep 'compatMatrix.rst'];

    % read in the file with the compatibility matrix
    persistent C;
    persistent compatMatrix;
    persistent testedOS;

    untestedFlag = false;
    if isempty(compatMatrix)
        C = {};
        compatMatrix = {};
        testedOS = {};
        fid = fopen(compatMatrixFile);
        lineCount = 0;
        while 1

            % save the previous line
            if lineCount > 0
                tline_prev = tline;
            end

            % read in a new line
            tline = fgetl(fid);

            if ~ischar(tline), break; end

            % skip a table line
            if ~isempty(strfind(tline, '+-------------------+')) || ~isempty(strfind(tline, '+====================+'))
                tline = '';
            end

            % if the line is not empty, read it in
            if length(tline) > 1
                if printLevel > 1
                    disp(tline);
                end

                % replace the markers
                tline = strrep(tline, '|x|', ':x:');
                tline = strrep(tline, '|white_check_mark|', ':white_check_mark:');
                tline = strrep(tline, '|warning|', ':warning:');

                if strcmp(tline(1), '|')
                    % split the line at the vertical bar
                    Cpart = strsplit(tline, '|');
                    Cpart = Cpart(2:end-1);

                    % convert incompatible flag
                    Cpart = strrep(Cpart, ':x:', '0');

                    % convert compatible flag
                    Cpart = strrep(Cpart, ':white_check_mark:', '1');

                    % convert untested flag
                    Cpart = strrep(Cpart, ':warning:', '2');

                    C{end+1} = strtrim(Cpart);
                else
                    if ~isempty(C)
                        compatMatrix{end+1} = C;
                        C = {};
                    end
                    if strcmp(tline(1:2), '~~')
                        testedOS{end+1} = strtrim(tline_prev);
                    end
                end
            end
            lineCount = lineCount + 1;
        end
        compatMatrix{end+1} = C;
    end
    % select the compatibility matrix based on the OS
    if isunix && ~ismac % linux
        tableNb = 1;
        % help system_dependent
        % system_dependent is an undocumented builtin function.
        resultVERS = system_dependent('getos');
        tmp = strsplit(testedOS{tableNb});
        if contains(lower(resultVERS), lower(tmp{2}))
            cMatrix = compatMatrix{tableNb};
        else
            untestedFlag = true;
            if printLevel > 0
                fprintf([' > isCompatible: The compatibility can only be verified on Linux ', tmp{2}, '.\n']);
            end
        end
    elseif ismac % macOS
        macpos = find(strncmp(testedOS,'macOS',5));
        macOSFound = false;
        [~, resultVERS] = unix('sw_vers');
        for tableNb = macpos %
            tmp = strsplit(testedOS{tableNb});
            if ~isempty(strfind(resultVERS, tmp{2}))
                cMatrix = compatMatrix{tableNb};
                macOSFound = true;
            end
        end
        if ~macOSFound
            untestedFlag = true;
            if printLevel > 0
                fprintf([' > isCompatible: The compatibility can only be verified on the following mac OS versions: ', strjoin(testedOS(macpos),', '), '.\n']);
            end
        end
    else % Windows
        resultVERS = system_dependent('getos');
        for tableNb = length(testedOS)-1:length(testedOS) % loop through the last 2 tables
            tmp = strsplit(testedOS{tableNb});
            if ~isempty(strfind(resultVERS, tmp{2}))
                cMatrix = compatMatrix{tableNb};
                break;
            else
                untestedFlag = true;
                if printLevel > 0 && tableNb == length(testedOS)
                    fprintf([' > isCompatible: The compatibility can only be verified on Windows 7 and Windows 10.\n']);
                end
            end
        end
    end

    if untestedFlag
        cMatrix = {};
    end

    % determine the version of MATLAB and the corresponding column
    versionMatlab = ['R' version('-release')];
    if ~isempty(cMatrix)
        compatMatlabVersions = cMatrix{1}(2:end);
        colIndexVersion = strmatch(versionMatlab, compatMatlabVersions);
    else
        colIndexVersion = [];
    end

    % any MATLAB version that is not explicitly supported yields a compatibility status of 2
    if isempty(colIndexVersion)
        compatibleStatus = 2;
        if printLevel > 0
            fprintf([' > isCompatible: The solver compatibility cannot be verified with MATLAB ', versionMatlab, '.\n']);
        end
    else
        % replace any underscores in the solvername
        solverNameAlias = strrep(upper(solverName), '_', '');

        % set the solver version to be checked
        if exist('specificSolverVersion', 'var') && ~isempty(specificSolverVersion)
            solverVersion = strrep(specificSolverVersion, '.', '');
        else
            solverVersion = getCobraSolverVersion(solverName, 0);
        end

        % check compatibility of solver
        for i = 1:length(cMatrix)
            % save the row of the compatibilitx matrix
            row = cMatrix{i};

            % determine the name of the solver
            solverNameRow = row{1};
            solverNameRow = upper(solverNameRow);
            solverNameRow = regexprep(solverNameRow, {'\ ', '\.'}, '');

            % find the correct row
            if strcmpi([solverNameAlias, solverVersion], solverNameRow)
                % retrieve the compatibility status from the table
                compatibilityBoolean = row{colIndexVersion + 1};

                % define output of solverVersion
                if ~isempty(solverVersion),
                    txtSolverVersion = [' (version ', solverVersion, ')'];
                else
                    txtSolverVersion = '';
                end

                % determine the compatibility status
                if strcmpi(compatibilityBoolean, '1')
                    compatibleStatus = 1;
                    if printLevel > 0
                        fprintf([' > isCompatible: ', lower(solverName), txtSolverVersion, ' is compatible and fully tested with MATLAB ', versionMatlab, ' on your type of operating system.\n']);
                    end
                elseif strcmpi(compatibilityBoolean, '0')
                    compatibleStatus = 0;
                    if printLevel > 0
                        fprintf([' > isCompatible: ', lower(solverName), txtSolverVersion, ' is NOT compatible with MATLAB ', versionMatlab, ' on your  type of operating system.\n']);
                    end
                else
                    compatibleStatus = 2;
                    if printLevel > 0
                        fprintf([' > isCompatible: The compatibility of ', upper(solverName), txtSolverVersion, ' is not fully tested and might not be compatible with MATLAB ', versionMatlab, ' on your  type of operating system.\n']);
                    end
                end
            end
        end
    end

    matver = split(version, '.');
    matver = str2double(strcat(char(matver(1)), '.', char(matver(2))));

    % special case: cplex_direct
    if strcmp(solverName, 'cplex_direct') && matver < 8.4
        compatibleStatus = 0;
        if printLevel > 0
            fprintf([' > ', 'cplex_direct is NOT compatible with MATLAB ', versionMatlab, ' on your  type of operating system. Try using the tomlab_cplex interface.\n']);
        end
    end
end


text:
function L_adapt = ReinhardBilateralFiltering(L, pAlpha, pPhi, pEpsilon)
%
%
%      L_adapt = ReinhardBilateralFiltering(L, pAlpha, pPhi, pEpsilon)
%
%
%       Input:
%           -L: input grayscale image
%           -pAlpha: value of exposure of the image
%           -pPhi: a parameter which controls the sharpening
%           -pEpsilon: smoothing threshold
%
%       Output:
%           -L_adapt: filtered image
%
%     Copyright (C) 2015  Francesco Banterle
%
%     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/>.
%

if(~exist('pAlpha', 'var'))
    pAlpha = ReinhardAlpha(L);
end

if(~exist('pPhi', 'var'))
    pPhi = 8;
end

if(~exist('pEpsilon', 'var'))
    pEpsilon = 0.05; %as in the original paper
end

sMax = 8;
tmp = ((2^pPhi) * pAlpha) / (sMax^2);
L_tmp = L ./ (L + tmp);
L_adapt = bilateralFilter(L_tmp, [], 0, 1.0, 1.6, pEpsilon / 2.0);
L_adapt = RemoveSpecials(L_adapt * tmp ./ (1 - L_adapt));

end

author: banterle


text:
%% Create MATLAB pool job
% configuration??????????II. (6)?Profile????
j = createMatlabPoolJob('configuration', 'hpcfig1');
% Directories and files that worker can access
set(j, 'FileDependencies', {'RunDanmIt2.m'})

%% Create new task in job
createTask(j, @RunDanmIt2);
% Submit a Job to the Job Queue
submit(j);
% Wait for a job to finish running
waitForState(j);
% Retrieve the Job's Results
results = getAllOutputArguments(j);
results{:}
% Remove job
%destroy(j);

repo: Aero-matlab


text:
function [ ival1, ival2,  val1, val2 ] = p08_hx ( option )

%*****************************************************************************80
%
%% P08_HX reports the constraint equation data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, integer VAL1, IVAL2, the indices of the two constrained
%    variables.
%
%    Output, real VAL1, VAL2, the values to which the two constrained
%    variables are to be set.
%
  ival1 = 4;
  val1 = 0.0;
  ival2 = 5;
  val2 = 0.0;

  return
end


text:
function printtp(S)
%PRINTTP Print LTI models of a TP model
%	PRINTTP(S)
%
%	S - core tensor of a TP model

n = size(S);
LTI = shiftdim(S, length(n)-2)


text:
function z = conv(x,y)

%Disciplined convex programming information for CONV:
%   Convolutions are a combination of both multiplications and additions.
%   Therefore, to insure that the no-product rule is satisfied, CONV(X,Y)
%   requires that either X or Y be constant. Further conditions may also
%   apply depending on the exact content of X and Y. For example, if the
%   elements of Y are all convex, then the convolution kernel X must be
%   real, and all of the elements must have the same sign.
%
%Disciplined geometric programming information for TIMES:
%   CONV(X,Y) requires that either X or Y be constant, and that the
%   non-constant term be log-convex or log-affine. Strictly speaking,
%   CONV(X,Y) where X and Y are both log-convex satisfies the DGP rulest,
%   but this version does not support that scenario.

narginchk(2,2);
sx = size(x);
sy = size(y);
if sum(sx~=1)>1 || sum(sy~=1)>1,

    error( 'Arguments must be vectors.' );

elseif any(sx==0) && any(sy==0),

    error( 'At least one argument must be non-empty.' );

elseif cvx_isconstant(x) || cvx_isconstant(y),

    sz = sy;
    sx = prod(sx);
    sy = prod(sy);
    nz = sx + sy - 1;
    sz(sz>1) = nz;
    if cvx_isconstant(x),
        [xi,xj,xv] = find(cvx_constant(x));
        yb         = cvx_basis(y);
    else
        [xi,xj,xv] = find(cvx_constant(y));
        yb         = cvx_basis(x);
    end
    [yi,yj,yv] = find(yb);
    nx = length(xv);
    ny = length(yv);
    ox = ones(1,nx);
    yi = yi(:,ox);
    xi = reshape( xi + xj - 2, 1, nx );
    yj = yj(:,ox) + xi(ones(ny,1),:);
    yv = yv * reshape( xv, 1, nx );
    z  = sparse( yi, yj, yv, size(yb,1), nz );
    z  = cvx( sz, z );
    if nnz( cvx_classify( z ) == 13 ),
        error( 'Disciplined convex programming error:\n   Illegal affine combination of convex/concave terms in convolution.' );
    end

else

    error( 'At least one argument must be constant.' );

end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [IPPEPoses,refinedPoses,HHat,U,Q] = featureBasedPlanePoseFromImage(planarTemplate,inputImage,K,kc,featureOpts,ransacOpts,IPPEOpts,plottingOpts)
%featureBasedPlanePoseFromImage: An example of using IPPE to solve a plane's pose given a frontal 'template' image and a single input image. This follows the following pipeline:
%
%1. Detect features in the template and input image
%2. Match features based on their descriptor similarity
%3. Use the matched features to robustly estimate the homography between the template image and the input image.
%4. From the inlier matches, use IPPE to estimate the camera's pose.
%5. (Optionally) Refine the camera's pose with Levenberg–Marquardt.
%
% The camera should be intrinsically calibrated a priori. We use a perspective camera model with lens distortion, which
% is the standard model and used in OpenCV and Bouguet's Matlab camera calibration toolbox
%(http://www.vision.caltech.edu/bouguetj/calib_doc/). This is parameterised by a 3x3 calibration matrix K and a 5x1 distortion vector
%kc. You can determine K and kc by running Bouguet's calibration toolbox using a checkerboard calibration target.
%
%Inputs
%
%planarTemplate: a structure holding the planar object used to compute the
%camera's pose. See makePlanarTemplate.m for a description of its fields.
%
%inputImage: a 2D unit8 input image (grayscale or rgb). This is an image of the plane viewed
%from the camera.
%
%K: The input image's 3x3 intrinsic matrix (5x1 double).
%
%kc: The input image's distortion parameters (5x1 double)
%
%Note that when using ASIFT, it is best to use an undistorted image (i.e.
%the effects of lens distortion are undone). This is because ASIFT does som
%match filtering based on epipolar geometry, and does not handle lens
%distortion.
%
%featureOpts: Options for feature detection and matching. See
%detectAndMatchFeatures.m for details.
%
%ransacOpts: Options for ransac. See basicHomographyRansac.m for details.
%detectAndMatchFeatures.m for details.
%
%IPPEOpts: Options for IPPE pose estimation. See perspectiveIPPE.m for details.
%
%plottingOpts: Options for plotting the fitted homography. This has two
%fields:
%    plottingOpts.doPlot (true or false): true if we want to visualise
%    results, false otherwise
%    plottingOpts.figId (positive integer): the figure number for plotting the results.
%
%Outputs:
%IPPEPoses: structure containing the two pose solutions from IPPE. See
%perspectiveIPPE.m for details.
%
%refinedPoses: structure containing the two refined pose solutions using
%Levenberg–Marquardt (same format as IPPEPoses). See
%perspectiveIPPE.m for details.
%
%HHat is the homography matrix outputted from ransac (3x3 double). Note
%that this is the mapping from the template's image to the input image's
%normalised pixel coordinates.
%
%poseResults.U is the set of inlier correspondences detected in the template image (in mm) (2xN double)
%
%poseResults.Q is the set of inlier correspondences detected in the input image (in normalised pixels) (2xN double)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of the IPPE package for fast plane-based pose
% estimation from the paper "Infinitesimal Plane-based Pose Estimation" by Toby Collins and Adrien Bartoli,
% published in the International Journal of Computer Vision, September
% 2014. A copy of the author's pre-print version can be found here:
%
% http://isit.u-clermont1.fr/~ab/Publications/Collins_Bartoli_IJCV14.pdf
%
% This package is free and covered by the BSD licence without any warranty. We hope you find this code useful and please cite our paper in your work:
% (c) Toby Collins 2015
%
%
%
%@article{
%year={2014},
%issn={0920-5691},
%journal={International Journal of Computer Vision},
%volume={109},
%number={3},
%doi={10.1007/s11263-014-0725-5},
%title={Infinitesimal Plane-Based Pose Estimation},
%url={http://dx.doi.org/10.1007/s11263-014-0725-5},
%publisher={Springer US},
%keywords={Plane; Pose; SfM; PnP; Homography},
%author={Collins, Toby and Bartoli, Adrien},
%pages={252-286},
%language={English}
%}
%
%
%basic argument checking:
assert(isa(inputImage,'uint8'));
assert(size(K,1)==3);
assert(size(K,2)==3);
assert(size(kc,1) == 5);
assert(size(kc,2) == 1);

%first detect and match features:
[p,q] = detectAndMatchFeatures(planarTemplate.rectifiedImage_g,planarTemplate.roi,...
    inputImage,[],featureOpts);
qNormalised = normaliseImagePoints(q,K,kc);

%now estimate the homography with RANSAC:
[HHat,inliers] = basicHomographyRansac(p,qNormalised,ransacOpts);

%lets plot the results:
if plottingOpts.doPlot
    visualiseHomography(planarTemplate.rectifiedImage_g,planarTemplate.roi,inputImage,HHat,K,kc,p,q,inliers,plottingOpts.figId);
end

%Now use IPPE to estimate the camera's pose.

%First get the points on the template image into mm, and reject outliers:
U = planarTemplate.templatImgPixelSize*p(:,inliers);
Q = qNormalised(:,inliers);

[IPPEPoses,refinedPoses] = perspectiveIPPE(U,Q,[],IPPEOpts);

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function response = Golden(n,p)

persistent mems
if n==1; mems=0; end

if p==0 || p==3 || n==1
    response = 'cooperate';
else
    response = 'defect';
end

if length(mems)>=4
    if isequal( mems(end-3:end), [5 0 5 0] ); response = 'cooperate'; end
end
mems(n)=p;


text:
function [V, dV] = variance(X, dim)
% Computes the variance of X in dimension dim, and its gradient

mu = mean(X, dim);
V = mean(bsxfun(@minus, X, mu).^2, dim);

n = size(X, dim);
dV = bsxfun(@plus, (2/n - 4/n) * mu, 2/n*X);

repo: CategoryShapes


text:
function value = beta_inc ( a, b, x )

%*****************************************************************************80
%
%% BETA_INC returns the value of the incomplete Beta function.
%
%  Discussion:
%
%    This calculation requires an iteration.  In some cases, the iteration
%    may not converge rapidly, or may become inaccurate.
%
%    BETA_INC(A,B,X)
%
%      =   Integral ( 0 <= T <= X ) T^(A-1) (1-T)^(B-1) dT
%        / Integral ( 0 <= T <= 1 ) T^(A-1) (1-T)^(B-1) dT
%
%      =   Integral ( 0 <= T <= X ) T^(A-1) (1-T)^(B-1) dT
%        / BETA(A,B)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    Original FORTRAN77 version by Majumder, Bhattacharjee.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Majumder and Bhattacharjee,
%    Algorithm AS63,
%    Applied Statistics,
%    1973, volume 22, number 3.
%
%  Parameters:
%
%    Input, A, B, the parameters of the function.
%    0.0D+00 < A,
%    0.0D+00 < B.
%
%    Input, real X, the argument of the function.
%    Normally, 0.0D+00 <= X <= 1.0.
%
%    Output, BETA_INC, the value of the function.
%
  it_max = 1000;
  tol = 1.0E-07;

  if ( a <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BETA_INC - Fatal error!\n' );
    fprintf ( 1, '  A <= 0.\n' );
    error ( 'BETA_INC - Fatal error!' );
  end

  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BETA_INC - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    error ( 'BETA_INC - Fatal error!' );
  end

  if ( x <= 0.0 )
    value = 0.0;
    return
  elseif ( 1.0 <= x )
    value = 1.0;
    return
  end
%
%  Change tail if necessary and determine S.
%
  psq = a + b;

  if ( a < ( a + b ) * x )
    xx = 1.0 - x;
    cx = x;
    pp = b;
    qq = a;
    indx = 1;
  else
    xx = x;
    cx = 1.0 - x;
    pp = a;
    qq = b;
    indx = 0;
  end

  term = 1.0;
  i = 1;
  value = 1.0;

  ns = floor ( qq + cx * ( a + b ) );
%
%  Use Soper's reduction formulas.
%
  rx = xx / cx;

  temp = qq - i;
  if ( ns == 0 )
    rx = xx;
  end

  it = 0;

  while ( 1 )

    it = it + 1;

    if ( it_max < it )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'BETA_INC - Fatal error!\n' );
      fprintf ( 1, '  Maximum number of iterations exceeded!\n' );
      fprintf ( 1, '  IT_MAX = %d\n', it_max );
      error ( 'BETA_INC - Fatal error!' );
    end

    term = term * temp * rx / ( pp + i );
    value = value + term;
    temp = abs ( term );

    if ( temp <= tol & temp <= tol * value )
      break
    end

    i = i + 1;
    ns = ns - 1;

    if ( 0 <= ns )
      temp = qq - i;
      if ( ns == 0 )
        rx = xx;
      end
    else
      temp = psq;
      psq = psq + 1.0;
    end

  end
%
%  Finish calculation.
%
  value = value * exp ( pp * log ( xx )  + ( qq - 1.0 ) * log ( cx ) ) ...
    / ( beta ( a, b ) * pp );

  if ( indx )
    value = 1.0 - value;
  end

  return
end


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

%%
clear; close all; clc;

%% Syntax
% |Pb=parbound(P,Ps);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% Test script for optimum mixing
% 17/Feb./2016
% Hideki Kawahara

clear all
close all

nSignal = 6;
baseValue = 100;
nData = 10000;
nIteration = 5000;

stdBest = zeros(nIteration, 1);
stdMean = zeros(nIteration, 1);
stdSD = zeros(nIteration, 1);
for kk = 1:nIteration
    sdList = 0.2 + 3 * rand(nSignal, 1);
    sampleData = randn(nData, nSignal) * diag(sdList) + baseValue;

    H = ones(nSignal, nSignal) * baseValue^2 * 2;
    for ii = 1:nSignal
        H(ii, ii) = H(ii, ii) + 2 * sdList(ii)^2;
    end;
    v = ones(nSignal, 1) * 2 * baseValue^2;
    a = inv(H) * v;
    stdBest(kk) = std(sampleData * a);
    stdMean(kk) = std(sampleData  *ones(nSignal, 1) / nSignal);
    wSD = (1.0 ./ sdList);
    wSD = wSD / sum(wSD);
    stdSD(kk) = std(sampleData * wSD);
end;
figure;plot(stdBest, stdMean, '.');grid on;
hold all
plot([0 100],[0 100]);
axis([0 max([stdBest; stdMean]) 0 max([stdBest; stdMean])]);
axis('square');
xlabel('optimized SD');
ylabel('SD of simple mean');
figure;plot(stdBest, stdSD, '.');grid on;
hold all
plot([0 100],[0 100]);
axis([0 max([stdBest; stdSD]) 0 max([stdBest; stdSD])]);
axis('square');
xlabel('optimized SD');
ylabel('SD of 1/SD mixing');

author: HidekiKawahara


text:
function b = isPointOnLine3d(point, line, varargin)
%ISPOINTONLINE3D Test if a 3D point belongs to a 3D line.
%
%   B = isPointOnLine3d(POINT, LINE)
%   with POINT being [xp yp zp], and LINE being [x0 y0 z0 dx dy dz].
%   Returns 1 if point lies on the line, 0 otherwise.
%
%   If POINT is an N-by-3 array of points, B is a N-by-1 array of booleans.
%
%   If LINE is a N-by-6 array of lines, B is a N-by-1 array of booleans.
%
%   B = isPointOnLine3d(POINT, LINE, TOL)
%   Specifies the tolerance used for testing location on 3D line.
%
%   See also:
%   lines3d, distancePointLine3d, linePosition3d, isPointOnLine
%

% ---------
% author : David Legland
% e-mail: david.legland@inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 31/10/2003.
%

%   HISTORY
%   17/12/2013 create from isPointOnLine

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% test if lines are colinear, using norm of the cross product
b = bsxfun(@rdivide, vectorNorm3d( ...
        crossProduct3d(bsxfun(@minus, line(:,1:3), point), line(:,4:6))), ...
        vectorNorm3d(line(:,4:6))) < tol;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function bow_map = sift(I, siftParam, vocab)
	paths = getPaths();

	% Sift Param
	% siftParam = struct('ds_sampling', 1, 'scales', 1.2, 'descriptor', 'opponentsift');

	%% Compute the SIFT Features
	% Generate three filenames
	ext = {'png', 'txt'};
	r = randsample(100000, length(ext), false);
	pid = getPID();
	for i = 1:length(ext),
		%f{i} = fullfile('/tmp', sprintf('sgupta-sift-%07d-%06d.%s', pid, r(i)), ext{i});
		f{i} = fullfile('/dev/shm', sprintf('sgupta-sift-%07d-%06d.%s', pid, r(i), ext{i}));
	end

	imwrite(im2uint8(I), f{1}, 'png');
	str = sprintf('%s %s --detector densesampling --ds_spacing %d --ds_scales %0.4f --descriptor %s --output %s\n', paths.siftLib, f{1}, siftParam.ds_sampling, siftParam.scales, siftParam.descriptor, f{2});
	system(str);

	info = imfinfo(f{1});
	feats = feat_txt2mat(f{2}, siftParam.scales);

	% Remove the temporary files.
	for i = 1:length(f)
		system(sprintf('rm %s &', f{i}));
	end

	%% Compute the SIFT bow here
	data_img = feats.data';
	feats.pos = double(feats.pos);
	pos = sub2ind([info.Height, info.Width], feats.pos(:,1), feats.pos(:,2));
	cnt_id = 1:numel(pos);
	im_pos = zeros([info.Height info.Width]);
	im_pos(pos) = cnt_id;
	clear feats;

	[drop, binsa] = min(vl_alldist(vocab, single(data_img)), [], 1) ;

	ind = (im_pos > 0);

	bow_map = zeros(size(im_pos));
	bow_map(ind) = binsa(im_pos(ind));
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function [B,R] = FindRotation(X, bit)

UX = randn(size(X,1),bit);
UX = sign(UX);
UX(UX<=0)=0;
UX = normalize(UX);

% optimization
for i=1:1
    % find rotation
    C = UX' * X;
    [UB,~,UA] = svds(double(C),bit);
    %[UB,~,UA] = lansvd(double(C),bit,'L');

    R = UA * UB';

    % find B
    Z = X*R;
    UX = ones(size(Z)).*0;
    for j=1:size(Z,1)
        [b] = findBestBinary(Z(j,:));
        UX(j,:) = b./sqrt(sum(b));
    end
end

B = UX;

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function newDwImg = dtiRawRohdeApplyIntensityCorrection(mc, phaseDir, origDwImg)
% Applies intensity correction given Rohode eddy/motion params in mc
%
%   newDwImg = dtiRawRohdeApplyIntensityCorrection(mc, phaseDir, origDwImg)
%
% mc is the 14 motion/eddy parameters, with translations in 1-3, rotations
% in 4-6. The remaining 8 elements are the eddy-current warping parameters
% (Rohde et al's 'c' parameters used to estimate their b_eddy).
%
% phaseDir specifies the phase-encoding direction, either a 1x3 with a '1'
% indicating the phase-encode dir, or a scalar value of 1, 2 or 3.
%
% Applies the Rohde intensity correction to origDwImg, which should be the
% ORIGINAL diffusion-weighted image volume for which the Rohde parameters
% in mc were estimated. We undo the motion correction to apply this, since
% it is specified in the motion-corrected image space.
%
%   Rohde, Barnett, Basser, Marenco and Pierpaoli (2004). Comprehensive
%   Approach for Correction of Motion and Distortion in Diffusion-Weighted
%   MRI. MRM 51:103-114.
%
% HISTORY:
%
% 2007.05.03 RFD wrote it.

if(numel(phaseDir)>1) phaseDir = find(phaseDir); end
if(length(phaseDir)>1 || phaseDir<1 || phaseDir>3) error('phaseDir specification error.'); end

sz = size(origDwImg);
% For our purposes here, ndgrid and meshgrid are equivalent- the only
% difference is which dim increases fastest, but this doesn't matter for a
% simple list of coords.
[X,Y,Z] = ndgrid([1:sz(1)],[1:sz(2)],[1:sz(3)]);
x = [X(:), Y(:), Z(:)];

% x and y are coordinate lists (names are same as used in Rohde et. al.
% paper). x is a list of coords in native image space, y is a list of the
% corresponsding coordinates in motion-corrected space.
motionMat = affineBuild(mc(1:3),mc(4:6));
y = mrAnatXformCoords(motionMat,x);

% Intensity correction (ic) is computed in motion-corrected space (y):
c = mc(7:14);
if(phaseDir==1)
    ic = 1 + c(1) + c(4)*y(:,2) + c(5)*y(:,3) + 2*c(7)*y(:,1) - 2*c(8)*y(:,1);
elseif(phaseDir==2)
    ic = 1 + c(2) + c(4)*y(:,1) + c(6)*y(:,3) - 2*(c(7)+c(8))*y(:,2);
else
    ic = 1 + c(3) + c(5)*y(:,1) + c(6)*y(:,2) + 4*c(8)*y(:,3);
end

% We want to apply ic (in motion-corrected space) to the original,
% uncorrected space. This is easier for the downstream processing pipeline,
% since we will likely reslice the original to a standard space and thus
% never actually generate a motion-corrected version. Otherwise, we would
% be forced to resample twice. So, we interpolate ic back to the
% uncorrected space x.
ic = myCinterp3(double(ic),[sz(1) sz(2)], sz(3), x, 1.0);
ic = reshape(ic,sz);
% Now apply it:
newDwImg = double(origDwImg).*ic;

return;

% Rohde et. al. only work out the intensity correction for phaseDir = 2.
% Here we solve for the other two options. Writing their eq. 12 out:
% c1*y1 + c2*y2 + c3*y3 + c4*y1*y2 + c5*y1*y3 + c6*y2*y3 + (c7*y1^2 - c7*y2^2) + (2*c8*y3^2 - c8*y1^2 - c8*y2^2)
%
% For phaseDir = 1:
% c1    + 0     + 0     + c4*y2    + c5*y3    + 0        + (2*c7*y1 - 0      ) + (0         - 0       - 2*c8*y1)
% c1 + c4*y2 + c5*y3 + 2*c7*y1 - 2*c8*y1
%
% For phaseDir = 2:
% NOTE: there is a typo in the Rohde paper where they solve for this
% derivative. In eq. 13, they present: c2 + c4*y1 + c6*y3 + 2*(c7+c8)*y2
% 0     + c2*y2 + 0     + c4*y1    + 0        + c6*y3    + (0       - 2*c7*y2) + (0         - 0       - 2*c8*y2)
% c2 + c4*y1 + c6*y3 - 2*c7*y2 - 2*c8*y2
% c2 + c4*y1 + c6*y3 - 2*(c7+c8)*y2
%
% For phaseDir = 3:
% 0     +       + c3*y3 + 0        + c5*y1    + c6*y2    + (0       - 0      ) + (4*c8*y3   - 0       - 0     )
% c3 + c5*y1 + c6*y2 + 4*c8*y3
%
% Or, simply:
% eq = 'c1*y1 + c2*y2 + c3*y3 + c4*y1*y2 + c5*y1*y3 + c6*y2*y3 + (c7*y1^2 - c7*y2^2) + (2*c8*y3^2 - c8*y1^2 - c8*y2^2)';
% maple(['diff(' eq ',y1)'])
% maple(['diff(' eq ',y2)'])
% maple(['diff(' eq ',y3)'])

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% OPENMAS EVENT HANDLER (OMAS_eventHandler.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function is designed to handle an event occuring between objectA and
% objectB. An EVENT object is created at the current time, based on the two
% objects and the enumerated event type.

% Author: James A. Douthwaite 06/10/17

function [EVENT] = OMAS_eventHandler(time_event,metaObjectA,metaObjectB,eventEnumeration)
% INPUTS:
% time_event         - The time the event are occurred
% metaobjectA        - The reference object (entity or meta)
% metaobjectB        - The detected/collision object (entity or meta)
% eventEnumeration   - The type of event indicated by the simulator

% OUTPUTS:
% EVENT              - The generated detection event structure

% INPUT HANDLING
% if ~exist('eventEnumeration','var')
%     error('[ERROR]\tEvent "type" unspecified.\n');
% end

% /////////////// GENERATE THE APPROPRIATE EVENT CLASS ////////////////////
% The simulation will generate one of the following events inreponse to the
% request. The result is an event object and the associated updated meta
% object.
switch eventEnumeration
% DETECTION EVENTS
    case eventType.detection
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Detection notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = detectionEvent(time_event,metaObjectA,metaObjectB,infoString);
% DETECTION-LOSS EVENTS
    case eventType.null_detection
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Detection-loss notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = detectionEvent(time_event,metaObjectA,metaObjectB,infoString,eventType.null_detection);
 % WARNING/NEAR-MISS EVENTS
    case eventType.warning
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Proximity warning notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = warningEvent(time_event,metaObjectA,metaObjectB,infoString);
% A WARNING/NEAR MISS CONDITION NULLIFICATION
    case eventType.null_warning
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Proximity warning-clear notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = warningEvent(time_event,metaObjectA,metaObjectB,infoString,eventType.null_warning);
% COLLISION/GEOMETRIC VIOLATION EVENTS
    case eventType.collision
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Collision notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = collisionEvent(time_event,metaObjectA,metaObjectB,infoString);
% COLLISION/GEOMETRIC VIOLATION NULLIFICATION EVENTS
    case eventType.null_collision
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Collision-clear notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = collisionEvent(time_event,metaObjectA,metaObjectB,infoString,eventType.null_collision);
% WAYPOINT ACHIEVED EVENTS
    case eventType.waypoint
        % GENERATE C-FRIENDLY NOTIFICATION STRING
        infoString = char(['Waypoint notification [',metaObjectA.name,':',metaObjectB.name,']']);
        EVENTobj = waypointEvent(time_event,metaObjectA,metaObjectB,infoString);
% WAYPOINT-LOSS EVENTS
    case eventType.null_waypoint
        % A WAYPOINT RESET EVENT (NON-WAYPOINT)
        infoString = sprintf('Waypoint-reset notification [%s:%s]',metaObjectA.name,metaObjectB.name);
        EVENTobj = waypointEvent(time_event,metaObjectA,metaObjectB,infoString,eventType.null_waypoint);
% EVENT UNKNOWN
    otherwise
        fprintf('[ERROR]\tEvent type not recognised.');
        EVENT = [];
        return
end
% OUTPUT EVENT STRUCTURE
EVENT = [];
% CONVERT THE EVENT OBJECT INTO A REPORT STRUCTURE
if ~isempty(EVENTobj)
    % GENERATE THE SIMULATION EVENT STRUCTURE
    [EVENT] = event2Struct(EVENTobj);
    % IF VERBOSE, GENERATE NOTIFICATION OF EVENT CREATION
%     if verbosity > 2
%         fprintf('[%s]\tevent created (ID:%d @t=%.2fs):\t%s\n',EVENTobj.name,EVENTobj.eventID,EVENTobj.time,EVENTobj.info);
%     end
end
end

% GENERATE EVENT RECORD
function [EVENT] = event2Struct(EVENTobject)
% This function converts the event object into a structure for the
% simulation log.
% INPUT
% EVENTobject - The EVENT class object
% OUTPUT
% EVENT       - The EVEN log structure

% Append EVENT CLASS TO STRUCTURE FOR REPORTING
EVENT = struct('eventID',EVENTobject.eventID,...
                  'time',EVENTobject.time,...
                  'type',EVENTobject.type,...
                  'name',EVENTobject.name,...
            'objectID_A',EVENTobject.objectID_A,...
                'name_A',EVENTobject.name_A,...
               'state_A',EVENTobject.state_A,...
            'objectID_B',EVENTobject.objectID_B,...
                'name_B',EVENTobject.name_B,...
               'state_B',EVENTobject.state_B,...
            'separation',EVENTobject.separation,...
                  'info',EVENTobject.info);
% clear EVENTobject
end

repo: OpenMAS


text:
function [hdr, nameFlag] = read_buffer_v1_header(headerfile)
% function [hdr, nameFlag] = read_buffer_v1_header(headerfile)
%
% This function reads a FCDC buffer header from a binary file.
%
% On return, nameFlag has one of the following values:
%   0 = No labels were generated (fMRI etc.)
%   1 = Fake labels were generated
%   2 = Got channel labels from chunk information

% (C) 2010 S. Klanke

type = {
  'char'
  'uint8'
  'uint16'
  'uint32'
  'uint64'
  'int8'
  'int16'
  'int32'
  'int64'
  'single'
  'double'
};

wordsize = {
  1 % 'char'
  1 % 'uint8'
  2 % 'uint16'
  4 % 'uint32'
  8 % 'uint64'
  1 % 'int8'
  2 % 'int16'
  4 % 'int32'
  8 % 'int64'
  4 % 'single'
  8 % 'double'
};

hdr = [];
hdr.nSamples    = 0;
%hdr.nEvents     = 0;
hdr.nSamplesPre = 0;
hdr.nTrials     = 1;
hdr.orig =[];

nameFlag = 0;

FA = fopen([headerfile '.txt'], 'r');
while ~feof(FA);
  s = fgetl(FA);
  indEq = find(s=='=');
  if numel(indEq)==1
    key = s(1:(indEq-1));
    val = s((indEq+1):end);
    switch key
      case 'fSample'
        hdr.Fs = str2num(val);
      case 'nSamples'
        hdr.nSamples = str2num(val);
      case 'nEvents'
%        hdr.nEvents = str2num(val);
      case 'nChans'
        hdr.nChans = str2num(val);
        hdr.label = cell(hdr.nChans,1);
      case 'dataType'
        if strcmp(val,'float32')
           hdr.orig.data_type = 'single';
        elseif strcmp(val,'float64')
           hdr.orig.data_type = 'double';
        else
          hdr.orig.data_type = val;
        end
      case 'version'
        hdr.orig.version = str2num(val);
      case 'endian'
        if strcmp(val,'big')
          hdr.orig.endianness='ieee-be';
        elseif strcmp(val,'little')
          hdr.orig.endianness='ieee-le';
        else
          error('Invalid value for key ''endian''.');
        end
    end
    continue;
  end
  indCol = find(s==':');
  if numel(indCol)>=1
    ind = str2num(s(1:(indEq-1)));
    name = s((indEq+1):end);
    if ~isempty(ind) && ind>=1 && ind<=hdr.nChans
      hdr.label{ind} = name;
      nameFlag = 2;
    else
      error('Invalid channel name definition - broken header file?');
    end
  end
end

fclose(FA);

F = fopen(headerfile,'rb',hdr.orig.endianness);
bin=[];
bin.nChans    = double(fread(F,1,'uint32'));
bin.nSamples  = double(fread(F,1,'uint32'));
bin.nEvents   = double(fread(F,1,'uint32'));
bin.Fs        = double(fread(F,1,'single'));
bin.data_type = fread(F,1,'uint32');
bin.bufsize   = fread(F,1,'uint32');

if hdr.nChans ~= bin.nChans
   error('Number of channels in binary header does not match ASCII definition');
end
if hdr.nSamples ~= bin.nSamples
   warning('Number of samples in binary header does not match ASCII definition');
end
%if hdr.nEvents ~= bin.nEvents
%   error('Number of samples in binary header does not match ASCII definition');
%end
if strcmp(hdr.orig.data_type, type{bin.data_type+1})
   hdr.orig.wordsize = wordsize{bin.data_type+1};
else
   error('Data type in binary header does not match ASCII definition');
end

if hdr.nSamples == 0
  datafile = [headerfile(1:end-6) 'samples'];
  FD = fopen(datafile, 'rb');
  fseek(FD,0,'eof');
  fileSize = ftell(FD);
  fclose(FD);

  n = fileSize / hdr.orig.wordsize / hdr.nChans;
  hdr.nSamples = floor(n);
  if hdr.nSamples ~= n
    warning('Size of ''samples'' is not a multiple of the size of one sample');
  end
end

% TODO: add chunk handling
while ~feof(F)
	typeAndSize = fread(F, 2, 'uint32');
	if numel(typeAndSize) < 2
		break
	end
	switch typeAndSize(1)
		case 1 % channel names
			% already dealt with, TODO: maybe check consistency with ASCII stuff
			dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
		case 2 % FT_CHUNK_CHANNEL_FLAGS
			% FIXME: ignored for now
			dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
		case 3 % FT_CHUNK_RESOLUTIONS
			if typeAndSize(2) == 8*hdr.nChans
				hdr.resolutions = fread(F, [hdr.nChans, 1], 'double');
			else
				warning('Invalid size of RESOLUTIONS chunk - skipping');
				dummy = fread(F, typeAndSize(2), 'uint8=>uint8');
			end
		case 4 % FT_CHUNK_ASCII_KEYVAL
			dummy = fread(F, typeAndSize(2), 'uint8');
			% FIXME: ignore for now
		case 5 % FT_CHUNK_NIFTI1
			hdr.orig.nifti_1 = fread(F, [1, typeAndSize(2)], 'uint8=>uint8');
			if typeAndSize(2) == 348
				hdr.nifti_1 = decode_nifti1(hdr.orig.nifti_1);
			else
				warning('Invalid size of NIFTI_1 chunk - skipping');
			end
		case 6 % FT_CHUNK_SIEMENS_AP = 6
			hdr.orig.siemensap = fread(F, typeAndSize(2), 'uint8=>uint8');
			if exist('sap2matlab')==3
				hdr.siemensap = sap2matlab(hdr.orig.siemensap);
			end
		case 7 % FT_CHUNK_CTF_RES4 = 7
			hdr.orig.ctf_res4 = fread(F, typeAndSize(2), 'uint8=>uint8');
			tmp_name = tempname;
			FT = fopen(tmp_name, 'wb');
			fwrite(FT, orig.ctf_res4, 'uint8');
			fclose(FT);
			R4F = read_ctf_res4(tmp_name);
			delete(tmp_name);
			% copy over the labels
			hdr.label = R4F.label;
			% copy over the 'original' header
			hdr.orig = R4F;
	end
end

fclose(F);

if nameFlag < 2 && hdr.nChans < 2000
  nameFlag = 1; % fake labels generated - these are unique
  for i=1:hdr.nChans
    hdr.label{i} = sprintf('%d', i);
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function [shape] = normalize_first_shape( Data, options )

shape = Data.shape;

%% calculating bounding box %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[cropmin,cropmax,offset,minshape,marginW,marginH] = bounding_box ( shape );

%% calculate scale factor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W_H = cropmax - cropmin;
wh1 = W_H(1);
wh2 = W_H(2);

CanvasSize = options.canvasSize;

scf = CanvasSize(1)/wh1;
if(scf*wh2 > CanvasSize(2))
    scf = CanvasSize(2)/wh2;
end

%% croping image (for debug only) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
debug = 0;

if debug
    img = imread(Data.img);
    cropImage   = img(cropmin(2):cropmax(2), cropmin(1):cropmax(1));
    scaleImage  = imresize(cropImage, scf);
end

%% scale shape and image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shape = shape - repmat((minshape - [marginW marginH] + offset) ...
    , size(shape, 1), 1);
shape = shape*scf;

if debug
    % Displaying image and feature points.
    figure(1);
    imshow(scaleImage);
    hold on;
    plot(shape(:, 1), shape(:, 2), 'g*');
    pause;
end

end

author: tntrung


text:
function qwv_test02 ( )

%*****************************************************************************80
%
%% QWV_TEST02 tests QWV for a Clenshaw-Curtis rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2014
%
%  Author:
%
%    John Burkardt
%
  a =  -1.0;
  b = +1.0;
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWV_TEST02\n' );
  fprintf ( 1, '  Use the Vandermonde procedure to compute the\n' );
  fprintf ( 1, '  quadrature weights for a Clenshaw-Curtis rule.\n' );
  fprintf ( 1, '  Order N = %d\n', n );
  fprintf ( 1, '  Interval is [%g,%g]\n', a, b );
%
%  Set the points.
%
  x = zeros ( n, 1 );

  for i = 1 : n

    theta =  ( n - i ) * pi / ( n - 1 );

    x(i) = ( ( 1 - cos ( theta ) ) * a   ...
           + ( 1 + cos ( theta ) ) * b ) ...
           /   2.0;

  end

  r8vec_print ( n, x, '  Abscissas:' );
%
%  Determine the corresponding weights.
%
  w = qwv ( n, a, b, x );

  r8vec_print ( n, w, '  Weights:' )

  return
end


text:
function [model,secretionRxnsAdded] = secretionProductGapfill(model,microbeID,database,inputDataFolder)
% This function adds exchange, transport and biosynthesis reactions for
% experimentally shown secreted metabolites according to data collected for
% the DEMETER pipeline.
%
% USAGE:
%
%   [model,secretionRxnsAdded] = secretionProductGapfill(model,microbeID,database,inputDataFolder)
%
% INPUTS
% model:                COBRA model structure
% microbeID:            ID of the reconstructed microbe that serves as the
%                       reconstruction name and to identify it in input tables
% database:             rBioNet reaction database containing min. 3 columns:
%                       Column 1: reaction abbreviation, Column 2: reaction
%                       name, Column 3: reaction formula.
% inputDataFolder:      Folder with input tables with experimental data
%                       and databases that inform the refinement process
%
% OUTPUTS
% model:                COBRA model structure with added pathways if applies
% secretionRxnsAdded:   Reactions added based on experimental data
%
% .. Author:
%       - Almut Heinken, 2019-2020
%       - Bronson R. Weston, 2022, introduced global additional secretion
%       reactions
global additionalSecretionRxns

% structure of lists of reactions to add per secretion product
% non-alphanumeric characters are removed from secretion product names in
% the structure; the script accounts for this when matching with the
% experimental data input
secretionRxns = struct();
% secretionRxns.Folate = {'EX_fol(e)','FOLte','EX_5mthf(e)','MTHFTe','EX_thf(e)','THFte'};
secretionRxns.Folate = {'EX_fol(e)','FOLte','EX_5mthf(e)','MTHFTe','EX_thf(e)','THFte','EX_4abz(e)','4ABZt2','DM_GCALD','DHPS','DHNPA','AKP1','GTPCI'};
secretionRxns.Thiamin = {'EX_thm(e)','THMte'};
% secretionRxns.Riboflavin = {'EX_ribflv(e)','RIBFLVt2r'};
secretionRxns.Riboflavin = {'EX_ribflv(e)','RIBFLVt2r','PMDPHT'};
secretionRxns.Niacin = {'EX_nac(e)','EX_ncam(e)','NACt2r','EX_ncam(e)','NCAMt2r'};
secretionRxns.Pyridoxine = {'EX_pydx(e)','EX_pydxn(e)','EX_pydam(e)','PYDXtr','PYDXNtr','PYDAMtr'};
% secretionRxns.Cobalamin = {'EX_cbl1(e)','CBl1te','EX_adocbl(e)','CBLTDe'};
secretionRxns.Cobalamin = {'EX_cbl1(e)','CBl1te','EX_adocbl(e)','CBLTDe','sink_dmbzid','CPC4MT','CPC5MT','LTHRK','SHCHD2','CYRDAR'};
% secretionRxns.Menaquinone = {'EX_mqn7(e)','MQN7te','EX_mqn8(e)','MK8t'};
secretionRxns.Menaquinone = {'EX_mqn7(e)','MQN7te','EX_mqn8(e)','MK8t','AMMQT8r','AMMQT72','DHNAOT7','DHNAOT4','DHNAOT','NCOAH','DHNAS','NPHS','SUCBZL','SUCBZS','2S6HCC','ICHORS','HETT','PREN','EX_2obut(e)','2OBUTt2r','EX_adn(e)','ADNCNT3tc','DHQTi'};
% secretionRxns.GABA = {'EX_4abut(e)','ABUTt2r'};
secretionRxns.GABA = {'EX_4abut(e)','ABUTt2r','GLUDC'};
secretionRxns.Biotin ={'EX_btn(e)','BTNT5r'};
secretionRxns.Cholate = {'EX_cholate(e)','BIACt1'};
secretionRxns.Chenodeoxycholate = {'EX_C02528(e)','BIACt2'};
secretionRxns.Deoxycholate = {'EX_dchac(e)','DCAT'};
% secretionRxns.Tyramine = {'EX_tym(e)','TYMt2r'};
% secretionRxns.Tryptamine = {'EX_trypta(e)', 'TRYPTAte'};
secretionRxns.Tyramine = {'EX_tym(e)','TYMt2r','TYRCBOX'};
secretionRxns.Tryptamine = {'EX_trypta(e)','TRYPTAte','LTDCL'};
secretionRxns.Dimethylamine = {'EX_dma(e)', 'DMAt2r'};
secretionRxns.TrimethylamineNoxide = {'EX_tmao(e)'};
secretionRxns.Trimethylamine = {'EX_tmao(e)', 'EX_tma(e)'};
secretionRxns.Spermine = {'EX_sprm(e)','SPRMTDe'};
secretionRxns.Spermidine = {'EX_spmd(e)','SPMDtex2'};
secretionRxns.Putrescine = {'EX_ptrc(e)','PTRCtex2'};
% secretionRxns.pCresol = {'EX_pcresol(e)','PCRESOLt2r'};
secretionRxns.pCresol = {'EX_pcresol(e)','PCRESOLt2r','4HPHACDC'};
secretionRxns.Ammonia = {'EX_nh4(e)','NH4tb'};
secretionRxns.Nitrogen = {'EX_n2(e)','N2t','N2OFO','EX_n2o(e)','N2Ot'};
secretionRxns.Methylamine = {'EX_mma(e)','MMAt2e'};
secretionRxns.Methanol = {'EX_meoh(e)','MEOHt2'};
secretionRxns.Lthreonine = {'EX_thr_L(e)','THRt2r'};
secretionRxns.Linoleicacid = {'EX_lnlc(e)','LNLCt'};
secretionRxns.Lithocholate = {'EX_HC02191(e)','LCAT'};
secretionRxns.Lglutamate = {'EX_glu_L(e)','GLUt2r'};
secretionRxns.Lglutamine = {'EX_gln_L(e)','GLNt2r'};
secretionRxns.Lalanine = {'EX_ala_L(e)','ALAt2r'};
secretionRxns.Indole3acetate = {'EX_ind3ac(e)','IND3ACt2r'};
% secretionRxns.Histamine = {'EX_hista(e)','HISTAt2'};
secretionRxns.Histamine = {'EX_hista(e)','HISTAt2','HISDC'};
secretionRxns.Peroxide = {'EX_h2o2(e)','H2O2t'};
secretionRxns.x_5Aminovalerate = {'EX_5aptn(e)','5APTNt2r'};
secretionRxns.x_2Oxobutyrate = {'EX_2obut(e)','2OBUTt2r'};
% secretionRxns.x_12Ethanediol = {'EX_12ethd(e)','12ETHDt'};
secretionRxns.x_12Ethanediol = {'EX_12ethd(e)','12ETHDt','LCAR2'};
% secretionRxns.Hydrogen = {'EX_h2(e)','H2td'};
secretionRxns.Hydrogen = {'EX_h2(e)','H2td','HYD4'};
% secretionRxns.Propionate = {'EX_ppa(e)','PPAtr'};
% secretionRxns.x_2Aminobutyrate = {'EX_C02356(e)','C02356t2r'};
secretionRxns.x_2Aminobutyrate = {'EX_C02356(e)','C02356t2r','RE2034C'};
% secretionRxns.x_13Propanediol = {'EX_13ppd(e)','13PPDt'};
secretionRxns.x_13Propanediol = {'EX_13ppd(e)','13PPDt','13PPDH','GLYCDH'};
% secretionRxns.x_12propanediol = {'EX_12ppd_S(e)','12PPDt'};
secretionRxns.x_12propanediol = {'EX_12ppd_S(e)','12PPDt','LCARS'};
% secretionRxns.Acetone = {'EX_acetone(e)','ACETONEt2'};
secretionRxns.Acetone = {'EX_acetone(e)','ACETONEt2','ADCi'};
% secretionRxns.Butylamine = {'EX_butam(e)','BUTAMt2r','EX_norval_L(e)','NORVALt2r','EX_M03134(e)','M03134t2r'};
secretionRxns.Butylamine = {'EX_butam(e)','BUTAMt2r','EX_norval_L(e)','NORVALt2r','EX_M03134(e)','M03134t2r','AKVALS','NORVALS','NORVALDC'};
% secretionRxns.Cadaverine = {'EX_15dap(e)','15DAPt'};
secretionRxns.Cadaverine = {'EX_15dap(e)','15DAPt','LYSDC'};
secretionRxns.Formaldehyde = {'EX_fald(e)','r1421'};
% secretionRxns.Urea = {'EX_urea(e)','UREAt'};
secretionRxns.Urea = {'EX_urea(e)','UREAt','ARGN'};
% secretionRxns.Propanol = {'EX_ppoh(e)','PPOHt2r','PPALt2r','EX_ppal(e)'};
secretionRxns.Propanol = {'EX_ppoh(e)','PPOHt2r','PPALt2r','EX_ppal(e)','ALCD3ir'};
% secretionRxns.Propanal = {'EX_ppal(e)','PPALt2r','EX_12ppd_S(e)','12PPDt'};
secretionRxns.Propanal = {'EX_ppal(e)','PPALt2r','EX_12ppd_S(e)','12PPDt','12PPDSDH'};
% secretionRxns.Phenylethylamine = {'EX_peamn(e)','PEAMNt2r'};
secretionRxns.Phenylethylamine = {'EX_peamn(e)','PEAMNt2r','PHYCBOXL'};
% secretionRxns.Isopropanol = {'EX_2ppoh(e)','2PPOHt2r','EX_acetone(e)','ACETONEt2'};
secretionRxns.Isopropanol = {'EX_2ppoh(e)','2PPOHt2r','EX_acetone(e)','ACETONEt2','ALCD20y'};
secretionRxns.Lmalate = {'EX_mal_L(e)','MALt2r'};
secretionRxns.Sulfide = {'EX_h2s(e)','H2St'};

%Add any additional secretion rxns if specified
if ~isempty(additionalSecretionRxns)
    addFields=fields(additionalSecretionRxns);
    for f=1:length(addFields)
        secretionRxns.(addFields{f})=additionalSecretionRxns.(addFields{f});
    end
end

% read in the secretion product data
secretionTable = readInputTableForPipeline([inputDataFolder filesep 'secretionProductTable.txt']);
secretionTable(:,find(strncmp(secretionTable(1,:),'Ref',3)))=[];

secretionGapfillAddConditional = {
    % secretion products, condition, add reaction(s)
    'Sulfide', '~any(ismember(model.rxns, {''CYSDS'',''CYSTS_H2S''}))', {'CYSDS'}
    'Trimethylamine', '~any(ismember(model.rxns, {''TMAOR1e'',''TMAOR2e''}))', {'TMAOt2r','TMAt2r','TMAOR1','TMAOR2'}
    };

% secretion product list from input table
% modify names to agree with structure
% can only contain alphabetic, numeric, or underscore characters and cannot
% start with a number or underscore.
% Any non-alphabetic, numeric, or underscore characters removed. Any names
% starting with a number changed to x_number.
products = secretionTable(1, 2:end);
numStart = find(~cellfun(@isempty, regexp(products, '^\d+')));% find names that start with numbers
products(numStart) = strcat('x_', products(numStart));
products = regexprep(products, '\W', '');% remove special characters

% microbe index in file
orgRow = find(strcmp(microbeID, secretionTable(:, 1)));

% find the secretion products for this microbe
if contains(version,'(R202') % for Matlab R2020a and newer
spCols = find(cell2mat(secretionTable(orgRow, 2:end)) == 1);
else
    spCols = find(str2double(secretionTable(orgRow, 2:end)) == 1);
end

% added rxns list
secretionRxnsAdded = {};

if ~isempty(spCols)
    secProds = products(spCols);
    for i = 1:length(secProds)
        % add rxns that are not already in model
        rxns2Add = setdiff(secretionRxns.(secProds{i}), model.rxns)
        if ~isempty(rxns2Add)
            for j = 1:length(rxns2Add)
                RxnForm = database.reactions(find(ismember(database.reactions(:, 1), rxns2Add{j})), 3);
                model = addReaction(model, rxns2Add{j}, 'reactionFormula', RxnForm{1, 1}, 'geneRule','secretionProductGapfill');
            end
            secretionRxnsAdded = union(secretionRxnsAdded, rxns2Add);
        end
        % add conditional reactions
        if any(ismember(secretionGapfillAddConditional(:, 1), secProds{i}))
            conditions = find(ismember(secretionGapfillAddConditional(:, 1), secProds{i}));
            for k = 1:length(conditions)
                if eval(secretionGapfillAddConditional{conditions(k), 2})
                    addRxns = secretionGapfillAddConditional{conditions(k), 3};
                    for j = 1:length(addRxns)
                        if ~any(ismember(model.rxns, addRxns{j}))
                            formula = database.reactions{ismember(database.reactions(:, 1), addRxns{j}), 3};
                            model = addReaction(model, addRxns{j}, 'reactionFormula', formula, 'geneRule', 'secretionProductGapfill');
                            secretionRxnsAdded = union(secretionRxnsAdded, addRxns{j});
                        end
                    end
                end
            end
        end
    end
end

end


text:
function [dat, W, A, r_values] = proc_cspoc(dat, maxmin_flag, varargin)
% PROC_CSPOC - Canonical Source Power Co-modulation Analysis (cSPoC)
%
% Optimizes spatial filters such that the power of the filtered signals are
% maximally positively (or negatively) correlated. The function can be used
% in two ways:
% (1) The data of a single subject are bandpass filtered in multiple
% frequency bands, using the function proc_filterbank
% (2) The data of multiple subjects are bandpass filtered in a single
% frequency band
%
%Synopsis:
% [DAT, SPOC_W, SPOC_A, R_VALUES]= proc_cspoc(DAT, MINMAX_FLAG, <OPT>)
%
%Arguments:
% DAT         - either a data structure of epoched data of one subject,
%               containing multiband filtered channels (obtained via
%               proc_filterbank), or a cell array, each cell containing a
%               data structure of epoched data that has been multiband
%               filtered
% MINMAX_FLAG - either 1 for maximizing or -1 for minimizing correlation
% OPT         - struct or property/value list of optional properties:
%  .nComponentPairs   - either the string 'all' or an integer, determining
%                       the number of components pairs to be returned,
%                       default: 'all'
%  .nRepeats          - number of re-starts per component pair, default: 10
%  .maxIter           - maximum number of optimizer iterations, default: 200
%  .averageOverEpochs - when optimizing the correlations, average the
%                       source envelopes within epochs, default: false
%
%Returns:
% DAT    - updated data structure
% SPOC_W  - SPOC projection matrix (spatial filters, in the columns)
% SPOC_A  - estimated mixing matrix (activation patterns, in the columns)
% LAMBDA - eigenvalue score of SPOC projections

props= {'nComponentPairs'       'all'       'STRING|DOUBLE[1]'
        'nRepeats'              10          'DOUBLE[1]'
        'averageOverEpochs'     0           'BOOLE'
        'maxIter'               200         'DOUBLE[1]'
       };

if nargin==0,
  dat = props; return
end

dat = misc_history(dat);
misc_checkType(dat, 'STRUCT(x clab y)|CELL');

opt = opt_proplistToStruct(varargin{:});
opt = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

%% contruct data cells
% cSPoC is performed on one multi bandpass filtered data set
if not(iscell(dat))
    % check that all clabs contain the string 'flt'
    if any(cellfun(@(x) isempty(x),cellfun(@(x) strfind(x,'flt'),dat.clab,'UniformOutput',0)))
        error('data are not band-pass filtered by proc_filterbank')
    end
    N = length(unique(str2mat(cellfun(@(x) x(end),dat.clab)')));
    X = cell(1,N);
    for ii = 1:N
        dat2 = proc_selectChannels(dat,sprintf('*flt%d',ii));
        X{ii} = dat2.x;
    end
% cSPoC is performed on multiple single bandpass filtered data sets
else
    N = length(dat);
    X = cell(1,N);
    for ii = 1:N
        X{ii} = dat{ii}.x;
    end
end

%% run cSPoC
opt = renameStructField(opt,'nComponentPairs','n_component_sets');
opt = renameStructField(opt,'nRepeats','n_repeats');
opt = renameStructField(opt,'averageOverEpochs','average_over_epochs');

[W,A,r_values] = cspoc(X,maxmin_flag,opt);

%% project the data onto cSPoC filters
datnew = [];
for ii = 1:N
    if not(iscell(dat))
        dat2 = proc_selectChannels(dat,sprintf('*flt%d',ii));
    else
        dat2 = dat{ii};
    end
    dat2 = proc_linearDerivation(dat2,W{ii},'prependix',sprintf('cspoc%d_',ii));
    datnew = proc_appendChannels(datnew,dat2);
end

dat = datnew;

author: bbci


text:
tdir = '/biac2/wandell2/data/templates/';

% load the zeroth-iteration template (built from a simple average
% of all the ac-pc aligned brains
%im = loadAnalyze(fullfile(tdir,'SIRL55-0.img'));

baseDir = '/biac2/wandell2/data/reading_longitude/dti/*0*';
[files,subCodes] = findSubjects(baseDir, '_dt6', {'es041113','tk040817'});
N = length(files);

% Load all the talscales
for(ii=[1:N])
  fprintf('Loading scale factors from %s (%d of %d)...\n',subCodes{ii}, ii, N);
  d = load(files{ii},'anat');
  sf(ii) = d.anat.talScale;
end
% To get the mean extent, we average the Talairach scale factors
% and then divide Talairach's actual dimensions (in mm) by these
% scales.
tal = mrAnatGetTalairachDists;
meanExtent.sac = abs(tal.sac)/mean([sf(:).sac]);
meanExtent.iac = abs(tal.iac)/mean([sf(:).iac]);
meanExtent.lac = abs(tal.lac)/mean([sf(:).lac]);
meanExtent.rac = abs(tal.rac)/mean([sf(:).rac]);
meanExtent.aac = abs(tal.aac)/mean([sf(:).aac]);
meanExtent.acpc = abs(tal.acpc)/mean([sf(:).acpc]);
meanExtent.ppc = abs(tal.ppc)/mean([sf(:).ppc]);
for(ii=[1:N])
  extent(ii).sac = abs(tal.sac)/sf(ii).sac;
  extent(ii).iac = abs(tal.iac)/sf(ii).iac;
  extent(ii).lac = abs(tal.lac)/sf(ii).lac;
  extent(ii).rac = abs(tal.rac)/sf(ii).rac;
  extent(ii).aac = abs(tal.aac)/sf(ii).aac;
  extent(ii).acpc = abs(tal.acpc)/sf(ii).acpc;
  extent(ii).ppc = abs(tal.ppc)/sf(ii).ppc;
  newScale(ii).sac = meanExtent.sac/extent(ii).sac;
  newScale(ii).iac = meanExtent.iac/extent(ii).iac;
  newScale(ii).lac = meanExtent.lac/extent(ii).lac;
  newScale(ii).rac = meanExtent.rac/extent(ii).rac;
  newScale(ii).aac = meanExtent.aac/extent(ii).aac;
  newScale(ii).acpc = meanExtent.acpc/extent(ii).acpc;
  newScale(ii).ppc = meanExtent.ppc/extent(ii).ppc;
  newScale(ii).pcReference = -meanExtent.acpc;
end

% Reslice using the talscale measurements.
d = load(files{1},'anat');
ac = round(mrAnatXformCoords(inv(d.anat.xformToAcPc),[0 0 0]));
imSkull = d.anat.img;
im = imSkull; im(~d.anat.brainMask) = 0;
sz = size(im);
%imSkull = mrAnatHistogramClip(double(imSkull), 0.5, 0.99);
proportionNonBrain = sum(~d.anat.brainMask(:))./prod(size(d.anat.brainMask));
im = mrAnatHistogramClip(double(im), proportionNonBrain, 0.99);
bb = [-d.anat.mmPerVox.*(ac-1); d.anat.mmPerVox.*(sz-ac)];
ts = newScale(1); ts.talScaleDir = 'tal2acpc';
% warp T1 to tal space
ts.outMat = inv(d.anat.xformToAcPc);
[im,newXform] = mrAnatResliceSpm(im,ts,bb,[1 1 1],[7 7 7 0 0 0],0);
im(isnan(im)) = 0;
im(im<0) = 0; im(im>1) = 1;
ax = zeros(sz(2), sz(1), N+1);
cr = zeros(sz(3), sz(1), N+1);
sg = zeros(sz(3), sz(2), N+1);
meanIm = im;
ax(:,:,1) = flipud(permute(squeeze(im(:,:,ac(3))),[2,1]));
cr(:,:,1) = flipud(permute(squeeze(im(:,ac(2),:)),[2,1]));
sg(:,:,1) = flipud(permute(squeeze(im(ac(1),:,:)),[2,1]));
for(ii=[2:N])
  fprintf('Processing %s (%d of %d)...\n',subCodes{ii}, ii, N);
  d = load(files{ii},'anat');
  % All coords should be the same.
  ac = round(mrAnatXformCoords(inv(d.anat.xformToAcPc),[0 0 0]));
  im = d.anat.img; im(~d.anat.brainMask) = 0;
  im = mrAnatHistogramClip(double(im), 0.4, 0.99);
  ts = newScale(ii); ts.talScaleDir = 'tal2acpc';
  ts.outMat = inv(d.anat.xformToAcPc);
  [im,newXform] = mrAnatResliceSpm(im,ts,bb,[1 1 1],[7 7 7 0 0 0],0);
  im(isnan(im)) = 0;
  im(im<0) = 0; im(im>1) = 1;
  meanIm = meanIm + im;
  ax(:,:,ii) = flipud(permute(squeeze(im(:,:,ac(3))),[2,1]));
  cr(:,:,ii) = flipud(permute(squeeze(im(:,ac(2),:)),[2,1]));
  sg(:,:,ii) = flipud(permute(squeeze(im(ac(1),:,:)),[2,1]));
end
meanIm = meanIm./N;
figure; image(makeMontage(uint8(meanIm*255)));
colormap(gray(256)); axis equal tight off;
ax(:,:,N+1) = mean(ax(:,:,[1:N]), 3);
cr(:,:,N+1) = mean(cr(:,:,[1:N]), 3);
sg(:,:,N+1) = mean(sg(:,:,[1:N]), 3);
ax = uint8(ax.*255+0.5);
cr = uint8(cr.*255+0.5);
sg = uint8(sg.*255+0.5);
figure;image(makeMontage(ax));colormap(gray(256));axis equal tight off;
figure;image(makeMontage(cr));colormap(gray(256));axis equal tight off;
figure;image(makeMontage(sg));colormap(gray(256));axis equal tight off;

tdir = '/snarp/u1/data/templates/';
imwrite(makeMontage(ax), gray(256), fullfile(tdir, 'SIRL55scaledToMean_ax.png'));
imwrite(makeMontage(cr), gray(256), fullfile(tdir, 'SIRL55scaledToMean_cr.png'));
imwrite(makeMontage(sg), gray(256), fullfile(tdir, 'SIRL55scaledToMean_sg.png'));

% Now build an atlas from the mean image
V.dat = int16(meanIm.*(2^15-1)+0.5);
notes = ['SIRL55ms: average of 55 mean-Tal-scaled brains. Created at ' datestr(now,31)];
mmPerVox = d.anat.mmPerVox;
origin = ac;
hdr = saveAnalyze(V.dat, fullfile(tdir,'SIRL55ms'), mmPerVox, notes, origin);
save(fullfile(tdir,'SIRL55ms_details.mat'), 'meanIm','notes','extent','meanExtent','origin','mmPerVox');


text:
function varargout = sound_field_mono_localwfs_sbl(X,Y,Z,xs,src,f,conf)
%SOUND_FIELD_MONO_LOCALWFS_SBL sound field of local WFS using spatial bandwith
%limitation
%
%   Usage: [P,x,y,z,x0] = sound_field_mono_localwfs_sbl(X,Y,Z,xs,src,f,conf)
%
%   Input options:
%       X           - x-axis / m; single value or [xmin,xmax] or nD-array
%       Y           - y-axis / m; single value or [ymin,ymax] or nD-array
%       Z           - z-axis / m; single value or [zmin,zmax] or nD-array
%       xs          - position of point source or direction of plane
%                     wave / m [1x3]
%       src         - source type of the virtual source
%                         'pw' - plane wave
%                         'ps' - point source
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output options:
%       P           - simulated sound field
%       x           - corresponding x values / m
%       y           - corresponding y values / m
%       z           - corresponding z values / m
%       x0          - secondary sources / m
%
%   SOUND_FIELD_MONO_LOCALWFS_SBL(X,Y,Z,xs,src,f,conf) simulates a monochromatic
%   sound field of the given source type (src) synthesized with local wave field
%   synthesis using spatial bandwidth limitation for the frequency f.
%
%   To plot the result use:
%   plot_sound_field(P,X,Y,Z,x0,conf);
%   or simple call the function without output argument:
%   sound_field_mono_localwfs_sbl(X,Y,Z,xs,src,f,conf)
%   For plotting you may also consider to display the result in dB, by setting
%   the following configuration option before:
%   conf.plot.usedB = true;
%
%   See also: plot_sound_field, driving_function_mono_localwfs_sbl

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 7;
nargmax = 7;
narginchk(nargmin,nargmax);
isargnumeric(X,Y,Z);
isargxs(xs);
isargchar(src);
isargscalar(f);
isargstruct(conf);

%% ===== Configuration ==================================================
if strcmp('2D',conf.dimension)
    greens_function = 'ls';
else
    greens_function = 'ps';
end

%% ===== Computation =====================================================
% Get secondary sources
x0 = secondary_source_positions(conf);
% Get driving signals
D = driving_function_mono_localwfs_sbl(x0,xs,src,f,conf);
% Calculate sound field
[varargout{1:min(nargout,4)}] = ...
    sound_field_mono(X,Y,Z,x0,greens_function,D,f,conf);
% Return secondary sources if desired
if nargout==5, varargout{5}=x0; end

author: sfstoolbox

repo: sfs-matlab


text:
function M = spherefactory(n, m)
% Returns a manifold struct to optimize over unit-norm vectors or matrices.
%
% function M = spherefactory(n)
% function M = spherefactory(n, m)
%
% Manifold of n-by-m real matrices of unit Frobenius norm.
% By default, m = 1, which corresponds to the unit sphere in R^n. The
% metric is such that the sphere is a Riemannian submanifold of the space
% of nxm matrices with the usual trace inner product, i.e., the usual
% metric.
%
% See also: obliquefactory spherecomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    if ~exist('m', 'var')
        m = 1;
    end

    if m == 1
        M.name = @() sprintf('Sphere S^%d', n-1);
    else
        M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
    end

    M.dim = @() n*m-1;

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) real(acos(x(:).'*y(:)));

    M.typicaldist = @() pi;

    M.proj = @(x, d) d - x*(x(:).'*d(:));

    M.tangent = M.proj;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
	M.egrad2rgrad = M.proj;

	M.ehess2rhess = @ehess2rhess;
	function rhess = ehess2rhess(x, egrad, ehess, u)
        rhess = M.proj(x, ehess) - (x(:)'*egrad(:))*u;
	end

    M.exp = @exponential;

    M.retr = @retraction;

    M.log = @logarithm;
    function v = logarithm(x1, x2)
        v = M.proj(x1, x2 - x1);
        di = M.dist(x1, x2);
        % If the two points are "far apart", correct the norm.
        if di > 1e-6
            nv = norm(v, 'fro');
            v = v * (di / nv);
        end
    end

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() random(n, m);

    M.randvec = @(x) randomvec(n, m, x);

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n, m);

    M.transp = @(x1, x2, d) M.proj(x2, d);

    M.pairmean = @pairmean;
    function y = pairmean(x1, x2)
        y = x1+x2;
        y = y / norm(y, 'fro');
    end

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
    M.vecmatareisometries = @() true;

end

% Exponential on the sphere
function y = exponential(x, d, t)

    if nargin == 2
        t = 1;
    end

    td = t*d;

    nrm_td = norm(td, 'fro');

    if nrm_td > 4.5e-8
        y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
    else
        % If the step is too small to accurately evaluate sin(x)/x,
        % then sin(x)/x is almost indistinguishable from 1.
        y = x + td;
        y = y / norm(y, 'fro');
    end

end

% Retraction on the sphere
function y = retraction(x, d, t)

    if nargin == 2
        t = 1;
    end

    y = x + t*d;
    y = y / norm(y, 'fro');

end

% Uniform random sampling on the sphere.
function x = random(n, m)

    x = randn(n, m);
    x = x/norm(x, 'fro');

end

% Random normalized tangent vector at x.
function d = randomvec(n, m, x)

    d = randn(n, m);
    d = d - x*(x(:).'*d(:));
    d = d / norm(d, 'fro');

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
% Max Jaderberg 2011

function json_response = web_feval(headers, config)
%WEB_FEVAL Returns a json object of the result of calling the funciton
%   This allows you to run any local matlab file. To be used with webserver.m.
%   HTTP POST to /web_feval.m with the following parameters:
%       func_path: a string which is the path to the .m file you want to
%       call
%       arguments: a json string of the arguments structure to pass to the
%       function you are calling
%
%   Returns a json object containing the result

    response.success = 'false';
    field_names = fieldnames(headers.Content);

%     URL decode the POST data
    for i=1:numel(field_names)
        headers.Content.(field_names{i}) = urldecode(headers.Content.(field_names{i}));
    end

    response.content = '';

    func_path_check = false;
    arguments_check = false;
    if size(field_names)
        if isfield(headers.Content, 'func_path')
            func_path_check = true;
        end
        if isfield(headers.Content, 'arguments')
            arguments_check = true;
        end
    end

    if ~func_path_check
        response.message = 'No function given as func_path POST parameter';
        json_response = mat2json(response);
        return
    end

    func_path = headers.Content.func_path;

    if arguments_check
        arguments = json2mat(headers.Content.arguments);
    else
        arguments = '';
    end

    response.result = run_dot_m(func_path, arguments);
    response.success = 'true';
    response.message = 'Successfully completed request';

    response.content.func_path = func_path;
    response.content.arguments = arguments;

    json_response = mat2json(response);

    return

end

author: jaderberg

repo: python-matlab-bridge