text:
% performs argument checks/transformations similar to those found in
% gp.m from GPML, but Hessian friendly

% Copyright (c) 2014--2015 Roman Garnett.

function gpstruct = gpstruct_check(gpstruct, x)

% Default gp values
defaultstruct = struct( ...
    'inf',{@exact_inference}, ...
    'cov', [], ...
    'lik',{@likGauss}, ...
    'mean', {@zero_mean}, ...
    'hyp', struct('cov',[],'lik',[],'mean',[]), ...
    'prior', [], ...
    'bounds', [] ...
    );

% Assign default values
for field = fields(defaultstruct)'
    if ~isfield(gpstruct,field{:}) || isempty(gpstruct.(field{:}))
        gpstruct.(field{:}) = defaultstruct.(field{:});
    end
end

% No default covariance function
if isempty(gpstruct.cov)
    error('gpstruct:missing_argument', ...
          'covariance function must be defined');
end

% Allow string/function handle input; convert to cell arrays if necessary
for field = {'inf', 'cov', 'lik', 'mean'}
    if ischar(gpstruct.(field{:})) || isa(gpstruct.(field{:}), 'function_handle')
        gpstruct.(field{:}) = {gpstruct.(field{:})};
    end
end

% Check struct of hyperparameters
if ~isfield(gpstruct,'hyp')
    error('gpstruct:missing_argument', ...
        'hyperparameter structure must be defined');
end

% If HYP is a cell array, convert to struct array
if iscell(gpstruct.hyp) && isstruct(gpstruct.hyp{1})
    hyptemp = struct('cov','lik','mean');
    for i = 1:length(gpstruct.hyp)
        hyptemp(i) = gpstruct.hyp{i};
    end
    gpstruct.hyp = hyptemp;
end

% Ensure all hyperparameter fields exist
for field = {'cov', 'lik', 'mean'}
    for i = 1:numel(gpstruct.hyp)
        if (~isfield(gpstruct.hyp(i), field{:}))
            gpstruct.hyp(i).(field{:}) = [];
        end
    end
end

% Check dimension of hyperparameter fields
D = size(x, 2);
shortfields = {'cov','lik','mean'};
longfields = {'covariance','likelihood','mean'};
for ifield = 1:numel(shortfields)
    field = shortfields{ifield};
    expression = feval(gpstruct.(field){:});
    for i = 1:numel(gpstruct.hyp)
        nparams = numel(gpstruct.hyp(i).(field));
        if (nparams ~= eval(expression))
            error('gpstruct:incorrect_specification', ...
                  'wrong number of %s hyperparameters (%i given, %s expected)', ...
                  longfields{ifield}, ...
                  nparams, ...
                  expression);
        end
    end
end

% Check prior and bounds
for field = {'cov', 'lik', 'mean'}
    % Ensure all prior fields exist
    if (~isfield(gpstruct.prior, field{:}))
        gpstruct.prior.(field{:}) = [];
    end
    % Ensure all bounds fields exist
    if (~isfield(gpstruct.bounds, field{:}))
        gpstruct.bounds.(field{:}) = [];
    end
    % Assign [-Inf;Inf] bounds by default
    num_entries = numel(gpstruct.hyp(1).(field{:}));
    for i = 1:num_entries
        if numel(gpstruct.prior.(field{:})) < i
            gpstruct.prior.(field{:}){i} = [];
        end
        if numel(gpstruct.bounds.(field{:})) < i || ...
                isempty(gpstruct.bounds.(field{:}){i})
            gpstruct.bounds.(field{:}){i} = [-Inf,Inf];
        end
    end
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [H,S]=create_CHS(A,B,C,N,M);
% CREATE_CHS Internal function to define matrices for MPC problem

% How many states
n=length(A);p=size(C,1);
[dummy,m]=size(B);

H=zeros(N*p,length(A));
S=zeros(N*p,N*m);

Acum=A;
for j=1:N,
  H(1+p*(j-1):p*j,:)=C*Acum;
  Acum=Acum*A;
end;

for j=1:N,
  Acum=eye(n);
  for k=j:-1:1,
    S(1+p*(j-1):p*j,1+m*(k-1):k*m)=C*Acum*B;
    Acum=Acum*A;
  end;
end;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [nav,stat]=decode_navb(nav,fid,opt,headinfo)
global glc gls
stat=1; NMAX=10000;
nav.eph=repmat(gls.eph,NMAX,1); nav.geph=repmat(gls.geph,NMAX,1);

switch headinfo.type
    case 'N',sys=headinfo.sys;
    case 'G',sys=glc.SYS_GLO;
    case 'L',sys=glc.SYS_GAL;
    case 'J',sys=glc.SYS_QZS;
    otherwise,stat=0;return;
end

while ~feof(fid)

    [eph,geph,type,fid,stat0]=decode_navb_data(opt,headinfo.ver,sys,fid);

    if stat0==1
        switch type
            case 1
                if nav.n+1>size(nav.eph,1)
                    nav.eph(nav.n+1:nav.n+NMAX,1)=repmat(gls.eph,NMAX,1);
                end
                nav.eph(nav.n+1)=eph;
                nav.n=nav.n+1;
            case 2
                if nav.ng+1>size(nav.geph,1)
                    nav.geph(nav.ng+1:nav.ng+NMAX,1)=repmat(gls.geph,NMAX,1);
                end
                nav.geph(nav.ng+1)=geph;
                nav.ng=nav.ng+1;
        end
    end
end
fclose(fid);

if nav.n<size(nav.eph,1)
    nav.eph(nav.n+1:end,:)=[];
end
if nav.ng<size(nav.geph,1)
    nav.geph(nav.ng+1:end,:)=[];
end

return

repo: GINav


text:
function determ = fourier_sine_determinant ( n )

%*****************************************************************************80
%
%% FOURIER_SINE_DETERMINANT returns the determinant of the FOURIER_SINE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  if ( mod ( n, 2 ) == 1 )
    determ = + 1.0;
  else
    determ = - 1.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pde = StokesZulehnerdata
%% STOKESDATA2 data for Stokes equations
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f', @f, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D);

    function z = f(p)
        x = p(:,1); y = p(:,2);
        z(:,1) = zeros(size(x));
        z(:,2) = 4*cos(x).*cos(y);
    end
    % exact velocity
    function z = exactu(p)
        x = p(:,1); y = p(:,2);
        z(:,1) = sin(x).*sin(y);
        z(:,2) = cos(x).*cos(y);
    end
    % exact pressure
    function z = exactp(p)
        x = p(:,1); y = p(:,2);
        z = 2*cos(x).*sin(y) - 2*sin(1)*(1-cos(1));
    end
    % Dirichlet boundary condition of velocity
    function z = g_D(p)
        z = exactu(p);
    end
end

repo: ifem


text:
function cdf = log_series_cdf ( x, a )

%*****************************************************************************80
%
%% LOG_SERIES_CDF evaluates the Logarithmic Series CDF.
%
%  Discussion:
%
%    Simple summation is used, with a recursion to generate successive
%    values of the PDF.
%
%    Thanks to Oscar van Vlijmen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 December 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the PDF.
%    0 < X
%
%    Input, real A, the parameter of the PDF.
%    0.0 < A < 1.0.
%
%    Output, real CDF, the value of the CDF.
%
  cdf = 0.0;

  for x2 = 1 : x

    if ( x2 == 1 )
      pdf = - a / log ( 1.0 - a );
    else
      pdf = ( x2 - 1 ) * a * pdf / x2;
    end

    cdf = cdf + pdf;

  end

  return
end


text:
function Stats = stat_getSecondOrderStats(PConn,alpha)
% Return the mean, standard deviation, and confidence intervals of a
% surrogate distribution
%
% Inputs:
%
%       PConn:     connectivity distribution returned from stat_bootstrap()
%       alpha:     significance level for empircal confidence intervals
%                  (e.g. alpha=0.95 for 95% confidence intervals)
% Outputs:
%
%       Stats:     Structure containing
%                  .mean
%                  .stdev
%                  .ci
%                  for each connectivity measure
%
% See Also: stat_surrogate(), stat_surrogateStats()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 6.8
%   Available at: http://www.sccn.ucsd.edu/wiki/SIFT
%
% Author: Tim Mullen, 2011, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

connmethods = hlp_getConnMethodNames(PConn);

for cnd=1:length(PConn)

    for m=1:length(connmethods)

        sz = size(PConn(cnd).(connmethods{m}));

        Stats(cnd).(connmethods{m}).mean  = mean(PConn(cnd).(connmethods{m}),length(sz));
        Stats(cnd).(connmethods{m}).stdev = std(PConn(cnd).(connmethods{m}),0,length(sz));
        Stats(cnd).(connmethods{m}).ci(1,:,:,:,:,:) = ...
            single(prctile(PConn(cnd).(connmethods{m}),(100*alpha)/2,length(sz)));     % lower ci
        Stats(cnd).(connmethods{m}).ci(2,:,:,:,:,:) = ...
            single(prctile(PConn(cnd).(connmethods{m}),100-100*alpha/2,length(sz)));   % upper ci

    end

end

repo: REST


text:
function [y,g] = objfun(x)
y = 2*x(1)^2+2*x(2)^2-2*x(1)*x(2)-4*x(1)-6*x(2);
g = [4*x(1)-2*x(2)-4;4*x(2)-2*x(1)-6];

repo: Optimization-Theory-and-Algorithm


text:
function [minr, maxr, minc, maxc, slcV, modelMask3M, planC] = ...
    getCropLimits(planC,mask3M,scanNum,cropS)
% getCropLimits.m
% Get limits of bounding box for various cropping options.
%
% AI 5/2/19
%--------------------------------------------------------------------------
%INPUTS:
% planC
% scan3M       : Scan array
% mask3M       : Mask
% cropS        : Dictionary of parameters for cropping
%                Supported methods: 'crop_fixed_amt','crop_to_bounding_box',
%                'crop_to_str', 'crop_around_center', 'none'
%--------------------------------------------------------------------------
% AI 5/2/19
% AI 7/23/19

%Get mask for model config
[modelMask3M, planC] = getMaskForModelConfig(planC,mask3M,scanNum,cropS);

%Compute bounding box
methodC = lower({cropS.method});
crop2DMethodsC = {'crop_to_bounding_box_2d','crop_to_str_2d',...
    'crop_pt_outline_2d'};%Suported 2D crop options
if length(methodC) == 1 && strcmp(methodC{1},'none')
    scanSizeV = size(getScanArray(scanNum,planC));
    minr = 1;
    maxr = scanSizeV(1);
    minc = 1;
    maxc = scanSizeV(2);
    slcV = 1:scanSizeV(3);

elseif any(ismember(methodC,crop2DMethodsC))

    slcV = find(sum(sum(modelMask3M))>0);
    numSlcs = length(slcV);
    minr = nan(numSlcs,1);
    minc = nan(numSlcs,1);
    maxr = nan(numSlcs,1);
    maxc = nan(numSlcs,1);
    for n = 1:numSlcs
        [minr(n), maxr(n), minc(n), maxc(n)] = compute_boundingbox(modelMask3M(:,:,slcV(n)));
    end
else
    [minr, maxr, minc, maxc] = compute_boundingbox(modelMask3M);
    slcV = find(sum(sum(modelMask3M))>0);
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function timeplot(axs_hand,mode)
%
% Utility Function: TIMEPLOT
%
% The purpose of this function is to handle the plotting of all the time
% response options

% Author: Craig Borghesani
% Date: 8/8/94
% Revised: 11/15/94
% Copyright (c) 1999, Prentice-Hall

% if no time response axis present, get out
if ~length(axs_hand), return; end

% obtain handle information
f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};
open_loop = ui_han(78);
clos_loop = ui_han(79);
unit_imp  = ui_han(80);
unit_step = ui_han(81);
unit_ramp = ui_han(82);
rise_time = ui_han(83);
dlay_time = ui_han(84);
peak_time = ui_han(85);
over_shot = ui_han(86);
setl_2per = ui_han(87);
setl_5per = ui_han(88);
cur_sys = get(ui_han(30),'userdata');
stat_bar = get(ui_han(43),'userdata');
range1 = 78:82;
range2 = 83:88;

% distribute handles
axs_data = get(axs_hand,'userdata');

if ~length(mode) | any(mode == [13,14,15]),

% determine various states of time response environment

 if ~length(mode), mode = 0; end
 ct = 2;
 for k = range1,
  if strcmp(get(ui_han(k),'checked'),'on'),
   mode(ct) = k-(range1(1)-1);
  else
   mode(ct) = 0;
  end
  ct = ct + 1;
 end

end

if length(mode) > 1,
 set(axs_data(2:length(axs_data)),'vis','off');
end

open_imp = axs_data(2);
clos_imp = axs_data(5);
open_step = axs_data(3);
clos_step = axs_data(6);
open_ramp = axs_data(4);
clos_ramp = axs_data(7);
tim_txt = axs_data(8:length(axs_data));

plant_mat = get(ui_han(3+cur_sys),'userdata');
for_mat  = get(ui_han(6+cur_sys),'userdata');
bac_mat  = get(ui_han(9+cur_sys),'userdata');

term_mat  = termjoin(plant_mat,for_mat,bac_mat);
[num,den] = termextr(term_mat);

if length(num)==1 & length(den) == 1,
 set(stat_bar,'string','No Time Response computed');
 return;
end

ln = length(num); ld = length(den);
if ln > ld,
 set(stat_bar,'string','You have more Zeros than Poles.  Time response not plotted.');
 return
end

if polyval(den,0) ~= 0,
 final = polyval(num,0)/polyval(den,0);
else
 final = inf;
end

[nump,denp] = termextr(plant_mat);
[numg,deng] = termextr(for_mat);
[numh,denh] = termextr(bac_mat);
numt = conv(conv(nump,numg),numh);
dent = conv(conv(denp,deng),denh);
lnt = length(numt); ldt = length(dent);
numt = [zeros(1,ldt-lnt),numt];
dent = [zeros(1,lnt-ldt),dent];
numcl = conv(conv(nump,numg),denh);
dencl = dent + numt;

finalcl = polyval(numcl,0)/polyval(dencl,0);

tvec = get(ui_han(40),'userdata');
imp_input = [1/(tvec(2)-tvec(1));zeros(length(tvec)-1,1)];
stp_input = ones(length(tvec),1);
rmp_input = tvec';

if term_mat(1,2)~=0, % pure delay time
 loc_t = find(tvec <= term_mat(1,2));
 len = length(loc_t);
 if len > 1,
  imp_input = [zeros(len,1);imp_input(1:(length(imp_input)-len))];
  stp_input = [zeros(len,1);stp_input(1:(length(stp_input)-len))];
  rmp_input = [zeros(len,1);rmp_input(1:(length(rmp_input)-len))];
 end
end

if any(mode==1), % Open-Loop

 if strcmp(get(open_loop,'checked'),'off') | length(mode)>1,
  set(open_loop,'checked','on');

  if strcmp(get(unit_imp,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop impulse response');
   imp_ol = lsim(num,den,imp_input,tvec);
   set(open_imp,'xdata',tvec,'ydata',imp_ol,'vis','on');
   if length(mode) > 1, mode(4) = 0; end

  end

  if strcmp(get(unit_step,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop step response');
   step_ol = lsim(num,den,stp_input,tvec);
   set(open_step,'xdata',tvec,'ydata',step_ol,'vis','on');
   if length(mode) > 1, mode(5) = 0; end

  end

  if strcmp(get(unit_ramp,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop ramp response');
   ramp_ol = lsim(num,den,rmp_input,tvec);
   set(open_ramp,'xdata',tvec,'ydata',ramp_ol,'vis','on');
   if length(mode) > 1, mode(6) = 0; end

  end

 else

  set(stat_bar,'string','Removing Open-Loop');
  set(open_loop,'checked','off');
  set([open_imp,open_step,open_ramp,tim_txt([1:2:9])],'vis','off');

 end

end

if any(mode==2), % Closed-Loop

 if strcmp(get(clos_loop,'checked'),'off') | length(mode)>1,
  set(clos_loop,'checked','on');

  if strcmp(get(unit_imp,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop impulse response');
   imp_cl = lsim(numcl,dencl,imp_input,tvec);
   set(clos_imp,'xdata',tvec,'ydata',imp_cl,'vis','on');
   if length(mode) > 1, mode(4) = 0; end

  end

  if strcmp(get(unit_step,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop step response');
   step_cl = lsim(numcl,dencl,stp_input,tvec);
   set(clos_step,'xdata',tvec,'ydata',step_cl,'vis','on');
   if length(mode) > 1, mode(5) = 0; end

  end

  if strcmp(get(unit_ramp,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop ramp response');
   ramp_cl = lsim(numcl,dencl,rmp_input,tvec);
   set(clos_ramp,'xdata',tvec,'ydata',ramp_cl,'vis','on');
   if length(mode) > 1, mode(6) = 0; end

  end

 else

  set(stat_bar,'string','Removing Closed-Loop');
  set(clos_loop,'checked','off');
  set([clos_imp,clos_step,clos_ramp,tim_txt([2:2:10])],'vis','off');

 end

end

if any(mode==3), % unit impulse

 if strcmp(get(unit_imp,'checked'),'off') | length(mode)>1,
  set(unit_imp,'checked','on');

  if strcmp(get(open_loop,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop impulse response');
   imp_ol = lsim(num,den,imp_input,tvec);
   set(open_imp,'xdata',tvec,'ydata',imp_ol,'vis','on');

  end

  if strcmp(get(clos_loop,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop impulse response');
   imp_cl = lsim(numcl,dencl,imp_input,tvec);
   set(clos_imp,'xdata',tvec,'ydata',imp_cl,'vis','on');

  end

 else

  set(stat_bar,'string','Removing impulse response');
  set(unit_imp,'checked','off');
  set([open_imp,clos_imp],'vis','off');

 end

end

if any(mode==4), % unit step

 if strcmp(get(unit_step,'checked'),'off') | length(mode)>1,
  set(unit_step,'checked','on');

  if strcmp(get(open_loop,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop step response');
   step_ol = lsim(num,den,stp_input,tvec);
   set(open_step,'xdata',tvec,'ydata',step_ol,'vis','on');

  end

  if strcmp(get(clos_loop,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop step response');
   step_cl = lsim(numcl,dencl,stp_input,tvec);
   set(clos_step,'xdata',tvec,'ydata',step_cl,'vis','on');

  end

  set(ui_han(range2),'enable','on');

 else

  set(stat_bar,'string','Removing step response');
  set(unit_step,'checked','off');
  set([open_step,clos_step,tim_txt],'vis','off');
  set(ui_han(range2),'enable','off');

 end

end

if any(mode==5), % unit ramp

 if strcmp(get(unit_ramp,'checked'),'off') | length(mode)>1,
  set(unit_ramp,'checked','on');

  if strcmp(get(open_loop,'checked'),'on'),

   set(stat_bar,'string','Computing open-loop ramp response');
   ramp_ol = lsim(num,den,rmp_input,tvec);
   set(open_ramp,'xdata',tvec,'ydata',ramp_ol,'vis','on');

  end

  if strcmp(get(clos_loop,'checked'),'on'),

   set(stat_bar,'string','Computing closed-loop ramp response');
   ramp_cl = lsim(numcl,dencl,rmp_input,tvec);
   set(clos_ramp,'xdata',tvec,'ydata',ramp_cl,'vis','on');

  end

 else

  set(stat_bar,'string','Removing ramp response');
  set(unit_ramp,'checked','off');
  set([open_ramp,clos_ramp],'vis','off');

 end

end

if any(mode < 6),

 pageview(1,axs_hand);
 ct = 1;
 mode = 0;
 for k = range2,
  if strcmp(get(ui_han(k),'checked'),'on'),
   mode(ct) = k-(range1(1)-1);
  else
   mode(ct) = 0;
  end
  ct = ct + 1;
 end

end

if any(mode==6) & strcmp(get(rise_time,'enable'),'on'), % rise time

 if strcmp(get(rise_time,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   set(stat_bar,'string','Computing open-loop rise time');
   set(rise_time,'checked','on');
   amp = get(open_step,'ydata');
   diff_amp = diff(amp);
   loc_neg = find(diff_amp < 0);
   if ~length(loc_neg), % overdamped system
    less10 = find(amp>=0.1*final);
    less90 = find(amp>=0.9*final);
    riset = tvec(less90(1)) - tvec(less10(1));
    risea = amp(less90(1));
    set(tim_txt(1),'pos',[riset,risea,0],...
            'string',['Rise time (10-90%) = ',num2str(riset)],'vis','on');
   else  % underdamped system
    less0 = 1;
    less100 = find(amp>=final);
    riset = tvec(less100(1)) - tvec(less0);
    risea = amp(less100(1));
    set(tim_txt(1),'pos',[riset,risea,0],...
            'string',['Rise time (0-100%) = ',num2str(riset)],'vis','on');
   end
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   set(stat_bar,'string','Computing closed-loop rise time');
   set(rise_time,'checked','on');
   amp = get(clos_step,'ydata');
   diff_amp = diff(amp);
   loc_neg = find(diff_amp < 0);
   if ~length(loc_neg), % overdamped system
    less10 = find(amp>=0.1*finalcl);
    less90 = find(amp>=0.9*finalcl);
    riset = tvec(less90(1)) - tvec(less10(1));
    risea = amp(less90(1));
    set(tim_txt(2),'pos',[riset,risea,0],...
          'string',['Rise time (10-90%) = ',num2str(riset)],'vis','on');
   else % underdamped system
    less0 = 1;
    less100 = find(amp>=finalcl);
    riset = tvec(less100(1)) - tvec(less0);
    risea = amp(less100(1));
    set(tim_txt(1),'pos',[riset,risea,0],...
            'string',['Rise time (0-100%) = ',num2str(riset)],'vis','on');
   end
  end

 else

  set(stat_bar,'string','Removing rise time');
  set(rise_time,'checked','off');
  set(tim_txt(1:2),'vis','off');
 end

end

if any(mode==7) & strcmp(get(dlay_time,'enable'),'on'), % delay time

 if strcmp(get(dlay_time,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   set(stat_bar,'string','Computing open-loop delay time');
   set(dlay_time,'checked','on');
   amp = get(open_step,'ydata');
   less50 = find(amp>=0.5*final);
   dlayt = tvec(less50(1));
   dlaya = amp(less50(1));
   set(tim_txt(3),'pos',[dlayt,dlaya,0],...
         'string',['Delay time = ',num2str(dlayt)],'vis','on');
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   set(stat_bar,'string','Computing closed-loop delay time');
   set(dlay_time,'checked','on');
   amp = get(clos_step,'ydata');
   less50 = find(amp>=0.5*finalcl);
   dlayt = tvec(less50(1));
   dlaya = amp(less50(1));
   set(tim_txt(4),'pos',[dlayt,dlaya,0],...
         'string',['Delay time = ',num2str(dlayt)],'vis','on');
  end

 else

  set(stat_bar,'string','Removing delay time');
  set(dlay_time,'checked','off');
  set(tim_txt(3:4),'vis','off');
 end

end

if any(mode==8) & strcmp(get(peak_time,'enable'),'on'), % peak time

 if strcmp(get(peak_time,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   amp = get(open_step,'ydata');
   diff_amp = diff(amp);
   loc_neg = find(diff_amp < 0);
   if length(loc_neg),
    set(stat_bar,'string','Computing open-loop peak time');
    set(peak_time,'checked','on');
    peakt = tvec(loc_neg(1)-1);
    peaka = amp(loc_neg(1)-1);
    set(tim_txt(5),'pos',[peakt,peaka,0],...
         'string',['Peak time = ',num2str(peakt)],'vis','on');
   else
    set(stat_bar,'string','Unable to determine open-loop peak time');
   end
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   amp = get(clos_step,'ydata');
   diff_amp = diff(amp);
   loc_neg = find(diff_amp < 0);
   if length(loc_neg),
    set(stat_bar,'string','Computing closed-loop peak time');
    set(peak_time,'checked','on');
    peakt = tvec(loc_neg(1)-1);
    peaka = amp(loc_neg(1)-1);
    set(tim_txt(6),'pos',[peakt,peaka,0],...
         'string',['Peak time = ',num2str(peakt)],'vis','on');
   else
    set(stat_bar,'string','Unable to determine closed-loop peak time');
   end
  end

 else

  set(stat_bar,'string','Removing peak time');
  set(peak_time,'checked','off');
  set(tim_txt(5:6),'vis','off');
 end

end

if any(mode==9) & strcmp(get(over_shot,'enable'),'on'), % percent overshoot

 if strcmp(get(over_shot,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   set(stat_bar,'string','Computing open-loop % overshoot');
   set(over_shot,'checked','on');
   amp = get(open_step,'ydata');
   [max_amp,k] = max(amp);
   perover = 100*(max_amp-final)/final;
   set(tim_txt(7),'pos',[tvec(k),max_amp,0],...
       'string',[sprintf('%0.2f',perover),'% overshoot'],'vis','on');
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   set(stat_bar,'string','Computing closed-loop % Overshoot');
   set(over_shot,'checked','on');
   amp = get(clos_step,'ydata');
   [max_amp,k] = max(amp);
   perover = 100*(max_amp-finalcl)/finalcl;
   set(tim_txt(8),'pos',[tvec(k),max_amp,0],...
       'string',[sprintf('%0.2f',perover),'% overshoot'],'vis','on');
  end

 else

  set(stat_bar,'string','Removing % overshoot');
  set(over_shot,'checked','off');
  set(tim_txt(7:8),'vis','off');
 end

end

if any(mode==10) & strcmp(get(setl_2per,'enable'),'on'), % 2% settling time

 if strcmp(get(setl_2per,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   set(stat_bar,'string','Computing open-loop 2% settling time');
   set(setl_2per,'checked','on');
   amp = get(open_step,'ydata');
   per2 = length(amp);
   while (amp(per2)>0.98*final & amp(per2)<1.02*final), per2=per2-1; end
   per2t = tvec(per2);
   per2a = amp(per2);
   set(tim_txt(9),'pos',[per2t,per2a,0],...
         'string',['2%=',sprintf('%0.2f',per2t)],'vis','on');
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   set(stat_bar,'string','Computing closed-loop 2% settling time');
   set(setl_2per,'checked','on');
   amp = get(clos_step,'ydata');
   per2 = length(amp);
   while (amp(per2)>0.98*finalcl & amp(per2)<1.02*finalcl), per2=per2-1; end
   per2t = tvec(per2);
   per2a = amp(per2);
   set(tim_txt(10),'pos',[per2t,per2a,0],...
         'string',['2%=',sprintf('%0.2f',per2t)],'vis','on');
  end

 else

  set(stat_bar,'string','Removing 2% settling time');
  set(setl_2per,'checked','off');
  set(tim_txt(9:10),'vis','off');
 end

end

if any(mode==11) & strcmp(get(setl_5per,'enable'),'on'), % 5% settling time

 if strcmp(get(setl_5per,'checked'),'off') | length(mode)>1,

  if strcmp(get(open_loop,'checked'),'on'),
   set(stat_bar,'string','Computing open-loop 5% settling time');
   set(setl_5per,'checked','on');
   amp = get(open_step,'ydata');
   per5 = length(amp);
   while (amp(per5)>0.95*final & amp(per5)<1.05*final), per5=per5-1; end
   per5t = tvec(per5);
   per5a = amp(per5);
   set(tim_txt(11),'pos',[per5t,per5a,0],...
         'string',['5%=',sprintf('%0.2f',per5t)],'vis','on');
  end

  if strcmp(get(clos_loop,'checked'),'on'),
   set(stat_bar,'string','Computing closed-loop 5% settling time');
   set(setl_5per,'checked','on');
   amp = get(clos_step,'ydata');
   per5 = length(amp);
   while (amp(per5)>0.95*finalcl & amp(per5)<1.05*finalcl), per5=per5-1; end
   per5t = tvec(per5);
   per5a = amp(per5);
   set(tim_txt(12),'pos',[per5t,per5a,0],...
         'string',['5%=',sprintf('%0.2f',per5t)],'vis','on');
  end

 else

  set(stat_bar,'string','Removing 5% settling time');
  set(setl_5per,'checked','off');
  set(tim_txt(11:12),'vis','off');
 end

end

author: Sable

repo: mcbench-benchmarks


text:
function M = fixedrankembeddedfactory(m, n, k)
% Manifold struct to optimize fixed-rank matrices w/ an embedded geometry.
%
% function M = fixedrankembeddedfactory(m, n, k)
%
% Manifold of m-by-n real matrices of fixed rank k. This follows the
% embedded geometry described in Bart Vandereycken's 2013 paper:
% "Low-rank matrix completion by Riemannian optimization".
%
% Paper link: http://arxiv.org/pdf/1209.3834.pdf
%
% A point X on the manifold is represented as a structure with three
% fields: U, S and V. The matrices U (mxk) and V (nxk) are orthonormal,
% while the matrix S (kxk) is any /diagonal/, full rank matrix.
% Following the SVD formalism, X = U*S*V'. Note that the diagonal entries
% of S are not constrained to be nonnegative.
%
% Tangent vectors are represented as a structure with three fields: Up, M
% and Vp. The matrices Up (mxk) and Vp (mxk) obey Up'*U = 0 and Vp'*V = 0.
% The matrix M (kxk) is arbitrary. Such a structure corresponds to the
% following tangent vector in the ambient space of mxn matrices:
%   Z = U*M*V' + Up*V' + U*Vp'
% where (U, S, V) is the current point and (Up, M, Vp) is the tangent
% vector at that point.
%
% Vectors in the ambient space are best represented as mxn matrices. If
% these are low-rank, they may also be represented as structures with
% U, S, V fields, such that Z = U*S*V'. Their are no resitrictions on what
% U, S and V are, as long as their product as indicated yields a real, mxn
% matrix.
%
% The chosen geometry yields a Riemannian submanifold of the embedding
% space R^(mxn) equipped with the usual trace (Frobenius) inner product.
%
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{vandereycken2013lowrank,
%       Title   = {Low-rank matrix completion by {Riemannian} optimization},
%       Author  = {Vandereycken, B.},
%       Journal = {SIAM Journal on Optimization},
%       Year    = {2013},
%       Number  = {2},
%       Pages   = {1214--1236},
%       Volume  = {23},
%       Doi     = {10.1137/110845768}
%     }
%
% See also: fixedrankfactory_2factors fixedrankfactory_3factors

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%	Feb. 20, 2014 (NB):
%       Added function tangent to work with checkgradient.
%
%   June 24, 2014 (NB):
%       A couple modifications following
%       Bart Vandereycken's feedback:
%       - The checksum (hash) was replaced for a faster alternative: it's a
%         bit less "safe" in that collisions could arise with higher
%         probability, but they're still very unlikely.
%       - The vector transport was changed.
%       The typical distance was also modified, hopefully giving the
%       trustregions method a better initial guess for the trust region
%       radius, but that should be tested for different cost functions too.
%
%    July 11, 2014 (NB):
%       Added ehess2rhess and tangent2ambient, supplied by Bart.
%
%    July 14, 2014 (NB):
%       Added vec, mat and vecmatareisometries so that hessianspectrum now
%       works with this geometry. Implemented the tangent function.
%       Made it clearer in the code and in the documentation in what format
%       ambient vectors may be supplied, and generalized some functions so
%       that they should now work with both accepted formats.
%       It is now clearly stated that for a point X represented as a
%       triplet (U, S, V), the matrix S needs to be diagonal.

    M.name = @() sprintf('Manifold of %dx%d matrices of rank %d', m, n, k);

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

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

    M.norm = @(x, d) sqrt(M.inner(x, d, d));

    M.dist = @(x, y) error('fixedrankembeddedfactory.dist not implemented yet.');

    M.typicaldist = @() M.dim();

    % Given Z in tangent vector format, projects the components Up and Vp
    % such that they satisfy the tangent space constraints up to numerical
    % errors. If Z was indeed a tangent vector at X, this should barely
    % affect Z (it would not at all if we had infinite numerical accuracy).
    M.tangent = @tangent;
    function Z = tangent(X, Z)
        Z.Up = Z.Up - X.U*(X.U'*Z.Up);
        Z.Vp = Z.Vp - X.V*(X.V'*Z.Vp);
    end

    % For a given ambient vector Z, applies it to a matrix W. If Z is given
    % as a matrix, this is straightfoward. If Z is given as a structure
    % with fields U, S, V such that Z = U*S*V', the product is executed
    % efficiently.
    function ZW = apply_ambient(Z, W)
        if ~isstruct(Z)
            ZW = Z*W;
        else
            ZW = Z.U*(Z.S*(Z.V'*W));
        end
    end

    % Same as apply_ambient, but applies Z' to W.
    function ZtW = apply_ambient_transpose(Z, W)
        if ~isstruct(Z)
            ZtW = Z'*W;
        else
            ZtW = Z.V*(Z.S'*(Z.U'*W));
        end
    end

    % Orthogonal projection of an ambient vector Z represented as an mxn
    % matrix or as a structure with fields U, S, V to the tangent space at
    % X, in a tangent vector structure format.
    M.proj = @projection;
    function Zproj = projection(X, Z)

        ZV = apply_ambient(Z, X.V);
        UtZV = X.U'*ZV;
        ZtU = apply_ambient_transpose(Z, X.U);

        Zproj.M = UtZV;
        Zproj.Up = ZV  - X.U*UtZV;
        Zproj.Vp = ZtU - X.V*UtZV';

    end

    M.egrad2rgrad = @projection;

    % Code supplied by Bart.
    % Given the Euclidean gradient at X and the Euclidean Hessian at X
    % along H, where egrad and ehess are vectors in the ambient space and H
    % is a tangent vector at X, returns the Riemannian Hessian at X along
    % H, which is a tangent vector.
    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)

        % Euclidean part
        rhess = projection(X, ehess);

        % Curvature part
        T = apply_ambient(egrad, H.Vp)/X.S;
        rhess.Up = rhess.Up + (T - X.U*(X.U'*T));
        T = apply_ambient_transpose(egrad, H.Up)/X.S;
        rhess.Vp = rhess.Vp + (T - X.V*(X.V'*T));

    end

    % Transforms a tangent vector Z represented as a structure (Up, M, Vp)
    % into a structure with fields (U, S, V) that represents that same
    % tangent vector in the ambient space of mxn matrices, as U*S*V'.
    % This matrix is equal to X.U*Z.M*X.V' + Z.Up*X.V' + X.U*Z.Vp'. The
    % latter is an mxn matrix, which could be too large to build
    % explicitly, and this is why we return a low-rank representation
    % instead. Note that there are no guarantees on U, S and V other than
    % that USV' is the desired matrix. In particular, U and V are not (in
    % general) orthonormal and S is not (in general) diagonal.
    % (In this implementation, S is identity, but this might change.)
    M.tangent2ambient = @tangent2ambient;
    function Zambient = tangent2ambient(X, Z)
        Zambient.U = [X.U*Z.M + Z.Up, X.U];
        Zambient.S = eye(2*k);
        Zambient.V = [X.V, Z.Vp];
    end

    % This retraction is second order, following general results from
    % Absil, Malick, "Projection-like retractions on matrix manifolds",
    % SIAM J. Optim., 22 (2012), pp. 135-158.
    M.retr = @retraction;
    function Y = retraction(X, Z, t)
        if nargin < 3
            t = 1.0;
        end

        % See personal notes June 28, 2012 (NB)
        [Qu, Ru] = qr(Z.Up, 0);
        [Qv, Rv] = qr(Z.Vp, 0);

        % Calling svds or svd should yield the same result, but BV
        % advocated svd is more robust, and it doesn't change the
        % asymptotic complexity to call svd then trim rather than call
        % svds. Also, apparently Matlab calls ARPACK in a suboptimal way
        % for svds in this scenario.
        % [Ut St Vt] = svds([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)], k);
        [Ut, St, Vt] = svd([X.S+t*Z.M , t*Rv' ; t*Ru , zeros(k)]);

        Y.U = [X.U Qu]*Ut(:, 1:k);
        Y.V = [X.V Qv]*Vt(:, 1:k);
        Y.S = St(1:k, 1:k) + eps*eye(k);

        % equivalent but very slow code
        % [U S V] = svds(X.U*X.S*X.V' + t*(X.U*Z.M*X.V' + Z.Up*X.V' + X.U*Z.Vp'), k);
        % Y.U = U; Y.V = V; Y.S = S;

    end

    M.exp = @exponential;
    function Y = exponential(X, Z, t)
        if nargin < 3
            t = 1.0;
        end
        Y = retraction(X, Z, t);
        warning('manopt:fixedrankembeddedfactory:exp', ...
               ['Exponential for fixed rank ' ...
                'manifold not implemented yet. Used retraction instead.']);
    end

    % Less safe but much faster checksum, June 24, 2014.
    % Older version right below.
    M.hash = @(X) ['z' hashmd5([sum(X.U(:)) ; sum(X.S(:)); sum(X.V(:)) ])];
    %M.hash = @(X) ['z' hashmd5([X.U(:) ; X.S(:) ; X.V(:)])];

    M.rand = @random;
    % Factors U and V live on Stiefel manifolds, hence we will reuse
    % their random generator.
    stiefelm = stiefelfactory(m, k);
    stiefeln = stiefelfactory(n, k);
    function X = random()
        X.U = stiefelm.rand();
        X.V = stiefeln.rand();
        X.S = diag(sort(rand(k, 1), 1, 'descend'));
    end

    % Generate a random tangent vector at X.
    % TODO: consider a possible imbalance between the three components Up,
    % Vp and M, when m, n and k are widely different (which is typical).
    M.randvec = @randomvec;
    function Z = randomvec(X)
        Z.Up = randn(m, k);
        Z.Vp = randn(n, k);
        Z.M  = randn(k);
        Z = tangent(X, Z);
        nrm = M.norm(X, Z);
        Z.Up = Z.Up / nrm;
        Z.Vp = Z.Vp / nrm;
        Z.M  = Z.M  / nrm;
    end

    M.lincomb = @lincomb;

    M.zerovec = @(X) struct('Up', zeros(m, k), 'M', zeros(k, k), ...
                                                        'Vp', zeros(n, k));

    % New vector transport on June 24, 2014 (as indicated by Bart)
    % Reference: Absil, Mahony, Sepulchre 2008 section 8.1.3:
    % For Riemannian submanifolds of a Euclidean space, it is acceptable to
    % transport simply by orthogonal projection of the tangent vector
    % translated in the ambient space.
    M.transp = @project_tangent;
    function Z2 = project_tangent(X1, X2, Z1)
        Z2 = projection(X2, tangent2ambient(X1, Z1));
    end

    M.vec = @vec;
    function Zvec = vec(X, Z)
        Zamb = tangent2ambient(X, Z);
        Zamb_mat = Zamb.U*Zamb.S*Zamb.V';
        Zvec = Zamb_mat(:);
    end
    M.mat = @(X, Zvec) projection(X, reshape(Zvec, [m, n]));
    M.vecmatareisometries = @() true;

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>

    if nargin == 3
        d.Up = a1*d1.Up;
        d.Vp = a1*d1.Vp;
        d.M  = a1*d1.M;
    elseif nargin == 5
        d.Up = a1*d1.Up + a2*d2.Up;
        d.Vp = a1*d1.Vp + a2*d2.Vp;
        d.M  = a1*d1.M  + a2*d2.M;
    else
        error('fixedrank.lincomb takes either 3 or 5 inputs.');
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function brs = ssr_brs(X,phi,x0,d,irs,conf)
%SSR_BRS binaural room scanning (BRS) set for use with the SoundScape Renderer
%
%   Usage: brs = ssr_brs(X,phi,x0,d,irs,conf)
%
%   Input parameters:
%       X       - listener position / m
%       phi     - azimuthal head orientation / rad
%       x0      - secondary sources
%       d       - corresponding driving signals
%       irs     - impulse response data set for the second sources
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       brs     - conf.N x 2*nangles matrix containing all impulse responses
%                 (2 channels) for every angle of the BRS set
%
%   SSR_BRS(X,phi,x0,d,irs,conf) prepares a BRS set for the given secondary
%   sources and its driving signals for the given listener position.
%   One way to use this BRS set is using the SoundScapeRenderer (SSR), see
%   http://spatialaudio.net/ssr/
%
%   See also: ssr_brs_wfs, ssr_brs_nfchoa, ir_generic

%*****************************************************************************
% 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 = 6;
nargmax = 6;
narginchk(nargmin,nargmax);
isargposition(X);
isargsecondarysource(x0);
isargmatrix(d);
isargscalar(phi);
isargstruct(conf);

%% ===== Configuration ===================================================
N = conf.N;                       % Target length of BRIR impulse responses
angles = rad(conf.ir.brsangles);  % Angles for the BRIRs
showprogress = conf.showprogress; % Progress bar

%% ===== Computation =====================================================
nangles = length(angles);
% Initial values
brs = zeros(N,2*nangles);
% Generate a BRS set for all given angles
for ii = 1:nangles
    % Progress bar
    if showprogress, progress_bar(ii,nangles); end
    % Compute BRIR for the desired driving signals
    brs(:,(ii-1)*2+1:ii*2) = ...
        ir_generic(X,angles(ii)+phi,x0,d,irs,conf);
end

author: sfstoolbox

repo: sfs-matlab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% match fvec to atlas by rotating the parameterization of fvec
% - object does not transform in the object space
% - alpha, beta, gamma: samples in the rotation space
% - gran = # of alphas to process together
% - res: base res R + step of hierarchy Hs + depth of hierarchy Hd + top N
%

function [best_fvec, best_angl] = match_param_hie(fvec,atlas,A,B,G,gran,res,max_d)

% base resolution R, step of hierarchy Hs, depth of hierarchy Hd, top N
R = res(1); Hs = res(2); Hd = res(3); N = res(4);
disp(sprintf('----- Base Res %d, Step %d, Depth %d -----',R,Hs,Hd));
alpha = A{1}; beta = B{1}; gamma = G{1};
[fvec, Aprm] = match_param(fvec,atlas,alpha,beta,gamma,gran,N,max_d);
fvec = fvec(:,:,1); Aprm = Aprm(1,:); % remove later on
best_angl = Aprm;
for i=1:Hd
    disp(sprintf('----- Res %d -----',R+Hs*i));
    alpha = A{i+1}; beta = B{i+1}; gamma = G{i+1};
    [fvec, Aprm] = match_param(fvec,atlas,alpha,beta,gamma,gran,1,max_d);
    best_angl(i+1,:) = Aprm;
end
best_fvec = fvec;

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function SDP = chordal_relax_point_cloud_registration(problem)
%% Apply a chordal sparse second-order relaxation to point cloud registration
%% Generate SDP relaxations with multiple smaller blocks
%% This relaxation is likely to be looser than the single block relaxation
%% Depending on multivariate polynomial package SPOT
%% Heng Yang, June 29, 2021

fprintf('\n===================================================================')
fprintf('\nApplying Chordal SDP relaxation to point cloud registration')
fprintf('\n===================================================================\n')
t0              = tic;

%% define POP variables
N               = problem.N;
cloudA          = problem.cloudA;
cloudB          = problem.cloudB;
noiseBoundSq    = problem.noiseBoundSq;
tBound          = problem.translationBound;
tBoundSq        = tBound^2;
barc2           = 1.0;

nrPrimalVars    = 9 + 3 + N;
p               = msspoly('p',nrPrimalVars);
r               = p(1:9);
R               = reshape(r,3,3);
col1            = r(1:3);
col2            = r(4:6);
col3            = r(7:9);
t               = p(10:12);
theta           = p(13:nrPrimalVars);
x               = [r;t];

%% compute the cost function
residuals = [];
for i = 1:N
    ai               = cloudA(:,i);
    bi               = cloudB(:,i);
    res              = bi'*bi + ai'*ai + (t'*t) - 2*(bi'*(R*ai)) - (2)*(bi'*t) + (2)*(t'*(R*ai));
    residuals        = [residuals; res / noiseBoundSq];
end
f_cost = [];
for i = 1:N
    f_cost           = [f_cost; (1+theta(i))/2 * residuals(i) + (1-theta(i))/2 * barc2];
end

%% Define the equality and inequality constraints
h_x = [1.0-col1'*col1;...
       1.0-col2'*col2;...
       1.0-col3'*col3;... % columns unit length
       col1'*col2;...
       col2'*col3;...
       col3'*col1;... % colums orthogonal
       cross(col1,col2) - col3;...
       cross(col2,col3) - col1;...
       cross(col3,col1) - col2]; % columns righthandedness

h_theta = [];
for i = 1:N
    h_theta =[h_theta; 1-theta(i)^2];
end

g_x = tBoundSq - t'*t; % Translation bounded

%% Formulate the chordal sparse second-order relaxation
%% the 0-th block [1;x] * [1;x]'
basis0          = [1;x];
n0              = length(basis0);
basis_x0        = 1;

pop0            = [mykron(basis_x0,h_x);...
                   mykron(basis0,basis0)];
[~,degmat,coef_all] = decomp(pop0);
coef_all            = coef_all';
dim_loc0        = length(basis_x0) * length(h_x);
n0delta         = triangle_number(n0);
nterms          = size(degmat,1);
m_mom0          = n0delta - nterms;

assert(m_mom0==0,'The zero-th blk should have 0 moment constraints.')

coef_mom    = coef_all(:,dim_loc0+1:end);
coef_mom    = coef_mom';
B           = {};
B_normalize = {};

for i = 1:nterms
    [row,~,~]   = find(coef_mom(:,i));
    SDP_coli    = floor((row-1)./n0) + 1;
    SDP_rowi    = mod(row-1,n0) + 1;
    nnz         = length(SDP_rowi);

    Bi          = sparse(SDP_rowi,SDP_coli,ones(nnz,1),n0,n0);
    B{end+1}    = Bi;
    B_normalize{end+1} = Bi/nnz;
end

coef_loc0       = coef_all(:,1:dim_loc0);
A0_local        = {};

for i = 1:dim_loc0
    [rowi,~,vi] = find(coef_loc0(:,i));
    Ai      = sparse(n0,n0);
    for j   = 1:length(rowi)
        Ai  = Ai + vi(j) * B_normalize{rowi(j)};
    end
    A0_local = [A0_local;{Ai}];
end
A0_0    = sparse([1],[1],[1],n0,n0);
% The first block satisfies A0(X0) = b0;
A0      = [{A0_0};A0_local];
b0      = sparse(1,1,1,length(A0),1);

%% the 1-N blocks [1;x;theta(i);theta(i)*x]' * [1;x;theta(i);theta(i)*x]
%% Since there is an inequality constraint too, it will generate 2*N blocks
Aall            = {};
A1all           = {};
ball            = [];
Call            = {};
A0append        = {};
for blkidx = 1:N
    basis       = [1;x;theta(blkidx);theta(blkidx)*x];
    n           = length(basis);
    basis_x     = monomials(theta(blkidx),1:2);
    basis_theta = monomials(x,0:2);
    basis_g     = [1;theta(blkidx)];

    out         = gen_chordal_subblk_pcr(...
                    basis,basis_x,h_x,basis_theta,h_theta(blkidx),g_x,basis_g,f_cost(blkidx));
    Acell       = out.A;
    A1          = out.A1;
    b           = out.b;
    C           = out.C;

    % add constraint that the top-left [1;x]*[1;x]' block is the same as
    % the 0-th block
    A0blk = {};
    for i = 1:n0
        for j = i:n0
            if i == j
                A0i  = sparse(i,j,-1,n0,n0);
                Ai   = sparse(i,j,1,n,n);
            else
                A0i  = sparse([i,j],[j,i],[-0.5,-0.5],n0,n0);
                Ai   = sparse([i,j],[j,i],[0.5,0.5],n,n);
            end
            A0blk    = [A0blk;{A0i}];
            Acell    = [Acell;{Ai}];
        end
    end

    ball             = [ball;b;sparse(n0delta,1)];
    A0append{end+1}  = A0blk;
    Aall{end+1}      = Acell;
    Call             = [Call;C];
    A1all{end+1}     = A1;
end

%% Convert to standard SDPT3 format
b           = [b0;ball];
blk         = cell(2*N+1,2);
blk{1,1}    = 's'; blk{1,2} = n0;
n1          = out.blk{2,2};
n1delta     = triangle_number(n1);
for i = 1:N
    blk{2*i,1} = 's';
    blk{2*i,2} = n;
    blk{2*i+1,1} = 's';
    blk{2*i+1,2} = n1;
end

A0t     = sparsesvec(blk(1,:),A0);
for i = 1:N
    A0t = [A0t,...
           sparse(n0delta,out.m),...
           sparsesvec(blk(1,:),A0append{i})];
end
ndelta  = triangle_number(n);
At    = {A0t};
for i = 1:N
    Ait = [sparse(ndelta,length(b0)),...
           sparse(ndelta,(i-1)*length(Aall{i})),...
           sparsesvec(blk(2*i,:),Aall{i}),...
           sparse(ndelta,(N-i)*length(Aall{i}))];

    A1it = [sparse(n1delta,length(b0)),...
            sparse(n1delta,(i-1)*length(Aall{i})),...
            sparse(n1delta,out.m_mom+out.m_loc),...
            sparsesvec(blk(2*i+1,:),A1all{i}),...
            sparse(n1delta,n0delta),...
            sparse(n1delta,(N-i)*length(Aall{i}))];
    At  = [At;{Ait};{A1it}];
end

SDP.blk = blk;
SDP.At  = At;
SDP.n   = n;
SDP.m   = length(b);
SDP.C   = [{sparse(n0,n0)};Call];
SDP.b   = b;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')

%% Convert to sedumi
fprintf('Convert to sedumi')
t0 = tic;

sK.s  = [n0];
for i = 1:N
    sK.s        = [sK.s,n,n1];
end

A0t     = sparsevec(blk(1,:),A0);
n0sq    = n0^2;
for i = 1:N
    A0t = [A0t,...
           sparse(n0sq,out.m),...
           sparsevec(blk(1,:),A0append{i})];
end

nsq     = n^2;
n1sq    = n1^2;
At      = {A0t};
for i = 1:N
    Ait = [sparse(nsq,length(b0)),...
           sparse(nsq,(i-1)*length(Aall{i})),...
           sparsevec(blk(2*i,:),Aall{i}),...
           sparse(nsq,(N-i)*length(Aall{i}))];

    A1it = [sparse(n1sq,length(b0)),...
            sparse(n1sq,(i-1)*length(Aall{i})),...
            sparse(n1sq,out.m_mom+out.m_loc),...
            sparsevec(blk(2*i+1,:),A1all{i}),...
            sparse(n1sq,n0delta),...
            sparse(n1sq,(N-i)*length(Aall{i}))];
    At  = [At;{Ait};{A1it}];
end

sdata.K     = sK;
sdata.At    = cat(1,At{:});
sdata.b     = b;

sc          = [];
for i = 1:length(SDP.C)
    sc      = [sc;sparsevec(blk(i,:),SDP.C(i))];
end
sdata.c     = sc;

SDP.sedumi   = sdata;

tf    = toc(t0);
fprintf('\nDone in %g seconds.\n',tf);
fprintf('===================================================================\n')

end

repo: CertifiablyRobustPerception


text:
function h = tc_sessionPlots(tcs, plotType, titleType);
% h = tc_sessionPlots(tcs, plotType, [titleType = 1]);
%
% produce a figure showing individual subjects time course or amplitude data
% each subject's data will be a subplot of the figure
% plotType (string):
% 'wholetc'         raw timecourse
% 'meanampsplustcs' left panel will show mean amps according to tc.params.ampType; right panel will show event-triggered average
% 'meantcs'         event trigerred averaging mean time course
% 'deconvolvedTcs'  deconvolved time courses
% 'amps'            Means amplitudes of each condition according to the
%                   value in the tc.params.ampType filed which can be
%                   one of the following: 'difference' 'betas' 'relamps' or 'deconvolved';
% 'meanamps'        estimated from the event-triggered average based on the
%                   baseline and peak periods
% 'betas'           one beta per condition (no beta for baseline)estimated
%                   from a glm on the tc
% 'relamps'         David Ress's method for amplitude estimatation
% 'deconvolvedamps' Etimate of the mean amps from the deconvolved time
%                   course based on the peak and baseline period
%
% titleType: 1 (default) uses tcs.sessionCode for chart titles (usually
% corresponds to the directory of the session); 2 uses tcs.description.
%
%
% ras, 02/2007.
% ras, 04/2007: was part of tc_acrossSessions, but broke off to be its own
% function, since it's pretty useful. Also allows cells of tc structs as
% inputs.
% remus 10/2007: added titleType flag.
% kgs 3/2008: added deconvolvedtc and deconvolvedamps options
%
h = [];

% set titleType parameters
if exist('titleType') & titleType == 2
    titleField = 'description';
else
    titleField = 'sessionCode';
end

% allow input to be a cell array
if iscell(tcs)
	% remove empty entries
	tcs = tcs( cellfind(tcs) );

	C = tcs;
	tcs = tcs{1};
	for i = 2:length(C)
		tcs(i) = mergeStructures(tcs(i-1), C{i});
	end
end

nSessions = length(tcs);
maxSubplotsPerFig = 16;  % if length(tcs)>this, make multiple figures
nFigs = ceil(nSessions / maxSubplotsPerFig);
subplotsPerFig = min(nSessions, maxSubplotsPerFig);

if nFigs > 1
    nRows = 5;
    nCols = 4;
else
    nRows = ceil( sqrt(subplotsPerFig) );
    nCols = ceil( subplotsPerFig / nRows );
end

% adjustment for plot types which display >1 axis per session
if isequal(lower(plotType), 'meanampsplustcs')
    nFigs = ceil(2 * nSessions / maxSubplotsPerFig); % 2 subplots per session
    subplotsPerFig = maxSubplotsPerFig / 2;
    nRows = nRows*2;
    nCols = nCols*2;
end

for s = 1:nSessions
    if mod(s, subplotsPerFig)==1
        h(end+1) = figure('Color', 'w', 'Position', [100 100 600 800], ...
                    'Name', sprintf('%s Session Plots', tcs(1).roi.name));
	end

	% get text for the title for this subplot
	titleText = tcs(s).params.(titleField);
	titleText( titleText=='_' ) = '-';  % TeX markup--avoid subscripts

    switch lower(plotType)
        case {'wholetc' 'wholetcs'}
            subplot(nRows, nCols, s);
            tc_plotWholeTc(tcs(s), gca);
            title(titleText, 'FontSize', 10);

        case 'meanampsplustcs'
            subplot(ceil(nRows/2), nCols, (2*s-1));
            tc_barMeanAmplitudes(tcs(s), gca, 0); % this will plot the event triggered average plots
            title(titleText, 'FontSize', 10);

            subplot(ceil(nRows/2), nCols, (2*s));
            tc_plotMeanTrials(tcs(s), gca);

			if isfield(tcs(s), 'SNR') & ~isempty(tcs(s), 'SNR')
				snrNums = tcs(s).condNums(tcs(s).params.snrConds);
				snrConds = sprintf('%i ', snrNums);
				title(sprintf('SNR: %3.2f (conds [%s])', tcs(s).SNR, snrConds));
			end

        case {'meantc' 'meantcs'}
            subplot(nRows, nCols, s);
            tc_plotMeanTrials(tcs(s), gca); % this will plot what is in the meantc field
            title(titleText, 'FontSize', 10);

        case {'deconvolved' 'deconvolvedTcs' 'deconvolvedtcs'}  % will plot the deconvolved time courses
                currTc=tcs(s);
     			currTc.params.glmHRF=0 % set glm flag for deconvolution
                % deconvolve
                 currTc = tc_applyGlm(currTc);
                 % add a blank condition at the beginning betas are zero
                % because this is the baseline to which the glm is
                 % estimated
                currTc.meanTcs = [zeros(currTc.glm.nh, 1) currTc.glm.betas];
                subplot(nRows, nCols, s);
                tc_plotMeanTrials(currTc, gca);
                title(titleText, 'FontSize', 10);

        case 'amps'
            subplot(nRows, nCols, s);
            %   0 - Use the value in the tc.params.ampType field.
            tc_barMeanAmplitudes(tcs(s), gca, 0);
            title(titleText, 'FontSize', 10);

        case 'meanamps'
            subplot(nRows, nCols, s);
            %  1 - Plot Peak/Baseline difference of event-triggered average
            tc_barMeanAmplitudes(tcs(s), gca, 1);
            title(titleText, 'FontSize', 10);

        case 'betas'
            subplot(nRows, nCols, s);
            %   2 - Plot GLM Beta values
            tc_barMeanAmplitudes(tcs(s), gca, 2);
            title(titleText, 'FontSize', 10);

        case {'relamps' 'dotproductamps' 'projamps'}
            subplot(nRows, nCols, s);
            %   3 - Plot Dot-Product Projection Amplitudes
            tc_barMeanAmplitudes(tcs(s), gca, 3);
            title(titleText, 'FontSize', 10);

        case { 'deconvolvedamps'  'deconvolvedAmps'}
            subplot(nRows, nCols, s);
            %   4 - Plot Peak/Baseline from deconvolved time courses
            tc_barMeanAmplitudes(tcs(s), gca, 4);
            title(titleText, 'FontSize', 10);

        otherwise, warning( sprintf('Unknown plot type %s', plotType) )
    end

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_zero_rc_all ( )

%*****************************************************************************80
%
%% TEST_ZERO_RC_ALL tests ZERO_RC on all test functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_ZERO_RC_ALL\n' );
  fprintf ( 1, '  Test the Brent ZERO_RC routine, which seeks\n' );
  fprintf ( 1, '  a root of a function F(X)\n' );
  fprintf ( 1, '  in an interval [A,B].\n' );

  machep = eps;
  t = machep;

  a = 1.0;
  b = 2.0;

  test_zero_rc_one ( a, b, machep, t, @f_01, ...
    'f_01(x) = sin ( x ) - x / 2' );

  a = 0.0;
  b = 1.0;

  test_zero_rc_one ( a, b, machep, t, @f_02, ...
    'f_02(x) = 2 * x - exp ( - x )' );

  a = -1.0;
  b =  0.5;

  test_zero_rc_one ( a, b, machep, t, @f_03, ...
    'f_03(x) = x * exp ( - x )' );

  a =  0.0001;
  b =  20.0;

  test_zero_rc_one ( a, b, machep, t, @f_04, ...
    'f_04(x) = exp ( x ) - 1 / ( 100 * x * x )' );

  a = -5.0;
  b =  2.0;

  test_zero_rc_one ( a, b, machep, t, @f_05, ...
    'f_05(x) = (x+3) * (x-1) * (x-1)' );

  return
end


text:
function [der_loc, der_dir] = msi_file_DER_loc(set_fid)
%function [der_loc, der_dir] = msi_file_DER_loc(set_fid, chans)
%	MSI_FILE_GET_LONG(set_fid, chans)	Skip to a keyword in a file FID
%	and return the single float value
%

key_count = msi_file_find_keyword(set_fid, 'MSI.DerChanCount:');
if (key_count == 1)
    DerChanCount = msi_file_get_long(set_fid, 'MSI.DerChanCount:');
end

frewind(set_fid);
keyword = 'MSI.Der_Position_Information.Begin:';

while feof(set_fid) == 0
    	line = fgetl(set_fid);
	matches = findstr(line, keyword);
	num = length(matches);
	if num > 0
	    fprintf(1, 'I found a %s\n', line);
	    for (ii=1:DerChanCount)
		line = fgetl(set_fid);
		[name, count, errmsg, nextindex] = sscanf(line, '%s\t', 1);
		[token, nextline] = strtok(line, '	');
		[bigresult, count, errmsg, nextindex] = sscanf(nextline, '\t%f\t%f\t%f\t%f\t%f\t%f', 6);
		loc(1) = bigresult(1);
		loc(2) = bigresult(2);
		loc(3) = bigresult(3);
		dir(1) = bigresult(4);
		dir(2) = bigresult(5);
		dir(3) = bigresult(6);
		%fprintf(1, 'chan: %s\n',  name);
		%fprintf(1, '  %f %f %f   ',  loc(1), loc(2), loc(3));
		%fprintf(1, '  %f %f %f \n',  dir(1), dir(2), dir(3));
		%fprintf(1, 'chan: %s (%f,%f,%f) (%f,%f,%f)\n',  name, loc[1], loc[2], loc[3], dir[1], dir[2], dir[3]);
		outloc(ii,1) = loc(1);
		outloc(ii,2) = loc(2);
		outloc(ii,3) = loc(3);

		outdir(ii,1) = dir(1);
		outdir(ii,2) = dir(2);
		outdir(ii,3) = dir(3);
	    end
der_loc = outloc;
der_dir = outdir;
	end
end

author: PatternRecognition

repo: OpenBMI


text:
function varargout = sltakeval(A)
%SLTAKEVAL Extracts the values in an array/cell array to output
%
% $ Syntax $
%   -[O1, O2, ..., On] = sltakeval(A)
%
% $ Description $
%   -[O1, O2, ..., On] = sltakeval(A) takes the values in the array
%    or cell array A and assigns them to output variables.
%    The number of elements in A should be equal to the number of
%    output arguments.
%
% $ History $
%   - Created by Dahua Lin, on Sep 1st, 2006
%

n = numel(A);
if nargout ~= n
    error('sltoolbox:sizmismatch', ...
        'The number of elements in A is not equal to the number of outputs');
end

if n > 0
    if iscell(A)
        if isequal(size(A), [1, n])
            varargout = A;
        else
            varargout = reshape(A, 1, n);
        end
    else
        varargout = cell(1, n);
        for i = 1 : n
            varargout{i} = A(i);
        end
    end
end

author: lmthang


text:
function qInv=quatInv(q)
%%QUATINV Find the inverse of a quaternion. The inverse of a quaternion
%         left-or-right multiplied by the quaternion (using quatMult) is
%         equal to one (a real number). As with real numbers, issues arise
%         taking the inverse of the zero quaternion. Note that the inverse
%         does not depend on the handedness of the quaternion algebra, even
%         though the quatMult function does.
%
%%INPUTS: q  A 4XN set of N quaternions, each whose inverse is desired,
%           where the first element in each column is the scalar part of
%           the quaternion (sometimes called q0 or q4) and the next three
%           elements are the (hypercomplex) vector part. That is, the
%           hypercomplex quaternion given by q(:,1) can be written in
%           hypercomplex, non-vector form as
%           q(1,1)+i*q(2,1)+j*q(3,1)+k*q(4,1), where i, j, and k are all
%           roots of -1.
%
%OUTPUTS: qInv The 4XN set of inverses of the quaternions.
%
%Properties of quaternions including inversion are described in [1]. When a
%quaternion has unit magnitude, its conjugate it also its inverse.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Quaternion." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/Quaternion.html
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

qC=quatConj(q);
%The multiplication of a quaternion with its conjugate produces scalar
%values given in the first entries of qMags
qMags=quatMult(q,qC);
qMags=qMags(1,:);

qInv=bsxfun(@rdivide,qC,qMags);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%HOMTRANS Apply a homogeneous transformation
%
% P2 = HOMTRANS(T, P) applies the homogeneous transformation T to the points
% stored columnwise in P.
%
% - If T is in SE(2) (3x3) and
%   - P is 2xN (2D points) they are considered Euclidean (R^2)
%   - P is 3xN (2D points) they are considered projective (P^2)
% - If T is in SE(3) (4x4) and
%   - P is 3xN (3D points) they are considered Euclidean (R^3)
%   - P is 4xN (3D points) they are considered projective (P^3)
%
% P2 and P have the same number of rows, ie. if Euclidean points are given
% then Euclidean points are returned, if projective points are given then
% projective points are returned.
%
% TP = HOMTRANS(T, T1) applies homogeneous transformation T to the
% homogeneous transformation T1, that is TP=T*T1.  If T1 is a 3-dimensional
% transformation then T is applied to each plane as defined by the first two
% dimensions, ie. if T is NxN and T1 is NxNxM then the result is NxNxM.
%
% Notes::
% - If T is a homogeneous transformation defining the pose of {B} with respect to {A},
%   then the points are defined with respect to frame {B} and are transformed to be
%   with respect to frame {A}.
%
% See also E2H, H2E, RTBPose.mtimes.

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

    if numrows(p) == numrows(T)
        if ndims(p) == 3
            % homtrans(T1, T2)
            pt = [];
            for i=1:size(p,3)
                pt = cat(3, pt, T*p(:,:,i));
            end
        else
            % points are in projective coordinates
            pt = T * p;
        end
    elseif (numrows(T)-numrows(p)) == 1
        % points are in Euclidean coordinates, promote to homogeneous
        pt = h2e( T * e2h(p) );
    else
        error('SMTB:homtrans:badarg', 'matrices and point data do not conform')
    end

author: petercorke

repo: spatialmath-matlab


text:
% *************************************************************************************
% * DESCRIPTION:                                                                      *
% * This script computes parts of the experiments performed in ref. [1]               *
% * ('fast' version).                                                                 *
% * --------------------------------------------------------------------------------- *
% * REFERENCE:                                                                        *
% * [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and         *
% *     MRI texture features for the prediction of lung metastases in soft-tissue     *
% *     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),         *
% *     5471-5496. doi:10.1088/0031-9155/60/14/5471                                   *
% * --------------------------------------------------------------------------------- *
% * AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>                            *
% * --------------------------------------------------------------------------------- *
% * HISTORY:                                                                          *
% * - Creation: May 2015                                                              *
% * --------------------------------------------------------------------------------- *
% * STATEMENT:                                                                        *
% * This file is part of <https://github.com/mvallieres/radiomics/>,                  *
% * a package providing MATLAB programming tools for radiomics analysis.              *
% * --> Copyright (C) 2015  Martin Vallieres                                          *
% *                                                                                *
% *   This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.           *
% *************************************************************************************

fprintf('\n')
help masterScriptExample_STS
fprintf('\n')
warning off

% INITIALIZATION
pathWORK = pwd;
roiNumb = load('contour_Mass'); roiNumb = struct2cell(roiNumb); roiNumb = roiNumb{1};
outcome = load('outcome'); outcome = struct2cell(outcome); outcome = outcome{1};
nPatient = numel(outcome);

% TEXTURE EXTRACTION PARAMETERS AND DEGREES OF FREEDOM
MRIinv_cell = {'NoInv','Inv'};
MRIweight_mat = [1/3,2/3];
R_mat = [2/3,3/2];
scale_cell = {4,5};
algo_cell = {'Equal','Lloyd'};
Ng_mat = [8,16];
paramSEP = {R_mat,scale_cell,algo_cell,Ng_mat};                           % NOTE: paramSEP and paramFUS must always be of the same format and size for SEPARATE
paramFUS = {MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat}; % and FUSED SCANS as presented here, with the same ordering of different extraction parameters.
baselineSEP = [2,2,1,1]; % Arbitrarily chosen for the example
baselineFUS = [1,1,2,2,1,1]; % Arbitrarily chosen for the example
freedomSEP = [1 1 1 1 ; 0 0 0 0 ; 0 1 0 1]; % Arbitrarily chosen for the example
freedomFUS = [1 1 1 1 1 1 ; 0 0 0 0 0 0 ;  0 1 0 1 0 1]; % Arbitrarily chosen for the example

% MULTIVARIABLE ANALYSIS PARAMETERS
nBoot = 100;
alpha = 0.5; delta = 0.5;
setSize = 25;
fSetName = {'PET','SEPARATE','FUSED'};
maxOrder = 10;

% ************** START COMPUTATION **************

% 1. READ DATA DOWNLOADED FROM THE TCIA WEBSITE (http://dx.doi.org/10.7937/K9/TCIA.2015.7GO2GSKS)
fprintf('\n\n******************* ORGANIZING AND PROCESSING DICOM DATA FROM TCIA WEBSITE *******************\n')
readAllDICOM_STS([pathWORK,'/Soft-tissue-Sarcoma'],nPatient)

% 2. COMPUTE NON-TEXTURE FEATURES
fprintf('\n\n*********************** COMPUTING NON-TEXTURE FEATURES ***********************\n')
calcAllNonTextureFeatures_STS(pathWORK,nPatient,roiNumb,outcome)

% 3. COMPUTE AND ORGANIZING ALL TEXTURE FEATURES
mkdir('TEXTURES'), fprintf('\n')
calcAllSeparateTextures_STS(pathWORK,nPatient,roiNumb,R_mat,scale_cell,algo_cell,Ng_mat)
organizeSeparateTextures_STS(pathWORK,nPatient,R_mat,scale_cell,algo_cell,Ng_mat)
calcAllFusedTextures_STS(pathWORK,nPatient,roiNumb,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
organizeFusedTextures_STS(pathWORK,nPatient,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)

% 5. PERFORM FEATURE SET REDUCTION
mkdir('FSET')
nonTextures = load('nonTextures'); nonTextures = struct2cell(nonTextures); nonTextures = nonTextures{1};
PET  = load('textures_PET');  PET  = struct2cell(PET);  PET  = PET{1};
T1   = load('textures_T1');   T1   = struct2cell(T1);   T1   = T1{1};
T2FS = load('textures_T2FS'); T2FS = struct2cell(T2FS); T2FS = T2FS{1};
PET_T1   = load('textures_PET_T1');   PET_T1   = struct2cell(PET_T1);   PET_T1   = PET_T1{1};
PET_T2FS = load('textures_PET_T2FS'); PET_T2FS = struct2cell(PET_T2FS); PET_T2FS = PET_T2FS{1};

tic
fprintf('\n\nFinding the path to the ''MINE.jar'' application on the system ... ')
pathMINE = findMINE('Linux');
fprintf('DONE\n')
toc

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''PET'' FEATURE SET ***********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{1},outcome,setSize,nonTextures,{PET},{'PET'},paramSEP,freedomSEP,baselineSEP,alpha,delta,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''SEPARATE'' FEATURE SET **********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{2},outcome,setSize,nonTextures,{PET,T1,T2FS},{'PET','T1','T2FS'},paramSEP,freedomSEP,baselineSEP,alpha,delta,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING FEATURE SET REDUCTION FOR ''FUSED'' FEATURE SET ***********************\n')
calcAllFeatureSets_STS(pathWORK,pathMINE,fSetName{3},outcome,setSize,nonTextures,{PET_T1,PET_T2FS},{'PET_T1','PET_T2FS'},paramFUS,freedomFUS,baselineFUS,alpha,delta,nBoot)

% 6. PERFORM FEATURE SET SELECTION
mkdir('MODELS')
fprintf('\n\nFinding the path to''fastAUC.cpp'' on the system --> COMPILATION ... ')
try
    compileFastAUC('Linux')
    fprintf('DONE\n')
catch
    fprintf('FAILED (AUC computations will be slower)\n')
end

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING FEATURE SELECTION FOR ''PET'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{1},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING FEATURE SELECTION  FOR ''SEPARATE'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{2},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING FEATURE SELECTION  FOR ''FUSED'' FEATURE SET ***********************\n')
computeAllModelChoice_STS(pathWORK,fSetName{3},outcome,freedomFUS,maxOrder,nBoot)

% 7. PERFORM PREDICTION PERFORMANCE ESTIMATION
mkdir('RESULTS'), fprintf('\n')

% For Feature set 1: PET
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''PET'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{1},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 2: PET, T1, T2FS
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''SEPARATE'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{2},outcome,freedomSEP,maxOrder,nBoot)

% For Feature set 3: PET_T1, PET_T2FS
fprintf('\n*********************** PERFORMING PREDICTION PERFORMANCE ESTIMATION FOR ''FUSED'' FEATURE SET ***********************\n')
computeAllPrediction_STS(pathWORK,fSetName{3},outcome,freedomFUS,maxOrder,nBoot)

% Finding the best combinations of model order and texture extraction parameter degree of freedom for all feature set types
groupExperiments_STS(pathWORK,fSetName{1},freedomSEP,maxOrder)
groupExperiments_STS(pathWORK,fSetName{2},freedomSEP,maxOrder)
groupExperiments_STS(pathWORK,fSetName{3},freedomFUS,maxOrder)

% 8. CHOICE OF BEST PARSIMONIOUS MODEL (requires user input)
fprintf('\n')
plotPredictionResults_STS([pathWORK,'/RESULTS'],fSetName,{'AUC632','Sensitivity632','Specificity632'},maxOrder)
while 1
    set = input(['\nWhich feature set provides the best parsimonious model? \n' ...
                 '--> For the ',fSetName{1},' feature set, type ''1'' and press ENTER \n' ...
                 '--> For the ',fSetName{2},' feature set, type ''2'' and press ENTER \n' ...
                 '--> For the ',fSetName{3},' feature set, type ''3'' and press ENTER \n' ...
                 'ANSWER: ']);
    fprintf('\n')
    if isnumeric(set) && (set == 1 || set == 2 || set == 3)
        break
    end
end
while 1
    order = input(['Which model order of the ',fSetName{set},' feature set provides the best parsimonious model? \n' ...
                   '--> Type a number between 1 to ',num2str(maxOrder),' and press ENTER \n' ...
                   'ANSWER: ']);
    fprintf('\n')
    if isnumeric(order) && order <= 10 && order >= 1
        break
    end
end
cd([pathWORK,'/RESULTS'])
results = load(['RESULTS_',fSetName{set},'_BEST']); results = struct2cell(results); results = results{1};
finalModel = results.(['Order',num2str(order)]);
cd(pathWORK), mkdir('FINAL_MODEL'), cd('FINAL_MODEL'), save('finalModel', 'finalModel')

% 9. COMPUTING THE LOGISTIC REGRESSION COEFFICIENTS AND BOOTSTRAP CONFIDENCE INTERVALS OF THE FINAL MODEL
fprintf('\n\nCOMPUTING THE LOGISTIC REGRESSION COEFFICIENTS OF THE FINAL MODEL ... ')
tic
[coeff,response,modelCI] = computeModelCoefficients(finalModel.Data,outcome,'IABR');
save('coeff','coeff'), save('response','response'), save('modelCI','modelCI')
fprintf('DONE\n')
toc

% 10. DISPLAYING THE FINAL MODEL AND CORRESPONDING PREDICTION PERFORMANCE ESTIMATION
plotSigmoidalResponse(response,outcome,modelCI,'LungMets')
fprintf(['\n\n\n --> THE FINAL MULTIVARIABLE MODEL IS:\n\n'...
         '               g(x) =               \n'])
for i = 1:order
    fprintf([num2str(coeff(i)),' X ',finalModel.Name{i},'\n'])
    fprintf('                    +               \n')
end
fprintf(['                   ',num2str(coeff(end)),'\n'])
fprintf('\nWITH CORRESPONDING PREDICTION PERFORMANCE ESTIMATION:\n')
fprintf(['AUC = ',num2str(roundsd(finalModel.AUC632,ceil(log10(finalModel.AUC632/roundsd(finalModel.SE_AUC632,1))))),' ± ',num2str(roundsd(finalModel.SE_AUC632,1)),'\n'])
fprintf(['Sensitivity = ',num2str(roundsd(finalModel.Sensitivity632,ceil(log10(finalModel.Sensitivity632/roundsd(finalModel.SE_Sensitivity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Sensitivity632,1)),'\n'])
fprintf(['Specificity = ',num2str(roundsd(finalModel.Specificity632,ceil(log10(finalModel.Specificity632/roundsd(finalModel.SE_Specificity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Specificity632,1)),'\n'])
fprintf('\n')
cd(pathWORK)

author: mvallieres

repo: radiomics


text:
function[]=makefigs_twodstats
%MAKEFIGS_TWODSTATS  Makes a sample figure for TWODSTATS.

if isempty(which('drifters.mat'))
    disp('Sorry, TWODSTATS can''t find DRIFTERS.MAT.')
    return
end

%This make take a few minutes...

load drifters,use drifters
cv=cellpair(u,v);
vindex(lat,lon,cv,1:10:length(lat),1);

tic;[mat,xmid,ymid]=twodstats(lon,lat,cellabs(cv),-180.5:180.5,-89.5:89.5);etime1=toc;

figure,jpcolor(xmid,ymid,mat)
axis([-180 180 -70 90]),latratio(30),topoplot
xlabel('Longitude'),ylabel('Latitude')
title('Mean Speed from the Global Surface Drifter Dataset')
caxis([8 60]),h=colorbar('EastOutside');colormap lansey
h.Label.String='Mean Speed (cm/s)';

%To print
if 0
    currentdir=pwd;
    cd([whichdir('jlab_license') '/figures'])
    print -dpng twodstats
    crop twodstats.png
    cd(currentdir)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
classdef latticeType < int32
% class representing the different Bravais lattices

  enumeration
    triclinic    (1)
    monoclinic   (2)
    orthorhombic (3)
    trigonal     (4)
    tetragonal   (5)
    hexagonal    (6)
    cubic        (7)
    none         (8)
  end

  methods

    function abg = defaultAngles(this)

      switch this
        case {'trigonal','hexagonal'}
          abg = [90 90 120] * degree;

        otherwise
          abg = [90 90 90] * degree;
      end

    end

    function out = isTriHex(this)

      out = this == latticeType.trigonal || this == latticeType.hexagonal;

    end

    function out = isEucledean(this)

      out = this == latticeType.orthorhombic || ...
        this == latticeType.tetragonal || this == latticeType.cubic;

    end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
% This subroutine assigns creates a grid with
% spacing dx,dy (in degrees). The size will
% be selected interactively. The bvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/95

% TODO delete this. I'm higly suspicious this plays a role in anything. Message2 was not fully defined -CGR
report_this_filefun(mfilename('fullpath'));

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 1.00;
    dy = 1.00 ;
    ni = 100;

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ wex+200 wey-200 450 250]);
    axis off

    % creates a dialog box to input grid parameters
    %
    freq_field=uicontrol('Style','edit',...
        'Position',[.60 .50 .22 .10],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.60 .40 .22 .10],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.60 .30 .22 .10],...
        'Units','normalized','String',num2str(dy),...
        'Callback','dy=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dy));');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback','close,sel =''ca'', bvalgrid',...
        'String','Go');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.74 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.42 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in x (dx) in deg:');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.32 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in y (dy) in deg:');

    txt1 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.53 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m,...
        'FontWeight','bold',...
        'String','Number of Events (Ni):');
    set(gcf,'visible','on');
    watchoff

end   % if nargin ==0

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'
    message2 = ['right corner (with same mouse button). Allow  '
        'some time to complete calculation of curves.  '];
    try
        zmap_message_center.set_message(' ',message2);
    end

    figure_w_normalized_uicontrolunits(map)
    [x0,y0]  = ginput(1);
    mark1 =    plot(x0,y0,'ro','era','normal');
    set(mark1,'MarkerSize',10,'LineWidth',2.0)
    [x1,y1]  = ginput(1);
    f = [x0 y0 ; x1 y0 ; x1 y1 ; x0 y1 ; x0 y0];

    fplo = plot(f(:,1),f(:,2),'r','era','normal');
    set(fplo,'LineWidth',2)
    gx = x0:dx:x1;
    gy = y0:dy:y1;
    itotal = length(gx) * length(gy);

    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    cumu = zeros(length(t0b:par1/365:teb)+2);
    ncu = length(cumu);
    cumuall = zeros(ncu,length(gx)*length(gy));
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    bvg = [];
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','b-value grid - percent done');;
    drawnow
    %
    % longitude  loop
    %
    for x =  x0:dx:x1
        i1 = i1+ 1;

        % latitude loop
        %
        for  y = y0:dy:y1
            allcount = allcount + 1.;
            i2 = i2+1;

            % calculate distance from center point and sort wrt distance
            l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2) ;
            [s,is] = sort(l);
            b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

            % take first ni points
            b = b(1:ni,:);      % new data per grid point (b) is sorted in distance

            % call the b-value function
            [bv, magco] =  bvalcalc(b);
            l = sort(l);
            bvg = [bvg ; bv magco x y l(ni)];
            waitbar(allcount/itotal)
        end  % for y0
        i2 = 0;
    end  % for x0

    save  bvalgrid.mat bvg gx gy ni dx dy

    close(wai)
    watchoff

    % plot the results
    %
    % old and re3 (initially ) is the b-value matrix
    re3 = reshape(bvg(:,1),length(gy),length(gx));
    r = reshape(bvg(:,5),length(gy),length(gx));
    old = re3;
    % old1 is the magnitude of completness matrx
    old1 = reshape(bvg(:,2),length(gy),length(gx));
    view_bva

end   % if nargin ==3

author: CelsoReyes


text:
% project the 3D points to generate 2D points according to the viewpoint
function xproj = projectp3d(x3d, object)

if isfield(object, 'viewpoint') == 1
    % project the 3D points
    viewpoint = object.viewpoint;
    a = viewpoint.azimuth*pi/180;
    e = viewpoint.elevation*pi/180;
    d = viewpoint.distance;
    f = viewpoint.focal;
    theta = viewpoint.theta*pi/180;
    principal = [viewpoint.px viewpoint.py];
    viewport = viewpoint.viewport;
else
    xproj = [];
    return;
end

if d == 0
    xproj = [];
    return;
end

% camera center
C = zeros(3,1);
C(1) = d*cos(e)*sin(a);
C(2) = -d*cos(e)*cos(a);
C(3) = d*sin(e);

% Rotate coordinate system by theta is equal to rotating the model by -theta.
a = -a;
e = -(pi/2-e);

% rotation matrix
Rz = [cos(a) -sin(a) 0; sin(a) cos(a) 0; 0 0 1];   %rotate by a
Rx = [1 0 0; 0 cos(e) -sin(e); 0 sin(e) cos(e)];   %rotate by e
R = Rx*Rz;

% perspective project matrix
% however, we set the viewport to 3000, which makes the camera similar to
% an affine-camera. Exploring a real perspective camera can be a future work.
M = viewport;
P = [M*f 0 0; 0 M*f 0; 0 0 -1] * [R -R*C];
% project
x = P*[x3d ones(size(x3d,1), 1)]';
dist = max(x(1,:))-min(x(1,:));
x(1,:) = x(1,:) ./ x(3,:);
x(2,:) = x(2,:) ./ x(3,:);
new_dist = max(x(1,:))-min(x(1,:));
%x(1,:) = x(1,:) / mean(x(3,:));
%x(2,:) = x(2,:) / mean(x(3,:));
xz = x(3,:)*(M*f)*new_dist/dist;
x = x(1:2,:);
% rotation matrix 2D
R2d = [cos(theta) -sin(theta); sin(theta) cos(theta)];
x = (R2d * x)';
% x = x';

% transform to image coordinates
x(:,2) = -1 * x(:,2);
x = x + repmat(principal, size(x,1), 1);
xz = (xz-mean(xz))';
xproj = [x xz];

repo: CategoryShapes


text:
function [z, r] = fisherz(r)
% Fisher's r to z' transform, and the inverse
%
% :Outputs:
%
%   **z:**
%        z = z', treating input r as correlation
%
%  **r:**
%        treating input r as a z' score

    z = .5 .* log((1+r) ./ (1-r));     % Fisher's r-to-z transform

    if nargout > 0
        r = (exp(2.*r) - 1) ./ (exp(2.*r) + 1);    % inverse
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
close all;
clear all;
clc

A=imread('lowcc.jpg');
A=rgb2gray(A);

p=0:255;
C=hist(double(A(:)),256);
cuf=zeros(size(C));

for i=1:length(C)
    if i==1
        cuf(i)=C(i);
    else
        cuf(i)=cuf(i-1)+C(i);
    end
end
%For Histogram Matching give reference histogram to C2 variable
C2=(cuf(end)/length(cuf))*ones(size(C));
cuf2=zeros(size(C));
for i=1:length(C2)
    if i==1
        cuf2(i)=C2(i);
    else
        cuf2(i)=cuf2(i-1)+C2(i);
    end
end
q=zeros(size(p));
for i=1:length(C)
    t=closeone(cuf2,cuf(i));
    q(i)=p(t);
end
B=interp1(p,q,single(A(:)));
B=reshape(B,size(A));
figure;
imshow(A);
t=min(min(B));
h=max(max(B));
p=size(B);
B=floor(((B-t)/h)*255);
B=uint8(B);
figure;
imshow(B);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function poly_cof = r8poly_shift ( scale, shift, n, poly_cof )

%*****************************************************************************80
%
%% R8POLY_SHIFT adjusts the coefficients of a polynomial for a new argument.
%
%  Discussion:
%
%    Assuming P(X) is a polynomial in the argument X, of the form:
%
%      P(X) =
%          C(N) * X**(N-1)
%        + ...
%        + C(2) * X
%        + C(1),
%
%    and that Z is related to X by the formula:
%
%      Z = SCALE * X + SHIFT
%
%    then this routine computes coefficients C for the polynomial Q(Z):
%
%      Q(Z) =
%          C(N) * Z**(N-1)
%        + ...
%        + C(2) * Z
%        + C(1)
%
%    so that:
%
%      Q(Z(X)) = P(X)
%
%  Example:
%
%    P(X) = 2 * X**2 - X + 6
%
%    Z = 2.0 * X + 3.0
%
%    Q(Z) = 0.5 *         Z**2 -  3.5 * Z + 12
%
%    Q(Z(X)) = 0.5 * ( 4.0 * X**2 + 12.0 * X +  9 )
%            - 3.5 * (               2.0 * X +  3 )
%                                            + 12
%
%            = 2.0         * X**2 -  1.0 * X +  6
%
%            = P(X)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real SHIFT, SCALE, the shift and scale applied to X,
%    so that Z = SCALE * X + SHIFT.
%
%    Input, integer N, the order of the polynomial.
%
%    Input, real POLY_COF(N), the coefficient array in terms of the X variable.
%
%    Output, real POLY_COF(N), the coefficient array in terms of the Z variable.
%
  for i = 1 : n
    poly_cof(i+1:n) = poly_cof(i+1:n) / scale;
  end

  for i = 1 : n
    for j = n-1 : -1 : i
      poly_cof(j) = poly_cof(j) - shift * poly_cof(j+1);
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s=cosmo_type(fn)
% print or return ASCII contents of a file
%
% Usages
%  - cosmo_type(fn);     prints the contents of file fn
%  - s=cosmo_type(fn);   return the contents of file fn
%
% See also: type
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    fid=fopen(fn);
    s=fread(fid,inf,'char=>char')';
    fclose(fid);

    if nargout==0
        fprintf('%s',s);
        clear s;
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [C, log] = SOFAmerge(A,B)
%SOFAmerge
%   [C, log] = SOFAmerge(A,B) merges the SOFA objects A and B to a single one, C.
%
%   A and B are structs containing the data and meta. A and B
%   must be of the same SOFA conventions.
%
%   C is a struct containing merged data
%   log contains log data

% SOFA API - function SOFAupdateDimensions
% Copyright (C) 2012-2013 Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.1 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: http://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Initial check
if ~strcmp(A.GLOBAL_SOFAConventions,B.GLOBAL_SOFAConventions)
	error('Both SOFA objects must use the same SOFA conventions');
end
if A.API.N~=B.API.N
	error('Data size of both SOFA objects must be the same');
end
Def = SOFAdefinitions;
log={''};
OC = SOFAgetConventions(A.GLOBAL_SOFAConventions,'a');

%% Update dimensions
A=SOFAupdateDimensions(A);
B=SOFAupdateDimensions(B);

%% create field names which have to be checked
C=rmfield(B,{'API','Data'});
if isfield(C,'PRIVATE'), C=rmfield(C,'PRIVATE'); end
Bf=fieldnames(C);

%% Copy and merge metadata variables
C=A;
for ii=1:size(Bf,1)
	if ~isfield(A,Bf{ii})
		C.(Bf{ii}) = B.(Bf{ii});	% field in B but not in A. Simple copy.
	else	% here we have a potential conflict and have to merge
		if strfind(Bf{ii},'_')	% is it an attribute?
			if strcmp(A.(Bf{ii}),B.(Bf{ii})),
				C.(Bf{ii}) = B.(Bf{ii});	% content the same, no conflict
			else
				switch Bf{ii}
					case 'GLOBAL_TimeCreated'	% use the oldest date
						dateNew=datenum(A.GLOBAL_TimeCreated,Def.dateFormat);
						if datenum(B.GLOBAL_TimeCreated,Def.dateFormat)<dateNew, dateNew=datenum(B.GLOBAL_TimeCreated,Def.dateFormat); end;
						C.(Bf{ii}) = datestr(dateNew,Def.dateFormat);
						log{end+1}=[Bf{ii} ' set to ' C.(Bf{ii})];
					case 'GLOBAL_TimeModified' % now
						C.(Bf{ii}) = datestr(now,Def.dateFormat);
						log{end+1}=[Bf{ii} ' updated'];
					otherwise
						C.(Bf{ii}) = [A.(Bf{ii}) '; ' B.(Bf{ii})]; % concatate [A; B]
						log{end+1}=[Bf{ii} ' merged'];
				end
			end
		else	% a variable
			if isfield(OC.API.Dimensions, Bf{ii})	% is a known variable?
				AExp=SOFAexpand(A,Bf{ii});
				BExp=SOFAexpand(B,Bf{ii});
				dim=strfind(AExp.API.Dimensions.(Bf{ii}),'M');
				if isempty(dim),
					error([Bf{ii} ' can not be merged because it does not depend on M']);
				end
				C.(Bf{ii})=cat(dim,AExp.(Bf{ii}),BExp.(Bf{ii}));
				log{end+1}=[Bf{ii} ' expanded and merged'];
			else	% user-defined variable, dimensions must be stated
				if ~isfield(A.API.Dimensions, Bf{ii})
					error(['Dimension missing for ' Bf{ii} ' in A.']); end
				if ~isfield(B.API.Dimensions, Bf{ii})
					error(['Dimension missing for ' Bf{ii} ' in B.']); end
				dim=strfind(A.API.Dimensions.(Bf{ii}),'M');
				if ~isempty(dim),
          C.(Bf{ii})=cat(dim,A.(Bf{ii}),B.(Bf{ii})); % depends on M, merge
          log{end+1}=[Bf{ii} ' merged'];
        else  % not M-dependend, must be identical in A and B
          if prod(A.(Bf{ii})==B.(Bf{ii}))==1,
            C.(Bf{ii})=A.(Bf{ii});
            log{end+1}=[Bf{ii} ' identical'];
          else
            error([Bf{ii} ' must depend on M or be equal in both objects.']);
          end
        end
			end
		end
	end
end

%% Copy and merge Data variables
Bf=fieldnames(B.Data);
for ii=1:size(Bf,1)
	if ~isfield(A.Data,Bf{ii})
		C.Data.(Bf{ii}) = B.Data.(Bf{ii});	% field in B but not in A. Simple copy.
	else	% here we have a potential conflict and have to merge
		if strfind(Bf{ii},'_')	% is it an attribute?
			if strcmp(A.Data.(Bf{ii}),B.Data.(Bf{ii})),
				C.Data.(Bf{ii}) = B.Data.(Bf{ii});	% content the same, no conflict
			else
				C.Data.(Bf{ii}) = [A.Data.(Bf{ii}) '; ' B.Data.(Bf{ii})]; % concatate [A; B]
				log{end+1}=['Data.' Bf{ii} ' merged'];
			end
		else	% a variable in Data
			if isfield(OC.API.Dimensions.Data,Bf{ii})	% is a known variable?
				dim=strfind(A.API.Dimensions.Data.(Bf{ii}),'M'); % is a matrix
				if ~isempty(dim),
          C.Data.(Bf{ii})=cat(dim,A.Data.(Bf{ii}),B.Data.(Bf{ii})); % depends on M, cat
          log{end+1}=['Data.' Bf{ii} ' merged'];
        else  % not M-dependend, must be identical in A and B
          if all(A.Data.(Bf{ii})==B.Data.(Bf{ii}))==1,
            C.Data.(Bf{ii})=A.Data.(Bf{ii});
            log{end+1}=['Data.' Bf{ii} ' identical'];
          else
            error(['Data.' Bf{ii} ' must depend on M or be equal in both objects.']);
          end
        end
			else	% user-defined variable, dimensions must be stated
				if ~isfield(A.API.Dimensions.Data, Bf{ii})
					error(['Dimension missing for Data.' Bf{ii} ' in A.']); end
				if ~isfield(B.API.Dimensions.Data, Bf{ii})
					error(['Dimension missing for Data.' Bf{ii} ' in B.']); end
				dim=strfind(A.API.Dimensions.Data.(Bf{ii}),'M');
				if ~isempty(dim),
          C.Data.(Bf{ii})=cat(dim,A.Data.(Bf{ii}),B.Data.(Bf{ii})); % depends on M, cat
          log{end+1}=['Data.' Bf{ii} ' merged'];
        else  % not M-dependend, must be identical in A and B
          if prod(A.Data.(Bf{ii})==B.Data.(Bf{ii}))==1,
            C.Data.(Bf{ii})=A.Data.(Bf{ii});
            log{end+1}=['Data.' Bf{ii} ' identical'];
          else
            error(['Data.' Bf{ii} ' must depend on M or be equal in both objects.']);
          end
        end
			end
		end
	end
end

%% Update the new dimensions and finish
C=SOFAcompact(C);
C=SOFAupdateDimensions(C);
if length(log)>1, log=log(2:end); else log={}; end;

%% Get the sizes of the dimension variables according the dimension variables in str
function vec=getdim(Obj,str)
	vec=arrayfun(@(f)(Obj.API.(f)),upper(str));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function C = cache_test(A,B)
  % CACHE_TEST Dummy program that computes C = A+B, but demonstrates how to use
  % md5 caching of function results on input parameters
  %
  % C = cache_test(A,B)
  %

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Check for cached result, do NOT edit variables until cache is checked,
  % your function code comes later. See below
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get a list of current variables in this scope, this is the input "state"
  variables = who;
  % get a temporary file's name
  tmpf = [tempname('.') '.mat'];
  % save the "state" to file, so we can get a md5 checksum
  save(tmpf,'-regexp',sprintf('^%s$|',variables{:}),'-ascii');
  % get md5 checksum on input "state", we append .cache.mat to the check sum
  % because we'll use the checksum as the cache file name
  [s,cache_name] = system(['md5sum ' tmpf ' | awk ''{printf "."$1".cache.mat"}''']);
  % clean up
  delete(tmpf);
  clear s tmpf variables;

  % If the checksum cache file exists then we've seen this input "state"
  % before, load in cached output "state"
  if(exist(cache_name,'file'))
    fprintf('Cache exists. Using cache...\n');
    % use cache
    load(cache_name);
  % Otherwise this is the first time we've seen this input "state", so we
  % execute the function as usual and save the output "state" to the cache
  else
    fprintf('First time. Creating cache...\n');

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Your function code goes here
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    C = A+B;

    % get list of variables present in this scope at finish of function code,
    % this is the output "state"
    variables = who;
    % save output "state" to file, using md5 checksum cache file name
    save(cache_name,'-regexp',sprintf('^%s$|',variables{:}));
  end
end

author: alecjacobson

repo: gptoolbox