text:
function [ucm_sigmoid, ucm_raw, ws_wt2] = contours2ucm_RGBD(pb_oriented, pb_wt, thr)
% changes:
% - outputs raw ucm also
%
% Creates Ultrametric Contour Map from oriented contours
%
% syntax:
%   [ucm] = contours2ucm(pb_oriented, fmt)
%
% description:
%   Computes UCM by considering
%   the mean pb value on the boundary between regions as dissimilarity.
%
% arguments:
%   pb_oriented: Oriented Probability of Boundary
%   fmt:         Output format. 'imageSize' (default) or 'doubleSize'
%
% output:
%   ucm:    Ultrametric Contour Map in double
%
% Pablo Arbelaez <arbelaez@eecs.berkeley.edu>
% August 2013

% create finest partition and transfer contour strength

[ws_wt] = create_finest_partition(pb_oriented);

if nargin>1,

    % align old and new signal
    pb_wt(2:end,2:end)=pb_wt(1:end-1,1:end-1);
    pb_wt(1,:)=pb_wt(2,:); pb_wt(:,1)=pb_wt(:,2);
    pb_wt = pb_wt.*(ws_wt>0);

    ws_wt(ws_wt>thr)=thr;
    ws_wt=ws_wt+pb_wt;

end

% prepare pb for ucm
ws_wt2 = double(super_contour_4c(ws_wt));
ws_wt2 = clean_watersheds(ws_wt2);
labels2 = bwlabel(ws_wt2 == 0, 8);
labels = labels2(2:2:end, 2:2:end) - 1; % labels begin at 0 in mex file.
ws_wt2(end+1, :) = ws_wt2(end, :);
ws_wt2(:, end+1) = ws_wt2(:, end);

% compute ucm with mean pb.
super_ucm = double(ucm_mean_pb(ws_wt2, labels));

ucm_raw = super_ucm;

% normalize output
ucm_sigmoid = apply_sigmoid(ucm_raw,0.55,10).*(ucm_raw>0);

function [ws_wt] = create_finest_partition(pb_oriented)

pb = max(pb_oriented,[],3);
ws = watershed(pb);
ws_bw = (ws == 0);

contours = fit_contour(double(ws_bw));
angles = zeros(numel(contours.edge_x_coords), 1);

for e = 1 : numel(contours.edge_x_coords)
    if contours.is_completion(e), continue; end
    v1 = contours.vertices(contours.edges(e, 1), :);
    v2 = contours.vertices(contours.edges(e, 2), :);

    if v1(2) == v2(2),
        ang = 90;
    else
        ang = atan((v1(1)-v2(1)) / (v1(2)-v2(2)));
    end
    angles(e) = ang*180/pi;
end

orient = zeros(numel(contours.edge_x_coords), 1);
orient((angles<-78.75) | (angles>=78.75)) = 1;
orient((angles<78.75) & (angles>=56.25)) = 2;
orient((angles<56.25) & (angles>=33.75)) = 3;
orient((angles<33.75) & (angles>=11.25)) = 4;
orient((angles<11.25) & (angles>=-11.25)) =5;
orient((angles<-11.25) & (angles>=-33.75)) = 6;
orient((angles<-33.75) & (angles>=-56.25)) = 7;
orient((angles<-56.25) & (angles>=-78.75)) = 8;

ws_wt = zeros(size(ws_bw));
for e = 1 : numel(contours.edge_x_coords)
    if contours.is_completion(e), continue; end
    for p = 1 : numel(contours.edge_x_coords{e}),
        ws_wt(contours.edge_x_coords{e}(p), contours.edge_y_coords{e}(p)) = ...
            max(pb_oriented(contours.edge_x_coords{e}(p), contours.edge_y_coords{e}(p), orient(e)), ws_wt(contours.edge_x_coords{e}(p), contours.edge_y_coords{e}(p)));
    end
    v1=contours.vertices(contours.edges(e,1),:);
    v2=contours.vertices(contours.edges(e,2),:);
    ws_wt(v1(1),v1(2))=max( pb_oriented(v1(1),v1(2), orient(e)),ws_wt(v1(1),v1(2)));
    ws_wt(v2(1),v2(2))=max( pb_oriented(v2(1),v2(2), orient(e)),ws_wt(v2(1),v2(2)));
end
ws_wt=double(ws_wt);
%%

function [pb2, V, H] = super_contour_4c(pb)

V = min(pb(1:end-1,:), pb(2:end,:));
H = min(pb(:,1:end-1), pb(:,2:end));

[tx, ty] = size(pb);
pb2 = zeros(2*tx, 2*ty);
pb2(1:2:end, 1:2:end) = pb;
pb2(1:2:end, 2:2:end-2) = H;
pb2(2:2:end-2, 1:2:end) = V;
pb2(end,:) = pb2(end-1, :);
pb2(:,end) = max(pb2(:,end), pb2(:,end-1));

%%

function [ws_clean] = clean_watersheds(ws)
% remove artifacts created by non-thin watersheds (2x2 blocks) that produce
% isolated pixels in super_contour

ws_clean = ws;

c = bwmorph(ws_clean == 0, 'clean', inf);

artifacts = ( c==0 & ws_clean==0 );
R = regionprops(bwlabel(artifacts), 'PixelList');

for r = 1 : numel(R),
    xc = R(r).PixelList(1,2);
    yc = R(r).PixelList(1,1);

    vec = [ max(ws_clean(xc-2, yc-1), ws_clean(xc-1, yc-2)) ...
        max(ws_clean(xc+2, yc-1), ws_clean(xc+1, yc-2)) ...
        max(ws_clean(xc+2, yc+1), ws_clean(xc+1, yc+2)) ...
        max(ws_clean(xc-2, yc+1), ws_clean(xc-1, yc+2)) ];

    [nd,id] = min(vec);
    switch id,
        case 1,
            if ws_clean(xc-2, yc-1) < ws_clean(xc-1, yc-2),
                ws_clean(xc, yc-1) = 0;
                ws_clean(xc-1, yc) = vec(1);
            else
                ws_clean(xc, yc-1) = vec(1);
                ws_clean(xc-1, yc) = 0;

            end
            ws_clean(xc-1, yc-1) = vec(1);
        case 2,
            if ws_clean(xc+2, yc-1) < ws_clean(xc+1, yc-2),
                ws_clean(xc, yc-1) = 0;
                ws_clean(xc+1, yc) = vec(2);
            else
                ws_clean(xc, yc-1) = vec(2);
                ws_clean(xc+1, yc) = 0;
            end
            ws_clean(xc+1, yc-1) = vec(2);

        case 3,
            if ws_clean(xc+2, yc+1) < ws_clean(xc+1, yc+2),
                ws_clean(xc, yc+1) = 0;
                ws_clean(xc+1, yc) = vec(3);
            else
                ws_clean(xc, yc+1) = vec(3);
                ws_clean(xc+1, yc) = 0;
            end
            ws_clean(xc+1, yc+1) = vec(3);
        case 4,
            if ws_clean(xc-2, yc+1) < ws_clean(xc-1, yc+2),
                ws_clean(xc, yc+1) = 0;
                ws_clean(xc-1, yc) = vec(4);
            else
                ws_clean(xc, yc+1) = vec(4);
                ws_clean(xc-1, yc) = 0;
            end
            ws_clean(xc-1, yc+1) = vec(4);
    end
end

%%
function cmap_proba = apply_sigmoid(cmap, thr, fq)
cmap_proba=cmap;
cmap_proba(:) = 1./(1+exp(-fq*(cmap(:)-thr)));

author: s-gupta

repo: rgbd

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


text:
% test code for tt_qlaplace_dd()
%
% September 22, 2010
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
% Look for details in the Preprint No. 75, 2010 of
% Max-Planck Institute for Mathematics in the Sciences
% Vladimir A. Kazeev and Boris N. Khoromskij
% On explicit QTT representation of Laplace operator and its inverse
% http://www.mis.mpg.de/publications/preprints/2010/prepr2010-75.html

% d is the only parameter
d=[3,4,5];
%

D=size(d,2);
tt=tt_qlaplace_dd(d);
tt=tt_mat_to_vec(tt);

full=nd_to_full(tt);

Z=zeros(2^sum(d));
for k=1 : D
	L=2*eye(2^d(k));
	for i=1 : 2^d(k)-1
		L(i,i+1)=-1;
		L(i+1,i)=-1;
	end
	for m=k-1 : -1: 1
		L=kron(eye(2^d(m)),L);
	end
	for m=k+1 : D
		L=kron(L,eye(2^d(m)));
	end
	Z=Z+L;
end
err=norm(full-Z,'fro');
fprintf('fro err = %e\n', err);

author: oseledets

repo: TT-Toolbox


text:
function [amp_map, ang_map] = detectChange(obj, t1, t2)
% This method returns an amplitude map and a angle map as a byproduct.
% When called in this form:
%   change_map = detectChange(obj, t1, t2)
% Only the amp_map will be caught and copied to change_map, while the
% ang_map deserted.
%
% Unlike many counterparts, this implementation perfoms (t2 - t1)
% insead of (t1 - t2) as I think the former one easier for
% interpretation in vectors
diffMap = double(t2) - double(t1);
amp_map = sqrt(sum(diffMap.^2, 3));
ang_map = (diffMap ./ amp_map);
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function rad = deg2rad(deg)
%DEG2RAD Convert angle from degrees to radians
%
%   Usage:
%   R = deg2rad(D)
%   convert an angle in degrees to an angle in radians.
%
%   Example
%   deg2rad(180)    % gives pi
%   ans =
%       3.1416
%   deg2rad(60)     % gives pi/3
%   ans =
%       1.0472
%
%   See Also
%   angles2d, rad2deg
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 09/12/2004.
%

rad = deg*pi/180;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function pagelayo
% PAGELAYO Set up graphical objects for Tutor environment.
%          PAGELAYO sets up all the axes, line, and text graphical objects
%          depending upon the desired arrangement.

% Author: Craig Borghesani
% Date: 10/17/94
% Revised: 11/13/94
% Copyright (c) 1999, Prentice-Hall

page_chil = [];

f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};

page_setup = get(ui_han(34),'userdata');
page_layout = get(ui_han(35),'userdata');
if length(page_layout),
 delete(page_layout(find(page_layout~=0)));
 page_chil = get(ui_han(90),'userdata');
 set(page_chil,'vis','off');
end
page_layout = zeros(4,8);

% locations of axes
scrn_size = get(0,'screensize');
fig_size = get(gcf,'position');
%fig_wdt = scrn_size(3)-10*2;
fig_wdt = fig_size(3);
frm_wdt = 190;
ers = 'xor';

%axs1_l = (frm_wdt+2+60)/fig_wdt;
axs1_l = 0.07;
axs1_b = 0.15;
%axs1_w = 0.97 - axs1_l;
axs1_w = 0.97 - (frm_wdt+2+60)/fig_wdt;
axs1_h = 0.68;

axs2_b = 0.58;
axs2_w = axs1_w/2 - 0.05;
axs2_l = axs1_l + axs2_w + 0.1;
%axs2_h = axs2_w;
axs2_h = 0.26;

if ~length(page_chil),
 for k = 1:8,
  page_chil(k) = uimenu(ui_han(90),'vis','off');
 end
 set(ui_han(90),'userdata',page_chil);
end

page_ct = 0;
for page_num = 1:8,
 lab_tit = [];
 plots = page_setup(:,page_num); plots(find(plots==0)) = [];

 if length(plots)==1,
  axs_pos = [axs1_l,axs1_b,axs1_w,axs1_h];
 elseif length(plots)==2,
  axs_pos = [axs1_l,axs2_b,axs1_w,axs2_h;
             axs1_l,axs1_b,axs1_w,axs2_h];
 elseif length(plots)==3,
  axs_pos = [axs1_l,axs2_b,axs1_w,axs2_h;
             axs1_l,axs1_b,axs2_w,axs2_h;
             axs2_l,axs1_b,axs2_w,axs2_h];
 elseif length(plots)==4,
  axs_pos = [axs1_l,axs2_b,axs2_w,axs2_h;
             axs2_l,axs2_b,axs2_w,axs2_h;
             axs1_l,axs1_b,axs2_w,axs2_h;
             axs2_l,axs1_b,axs2_w,axs2_h];
 else
  break;
 end

 for k = 1:length(plots),
  if any(plots(k)==[2,3,6,7]),
   xscale = 'log';
  else
   xscale = 'linear';
  end
  if plots(k) == 1, % Nyquist
   xlab = 'Real';
   ylab = 'Imaginary';
   titl = 'Nyquist';
   xlim = [-1,1]; ylim = [-1,1];
  elseif plots(k) == 2, % Bode (Magnitude)
   xlab = 'Frequency (rad/sec)';
   ylab = 'Magnitude (dB)';
   titl = 'Bode (Magnitude)';
   xlim = [0.01,1000]; ylim = [-20,20];
  elseif plots(k) == 3, % Bode (Phase)
   xlab = 'Frequency (rad/sec)';
   ylab = 'Phase (degrees)';
   titl = 'Bode (Phase)';
   xlim = [0.01,1000]; ylim = [-360,0];
  elseif plots(k) == 4, % Nichols
   xlab = 'Phase (degrees)';
   ylab = 'Magnitude (dB)';
   titl = 'Nichols';
   xlim = [-360,0]; ylim = [-20,20];
  elseif plots(k) == 5, % Root Locus
   xlab = 'Real';
   ylab = 'Imaginary';
   titl = 'Root Locus';
   xlim = [-1,1]; ylim = [-1,1];
  elseif plots(k) == 6, % Gain (Magnitude)
   xlab = 'K';
   ylab = 'Magnitude';
   titl = 'Gain (Magnitude)';
   xlim = [0.01,1000]; ylim = [0.01,10];
  elseif plots(k) == 7, % Gain (Phase)
   xlab = 'K';
   ylab = 'Phase (degrees)';
   titl = 'Gain (Phase)';
   xlim = [0.01,1000]; ylim = [0,360];
  elseif plots(k) == 8, % Time Response
   xlab = 'Time (seconds)';
   ylab = 'Amplitude';
   titl = 'Time Response';
   xlim = [0,10]; ylim = [0,5];
  end

  axs = axes('vis','off','box','on','pos',axs_pos(k,:),...
             'xscale',xscale,'xgrid','on','ygrid','on',...
             'xlim',xlim,'ylim',ylim,...
             'drawmode','fast',...
             'tag',['1',int2str(plots(k)),'20']);
%             'buttondownfcn',['pagemous(1,',int2str(plots(k)),',10)']);
  XLabelHandle = get(axs,'xlabel');
  YLabelHandle = get(axs,'ylabel');
  TitleHandle  = get(axs,'title');
  set(XLabelHandle,'string',xlab);
  set(YLabelHandle,'string',ylab);
  set(TitleHandle,'string',titl);

  if any(plots(k)==[1,5]),
%   set(axs,'dataaspectratio',[1,1,1]);
   axis equal;
  end

  if any(plots(k)==[1,2,3,4]), % frequency response plots
   clr = 'g';
   lin = [];
   for k2 = 1:15,
    if k2 > 6, clr = 'r'; end
%    if k2 > 3, ers = 'norm'; end
    lin(k2) = line('xdata',0,'ydata',0,'color',clr,'erase',ers,...
                   'vis','off',...
                   'tag',['2',int2str(plots(k)),int2str(k2)]);
%     'buttondownfcn',['pagemous(2,',int2str(plots(k)),',',int2str(k2),')']);
   end
   set(axs,'userdata',[plots(k),lin]);

  elseif any(plots(k) == [5,6,7]), % root locus, gain plots

   set(axs,'userdata',plots(k),'nextplot','add');

   if plots(k) == 6,
    set(axs,'yscale','log');
   end

  else % time response
   clr = ['r';'g';'b';'y';'m';'c'];
   lin = [];
   for k2 = 1:6,
    lin(k2) = line('xdata',0,'ydata',0,'color',clr(k2),...
                   'vis','off',...
                   'erase',ers,...
                   'tag',['1',int2str(plots(k)),int2str(20+k2)]);
%     'buttondownfcn',['pagemous(1,',int2str(plots(k)),',10)']);
   end
   txt = [];
   for k2 = 1:12,
    txt(k2) = text('pos',[0,0,0],'vis','off',...
     'tag',['1',int2str(plots(k)),int2str(26+k2)],'color','k',...
     'erase',ers);
%     'buttondownfcn',['pagemous(1,',int2str(plots(k)),',10)']);
   end
   set(axs,'userdata',[plots(k),lin,txt]);

  end

  page_layout(k,page_num)=axs;
  lab_tit = [lab_tit,titl,', '];

 end

 if length(lab_tit),
  lab_tit((length(lab_tit)-1):length(lab_tit))=[];
  set(page_chil(page_num),'label',lab_tit,'vis','on',...
                   'callback',['pagesele(',int2str(page_num),',14)']);
 end
 page_ct = page_ct + 1;
end

if page_ct < 2,
 set(ui_han(90),'enable','off');
else
 set(ui_han(90),'enable','on');
end

set(ui_han(35),'userdata',page_layout);

author: Sable

repo: mcbench-benchmarks


text:
function roomhyp = sample_roomtype4(n, lines, vp, imgwidth, imgheight)

% 4 ways to create hypotheses
% roomtype2 + 3,top,right
% roomtype2 + 2,top
% roomtype7 + 3,bot,right
% roomtype7 + 2,bot

lc = [lines.lineclass]; % 0,1,2,3
lr = [lines.leftorright]; % 1 if right, -1 if left, 0 if neither
tb = [lines.above_horizon]; % 1 if above, -1 if below, 0 if neither

twotop = find((lc==2) & (tb==1));
twobot = find((lc==2) & (tb==-1));
threetopright = find((lc==3) & (tb==1) & (lr==1));
threebotright = find((lc==3) & (tb==-1) & (lr==1));

linecode{1} = twotop;
linecode{2} = twobot;
linecode{3} = threetopright;
linecode{4} = threebotright;

% recipe: [roomtypepart linecode oldcornerid newcornerid]
recipe(1,:) = [2 3 4 2];
recipe(2,:) = [2 1 4 2];
recipe(3,:) = [7 4 2 4];
recipe(4,:) = [7 2 2 4];

% roomhyp = [];
count = 0;
num_continue_without_progress = 0;
while count < n && num_continue_without_progress < 30
% pick recipe
    cur_recipe = recipe(randsample(size(recipe,1), 1), :);

    % sample partial room
    if cur_recipe(1) == 2
        roomhyppart = sample_roomtype2(1, lines, vp, imgwidth, imgheight);
    elseif cur_recipe(1) == 7
        roomhyppart = sample_roomtype7(1, lines, vp, imgwidth, imgheight);
    end

    if isempty(roomhyppart)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % sample lines
    l1 = linecode{cur_recipe(2)};
    if length(l1)<1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end
    ls1 = l1(randsample(length(l1), 1));

    % new corner
    oldpoint = roomhyppart.corner(cur_recipe(3)).pt;
    [newpoint degen] = line_intersect(vp{1},oldpoint,...
        lines(ls1).point1, lines(ls1).point2);
    if degen==1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end
    % check if inside img
    MARGIN = 5;
    if ~is_in_image(newpoint, imgwidth, imgheight, MARGIN)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % TODO: check location of new line sample
    %

    % add to hypothesis
    count = count + 1;
    num_continue_without_progress = 0;
    roomhyp(count) = roomhyppart;
    roomhyp(count).corner(cur_recipe(4)).pt = newpoint;
    roomhyp(count).type = 4;

%     global img;
%     disp_vanish(img, lines(ls1), vp);
%     plot(oldpoint(1), oldpoint(2), 'rx', 'MarkerSize',10, 'LineWidth',2);
%     plot(newpoint(1), newpoint(2), 'bx', 'MarkerSize',10, 'LineWidth',2);
%     pause;
%     close;
end

if ~exist('roomhyp','var')
    roomhyp = [];
end

author: zouchuhang

repo: LayoutNet


text:
function [ x, r, z, p, q, job ] = cg_rc ( n, b, x, r, z, p, q, job )

%*****************************************************************************80
%
%% CG_RC is a reverse communication conjugate gradient routine.
%
%  Discussion:
%
%    This routine seeks a solution of the linear system A*x=b
%    where b is a given right hand side vector, A is an n by n
%    symmetric positive definite matrix, and x is an unknown vector
%    to be determined.
%
%    Under the assumptions that the matrix A is large and sparse,
%    the conjugate gradient method may provide a solution when
%    a direct approach would be impractical because of excessive
%    requirements of storage or even of time.
%
%    The conjugate gradient method presented here does not require the
%    user to store the matrix A in a particular way.  Instead, it only
%    supposes that the user has a way of calculating
%      y = alpha * A * x + b * y
%    and of solving the preconditioned linear system
%      M * x = b
%    where M is some preconditioning matrix, which might be merely
%    the identity matrix, or a diagonal matrix containing the
%    diagonal entries of A.
%
%    This routine was extracted from the "templates" package.
%    There, it was not intended for direct access by a user;
%    instead, a higher routine called "cg()" was called once by
%    the user.  The cg() routine then made repeated calls to
%    cgrevcom() before returning the result to the user.
%
%    The reverse communication feature of cgrevcom() makes it, by itself,
%    a very powerful function.  It allows the user to handle issues of
%    storage and implementation that would otherwise have to be
%    mediated in a fixed way by the function argument list.  Therefore,
%    this version of cgrecom() has been extracted from the templates
%    library and documented as a stand-alone procedure.
%
%    The user sets the value of JOB to 1 before the first call,
%    indicating the beginning of the computation, and to the value of
%    2 thereafter, indicating a continuation call.
%    The output value of JOB is set by cgrevcom(), which
%    will return with an output value of JOB that requests a particular
%    new action from the user.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%    Charles Romine, Henk van der Vorst,
%    Templates for the Solution of Linear Systems:
%    Building Blocks for Iterative Methods,
%    SIAM, 1994,
%    ISBN: 0898714710,
%    LC: QA297.8.T45.
%
%  Parameters:
%
%    Input, integer N, the dimension of the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, real X(N).  On first call, the user
%    should store an initial guess for the solution in X.
%
%    Input, real R(N), Z(N), P(N), Q(N), work arrays.  The user should
%    create each of these before the first call, using the zeros() command.
%    On subsequent calls, the user may be asked to assign a value to one
%    of these vectors.
%
%    Input, integer JOB, communicates the task to be done.
%    The user needs to set the input value of JOB to 1, before the first call,
%    and then to 2 for every subsequent call for the given problem.
%
%    Output, real X(N), the current solution estimate.
%    Each time JOB is returned as 4, X has been updated.
%
%    Output, real R(N), Z(N), P(N), Q(N), work arrays.  Depending on the
%    output value of JOB, the user may be asked to carry out a computation
%    involving some of these vectors.
%
%    Output, integer JOB, communicates the task to be done.
%    * JOB = 1, compute Q = A * P;
%    * JOB = 2: solve M*Z=R, where M is the preconditioning matrix;
%    * JOB = 3: compute R = R - A * X;
%    * JOB = 4: check the residual R for convergence.
%               If satisfactory, terminate the iteration.
%               If too many iterations were taken, terminate the iteration.
%

%
%  Local variables should be SAVED between calls.
%
  persistent iter
  persistent rho
  persistent rho_old
  persistent rlbl
%
%  Initialization.
%  Ask the user to compute the initial residual.
%
  if ( job == 1 )

    r(1:n) = b(1:n);

    job = 3;
    rlbl = 2;
%
%  Begin first conjugate gradient loop.
%  Ask the user for a preconditioner solve.
%
  elseif ( rlbl == 2 )

    iter = 1;

    job = 2;
    rlbl = 3;
%
%  Compute the direction.
%  Ask the user to compute ALPHA.
%  Save A*P to Q.
%
  elseif ( rlbl == 3 )

    rho = r(:)' * z(:);

    if ( 1 < iter )
      beta = rho / rho_old;
      z(1:n) = z(1:n) + beta * p(1:n);
    end

    p(1:n) = z(1:n);

    job = 1;
    rlbl = 4;
%
%  Compute current solution vector.
%  Ask the user to check the stopping criterion.
%
  elseif ( rlbl == 4 )

    pdotq = p(:)' * q(:);
    alpha = rho / pdotq;
    x(1:n) = x(1:n) + alpha * p(1:n);
    r(1:n) = r(1:n) - alpha * q(1:n);

    job = 4;
    rlbl = 5;
%
%  Begin the next step.
%  Ask for a preconditioner solve.
%
  elseif ( rlbl == 5 )

    rho_old = rho;
    iter = iter + 1;

    job = 2;
    rlbl = 3;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, f, s, p, cdf ] = negative_binomial_cdf_values ( n_data )

%*****************************************************************************80
%
%% NEGATIVE_BINOMIAL_CDF_VALUES returns values of the negative binomial CDF.
%
%  Discussion:
%
%    Assume that a coin has a probability P of coming up heads on
%    any one trial.  Suppose that we plan to flip the coin until we
%    achieve a total of S heads.  If we let F represent the number of
%    tails that occur in this process, then the value of F satisfies
%    a negative binomial PDF:
%
%      PDF(F,S,P) = Choose ( F from F+S-1 ) * P**S * (1-P)**F
%
%    The negative binomial CDF is the probability that there are F or
%    fewer failures upon the attainment of the S-th success.  Thus,
%
%      CDF(F,S,P) = sum ( 0 <= G <= F ) PDF(G,S,P)
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`DiscreteDistributions`]
%      dist = NegativeBinomialDistribution [ s, p ]
%      CDF [ dist, f ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    F C Powell,
%    Statistical Tables for Sociology, Biology and Physical Sciences,
%    Cambridge University Press, 1982.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer F, the maximum number of failures.
%
%    Output, integer S, the number of successes.
%
%    Output, real P, the probability of a success on one trial.
%
%    Output, real CDF, the probability of at most F failures
%    before the S-th success.
%
  n_max = 27;

  cdf_vec = [ ...
     0.6367187500000000E+00, ...
     0.3632812500000000E+00, ...
     0.1445312500000000E+00, ...
     0.5000000000000000E+00, ...
     0.2265625000000000E+00, ...
     0.6250000000000000E-01, ...
     0.3437500000000000E+00, ...
     0.1093750000000000E+00, ...
     0.1562500000000000E-01, ...
     0.1792000000000000E+00, ...
     0.4096000000000000E-01, ...
     0.4096000000000000E-02, ...
     0.7047000000000000E-01, ...
     0.1093500000000000E-01, ...
     0.7290000000000000E-03, ...
     0.9861587127990000E+00, ...
     0.9149749500510000E+00, ...
     0.7471846521450000E+00, ...
     0.8499053647030009E+00, ...
     0.5497160941090026E+00, ...
     0.2662040052146710E+00, ...
     0.6513215599000000E+00, ...
     0.2639010709000000E+00, ...
     0.7019082640000000E-01, ...
     0.1000000000000000E+01, ...
     0.1990000000000000E-01, ...
     0.1000000000000000E-03 ];

  f_vec = [ ...
     4,  3,  2, ...
     3,  2,  1, ...
     2,  1,  0, ...
     2,  1,  0, ...
     2,  1,  0, ...
    11, 10,  9, ...
    17, 16, 15, ...
     9,  8,  7, ...
     2,  1,  0 ];

  p_vec = [ ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.40E+00, ...
     0.40E+00, ...
     0.40E+00, ...
     0.30E+00, ...
     0.30E+00, ...
     0.30E+00, ...
     0.30E+00, ...
     0.30E+00, ...
     0.30E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E-01, ...
     0.10E-01, ...
     0.10E-01 ];

  s_vec = [ ...
    4, 5, 6, ...
    4, 5, 6, ...
    4, 5, 6, ...
    4, 5, 6, ...
    4, 5, 6, ...
    1, 2, 3, ...
    1, 2, 3, ...
    1, 2, 3, ...
    0, 1, 2 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    f = 0;
    s = 0;
    p = 0.0;
    cdf = 0.0;
  else
    f = f_vec(n_data);
    s = s_vec(n_data);
    p = p_vec(n_data);
    cdf = cdf_vec(n_data);
  end

  return
end


text:
function drawIdpLin(MapFig,Lmk,color,MapOpt)

% DRAWIDPLIN  Draw inverse-depth line landmark in MapFig.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

% get the lmk from the Map
r   = Lmk.state.r;       % range in Map
idl = Map.x(r);          % mean
IDL = Map.P(r,r);        % covariances matrix
t   = [Lmk.par.endp.t]'; % abscissas of endpoints, t = [t1;t2]

% extract two endpoints - means and covariance
[e1,e2,E1_idl,E2_idl] = idpLinEndpoints(idl,t(1),t(2)); % means and Jacobians
E1 = E1_idl*IDL*E1_idl'; % covariances
E2 = E2_idl*IDL*E2_idl';

% draw the mean:
drawSeg(MapFig.Lmk(Lmk.lmk).mean,[e1;e2],color.mean)

% draw the covariance ellipses
if MapOpt.showEllip
    drawEllipse(MapFig.Lmk(Lmk.lmk).ellipse(1), e1, E1, color.ellip)
    drawEllipse(MapFig.Lmk(Lmk.lmk).ellipse(2), e2, E2, color.ellip)
end

% draw the label
if MapOpt.showLmkId
    e = e2-e1;
    n = null([e e e]);
    n = n(:,2)*sign(n(3,2)); % Inverse depth line's normal vector
    posOffset = 0.2*n;     % label orthogonally out of the line.
    drawLabel(MapFig.Lmk(Lmk.lmk).label,0.5*(e1+e2) + posOffset,num2str(Lmk.id))
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function varargout=estnoisem(varargin)
%
% For calling details please see v_estnoisem.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_estnoisem.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: estnoisem.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_estnoisem(varargin{:});
else
    v_estnoisem(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
classdef TestAccumulateProduct
    %TestAccumulateProduct

    methods (Static)
        function test_1
            sz = [10 20];
            dst = zeros(sz);
            for i=1:5
                dst = cv.accumulateProduct(rand(sz), rand(sz), dst);
            end
            validateattributes(dst, {'double'}, {'size',sz});
        end

        function test_2
            sz = [10 20];
            dst = zeros(sz, 'single');
            for i=1:5
                dst = cv.accumulateProduct(randi([0 255],sz,'uint8'), ...
                    randi([0 255],sz,'uint8'), dst);
            end
            validateattributes(dst, {'single'}, {'size',sz});
        end

        function test_3
            sz = [10 20];
            dst = zeros(sz);
            for i=1:5
                src = rand(sz);
                mask = (rand(sz) > 0.5);
                dst = cv.accumulateProduct(src, src, dst, 'Mask',mask);
            end
        end

        function test_error_argnum
            try
                cv.accumulateProduct();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function z = cumint3(x,y)
%
%   z = cumint3(x,y)
%
% This returns a vector z the same size as x and y
% which is the cumulative integral of y with respect
% to x, with the lower integration limit set to x(1) and
% the upper limit ranging from x(1) to x(n). The
% successive intervals in x need not be of equal lengths,
% though none should be of zero length. A third order
% approximation is made so that for up to cubic polynomials
% the values will be exact except for rounding.
% x and y must be column vectors of the same length
% and have at least four elements. Note that with only
% z = b.*g below, this would be computing matlab's
% 'cumtrapz' trapezoidal integration, the rest of the
% expression in z serving to carry out the additional
% third order approximation.  RAS - 03/11/08

% Test arguments
[n,m] = size(x);
if any([n,m]~=size(y)) | m~=1
 error('x and y must be column vectors of equal length.')
end
if n<4
 error('There must be at least four points.')
end

% Compute third order integral
xe = [x(4);x;x(n-3)]; ye = [y(4);y;y(n-3)]; % Provide for endpoints
x0 = xe(1:n-1); x1 = xe(2:n); x2 = xe(3:n+1); x3 = xe(4:n+2);
y0 = ye(1:n-1); y1 = ye(2:n); y2 = ye(3:n+1); y3 = ye(4:n+2);
a = x1-x0; b = x2-x1; c = x3-x2;
d = y1-y0; e = y2-y1; f = y3-y2;
g = (y1+y2)/2;

% Each z value will be the integral from x1 to x2 of a cubic
% polynomial running through (x0,y0),(x1,y1),(x2,y2),(x3,y3).
z = b.*g+1/12*b.^2.*(+c.*b.*(2*c+b).*(c+b).*d ...
    -a.*c.*(c-a).*(2*c+2*a+3*b).*e-a.*b.*(2*a+b).*(a+b).*f) ...
    ./(a.*c.*(a+b).*(c+b).*(c+a+b));

% Obtain cumulative integral values
z = [0;cumsum(z)];

return


text:
%kkreplace 'Replace All Occurrences of X with Y'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kreplace.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input file'
% OutputFile: o ' Output ', required: 'Output file'
%
% Example: o = kkreplace(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kreplace - Replace All Occurrences of X with Y
%
%  DESCRIPTION
%
% The \fBReplace Value" operator, kreplace, replaces each occurrence of
% the specified old value (X) with a new value (Y).
% The operations occur on the map segment if
% one is present in the source object otherwise the operations occur on the
% value segment.
%
% The optional tolerance values allow the value that is being
% replaced to fall within a given range (X - Lower tolerance to X +
% Upper tolerance).
%
%
%  "Data Type" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkreplace(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkreplace(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kreplace"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function kappa = evalkappa(lambda,varargin)
% eigenvalues of orientation tensor to bingham distribution parameters
%
% Options
%  approximated - approximated solution of kappas
%  precision    - precision of solvus
%  iteration    - number of newton interations
%
% See also
% orientation/mean EBSD/mean BinghamODF

% only approximated solution of kappas' which hold true only for large
% values of lambda see: F. Bachmann, R. Hielscher, P. Jupp, W. Pantleon, H.
% Schaeben, E.Wegert: Inferential statistics of electron backscatter
% diffraction data from within individual crystalline grains, Journal of
% applied Crystallogrpahy, 2010, 34, eq. 38.

lambda = sort(lambda);
kappa = -1./(2.*lambda);
kappa(4) = 0;
kappa = kappa - kappa(1);

% kappa = sort(1./((2.*lambda-1).*lambda));
% kappa = (kappa-min(kappa))/2;

end

author: mtex-toolbox

repo: mtex

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


text:
function run_mimic_imp_annotation
% run_mimic_imp_annotation is designed to annotate the quality and breaths
% of signals in the impedance pneumography signal quality index (SQI) dataset
% (extracted from the MIMIC III database). It was used to evaluate the
% impedance pneumography SQI (using only the breath annotation functionality).
%
%               run_mimic_imp_annotation
%
%	Inputs:
%       mimic_imp_sqi_data.mat - a file containing curated data from the MIMIC database.
%        - Specify the path to this file in the "setup_universal_params" function below.
%
%	Outputs:
%       Matlab files containing the annotations for each subject, called
%          "ImP_SQI#_##_an.mat" (where #_## is the subject number, followed
%          by the initials of the annotator)
%
%   Further Information:
%       This version of the run_mimic_imp_annotation is provided to facilitate
%       replication of the analysis reported in:
%           Charlton P. H. et al., "An impedance pneumography signal
%           quality index for respiratory rate monitoring: design,
%           assessment and application", [under review]
%       Further information on this study can be obtained at:
%           http://peterhcharlton.github.io/RRest/imp_sqi.html
%       In addition, further information on RRest, a toolbox of respiratory
%       algorithms which can be used with this dataset, can be obtained at:
%           http://peterhcharlton.github.io/RRest/index.html
%
%   Comments, Questions, Criticisms, Feedback, Contributions:
%       See: http://peterhcharlton.github.io/RRest/contributions.html
%
%   Version:
%       v.0.1 - accompanying peer review, 6th Aug 2020 by Peter H Charlton
%
%   Source:
%       This script has been adapted from 'MIMICII_data_importer.m', which
%       is one of the scripts in the RRest toolbox, available at:
%           https://github.com/peterhcharlton/RRest
%
%   Licence:
%       This program is available under the GNU public license. This
%       program is free software: you can redistribute it and/or modify
%       it under the terms of the GNU General Public License as published by
%       the Free Software Foundation, either version 3 of the License, or
%       (at your option) any later version. This program is distributed in
%       the hope that it will be useful, but WITHOUT ANY WARRANTY; without
%       even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%       PARTICULAR PURPOSE.  See the GNU General Public License for more
%       details: <http://www.gnu.org/licenses/>.
%

%% Setup
% ~~~ This function needs adjusting by the individual user as it contains paths
% specific to their computer ~~~
up = setup_universal_params;

%% Initialisation

% Load data
load(up.paths.data_file)
subjects = 1:length(data);   % create a list of subjects which are to be analysed.

% mode
up.mode = 2;        % select 1 if annotating high or low quality, 2 if annotating breaths
if up.mode == 1, fprintf('\n *** Mode selected: quality annotation *** '), elseif up.mode == 2, fprintf('\n *** Mode selected: Breath annotation *** '), end

fprintf(['\n\n\n\n\n~~~~~~~  Starting annotations  ~~~~~~~\n' ...
    'Instructions: Use the left mouse key to\n' ...
    'annotate a peak, and use the right mouse\n' ...
    'key to annotate a trough. To delete the\n' ...
    'most recent annotation, hold down SHIFT\n' ...
    'and press a mouse key.\n' ...
    'Please note that any previous annotations\n' ...
    'will be kept as long as the results files\n' ...
    'are in the correct location.\n\n']);

operator=input('(To enter review mode, type  review  )\n\nOtherwise, please type your first and last initial,\nand press enter:\n','s');
%operator = 'PC';
if strcmp(operator, 'review') == 1
    paw_verification_script_checker(subjects, data, up);
    return
end

question=input('\nWould you like to start from the \nfirst subject? (y or n, followed by enter)\n', 's');
%question = 'y';
if strcmp(question, 'n') == 1
    start_subject = str2num(input('\nWhich subject?\n', 's'));
else
    start_subject = 1;
end

%% Cycle through different subjects and periods
start_subject_el = find(subjects == start_subject);
% load data on which windows were high quality if annotating breaths
if up.mode == 2
    load(up.paths.data_file2)
    rel_col = find(strcmp(data_mat.all.response_header, 'subj'));
    data_mat.all.subj = data_mat.all.response_data(:,rel_col);
    rel_col = find(strcmp(data_mat.all.response_header, 'win_no'));
    data_mat.all.win_no = data_mat.all.response_data(:,rel_col);
end
if isempty(start_subject_el)
    fprintf('\nThis subject doesn''t exist. Try all over again.');
else
    subjects = subjects(start_subject_el:end);
    for s = subjects(1:end)

        %% Identify wins deemed to be high qual
        if up.mode == 2 % i.e. identifying RRs
            rel_col1 = find(strcmp(data_mat.all.response_header, 'sqi_novel'));
            rel_els1 = data_mat.all.subj==s;
            rel_els2 = data_mat.all.response_data(:, rel_col1) == 1;
            high_qual_wins = data_mat.all.win_no(rel_els1 & rel_els2);
            if length(high_qual_wins) > 5
                high_qual_wins = high_qual_wins(1:5);
            end
        end

        %% Load Impedance
        imp.fs = data(s).ref.resp_sig.imp.fs;
        imp.v = data(s).ref.resp_sig.imp.v;
        imp.t = [0:(length(imp.v)-1)]/imp.fs;

        %% Process Impedance
        imp = lpf_to_exclude_resp(imp, up);
        imp.v = -1*imp.v;

        %% Set up file for keeping annotations in:
        % see if this subject has already been annotated:
        if up.mode == 1
            qual.t = []; qual.v = [];
            save(up.paths.annotations_file, 'qual', 'up')
        else
            pk_anns.t=[]; pk_anns.v=[];
            save(up.paths.annotations_file, 'pk_anns', 'up')
        end

        %% plot waveforms and annotations (for further annotation)
        duration_of_signal = imp.t(end) - imp.t(1);
        NUMWINS = floor(duration_of_signal / up.win_length); clear duration_of_signal
        edge_time = up.edge_prop*up.win_length;
        for win = 1:NUMWINS
            % skip if annotating breaths and this window was deemed to be
            % low quality
            if up.mode == 2 && ~sum(win==high_qual_wins)
                continue
            end

            % Find window timings
            starttime = imp.t(1)+(win-1)*up.win_length;
            endtime = starttime + up.win_length;

            % Skip if this window has already been annotated
            savename = [up.paths.annotations_folder, 'ImP_SQI' num2str(s), '_' operator '_an.mat'];
            if exist(savename, 'file')
                temp = load(savename);
                rel_els = temp.pk_anns.t>=starttime & temp.pk_anns.t<=endtime;
                min_diff = min(diff(temp.pk_anns.t(rel_els)));
                if min_diff < 1
                    error(['Check window ', num2str(win), ', subj ', num2str(s)])
                end
                if sum(rel_els)>0
                    continue
                end
            end

            % Setup plot
            rel_els = find(imp.t >= starttime & imp.t < endtime);
            grey_els = find(imp.t >= (starttime - edge_time) & imp.t < (endtime + edge_time));
            scrsz = get(0,'ScreenSize');
            load(up.paths.annotations_file)

            % Plot
            h=figure('OuterPosition',[1 1 scrsz(3) scrsz(4)]);
            %save('plot_vars.mat', 'h', 'imp', 'Re', 'grey_paw_els', 'rel_paw_els', 'rel_Re_els', 'grey_Re_els', 'period', 's', 'subjects')
            [axis_h, h2] = annotation_plotter2(up.paths.annotations_file, h, imp, grey_els, rel_els, s, subjects, up, win, NUMWINS,starttime);

            % use pointer to identify breaths
            set(gcf, 'Pointer', 'cross');
            if up.mode == 1
                % Construct a questdlg with three options
                choice = MFquestdlg([0.1,0.1],'High quality?', ...
                    '', ...
                    'Yes','No','No');
                % Handle response
                switch choice
                    case 'Yes'
                        load(up.paths.annotations_file);
                        [qual.t,inds]=unique([qual.t; starttime]);
                        qual.v=[qual.v; 1]; qual.v = qual.v(inds);
                        [qual.t, inds] = sort(qual.t);
                        qual.v = qual.v(inds);
                        save(up.paths.annotations_file, 'qual');
                    case 'No'
                        load(up.paths.annotations_file);
                        [qual.t,inds]=unique([qual.t; starttime]);
                        qual.v=[qual.v; 0]; qual.v = qual.v(inds);
                        [qual.t, inds] = sort(qual.t);
                        qual.v = qual.v(inds);
                        save(up.paths.annotations_file, 'qual');
                end
            else
                set(h2,'ButtonDownFcn',{@Click_CallBack2 gca axis_h h2 up up.paths.annotations_file starttime});
                pause
            end
            close(h);

        end

        if up.DOSAVE
            fprintf(['~~~~~~ Saving subject ' num2str(s) '  ~~~~~~\n']);
            if exist(savename, 'file')
                old_data = load(savename);
            end
            new_data.up = up;
            temp = load(up.paths.annotations_file);
            % If there are annotations in both the old file and the current annotations file
            if exist('old_data', 'var') && sum(strcmp(fieldnames(old_data), 'pk_anns')) && ~isempty(old_data.pk_anns.t) ...
                    && sum(strcmp(fieldnames(temp), 'pk_anns')) && ~isempty(temp.pk_anns.t)
                all_t = [temp.pk_anns.t; old_data.pk_anns.t];
                all_v = [temp.pk_anns.v; old_data.pk_anns.v];
                [~, order] = sort(all_t);
                new_data.pk_anns.t = all_t(order);
                new_data.pk_anns.v = all_v(order);
                clear order all_t all_v
            % If there are annotations in the old file
            elseif exist('old_data', 'var') && sum(strcmp(fieldnames(old_data), 'pk_anns')) && ~isempty(old_data.pk_anns.t)
                new_data.pk_anns = old_data.pk_anns;
            % If there are annotations in the current annotations file
            elseif sum(strcmp(fieldnames(temp), 'pk_anns')) && ~isempty(temp.pk_anns.t)
                new_data.pk_anns = temp.pk_anns;
            else
                [new_data.pk_anns.t, new_data.pk_anns.v] =deal([]);
            end
            if exist('qual', 'var') && ~isempty(qual.t)
                new_data.qual = qual;
            elseif sum(strcmp(fieldnames(temp), 'qual'))  && ~isempty(temp.qual.t)
                new_data.qual = temp.qual;
            elseif exist('old_data', 'var') && sum(strcmp(fieldnames(old_data), 'qual'))  && ~isempty(old_data.qual.t)
                new_data.qual = old_data.qual;
            else
                new_data.qual.t = []; qual.v = [];
            end
            pk_anns = new_data.pk_anns;
            qual = new_data.qual;
            up = new_data.up;
            save(savename, 'pk_anns', 'qual', 'up')
            clear old_data new_data temp pk_anns qual
        end

        clear pk_anns qual NUMWINS pawmean STUDYID

    end

end

end

function up = setup_universal_params

close all

%% folderpaths
up.paths.root_data_folder = '/Users/petercharlton/Documents/Data/ImP_SQI_mimic_dataset_test/';
up.paths.data_file = [up.paths.root_data_folder, 'mimic_imp_sqi_data.mat'];
up.paths.data_file2 = [up.paths.root_data_folder, 'data_mat.mat'];
up.paths.resultspath = [up.paths.root_data_folder, 'SQI_development', filesep];
up.paths.annotations_folder = [up.paths.root_data_folder, '2019_annotations', filesep];
up.paths.annotations_file = [up.paths.root_data_folder, '2019_annotations', filesep, 'temp.mat'];
up.paths.dlg = '/Users/petercharlton/Google Drive/Work/Projects/PhD/Github/phd/Tools/Other Scripts/MFquestdlg/';
addpath(up.paths.dlg)
if ~exist(up.paths.annotations_folder, 'dir')
    mkdir(up.paths.annotations_folder)
end

%% waveform specifications
up.spikefs = 500;     % It is assumed in vortal_spike_conversion that all the spike waveforms are exported at this freq
up.reg = 1;
up.paworig = 2;
up.pawresample = 3;
up.resp = 4;
up.FS = [4, up.spikefs, 50, 62.5];

%% Operational specifications
up.DOSAVE = 1;                  % Save plots of filter characteristics?
up.DOPLOT = 0;
up.viewing_win_length = 32;     % duration of window to be viewed at once (in secs)
up.edge_prop = 0.3;             % proportion of window to show either side of the one currently being viewed.

%% Filtering specifications
up.filt_prop = 2;              % proportion of a second for which to apply the filter. (i.e. this * fs = order)
up.LP_CUTOFF = 60;                 % in rpm
up.time_offset = 1.1;          % The time (secs) by which the resp signal is delayed after the Paw signal.

%% Analysis specifications:
up.win_length = 32;            % window length for analysis in secs
up.true_BR_lims = 2;           % number of bpm either side at which ref BR estimates are considered reasonable.
up.mingap = 5;                 % If there is a gap in data of duration greater than this no. of secs, then this is considered to be a gap (anything less is ignored), for both numerics and waveforms

% duration of Tukey window taper in secs
up.paramSet.tukey_win_duration_taper = 2;

% Eliminate HFs (above resp freqs)
up.paramSet.elim_hf.Fpass = 1.2;  % in Hz
up.paramSet.elim_hf.Fstop = 0.8899;  % in Hz     (1.2 and 0.8899 provide a -3dB cutoff of 1 Hz)
up.paramSet.elim_hf.Dpass = 0.057501127785;
up.paramSet.elim_hf.Dstop = 0.01;

end

function [axis_h, h2] = annotation_plotter2(annotations_filepath, h, imp, grey_els, rel_els, s, subjects, up, win, NUMWINS, starttime)

%% ------------    annotation_plotter  ------------------------------
%
% Created by: Peter Charlton (17/10/2012)
%
% Purpose: To plot the Paw and impedance waveforms, and the annotated
%           breaths.
%
% Context: Some safety nets have been added to make it run ok even if the
% impedance signal isn't recorded.
%
% Inputs:
%
% Outputs:
%
% Files required:
%
% Called by:            Run_paw_verification
%
% ------------------------------------------------------------------------

ftsize = 16; lwidth = 2;
init_time = imp.t(rel_els(1));

load(annotations_filepath)
mean_imp = mean(imp.v(rel_els));
h3 = plot(imp.t(rel_els) - init_time, zeros(1,length(rel_els)), 'k', 'LineWidth', lwidth); hold on    % plot zero line on Paw plot
set(h3, 'HandleVisibility', 'off')
h1 = plot(imp.t(grey_els) - init_time, imp.v(grey_els) - mean_imp, 'Color', [0.5 0.5 0.5], 'LineWidth', lwidth); hold on,    % plot Paw (grey)
set(h1, 'HandleVisibility', 'off')
h2 = plot(imp.t(rel_els) - init_time, imp.v(rel_els) - mean_imp, 'LineWidth', lwidth);    % plot Paw
set(h2, 'HandleVisibility', 'off')
xlim([imp.t(grey_els(1)) - init_time, imp.t(grey_els(end)) - init_time])

%% Calculate and plot mix signal
if up.mode == 2
    % plot annotations which are relevant to this plot:
    %breathels = find(pk_anns.t>=imp.t(grey_els(1)) & pk_anns.t<=imp.t(grey_els(end)));
    %plot(pk_anns.t(breathels) - init_time,pk_anns.v(breathels)-mean_imp, 'ro','LineWidth',4)
    breathels = find(pk_anns.t>=imp.t(rel_els(1)) & pk_anns.t<=imp.t(rel_els(end)));
    plot(pk_anns.t(breathels) - init_time,pk_anns.v(breathels), 'ro','LineWidth',4)

    title(['Please annotate the peaks (red)  -  ' num2str(win) ' of ' num2str(NUMWINS) ' windows for subject ' num2str(s) ' of ' num2str(subjects(end))], 'FontSize', ftsize)
    ylabel('Impedance', 'FontSize', ftsize), xlabel('Time [s]', 'FontSize', ftsize);
    ylabel('Imepdance (above line = inhalation)', 'FontSize', ftsize)
    set(gca,'XTick',((imp.t(rel_els(1)) - init_time):5:ceil((imp.t(rel_els(end)) - init_time))), 'YTick', [])
    set(gca, 'FontSize', ftsize)
else
    title(['Please annotate high quality (left) or low quality (right)  -  ' num2str(win) ' of ' num2str(NUMWINS) ' windows for subject ' num2str(s) ' of ' num2str(subjects(end))], 'FontSize', ftsize)
    ylabel('Impedance', 'FontSize', ftsize), xlabel('Time [s]', 'FontSize', ftsize);
    ylabel('Imepdance (above line = inhalation)', 'FontSize', ftsize)
    set(gca,'XTick',((imp.t(rel_els(1)) - init_time):5:ceil((imp.t(rel_els(end)) - init_time))), 'YTick', [])
    set(gca, 'FontSize', ftsize)
    if find(qual.t == starttime)
        rel_qual = qual.v(qual.t == starttime);
        if rel_qual
            annotation('textbox',[0.75, 0.2, 0.1,0.1],'String','High Quality', 'Color', 'b', 'FontSize', ftsize*3, 'LineStyle', 'None')
        else
            annotation('textbox',[0.75, 0.2, 0.1,0.1],'String','Low Quality', 'Color', 'r', 'FontSize', ftsize*3, 'LineStyle', 'None')
        end
    end
    if exist('pk_anns', 'var')
        % plot annotations which are relevant to this plot:
        breathels = find(pk_anns.t>=imp.t(rel_els(1)) & pk_anns.t<=imp.t(rel_els(end)));
        plot(pk_anns.t(breathels) - init_time,pk_anns.v(breathels), 'ro','LineWidth',4)
    end
end

axis_h = gca;

end

function [pk_anns, TRS]= Click_CallBack2(h,e, a, axis_h, h2, up, annotations_filepath, start_time)

switch get(ancestor(a,'figure'),'SelectionType')

    case 'normal' %left click
        point = get(a,'CurrentPoint');
        load(annotations_filepath);
        [pk_anns.t,inds]=unique([pk_anns.t; point(1,1)+start_time]);
        pk_anns.v=[pk_anns.v; point(1,2)]; pk_anns.v = pk_anns.v(inds);
        [pk_anns.t, inds] = sort(pk_anns.t);
        pk_anns.v = pk_anns.v(inds);
        save(annotations_filepath, 'pk_anns');

    case 'alt'  % right click
        point = get(a,'CurrentPoint');
        load(annotations_filepath);
        [pk_anns.t,inds]=unique([pk_anns.t; point(1,1)]);
        pk_anns.v=[pk_anns.v; point(1,2)]; pk_anns.v = pk_anns.v(inds);
        [pk_anns.t, inds] = sort(pk_anns.t);
        pk_anns.v = pk_anns.v(inds);
        save(annotations_filepath, 'pk_anns');

    case 'extend' % right click whilst holding down shift
        load(annotations_filepath)
        pk_anns.t = pk_anns.t(1:(end-1));
        pk_anns.v = pk_anns.v(1:(end-1));
        save(annotations_filepath, 'pk_anns');

end

cla(axis_h)
plot(pk_anns.t-start_time,pk_anns.v, 'ro','LineWidth',4)

end

function paw_verification_script_checker(subjects, data, up)
% this function is used to check that the results can indeed be plotted
% back on the original waveform:
%% initialisation
fprintf('\n\n\n\n\n~~~~~~~  Reviewing Annotations  ~~~~~~~\n');

%% find out the file details
operator=input('Please type the first and last initial,\nof the operator whose annotations\nyou''d like to review, and press enter:\n','s');
question=input('\nWould you like to start from the \nfirst subject? (y or n, followed by enter)\n', 's');
if strcmp(question, 'n') == 1
    start_subject = input('\nWhich subject?\n', 's');
    start_subject = str2double(start_subject);
else
    start_subject = 1;
end

start_subject_el = find(subjects == start_subject);
if isempty(start_subject_el)
    fprintf('\nThis subject doesn''t exist. Try all over again.');
else
    subjects = subjects(start_subject_el:end);
    for s = subjects(1:end)

        %% Load Impedance
        imp.fs = data(s).ref.resp_sig.imp.fs;
        imp.v = data(s).ref.resp_sig.imp.v;
        imp.t = [0:(length(imp.v)-1)]/imp.fs;

        %% Process Impedance
        imp = lpf_to_exclude_resp(imp, up);
        imp.v = -1*imp.v;

        %% Load annotations
        annotations_filepath = [up.paths.annotations, 'RRLISTEN' num2str(s), '_' operator '-imp_an.mat'];
        up_copy = up;
        load(annotations_filepath);
        up = up_copy; clear up_copy

        %% plot waveforms and annotations
        duration_of_signal = imp.t(end) - imp.t(1);
        NUMWINS = floor(duration_of_signal / up.win_length); clear duration_of_signal
        edge_time = up.edge_prop*up.win_length;
        for win = 1:NUMWINS
            % Setup
            starttime = imp.t(1)+edge_time+(win-1)*up.win_length;
            endtime = starttime + up.win_length;
            rel_els = find(imp.t >= starttime & imp.t < endtime);
            grey_els = find(imp.t >= (starttime - edge_time) & imp.t < (endtime + edge_time));
            scrsz = get(0,'ScreenSize');

            % Plot
            h=figure('OuterPosition',[1 1 scrsz(3) scrsz(4)]);
            [axis_h, h2] = annotation_plotter2(annotations_filepath, h, imp, grey_els, rel_els, s, subjects, up, win, NUMWINS, starttime);

            % close plot after a pause
            pause
            close all
        end
    end
end

end

function imp = lpf_to_exclude_resp(imp, up)

imp.v(isnan(imp.v)) = median(imp.v(~isnan(imp.v)));

%% Window signal to reduce edge effects
duration_of_signal = imp.t(end) - imp.t(1);
prop_of_win_in_outer_regions = 2*up.paramSet.tukey_win_duration_taper/duration_of_signal;
tukey_win = tukeywin(length(imp.v), prop_of_win_in_outer_regions);
d_s_win = imp;    % copy time and fs
d_s_win.v = detrend(imp.v(:)).*tukey_win(:);

%% LPF to remove freqs above resp
respWave.t = d_s_win.t;
respWave.v = lp_filter_signal_to_remove_freqs_above_resp(d_s_win.v, d_s_win.fs, up);
respWave.fs = d_s_win.fs;
imp = respWave;
end

function s_filt = lp_filter_signal_to_remove_freqs_above_resp(s, Fs, up)
%% Filter pre-processed signal to remove freqs above resp

% parameters for the low-pass filter to be used
flag  = 'scale';
Dpass = up.paramSet.elim_hf.Dpass;
Dstop = up.paramSet.elim_hf.Dstop;
Fstop = up.paramSet.elim_hf.Fstop;
Fpass = up.paramSet.elim_hf.Fpass;

% create filter
[N,Wn,BETA,TYPE] = kaiserord([Fstop Fpass]/(Fs/2), [1 0], [Dstop Dpass]);
b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
AMfilter = dfilt.dffir(b);

%% Check frequency response
% Gives a -3 dB cutoff at ? Hz, using:
% freqz(AMfilter.Numerator)
% norm_cutoff_freq = 0.3998;    % insert freq here from plot
% cutoff_freq = norm_cutoff_freq*(Fs/2);

% Prepare signal
s_dt=detrend(s);
s_filt = filtfilt(AMfilter.numerator, 1, s_dt);
end

author: peterhcharlton

repo: RRest


text:
function linpack_s_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests SQRDC and SQRSL.
%
%  Discussion:
%
%    SQRDC and SQRSL compute the QR factorization, and use it
%    to solve linear systems.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  p = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  For a general matrix,\n' );
  fprintf ( 1, '  SQRDC computes the QR decomposition of a\n' );
  fprintf ( 1, '  matrix, but does not return Q and R explicitly.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Show how Q and R can be recovered using SQRSL.\n' );
%
%  Set the matrix A.
%
  a = [ ...
    1.0, 1.0, 0.0; ...
    1.0, 0.0, 1.0; ...
    0.0, 1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The original matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : p
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Decompose the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Decompose the matrix.\n' );

  job = 0;
  ipvt(1:p) = 0;

  [ a, qraux, ipvt ] = sqrdc ( a, lda, n, p, ipvt, job );
%
%  Print out what SQRDC has stored in A...
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The packed matrix A which describes Q and R:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : p
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  ...and in QRAUX.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The QRAUX vector, containing some additional\n' );
  fprintf ( 1, '  information defining Q:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %14f', qraux(i) );
  end
  fprintf ( 1, '\n' );
%
%  Print out the resulting R factor.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The R factor:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : p

      if ( j < i )
        r(i,j) = 0.0;
      else
        r(i,j) = a(i,j);
      end

    end

    for j = 1 : p
      fprintf ( 1, '  %14f', r(i,j) );
    end
    fprintf ( 1, '\n' );

  end
%
%  Call SQRSL to extract the information about the Q matrix.
%  We do this, essentially, by asking SQRSL to tell us the
%  value of Q*Y, where Y is a column of the identity matrix.
%
  job = 10000;

  for i = 1 : n
%
%  Set the vector Y.
%
    y(1:n) = 0.0;

    y(i) = 1.0;
%
%  Ask SQRSL to tell us what Q*Y is.
%
    [ qy, qty, b, rsd, xb, info ] = sqrsl ( a, lda, n, p, qraux, y, job );

    if ( info ~= 0 )
      fprintf ( 1, '  Error!  SQRSL returns INFO = %d\n', info );
      return
    end
%
%  Copy QY into the appropriate column of Q.
%
    q(1:n,i) = qy(1:n);

  end
%
%  Now print out the Q matrix we have extracted.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Q factor:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', q(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Compute Q*R to verify that it equals A.
%
  b(1:n,1:p) = q(1:n,1:n) * r(1:n,1:p);
%
%  Print the result.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The product Q * R:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : p
      fprintf ( 1, '  %14f', b(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function  [fx] = f_wsls(x,P,u,in)
% pseudo-evolution function for a "win/stay - lose/switch" strategy
% function  [fx,dfdx,dfdP] = f_wsls(x,P,u,in)
% The "win/stay - lose/switch" strategy is encoded in terms of the
% evolution of pseudo q-values, which swap sign depending upon the
% feedback. Note: the feedback should unambiguously code for "win" or
% "lose", i.e. in terms of a positive (resp. negative) reward.
% IN:
%   - x : pseudo q-values (1: stay, -1:switch)
%   - P : [useless]
%   - u : u(1)=previous action, u(2)=feedback
%   - in : [useless]
% OUT:
%   - fx: evolved pseudo q-values (2x1)
%   - dfdx/dfdP: [useless]

r = u(2);
if u(1)==1
    fx = sign(r)*[1;-1];
else
    fx = sign(r)*[-1;1];
end

author: MBB-team

repo: VBA-toolbox


text:
function ma_playMovieTimerCallback(timer, evt, hPlayBtn, hFig, mAxes, maData, frameRate, celBodyData, handles)

    userData = get(timer,'UserData');
    stateLogTime = userData(1);
    prevBody = userData(2);
    orbitToPlot = userData(3);

    stateLog = maData.stateLog;

    if(isnan(stateLogTime))
        stateLogTime = get(handles.movieFrameSlider,'Value');
    end

    if(get(hPlayBtn,'Value')==0 && ~isempty(evt))
        stop(timer); drawnow;

        return;
    end

    timeWarpMultiplier = get(handles.warpRateCombo,'UserData');
    [curBody] = ma_plotFrame(hFig, mAxes, maData, stateLog, stateLogTime, prevBody, orbitToPlot, celBodyData, handles);
    curStateLogTime = ma_getStateLogTimeToPlot(stateLogTime, frameRate, timeWarpMultiplier);

    if(prevBody~=curBody)
        orbitToPlot = orbitToPlot + 1;
    end
    prevBody = curBody;

    ma_setAnimatorSliderBounds(handles, stateLog);
    if(curStateLogTime > max(stateLog(:,1)) && stateLogTime < max(stateLog(:,1)))
        curStateLogTime = max(stateLog(:,1));
    end

    if(curStateLogTime > max(stateLog(:,1)))
        stop(timer);
        set(hPlayBtn,'Value',0);
        set(handles.movieFrameSlider,'Value',max(stateLog(:,1)));

        writeVideo();

        ma_handlePlayButton(handles, frameRate);
        set(handles.movieFrameSlider,'Value',get(handles.movieFrameSlider,'Min'));
        cb = get(handles.movieFrameSlider,'Callback');
        cb(handles.movieFrameSlider,[]);
        return;
    else
        set(handles.movieFrameSlider,'Value',curStateLogTime);
        set(handles.sliderUTText,'String',fullAccNum2Str(curStateLogTime));
    end
    stateLogTime = curStateLogTime;
    set(timer,'UserData',[stateLogTime, prevBody, orbitToPlot]);
end

function writeVideo()
    global GLOBAL_VideoWriter GLOBAL_isRecording;

    if(GLOBAL_isRecording)
        close(GLOBAL_VideoWriter);

        helpMsg = sprintf('Your recording has completed and is available at the following location:\n\n%s',GLOBAL_VideoWriter.Filename);
        h = msgbox(helpMsg,'Recording Complete','help','modal');
        uiwait(h);

        GLOBAL_isRecording = false;
        GLOBAL_VideoWriter = [];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Wgrid = Fdirs2grid(W, aziRes, polarRes, CLOSED)
%FDIRS2GRID Replicate vector function values on a regular grid
%
%   Fdirs2grid takes a vector of values of function evaluated at a
%   spherical grid with the grid2dirs function, and convert it back to a
%   2D grid. Useful for plotting with functions such as surf, or for
%   numerical integration numerically spherical functions.
%
%   W:  column vector of function values evaluated at each grid direction,
%       with the direction ordering given by grid2dirs. If W is a matrix
%       then each column is considered as a separate function to be
%       converted
%   aziRes: azimuth resolution of the grid in degrees (should be the same
%           as the one used in the grid2dirs function
%   polarRes:   inclination resolution of the grid in degrees (should be
%               the same as the one used in the grid2dirs function
%   CLOSED: {0,1} if true then the returned matrix replicates the first
%           column of function values at 0deg azimuth also at 360deg,
%           useful for 3D plotting so that the shape does not have a
%           hole in the end (see plotSphFunction test script)
%
%   Wgrid:  if W is a vector then Wgrid is the 2D matrix of the function
%           values replicated on the grid points. If W is a matrix, then
%           Wgrid is a 3D matrix with one grid per column of W, on the 3rd
%           dimension.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 10/10/2013
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if mod(360, aziRes) ~= 0 || mod(180, polarRes) ~= 0
    error('azimuth or elevation resolution should divide exactly 360 and 180deg')
end

if nargin<4
    CLOSED = 0;
end

Nphi = 360/aziRes;
Ntheta = 180/polarRes+1;

Nf = size(W, 2);
Wgrid = zeros(Nphi, Ntheta, Nf);
for i = 1:Nf

    Wgrid(:, 2:end-1, i) = reshape(W(2:end-1, i), Nphi, Ntheta-2);
    Wgrid(:, 1, i) = ones(Nphi, 1) * W(1, i);
    Wgrid(:, end, i) = ones(Nphi, 1) * W(end, i);
end

if Nf~=1
    Wgrid = permute(Wgrid, [2 1 3]);
else
    Wgrid = Wgrid.';
end

if CLOSED
    Wgrid = horzcat(Wgrid, Wgrid(:,1,:));
end

author: polarch

repo: Spherical-Harmonic-Transform


text:
%TEST_SmoothAbs
%
% This script plots the output of SmoothAbs for a few values of alpha, over
% the domain [-1,1]
%
% Written by Matthew Kelly
% October 2013
% Cornell University
%

t = linspace(-1,1,1000);
alpha = [0.02,0.1,0.5];
N=length(alpha);

figure(104); clf; hold on;
for i=1:N
    x1 = SmoothAbs(t,alpha(i));

    subplot(N,1,i); hold on;
    plot(t,x1,'b-','LineWidth',2)
    plot(t,abs(t),'k:')
    title(['Alpha = ' num2str(alpha(i))],'FontSize',14);
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [nodes, edges2] = grRemoveEdge(nodes, edges, edge)
%GRREMOVEEDGE Remove an edge in a graph.
%
%   [NODES2 EDGES2] = grRemoveEdge(NODES, EDGES, EDGE2REMOVE)
%   Remove an edge in the edges list, and return the modified graph.
%   The NODES array is not modified.

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-08-13
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

dim = size(edges);
edges2 = zeros(dim(1)-1, 2);
edges2(1:edge-1, :) = edges(1:edge-1, :);
edges2(edge:dim(1)-1, :) = edges(edge+1:dim(1), :);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [dat, mask] = apply_mask(dat, mask, varargin)
% Apply a mask image (image filename or fmri_mask_image object) to an image_vector object
% stored in dat.
%
% This can be used to:
%  - Mask an image_vector or fmri_data object with a mask
%  - Obtain "pattern expression" for a weight map (entered as the
%    mask, here) in a series of images stored in dat.
%
% The mask or weight map does not have to be in the same space as the dat;
% it will be resampled to the space of the data in dat.
%
% To extract pattern expression values for each ROI within a mask use extract_roi_averages()
%
% :Optional Inputs:
%
%   **pattern_expression:**
%        calculate and return the dot product of each
%        image in dat and the values in the mask.  This is useful if comparing
%        expression values that are comprised of different datasets or differing
%        number of voxels.
%
%   **correlation:**
%        calculate the pearson correlation coefficient of each
%        image in dat and the values in the mask.
%
%   **norm_mask:**
%        normalize the mask weights by L2 norm, for patt expression
%        only.
%
%   **ignore_missing:**
%        use with pattern expression only. Ignore weights on voxels
%        with zero values in test image. If this is not entered, the function will
%        check for these values and give a warning.
%
%   **invert:**
%        Invert the mask so that out-of-mask voxels are now in (using
%        the mask as an 'exclude mask' rather than an include-mask. If pattern
%        expression is requested, the behavior is different, and it inverts the
%        sign of in-mask pattern weights.
%
%   **cosine_similarity:**
%        use with pattern expression only. scales expression by product of
%        l2 norms (norm(mask)*norm(dat))
%
% :Examples:
% ::
%
%     [dat, mask] = apply_mask(dat, mask)
%     [dat, mask] = apply_mask(dat, mask image name)
%     [dat, mask] = apply_mask(dat, mask image vector object)
%     [pattern_exp_values] = apply_mask(dat, weight map image, 'pattern_expression', 'ignore_missing')
%     [pattern_exp_values] = apply_mask(dat, weight map image, 'pattern_expression', 'ignore_missing','correlation')
%
%
% :See also:
%
% extract_roi_averages, to get individual region averages / local pattern expression
% apply_nps, which does whole-pattern and local regional expression
%
% ..
%    Notes:
%    Last modified: 10/30/11 to add support for masks that are weight maps
%    12/15/13:  Luke Chang - added correlation option for pattern-expression
%    5/10/2016: Phil Kragel - added cosine similarity
%    8/2/2018: Tor Wager - if mask is a statistic_image object, apply .sig field
%
% ..

dopatternexpression = 0; % set options
donorm = 0;
doinvert = 0;

% Now handled by canlab_pattern_similarity
% doignoremissing = 0;
% docorr = 0; %run correlation instead of dot-product for pattern expression
% docosine = 0;

if any(strcmp(varargin, 'pattern_expression'))
    dopatternexpression = 1;

    % Now handled by canlab_pattern_similarity
    %     if any(strcmp(varargin, 'ignore_missing'))
    %         doignoremissing = 1;
    %     end
    %
    %     if any(strcmp(varargin, 'cosine_similarity'))
    %         docosine = 1;
    %     end
    %
    %
    %     if any(strcmp(varargin, 'correlation')) % run correlation instead of dot-product
    %         docorr = 1;
    %     end

end

if any(strcmp(varargin, 'invert'))
    doinvert = 1;
end

if any(strcmp(varargin, 'norm_mask')) % only good for pattern expression
    donorm = 1;
end

% create mask_image object if we have a filename
if ischar(mask)
    mask = fmri_mask_image(mask);
end

isdiff = compare_space(dat, mask);

if isdiff == 1 || isdiff == 2 % diff space, not just diff voxels
    % == 3 is ok, diff non-empty voxels

    % Both work, but resample_space does not require going back to original
    % images on disk.
    %mask = resample_to_image_space(mask, dat);
    mask = resample_space(mask, dat);

    % tor added may 1 - removed voxels was not legal otherwise
    %mask.removed_voxels = mask.removed_voxels(mask.volInfo.wh_inmask);
    % resample_space is not *always* returning legal sizes for removed
    % vox? maybe this was updated to be legal

    if length(mask.removed_voxels) == mask.volInfo.nvox
        disp('Warning: resample_space returned illegal length for removed voxels. Fixing...');
        mask.removed_voxels = mask.removed_voxels(mask.volInfo.wh_inmask);
    end

end

% dat = remove_empty(dat); Tor edited 9/10/17 for computational efficiency

% nonemptydat: Logical index of voxels with valid data, in in-mask space
nonemptydat = get_nonempty_voxels(dat);

dat = replace_empty(dat);

% Check/remove NaNs. This could be done in-object...
mask.dat(isnan(mask.dat)) = 0;

% Replace if necessary
mask = replace_empty(mask);

% if mask is a statistic_image object, apply .sig field
if isa(mask, 'statistic_image')

    mask.dat(~mask.sig) = 0;  % for masking purposes

end

if doinvert

    if dopatternexpression
        mask.dat = -mask.dat;
        disp('Reversing sign of mask values.');

    else
        disp('Inverting mask: Out-of-mask areas are now in.');
        maskdat = reconstruct_image(mask);
        maskdat = double(~maskdat);
        mask = rebuild_volinfo_from_dat(mask, maskdat(:));
        mask = resample_space(mask, dat);
    end

end

% save which are in mask, but do not replace with logical, because mask may
% have weights we want to preserve
inmaskdat = logical(mask.dat);

% Remove out-of-mask voxels
% ---------------------------------------------------

% mask.dat has full list of voxels
% need to find vox in both mask and original data mask

if size(mask.volInfo.image_indx, 1) == size(dat.volInfo.image_indx, 1)
    n = size(mask.volInfo.image_indx, 1);

    if size(nonemptydat, 1) ~= n % should be all vox OR non-empty vox
        nonemptydat = zeroinsert(~dat.volInfo.image_indx, nonemptydat);
    end

    if size(inmaskdat, 1) ~= n
        inmaskdat = zeroinsert(~mask.volInfo.image_indx, inmaskdat);
    end

    inboth = inmaskdat & nonemptydat;

    % List in space of in-mask voxels in dat object.
    % Remove these from the dat object
    to_remove = ~inboth(dat.volInfo.wh_inmask);

    to_remove_mask = ~inboth(mask.volInfo.wh_inmask);

elseif size(mask.dat, 1) == size(dat.volInfo.image_indx, 1)

    % mask vox are same as total image vox
    nonemptydat = zeroinsert(~dat.volInfo.image_indx, nonemptydat);
    inboth = inmaskdat & dat.volInfo.image_indx & nonemptydat;

    % List in space of in-mask voxels in dat object.
    to_remove = ~inboth(dat.volInfo.wh_inmask);

    to_remove_mask = ~inboth(mask.volInfo.wh_inmask);

elseif size(mask.dat, 1) == size(dat.volInfo.wh_inmask, 1)
    % mask vox are same as in-mask voxels in dat
    inboth = inmaskdat & dat.volInfo.image_indx(dat.volInfo.wh_inmask) & nonemptydat;

    % List in space of in-mask voxels in .dat field.
    to_remove = ~inboth;

    to_remove_mask = ~inboth;

else
    fprintf('Sizes do not match!  Likely bug in resample_to_image_space.\n')
    fprintf('Vox in mask: %3.0f\n', size(mask.dat, 1))
    fprintf('Vox in dat - image volume: %3.0f\n', size(dat.volInfo.image_indx, 1));
    fprintf('Vox in dat - image in-mask area: %3.0f\n', size(dat.volInfo.wh_inmask, 1));
    disp('Stopping to debug');
    keyboard
end

dat = remove_empty(dat, to_remove);
mask = remove_empty(mask, to_remove_mask);

if dopatternexpression
    %mask = replace_empty(mask); % need for weights to match

    weights = double(mask.dat); % force double b/c of matlab instabilities
    dat.dat = double(dat.dat); % force double b/c of matlab instabilities

    if donorm, weights = weights ./ norm(weights); end

    % Pass similarity metric in to canlab_pattern_similarity
    similarity_output = canlab_pattern_similarity(dat.dat, weights, varargin{:});

    dat = similarity_output;

    % CHECK and weight
    % ---------------------------------------------------------

    % Pattern weights can have zeros, which may be valid values in voxels,
    % i.e., with binary masks
    % Images with values of zero or NaN are considered out-of-mask, as they
    % are not valid values. These should be excluded from both image and
    % mask when calculating similarity.
    % Thus, there is an asymmetry between pattern mask and image data
    % in considering which voxels to use.
    % Otherwise, all dot product and similarity metrics are standard.

    % We also need to calculate bad values on an image-by-image basis, not
    % relying on remove_empty to exclude voxels with ineligible values
    % across the entire set of images.

    % End check and weight ---------------------------------------------------------

end % Pattern expression

% Handle special types
% So far, atlas

if isa(dat, 'atlas')

    dat = check_properties(dat, 'compress_index');

end

end  % Main function

function nonemptydat = get_nonempty_voxels(dat)
empty_voxels = all(dat.dat' == 0 | isnan(dat.dat'), 1)';

if size(empty_voxels, 1) == dat.volInfo.n_inmask
    % .dat is full in-mask length, we have not called remove_voxels or there are none to remove
    nonemptydat = ~empty_voxels;

elseif length(dat.removed_voxels) == dat.volInfo.n_inmask
    % .dat is not in-mask length, and we have .removed_voxels defined
    nonemptydat = false(dat.volInfo.n_inmask, 1);
    nonemptydat(~dat.removed_voxels) = true;

    % additional: we could have invalid voxels that have been changed/added since
    % remove_empty was last called.
    %     dat.removed_voxels(dat.removed_voxels) =
    %     nonemptydat(empty_voxels
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function hsp= gridutil_getSubplots(chans)
%GRIDUTIL_GETSUBPLOTS - Returns the handles of the subplots of a grid plot
% that correspond to given channels
%
%Synposis:
% HSP= gridutil_getSubplots(<CHANS>)
%
%Input:
% CHANS: channel indices, [] means all, default []

if ~exist('chans','var'), chans=[]; end
if ~isempty(chans) && ~iscell(chans),
  chans= {chans};
end

if ~isempty(chans) && ischar(chans{1}) && ...
      strcmp(chans{1},'plus'),
  search_type= '^ERP';   % starting with 'ERP'
  chans= chans(2:end);
else
  search_type= '^ERP$';  % only the word 'ERP'
end

hc= get(gcf, 'children');
isERPplot= zeros(size(hc));

for ih= 1:length(hc),
  ud= get(hc(ih), 'userData');
  if isstruct(ud) && isfield(ud,'type') && ischar(ud.type) && ...
        ~isempty(regexp(ud.type,search_type))
    if isempty(chans) || ~isempty(util_chanind(ud.chan,chans)),
      isERPplot(ih)= 1;
    end
  end
end

hsp= hc(find(isERPplot))';

author: bbci


text:
function test_suite = test_survival_coxph

%   Run specific demo and save values for comparison.
%
%   See also
%     TEST_ALL, DEMO_SURVIVAL_COXPH

% Copyright (c) 2011-2012 Ville Tolvanen

initTestSuite;

function testDemo
% Set random number stream so that failing isn't because randomness. Run
% demo & save test values.
prevstream=setrandstream(0);

disp('Running: demo_survival_coxph')
demo_survival_coxph;
path = which('test_survival_coxph.m');
path = strrep(path,'test_survival_coxph.m', 'testValues');
if ~(exist(path, 'dir') == 7)
    mkdir(path)
end
path = strcat(path, '/testSurvival_coxph');
save(path, 'Ef1', 'Ef2', 'Varf1', 'Varf2');

% Set back initial random stream
setrandstream(prevstream);
drawnow;clear;close all

% Compare test values to real values.

function testPredictionsCoxph
values.real = load('realValuesSurvival_coxph', 'Ef1', 'Varf1', 'Ef2', 'Varf2');
values.test = load(strrep(which('test_survival_coxph.m'), 'test_survival_coxph.m', 'testValues/testSurvival_coxph'), 'Ef1', 'Varf1', 'Ef2', 'Varf2');
assertElementsAlmostEqual(values.real.Ef1, values.test.Ef1, 'absolute', 0.10);
assertElementsAlmostEqual(values.real.Ef2, values.test.Ef2, 'absolute', 0.10);
assertElementsAlmostEqual(values.real.Varf1, values.test.Varf1, 'absolute', 0.10);
assertElementsAlmostEqual(values.real.Varf2, values.test.Varf2, 'absolute', 0.10);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x_opt = coordinate_search ( x, function_handle, flag )

%*****************************************************************************80
%
%% COORDINATE_SEARCH carries out a direct search minimization algorithm.
%
%  Discussion:
%
%    This function implements a direct search algorithm, to seek a
%    multidimensional point X_OPT which minimizes a scalar function F(X).
%
%    The function is currently set up for coordinate search by
%    using a default stencil along coordinate axes.  The stencil and the
%    initial size of the stencil (delta) should probably be input parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2009
%
%  Author:
%
%    Jeff Borggaard
%
%  Parameters:
%
%    Input, real X(N_VAR), a starting estimate for the minimizer.
%
%    Input, function handle FUNCTION_HANDLE ( X ), ...
%
%    Optional input, integer FLAG, selects features of the algorithm.
%    If FLAG is not input, it defaults to the value 0.
%    0, do not include graphical displays.
%    nonzero, display the status of the search graphically.
%    3, use the water stick stencil.
%
%    Output, real X_OPT(N_VAR), the estimated minimizer.
%
%  Local Parameters:
%
%    Local, real DELTA, the initial size of the stencil.
%
%    Local, real TOLERANCE, a tolerance for the variation between the maximum
%    and minimum values of F in the stencil.
%
%    Local, integer MAX_FEVAL, the maximum number of function evaluations.
%
  delta = 1.0;
  tolerance = 1.0E-06;
  max_feval = 250;
  verbose = 0;
%
%  The default value of FLAG is 0.
%
  if ( nargin <= 2 )
    flag = 0;
  end

  x = x(:)';
  n_var = length(x);
%
%  Define the search stencil
%
  if ( flag ~= 3 )

    for i = 1 : n_var
      stencil(2*i-1).v    = zeros(size(x));
      stencil(2*i-1).v(i) = 1;
      stencil(2*i  ).v    = zeros(size(x));
      stencil(2*i  ).v(i) =-1;
    end

  else

    stencil(1).v = [   1        0    ];
    stencil(2).v = [ -1/2  sqrt(3)/2 ];
    stencil(3).v = [ -1/2 -sqrt(3)/2 ];

  end

  n_sten = length ( stencil );
%
%  Compute the function value at the initial point.
%
  f = feval ( function_handle, x );
  n_feval = 1;

  if ( flag )
    xp = linspace ( -5.0, 5.0, 101 );
    yp = xp;
    for i=1:101
      for j=1:101
        fp(j,i) = feval(function_handle,[xp(i),yp(j)]);
      end
    end

    figure(27)
    hold on
    contour(xp,yp,fp,linspace(0,200,25))
  end
%
%  Begin the Direct Search iteration.
%
  while ( 1 )

    x_s = compute_points ( x, stencil, delta );
%
%  Plot the current search pattern.
%
    if ( flag )
      for i=1:n_sten
        plot([x(1) x_s(i,1)],[x(2) x_s(i,2)],'r')
        plot(x_s(i,1),x_s(i,2),'r*')
      end
    end

    f_s = evaluate ( x_s, function_handle );
    n_feval = n_feval + n_sten;
%
% "erase" the current search pattern after it has been considered
%
    if ( flag )

      for i = 1 : n_sten
        plot ( [x(1) x_s(i,1)],[x(2) x_s(i,2)],'b')
        plot ( x_s(i,1),x_s(i,2),'b*')
      end
    end
%
%  Find the element in the pattern with the lowest function value.
%
    [ fc_min, i_min ] = min ( f_s );
    [ fc_max        ] = max ( f_s );
%
%  If we found a lower value, move to that location.
%
    if ( fc_min < f )

      xb = x_s(i_min,:);
      shift = xb - x;

      x = xb;
      f = fc_min;

      for i = 1 : n_sten
        x_s(i,:) = x_s(i,:) + shift;
      end
%
%  Otherwise, shrink the pattern about the center and recalculate the stencil.
%
    else
      delta = delta / 2.0;
      x_s = compute_points ( x, stencil, delta );
    end
%
%  Test whether we should continue.
%
    converged = ( fc_max - fc_min < tolerance ) || ( delta < tolerance );

    if ( converged )
      break
    end

    diverged = ( max_feval < n_feval );

    if ( diverged )
      break
    end

  end

  if ( verbose )
    fprintf ( 'The best point x^* was: %d %d\n', x(1,:) );
    fprintf ( 'f(x^*) = %d\n', f(1) );
  end

  x_opt = x(1,:);

  if ( verbose )
    fprintf ( 'The algorithm terminated after %d function evaluations\n',...
      n_feval);
  end

  if ( diverged )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Convergence was not achieved.\n' );
    fprintf ( 1, '  The iteration limit of %d was exceeded.\n', max_feval );
  end

  return
end
function f = evaluate ( x, function_handle )

%*****************************************************************************80
%
%% EVALUATE evaluates the function at the stencil points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2009
%
%  Author:
%
%    Jeff Borggaard
%
%  Parameters:
%
%    Input, real X(N_STEN,N_VAR), a set of N_STEN points, each of
%    dimension N_VAR.
%
%    Input, f = FUNCTION_HANDLE ( x ), the handle of a function which
%    accepts as input a vector X and returns the scalar function value F.
%
%    Output, real F(N_STEN), the function evaluated at each of the points.
%
  [ n_sten, n_var ] = size ( x );

  f = zeros ( 1, n_sten );

  for i = 1 : n_sten
    f(i) = feval ( function_handle, x(i,:) );
  end

  return
end
function x_s = compute_points ( x, stencil, delta )

%*****************************************************************************80
%
%% COMPUTE_POINTS computes all of the points in the stencil.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2009
%
%  Author:
%
%    Jeff Borggaard
%
%  Parameters:
%
%    Input, real X(N_VAR), the center of the stencil.
%
%    Input, real STENCIL, a structure describing the stencil.
%
%    Input, real DELTA, determines the geometric size of the stencil.
%
%    Output, real X_S(N_STEN,N_VAR), the stencil points.
%
  n_var = length(x);
  n_sten = length(stencil);

  x_s = zeros ( n_sten, n_var );

  for i = 1 : n_sten
    x_s(i,:) = x + delta * stencil(i).v;
  end

  return
end


text:
function c = mpower(a,b)
% MPOWER implements c = a^b, where either a or b is an adiff object.
% For adiff objects, a^b is taken to be equivalent to a.^b

c=power(a,b);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function y = normalizeAngle(u)
%#codegen
% normalizes angle to -pi to pi
y = u;

y(u>pi) = mod(u(u>pi) + pi, 2*pi) - pi;
y(u<-pi) = -(mod(-(u(u<-pi) - pi), 2*pi) - pi);

author: TUMFTM