text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: September 9th, 2016
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Your_Plates_and_Input_Files()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  256;       % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
dx = Lx/Nx;      % Grid spatial resolution

%
% Immersed Structure Geometric / Dynamic Parameters %
%
ds = 0.5*dx;              % Lagrangian Pt. Spacing (2x resolution of Eulerian grid)
struct_name = 'plates';   % Name for .vertex, .spring, etc files. (must match what's in 'input2d')

% Call function to construct geometry
[xLag,yLag,dist] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,Lx);

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Lx]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
k_Spring = 2.5e4;                    % Spring stiffness (does not need to be equal for all springs)
ds_Plates = dist;                    % Spring resting length (does not need to be equal for all springs)
print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Plates,ds,struct_name);

% Prints .beam file!
% k_Beam = 0.5;                      % Beam Stiffness (does not need to be equal for all beams)
% C = compute_Curvatures(xLag,yLag)  % Computes curvature of initial configuration
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
%k_Target = 1e7;
%print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag); % Total # of Lag. Pts

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called 'struct_name'.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Plates,ds,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N-2 + N/2 );    % Print # of springs

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
            elseif s > N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
            end
    end

    % SPRINGS ACROSS PLATES
    for s=1:N/2
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+N/2, k_Spring, ds_Plates);
    end

    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called 'struct_name'.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s) );
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C(s) );
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes "curvature" of starting configuration
%
% NOTE: not curvature in the traditional geometric sense, in the 'discrete'
% sense through cross product.
%
% NOTE: assumes a CLOSED structure
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = compute_Curvatures(xLag,yLag)

N = length(xLag);
C = zeros( N );

%Note: needs to be done same order as you print .beam file!
for i=1:N

   % Pts Xp -> Xq -> Xr (same as beam force calc.)

   if ( (i > 1) && (i < N) )

        Xp = xLag(i-1); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(i-1); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==1)

        Xp = xLag(N); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(N); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==N)

        Xp = xLag(N-1); Xq = xLag(N); Xr = xLag(1);
        Yp = yLag(N-1); Yq = yLag(N); Yr = yLag(1);

   end

   C(i) = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp); %Cross product btwn vectors

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag,dist] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,L)

% ds: Lagrangian pt. spacing
% Nx: Eulerian grid resolution
% L:  Length of computational domain

yLag = 0:ds:L/5;             % Make plate of length L/10
xLag = zeros(1,length(yLag)); % Make corresponding xPts for Vertical Plate

% Translate points
yLag = yLag + (L/2-L/10);  % Translate them symmetrically

% Make sides
xLag_Left = xLag + (L/2-L/10);
xLag_Right = xLag + (L/2+L/10);

% Combine Pts Into ONE Vector
xLag = [xLag_Left xLag_Right];
yLag = [yLag yLag];

% Distance between plates
dist = (L/2+L/10) - (L/2-L/10);

% Plot the Geometry
% plot(xLag,yLag,'*'); hold on;
% axis([0 L 0 L]);

author: nickabattista


text:
function linpack_d_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests DGBFA and DGBDI.
%
%  Discussion:
%
%    Matrix A is ( 2 -1  0  0  0)
%                (-1  2 -1  0  0)
%                ( 0 -1  2 -1  0)
%                ( 0  0 -1  2 -1)
%                ( 0  0  0 -1  2)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    N is the number of equations.
%
%    ML is the number of subdiagonals,
%    MU the number of superdiagonals.
%
%    LDA is the leading dimension of the array used to store the
%    matrix, which must be at least 2*ML+MU+1.
%
  n_max = 128;
  ml = 1;
  mu = 1;
  lda = 2*ml+mu+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  DGBFA factors the matrix,\n' );
  fprintf ( 1, '  DGBDI computes the determinant as\n' );
  fprintf ( 1, '    det = MANTISSA * 10^EXPONENT\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Find the determinant of the -1,2,-1 matrix\n' );
  fprintf ( 1, '  for N = 2, 4, 8, 16, 32, 64, 128.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  (For this matrix, det ( A ) = N + 1.)\n' );
%
%  Set the matrix A.
%
  m = ml + mu + 1;
  fprintf ( 1, '  The bandwidth of the matrix is %d\n', m );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       N    Mantissa       Exponent\n' );
  fprintf ( 1, '\n' );

  n = 1;

  for n_log = 1 : 7

    n = 2 * n;

    a(1:lda,1:n) = 0.0;

    for j = 1 : n
      a(m-1,j) = -1.0;
      a(m,j) =    2.0;
      a(m+1,j) = -1.0;
    end

    [ a, ipivot, info ] = dgbfa ( a, lda, n, ml, mu );

    if ( info ~= 0 )
      fprintf ( 1, '  Error!  DGBFA returns INFO = %d\n', info );
      return
    end

    det = dgbdi ( a, lda, n, ml, mu, ipivot );

    fprintf ( 1, '  %6d  %14f  %14f\n', n, det(1), det(2) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = in_fread_avg(sFile, sfid, SamplesBounds)
% IN_FREAD_AVG:  Read an epoch from a Neuroscan .avg file (averaged file).
%
% USAGE:  F = in_fread_eeg(sFile, sfid, SamplesBounds) : Read selected times
%         F = in_fread_eeg(sFile, sfid)                : Read all file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, 2009-2011

%% ===== PARSE INPUTS =====
if (nargin < 3)
    SamplesBounds = [];
end
nTime = sFile.header.data.pnts;
nChannels = sFile.header.data.nchannels;

%% ===== READ DATA =====
% Initialize data structure
F = zeros(nChannels, nTime);
% Position at the beginning of the data block
fseek(sfid, double(sFile.header.data.datapos), 'bof');
% Read averaged data [nChannels x nTime]
for i = 1:nChannels
    % Read the channel data
    unused_header = fread(sfid, 5, 'char');
    F(i,:) = fread(sfid, nTime, 'float')';
end
% Keep only specific time indices
if ~isempty(SamplesBounds)
    iTimes = (SamplesBounds(1):SamplesBounds(2)) - round(sFile.prop.times(1) .* sFile.prop.sfreq) + 1;
    F = F(:,iTimes);
end
% Calibrate data
F = neuroscan_apply_calibration(F, sFile.header);

author: brainstorm-tools


text:
function [w,run] = train_dual(x,w,lambda,priority,mod)

% Written by Thomas P Minka

flops(0);
% kernel matrix
c = x'*x;
dc = diag(c);
v = 1/lambda;
% for flop counting, multiplication by lambda is free

[d,n] = size(x);
% this assumes w = 0
alpha = repmat(1e-4,n,1);
% Keerthi-type initialization
alpha = repmat(1/n,n,1);
if nargin < 4
  priority = 1;
end
if nargin < 5
  mod = 0;
end

count = 1;
for iter = 1:2000
  old_alpha = alpha;
  % annealing for expt2
  %v = v + 1e2
  if ~priority
    if mod
      for i = 1:n
	theta = log(alpha(i)./(1-alpha(i)));
	g = v*(c(i,:)*alpha) + theta;
	a = alpha(i)*(1-alpha(i));
	h = v*dc(i) + 1/a;
	%h = h*a*a;
	%h = h*a + g*(1-2*alpha(i));
	if 1
	  a = alpha(i);
	  c2 = ((1+log(a))*a*h-g)/((1+log(a))*a/(1-a) + 1+log(1-a));
	  c1 = a*(h - c2/(1-a));
	  as = linspace(eps,1-eps,100);
	  f = [];
	  f2 = [];
	  f0 = 0.5*v*(alpha'*c*alpha) + ...
	      sum(alpha.*log(alpha)) + sum((1-alpha).*log(1-alpha));
	  for k = 1:length(as)
	    alpha(i) = as(k);
	    f(k) = 0.5*v*(alpha'*c*alpha) + ...
		sum(alpha.*log(alpha)) + sum((1-alpha).*log(1-alpha));
	    f2(k) = c1*as(k)*log(as(k))+c2*(1-as(k))*log(1-as(k));
	    f2(k) = f2(k) - (c1*a*log(a)+c2*(1-a)*log(1-a)) + f0;
	  end
	  if 0
	    plot(as,f,as,f2)
	    drawnow
	    ax = axis;
	    line([a a],[ax(3) ax(4)],'Color','r')
	    axis_pct
	    pause
	  end
	end
	theta = theta - g/h;
	alpha(i) = 1/(1+exp(-theta));
      end
      flops(flops + n*(6+2+flops_exp));
    else
      for i = 1:n
	g = v*(c(i,:)*alpha) + log(alpha(i)./(1-alpha(i)));
	h = v*dc(i) + 1/alpha(i)/(1-alpha(i));
	alpha(i) = alpha(i) - g./h;
	if alpha(i) < eps
	  alpha(i) = eps;
	elseif alpha(i) > 1-eps
	  alpha(i) = 1-eps;
	end
      end
    end
    flops(flops + n*(flops_mul(1,n,1)+3+flops_exp + 4 + 2));
  else
    % incremental algorithm
    if iter == 1
      g = v*(c*alpha) + log(alpha./(1-alpha));
      flops(flops + flops_mul(c,alpha)+n*(flops_exp+3));
    end
    for j = 1:n
      [dummy,i] = max(abs(g));
      %disp(['i=' num2str(i) ' g=' num2str(g(i))])
      %i = j;
      o_alpha = alpha(i);
      h = v*dc(i) + 1/alpha(i)/(1-alpha(i));
      if mod
	theta = log(alpha(i)./(1-alpha(i)));
	a = alpha(i)*(1-alpha(i));
	gt = g(i)*a;
	h = h*a*a + gt*(1-2*alpha(i));
	theta = theta - g/h;
	alpha(i) = 1/(1+exp(-theta));
      else
	alpha(i) = alpha(i) - g(i)./h;
      end
      if alpha(i) < eps
	alpha(i) = eps;
      elseif alpha(i) > 1-eps
	alpha(i) = 1-eps;
      end
      % update all g(i)
      da = alpha(i) - o_alpha;  % no cost
      dg = da*c(:,i)*v;
      dg(i) = dg(i) + log(alpha(i)/(1-alpha(i))) - ...
	  log(o_alpha/(1-o_alpha));
      % no cost for second log
      g = g + dg;
    end
    flops(flops + n*(2*n-1 + 4 + 2 + 2) + n*(n + flops_exp+4 + n));
  end

  % computations here don't count
  w = v*x*alpha;
  run.w(:,count) = w;
  run.flops(count) = flops;
  run.e(count) = logProb(x,w) -0.5/v*w'*w;
  e2(count) = 0.5/v*w'*w + sum(alpha.*log(alpha)) + sum((1-alpha).*log(1-alpha));
  count = count + 1;
  if rem(count,200) == 0
    fprintf('Kernel count %d\n', count)
  end

  if max(abs(alpha - old_alpha)) < 1e-8
    break
  end
end
if iter == 2000
  warning('not enough iters')
end
if nargout > 1
  figure(2)
  % e should go up, e2 go down
  plot(1:length(run.e), run.e, 1:length(e2), e2)
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
% SLIDER - add slider to a figure
%
% Usage:
%   >>  slider( handler, horiz, vert, horizmag, vertmag);
%
% Inputs:
%   handler  - figure handler (for the current figure, use gcf)
%   horiz    - [0|1] add a horizontal slider
%   vert     - [0|1] add a horizontal slider
%   horizmag - magnify the width of the figure before adding the slider.
%              Default is 1.
%   vertmag  - magnify the height of the figure before adding the slider.
%              Default is 1.
%   allowsup - [0|1] allow suppression of slider by the 'x' button.
%              Default is 1.
%
% Note:
%   clicking on the 'x' the right corner restores the original setting
%
% Example: figure; plot(1:10); slider(gcf, 1, 1, 2, 2);
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function slider( handler, horiz, vert, horizmag, vertmag, allowsup);

if nargin < 2
	help slider;
	return;
end
if nargin < 3
	vert = 0;
end
if nargin < 4
	horizmag = 1;
end
if nargin < 5
	vertmag = 1;
end
if nargin < 6
	allowsup = 1;
end

pos   = get(gcf, 'position');
width  = 5/pos(3)*400;
height = 5/pos(4)*400;

pos = get(gca,'position'); % plot relative to current axes
q = [0.13 0.11 0 0];
s = [0.0077 0.0081 0.0077 0.0081];
if vert
	h = uicontrol('Parent',handler, ...
 	'style', 'slider', 'Units','Normalized', 'userdata', 1, 'value', 1, ...
	'Position', [113-width -9 width 119].*s+q, 'tag', 'winslider', ...
	'string','vertslider', 'callback', ...
	['if isempty(gcbf), fig = gcf; else fig = gcbf; end;' ...
	 'h = findobj(''parent'', fig);' ...
	 'h2 = findobj(''parent'', fig, ''tag'', ''winslider'');' ...
	 'h = setdiff_bc( h, h2);' ...
	 'curobj = findobj(''parent'', fig, ''string'', ''vertslider'');' ...
	 'shift = get(curobj, ''userdata'') - get(curobj, ''value'');' ...
	 'set( curobj, ''userdata'', get(curobj, ''value''));' ...
	 'for i = 1:length(h),' ...
	 '   curpos = get( h(i), ''position'');' ...
	 '   set( h(i), ''position'', [ curpos(1) curpos(2)+' num2str(vertmag-1) '*shift curpos(3:end)]);' ...
	 'end;' ...
	 'clear h2 h shift i curpos fig curobj;'] );
end

if horiz
	hz = uicontrol('Parent',handler, ...
 	'style', 'slider', 'Units','Normalized', 'userdata', 1, 'value', 1, ...
	'Position', [-17 -19+height 125 height].*s+q, 'tag', 'winslider', ...
	'string','horizslider', 'callback', ...
	['if isempty(gcbf), fig = gcf; else fig = gcbf; end;' ...
	 'h = findobj(''parent'', fig);' ...
	 'h2 = findobj(''parent'', fig, ''tag'', ''winslider'');' ...
	 'h = setdiff_bc( h, h2);' ...
	 'curobj = findobj(''parent'', fig, ''string'', ''horizslider'');' ...
	 'shift = get(curobj, ''userdata'') - get(curobj, ''value'');' ...
	 'set( curobj, ''userdata'', get(curobj, ''value''));' ...
	 'for i = 1:length(h),' ...
	 '   curpos = get( h(i), ''position'');' ...
	 '   set( h(i), ''position'', [ curpos(1)+' num2str(horizmag-1) '*shift curpos(2:end)]);' ...
	 'end;' ...
	 'clear h2 h shift i curpos fig curobj;'] );
end

% button to remove the slider
% ---------------------------
but = uicontrol( 'style', 'pushbutton', 'Units','Normalized', ...
		   'string', 'x', 'position', [113-width -19+height width height].*s+q, 'tag', 'winslider', 'callback', ...
		  ['hx = findobj(''parent'', gcbf, ''tag'', ''winslider'');' ... % put slider to their extremities
		   'hx = setdiff_bc(hx, gcbo);' ...
		   'set(hx, ''value'', 1);' ...
		   'eval(get(hx(1), ''callback''));' ...
		   'if length(hx) >1, eval(get(hx(2), ''callback'')); end;' ...
		   'h = findobj(''parent'', gcbf);' ... % recompute positions
		   'for i = 1:length(h),' ...
		   '   curpos = get( h(i), ''position'');' ...
		   '   set( h(i), ''position'', [(curpos(1)+(' num2str(horizmag) '-1))/' num2str(horizmag) '  (curpos(2)+(' num2str(vertmag) '-1))/' num2str(vertmag) ' curpos(3)/' num2str(horizmag) ' curpos(4)/' num2str(vertmag) ']);' ...
		   'end;' ...
		   'clear h hx curpos;' ...
		   'delete( findobj(''parent'', gcbf, ''tag'', ''winslider'') );' ]);

if ~allowsup
	set(but, 'enable', 'off');
end

% magnify object in the window
% ----------------------------
h = findobj('parent', handler);
set( h, 'units', 'normalized');
h2 = findobj('parent', handler, 'tag', 'winslider');
h = setdiff_bc( h, h2);
for i = 1:length(h)
	curpos = get( h(i), 'position');
	set( h(i), 'position', [curpos(1)*horizmag-(horizmag-1)  curpos(2)*vertmag-(vertmag-1) curpos(3)*horizmag curpos(4)*vertmag]);
end

if horiz
	% set the horizontal axis to 0
	% ----------------------------
	set(hz, 'value', 0);
	eval(get(hz, 'callback'));
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef Drawable < handle

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
        surfaces
        parent
    end

    properties (Access = protected)
        pose
    end

    methods
        function obj = Drawable(parent, pose)
           obj.pose = pose;
           obj.parent = parent;
           obj.surfaces = mcodekit.list.dl_list();
        end
    end

    methods (Access = protected)

        function surface = add_surface(obj, geometry, color)
            surface = obj.add_surface_with_depth(geometry, color, 1);
        end

        function surface = add_surface_with_depth(obj, geometry, color, depth)
            surface_g = geometry;
            T = obj.pose.get_transformation_matrix();
            geometry_t = geometry*T';
            geometry_t(:,3) = depth;
            surface_h = patch('Parent', obj.parent, ...
                              'Vertices', geometry_t, ...
                              'Faces', 1:size(geometry,1), ...
                              'FaceColor', 'flat', ...
                              'FaceVertexCData', color);
            surface = simiam.ui.Surface2D(surface_h, surface_g);
            surface.set_surface_depth(depth);
            surface.transform_surface(T);
            obj.surfaces.append_key(surface);
        end

        function draw_surfaces(obj)
            T = obj.pose.get_transformation_matrix();

            token_k = obj.surfaces.head_;
            while(~isempty(token_k))
                token_k.key_.transform_surface(T);
                token_k = token_k.next_;
            end
        end

%         function transform_surfaces(obj, T)
%             i = obj.surfaces.get_iterator();
%             while(i.has_next())
%                 surface = i.next();
%                 surface.transform(T);
%             end
%         end

        function update_pose(obj, pose)
            [x, y, theta] = pose.unpack();
            obj.pose.set_pose([x, y, theta]);
%             obj.transform_surfaces(pose.get_transformation_matrix());
            obj.draw_surfaces();
        end

    end
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
disp('1 end ')
TheImportantValue='';
if (1>2)
	end
disp('2')
% blah blah end FF
disp('3')
disp('4')

a=['a' 'b' 'asd'...
    'c' 'd'];

if (1)
    while(3<2)
        disp('!');
    end
    for i=0:100
        fprintf(num2str(i));
    end
    disp('!!!')
end

disp('5')
disp('6')

author: SPTHvx

repo: SPTH

save_path: github-repos/MATLAB/SPTHvx-SPTH


text:
%% Selecting Grains
%
%%
% In this section we discuss how to select grains by properties. We start
% our discussion by reconstructing the grain structure from a sample EBSD
% data set.

% load sample EBSD data set
mtexdata forsterite silent

% restrict it to a subregion of interest.
ebsd = ebsd(inpolygon(ebsd,[5 2 10 5]*10^3));

% remove all not indexed pixels
ebsd = ebsd('indexed');

% reconstruct grains
[grains, ebsd.grainId] = calcGrains(ebsd,'angle',5*degree);

% smooth them
grains = smooth(grains,5);

% plot the orientation data of the Forsterite phase
plot(ebsd('fo'),ebsd('fo').orientations)

% plot the grain boundary on top of it
hold on
plot(grains.boundary,'lineWidth',2)
hold off

%% Selecting grains by mouse
% The most easiest way to select a grain is by using the mouse and the
% command <grain2d.selectInteractive.html selectInteractive> which allows
% you to select an arbitrary amount of grains. The index of the selected
% grains appear as the global variable |indSelected| in your workspace

selectInteractive(grains,'lineColor','gold')

% this simulates a mouse click
pause(0.1)
simulateClick(9000,3500)
pause(0.1)

global indSelected;
grains(indSelected)

hold on
plot(grains(indSelected).boundary,'lineWidth',4,'lineColor','gold')
hold off

%% Indexing by orientation or position
% One can also to select a grain by spatial coordinates without user
% interaction. This is done using the syntax |grains(x,y)|, i.e.,

x = 12000; y = 4000;

hold on
plot(grains(x,y).boundary,'linewidth',4,'linecolor','blue')

plot(x,y,'marker','s','markerfacecolor','k',...
  'markersize',10,'markeredgecolor','w')
hold off

%%
% Alternatively one can also select all grains with a certain orientation.
% Lets find all grains with a similar orientation as the one marked in
% gold. As threshold we shall use 20 degree

% select grains by orientation
grains_selected = grains.findByOrientation(grains(indSelected).meanOrientation,20*degree)

hold on
plot(grains_selected.boundary,'linewidth',4,'linecolor','gold')
hold off

%% Indexing by a Property
% In order the generalize the above concept lets remember that the variable
% |grains| is essentially a large vector of grains. Thus when applying a
% function like <grain2d.area.html area> to this variable we obtain a
% vector of the same lenght with numbers representing the area of each
% grain

grain_area = grains.area;

%%
% As a first rather simple application we could colorize the grains
% according to their area, i.e., according to the numbers stored in
% |grain_area|

plot(grains,grain_area)

%%
% As a second application, we can ask for the largest grain within our data
% set. The maximum value and its position within a vector are found by the
% Matlab command |max|.

[max_area,max_id] = max(grain_area)

%%
% The number |max_id| is the position of the grain with a maximum area within
% the variable |grains|. We can access this specific grain by direct
% indexing

grains(max_id)

%%
% and so we can plot it

hold on
plot(grains(max_id).boundary,'linecolor','red','linewidth',4)
hold off

%%
% Note that this way of addressing individual grains can be generalized to
% many grains. E.g. assume we are interested in the largest 5 grains. Then
% we can sort the vector |grain_area| and take the indices of the 5 largest
% grains.

[sorted_area,sorted_id] = sort(grain_area,'descend');

large_grain_id = sorted_id(2:5);

hold on
plot(grains(large_grain_id).boundary,'linecolor','Orange','linewidth',4)
hold off

%% Indexing by a Condition
% By the same syntax as above we can also single out grains that satisfy a
% certain condition. I.e., to access are grains that are at least one
% quarter as large as the largest grain we can do

condition = grain_area > max_area/4;

hold on
plot(grains(condition).boundary,'linecolor','Yellow','linewidth',4)
hold off

%%
% This is a very powerful way of accessing grains as the condition can be
% build up using any grain property. As an example let us consider the
% phase. The phase of the first five grains we get by

grains(1:5).phase

%%
% Now we can access or grains of the first phase Forsterite by the
% condition

condition = grains.phase == 1;
plot(grains(condition))

%%
% To make the above more directly you can use the mineral name for indexing

grains('forsterite')

%%
% Logical indexing allows also for more complex queries, e.g. selecting all
% grains perimeter larger than 6000 and at least 600 measurements within

condition = grains.perimeter>6000 & grains.grainSize >= 600;

selected_grains = grains(condition)

plot(selected_grains)

%% The grainId and how to select EBSD inside specific grains
%
% Besides, the list of grains the command <EBSD.calcGrains.html calcGrains>
% returns also two other output arguments.

plot(grains)
largeGrains = grains(grains.grainSize > 50);

text(largeGrains,largeGrains.id)

%%
% The second output argument grainId is a list with the same size as the
% EBSD measurements that stores for each measurement the corresponding
% grainId. The above syntax stores this list directly inside the ebsd
% variable. This enables MTEX to select EBSD data by grains. The following
% command returns all the EBSD data that belong to grain number 33.

ebsd(grains(33))

%%
% and is equivalent to the command

ebsd(ebsd.grainId == 33)

%%
% The following picture plots the largest grains together with its
% individual orientation measurements.

plot(ebsd(grains(max_id)),ebsd(grains(max_id)).orientations)
hold on
plot(grains(max_id).boundary,'lineWidth',2)
hold off

%% Boundary grains
% Sometimes it is desirable to remove all boundary grains as they might
% distort grain statistics. To do so one should remember that each grain
% boundary has a property |grainId| which stores the ids of the neigbouring
% grains. In the case of an outer grain boundary, one of the neighbouring
% grains has the id zero. We can filter out all these boundary segments by

% ids of the outer boundary segment
outerBoundary_id = any(grains.boundary.grainId==0,2);

% plot the outer boundary segments
plot(grains)
hold on
plot(grains.boundary(outerBoundary_id),'linecolor','red','linewidth',2)
hold off

%%
% Now |grains.boundary(outerBoundary_id).grainId| is a list of grain ids
% where the first column is zero, indicating the outer boundary, and the
% second column contains the id of the boundary grain. Hence, it remains to
% remove all grains with these ids.

% next we compute the corresponding grain_id
grain_id = grains.boundary(outerBoundary_id).grainId;

% remove all zeros
grain_id(grain_id==0) = [];

% and plot the boundary grains
plot(grains(grain_id))

%%
% finally, we could remove the boundary grains by
%
%   grains(grain_id) = []
%
% However, boundary grains can be selected more easily be the command
% |<grain2d/isBoundary.html isBoundary>|.

plot(grains(~grains.isBoundary))

author: mtex-toolbox

repo: mtex

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


text:
close all

% % setup
for nr= 1%1:nRuns
    figure
    plot(xR_true(1,:,nr), xR_true(2,:,nr), 'b-.','Linewidth',1), hold on
    plot(xL_true(1,:,nr), xL_true(2,:,nr), 'b*','Linewidth',1)
    title('Robot trajectory and landmarks')
    legend('True Trajectory','True Landmarks')
    xlabel('x (m)','FontWeight','bold'), ylabel('y (m)','FontWeight','bold')
    axis equal
    print('-dpng','setup_sim')
end

% % robot % %

start= 1;  incr= 5;

% Robot NEES
figure, hold on
plot([start:incr:nSteps],neesR_avg_id(start:incr:end)','g:','Linewidth',3);
plot([start:incr:nSteps],neesR_avg_std(start:incr:end)','b-o','Linewidth',1);
plot([start:incr:nSteps],neesR_avg_fej(start:incr:end)','m-.','Linewidth',2);
plot([start:incr:nSteps],neesR_avg_ocekf_1(start:incr:end)','k--','Linewidth',2);

xlabel('Time (sec)','FontWeight','bold'), ylabel('Robot pose NEES','FontWeight','bold')
legend( 'Ideal-EKF','Std-EKF','FEJ-EKF','OC-EKF')%,'OC-EKF3','Robocentric')
print('-dpng','slam_robot_nees_oc3_sim')

% Robot RMSE
figure
subplot(2,1,1), hold on
plot([start:incr:nSteps],rmsRp_avg_id(start:incr:end)','g:','Linewidth',3);
plot([start:incr:nSteps],rmsRp_avg_std(start:incr:end)','b-o','Linewidth',1);
plot([start:incr:nSteps],rmsRp_avg_fej(start:incr:end)','m-.','Linewidth',2);
plot([start:incr:nSteps],rmsRp_avg_ocekf_1(start:incr:end)','k--','Linewidth',2);
ylabel('Position RMSE (m)','FontWeight','bold')
legend( 'Ideal-EKF','Std-EKF','FEJ-EKF','OC-EKF')%,'OC-EKF3','Robocentric') %#ok<*LEGINTPAR>
subplot(2,1,2), hold on
plot([start:incr:nSteps],rmsRth_avg_id(start:incr:end)','g:','Linewidth',3);
plot([start:incr:nSteps],rmsRth_avg_std(start:incr:end)','b-o','Linewidth',1);
plot([start:incr:nSteps],rmsRth_avg_fej(start:incr:end)','m-.','Linewidth',2);
plot([start:incr:nSteps],rmsRth_avg_ocekf_1(start:incr:end)','k--','Linewidth',2);
xlabel('Time (sec)','FontWeight','bold'),
ylabel('Heading RMSE (rad)','FontWeight','bold')
print('-dpng','slam_robot_rms_oc3_sim')

NEES_Robot = [mean(neesR_avg_id(1:end)),mean(neesR_avg_std(1:end)), mean(neesR_avg_fej(1:end)), mean(neesR_avg_ocekf_1(1:end)) ]
RMS_Position = [mean(rmsRp_avg_id(1:end)),mean(rmsRp_avg_std(1:end)), mean(rmsRp_avg_fej(1:end)),   mean(rmsRp_avg_ocekf_1(1:end)) ]
RMS_Heading = [mean(rmsRth_avg_id(1:end)),mean(rmsRth_avg_std(1:end)), mean(rmsRth_avg_fej(1:end)), mean(rmsRth_avg_ocekf_1(1:end))]

% % landmarks % %

% NEES
figure
bar([neesL_avg_id,...
    neesL_avg_std,...
    neesL_avg_fej,...
    neesL_avg_ocekf_1 ]);
ylabel('Avg. Landmark NEES','FontWeight','bold')
print('-dpng','lm_nees')

% RMSE
figure
bar([rmsL_avg_id,...
    rmsL_avg_std,...
    rmsL_avg_fej,...
    rmsL_avg_ocekf_1  ]);
ylabel('Avg. Landmark RMSE','FontWeight','bold')
print('-dpng','lm_rms')

NEES_L = [neesL_avg_id,neesL_avg_std,neesL_avg_fej, neesL_avg_ocekf_1]
RMS_L = [rmsL_avg_id,rmsL_avg_std,rmsL_avg_fej, rmsL_avg_ocekf_1 ]

author: rpng

repo: ocekf-slam

save_path: github-repos/MATLAB/rpng-ocekf-slam


text:
function I_s = segImage(I,S)

    S = single(S);
    [cx,cy] = gradient(S);
    ccc = (abs(cx)+abs(cy))~=0;
    ccc = uint8(ccc)*255;
    I_s = I;
    I_s(:,:,1) = max(I_s(:,:,1),ccc);
    I_s(:,:,2) = min(I_s(:,:,2),255-ccc);
    I_s(:,:,3) = min(I_s(:,:,3),255-ccc);

author: andrewssobral

repo: lrslibrary


text:
function analyzeMgPipeResults(infoFilePath,resPath,varargin)
% This function takes simulation results generated by mgPipe as the input
% and determines which computed fluxes and reaction abundances are
% significantly different between groups. Also creates violin plots of the
% simulation results. Requires a file with sample information (e.g.,
% disease group, age) for the microbiome models that were generated and
% interrogated through mgPipe.
%
% USAGE:
%
%    analyzeMgPipeResults(infoFilePath,resPath,varargin)
%
% REQUIRED INPUTS:
% infoFilePath:         Path to text file or spreadsheet with information
%                       on analyzed samples including group classification
%                       with sample IDs as rows
% resPath:              char with path of directory where simulation
%                       results are saved
%
% OPTIONAL INPUTS
% statPath:             char with path of directory where results of
%                       statistical analysis are saved
% violinPath:           char with path of directory where violin plots are
%                       saved
% sampleGroupHeaders    list of one or more column headers in file with the
%                       sample information that should be analyzed
%                       (e.g., disease status, age). If not provided, the
%                       second column will be used.
%
% .. Author: Almut Heinken, 12/2020

parser = inputParser();
parser.addRequired('infoFilePath', @ischar);
parser.addRequired('resPath', @ischar);
parser.addParameter('statPath', [pwd filesep 'Statistics'], @ischar);
parser.addParameter('violinPath', [pwd filesep 'ViolinPlots'], @ischar);
parser.addParameter('sampleGroupHeaders', '', @iscellstr);

parser.parse(infoFilePath, resPath, varargin{:});

infoFilePath = parser.Results.infoFilePath;
resPath = parser.Results.resPath;
statPath = parser.Results.statPath;
violinPath = parser.Results.violinPath;
sampleGroupHeaders = parser.Results.sampleGroupHeaders;

% create the folders
mkdir(statPath)
mkdir(violinPath)

% Read in the file with sample information
infoFile = readInputTableForPipeline(infoFilePath);

% get all spreadsheet files in results folder
dInfo = dir(resPath);
fileList={dInfo.name};
fileList=fileList';
fileList(~contains(fileList(:,1),{'.csv','.txt'}))=[];
fileList(contains(fileList(:,1),{'ModelStat'}))=[];
fileList(contains(fileList(:,1),{'GrowthRates'}))=[];

% analyze data in spreadsheets
for i=1:length(fileList)
    sampleData = readInputTableForPipeline([resPath filesep fileList{i}]);

    % merge columns for shadow price results
    if strcmp(sampleData{1,2},'Source')
        for j=2:size(sampleData,1)
            sampleData{j,1}=[sampleData{j,1} '_' sampleData{j,2}];
        end
        sampleData(:,2)=[];
    end
    if strcmp(sampleData{1,3},'Source')
        for j=2:size(sampleData,1)
            sampleData{j,2}=[sampleData{j,2} '_' sampleData{j,3}];
        end
        sampleData(:,2:3)=[];
    end
    if strcmp(sampleData{1,2},'Objective')
        for j=2:size(sampleData,1)
            sampleData{j,1}=[sampleData{j,1} '_' sampleData{j,2}];
        end
        sampleData(:,2)=[];
    end

    sampleData(1,2:end)=strrep(sampleData(1,2:end),'microbiota_model_samp_','');
    sampleData(1,2:end)=strrep(sampleData(1,2:end),'microbiota_model_diet_','');
    sampleData(1,2:end)=strrep(sampleData(1,2:end),'microbiota_model_pDiet_','');
    sampleData(1,2:end)=strrep(sampleData(1,2:end),'microbiota_model_rich_','');
    sampleData(1,2:end)=strrep(sampleData(1,2:end),'host_microbiota_model_samp_','');

    % remove entries not in data
    [C,IA]=intersect(infoFile(2:end,1),sampleData(1,2:end));
    if length(C)<length(sampleData(1,2:end))
        error('Some sample IDs are not found in the file with sample information!')
    end

    % perform statistical analysis
    if ~isempty(sampleGroupHeaders)
        for j=1:length(sampleGroupHeaders)
            [Statistics,significantFeatures] = performStatisticalAnalysis(sampleData,infoFile,'stratification',sampleGroupHeaders{j});

            % Print the results as a text file
            filename = strrep(fileList{i},'.csv','');
            filename = strrep(filename,'.txt','');
            writetable(cell2table(Statistics),[statPath filesep filename '_' sampleGroupHeaders{j} '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
            if size(significantFeatures,2)>1
                writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
            end

            % create violin plots
            currentDir=pwd;
            cd(violinPath)

            % create violin plots for net uptake and secretion files
            if any(contains(fileList{i,1},{'net_uptake_fluxes.csv','net_secretion_fluxes.csv'}))
                makeViolinPlots(sampleData, infoFile, 'stratification',sampleGroupHeaders{j}, 'plottedFeature', filename, 'unit', 'mmol/person/day')
            end
            cd(currentDir)
        end
    else
        [Statistics,significantFeatures] = performStatisticalAnalysis(sampleData,infoFile);

        % Print the results as a text file
        filename = strrep(fileList{i},'.csv','');
        filename = strrep(filename,'.txt','');
        writetable(cell2table(Statistics),[statPath filesep filename '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
        if size(significantFeatures,2)>1
            writetable(cell2table(significantFeatures),[statPath filesep filename '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
        end

        % create violin plots
        currentDir=pwd;
        cd(violinPath)

        % create violin plots for net uptake and secretion files
        if any(strcmp(fileList{i,1},{'net_uptake_fluxes.csv','net_secretion_fluxes.csv','Objectives.txt'}))
            makeViolinPlots(sampleData, infoFile, 'plottedFeature', filename, 'unit', 'mmol/person/day')
        end
        cd(currentDir)
    end

end

author: opencobra

repo: cobratoolbox


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

%
% Rotate around x, y, z (counterclockwise when looking towards the origin)
%

function R = rotate_mat(x, y, z)
% It is assumed that this matrix is muplied from the hind, rotating objects
% along x-axis, y-axis, and z-axis in this order.

Rx = [      1       0       0; ...
            0  cos(x) -sin(x); ...
            0  sin(x)  cos(x)];

Ry = [ cos(y)       0  sin(y); ...
            0       1       0; ...
      -sin(y)       0  cos(y)];

Rz = [ cos(z) -sin(z)       0; ...
       sin(z)  cos(z)       0; ...
            0      0        1];

R = Rz*Ry*Rx;

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Modelling Power Law Absorption Example
%
% This example describes the characteristics of the absorption and
% dispersion encapsulated by the k-Wave simulation functions
%
% For a more detailed discussion of the absorption model used in k-Wave,
% see Treeby, B. E. and Cox, B. T., "Modeling power law absorption and
% dispersion for acoustic propagation using the fractional Laplacian," J.
% Acoust. Soc. Am., vol. 127, no. 5, pp. 2741-2748, 2010.
%
% author: Bradley Treeby
% date: 4th February 2011
% last update: 24th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% modify this parameter to run the different examples
example_number = 1;
% 1: Simulation using kspaceSecondOrder
% 2: Simulation using kspaceFirstOrder1D with default CFL
% 3: Simulation using kspaceFirstOrder1D with CFL = 0.05

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 1024;      % number of grid points in the x (row) direction
x = 12.8e-3;    % grid size in the x direction [m]
dx = x/Nx;      % grid point spacing in the x direction [m]
kgrid = makeGrid(Nx, dx);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]

% define time array
t_end = 4e-6;
if example_number < 3
    [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end);
else
    [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, 0.05, t_end);
end

% create a spatial delta pulse
source_pos = Nx/4;              % [grid points]
source.p0 = zeros(Nx, 1);
source.p0(source_pos) = 1;

% define the sensor positions
source_sensor_dist = 0.5e-3;    % [m]
sensor_sensor_dist = 1e-3;      % [m]
sensor_pos_1 = source_pos + round(source_sensor_dist/dx);
sensor_pos_2 = source_pos + round((source_sensor_dist + sensor_sensor_dist)/dx);

% calculate discrete distance between the sensor positions
d = (sensor_pos_2 - sensor_pos_1)*dx;   % [m]
d_cm = d*100;

% index where the relative dispersion is defined
f_index = 30;

% create a Binary sensor mask
sensor.mask = zeros(Nx, 1);
sensor.mask(sensor_pos_1) = 1;
sensor.mask(sensor_pos_2) = 1;

% preallocate the storage variables
attenuation = zeros(3,  floor(length(kgrid.t_array)/2) + 1);
attenuation_th = zeros(3, floor(length(kgrid.t_array)/2) + 1);
cp = zeros(3, floor(length(kgrid.t_array)/2) + 1);
cp_kk = zeros(3, floor(length(kgrid.t_array)/2) + 1);

for loop = 1:3

    % define the absorption properties of the propagation medium
    switch loop
        case 1
            medium.alpha_coeff = 0.5;
            medium.alpha_power = 1.1;
        case 2
            medium.alpha_coeff = 0.25;
            medium.alpha_power = 1.5;
        case 3
            medium.alpha_coeff = 0.1;
            medium.alpha_power = 1.9;
    end

    % run the simulation without visualisation
    if example_number == 1
        sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'PlotSim', false);
    else
        sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotSim', false);
    end

    % calculate the amplitude and phase spectrum at the two sensor
    % positions
    [f, as1, ps1] = spect(sensor_data(1, :), 1/dt);
    [f, as2, ps2] = spect(sensor_data(2, :), 1/dt);

    % calculate the attenuation from the amplitude spectrums
    attenuation(loop, :) = -20*log10(as2./as1)./d_cm;

    % calculate the corresponding theoretical attenuation in dB/cm
    attenuation_th(loop, :) = medium.alpha_coeff.*(f./1e6).^medium.alpha_power;

    % calculate the dispersion (dependence of the sound speed on frequency)
    % from the phase spectrums
    cp(loop, :) = 2*pi.*f.*d./(unwrap(ps1) - unwrap(ps2));

    % calculate the corresponding theoretical dispersion using the
    % Kramers-Kronig relation for power law absorption
    cp_kk(loop, :) = powerLawKramersKronig(2*pi*f, 2*pi*f(f_index), cp(loop, f_index), db2neper(medium.alpha_coeff, medium.alpha_power), medium.alpha_power);
end

% =========================================================================
% VISUALISATION
% =========================================================================

% set downsampling factor so there is sufficient space between plot markers
ds = 8;

% plot the attenuation
figure;
f_max = 50;
plot(f(1:ds:end)./1e6, attenuation(:, 1:ds:end), 'ko', f./1e6, attenuation_th, 'k-');
set(gca, 'XLim', [0 f_max]);
box on;
xlabel('Frequency [MHz]');
ylabel('\alpha [dB/cm]');

% label the plots
text(40, 160, 'y = 1.9');
text(40, 90, 'y = 1.5');
text(40, 40, 'y = 1.1');

% plot the dispersion
figure
plot(f(1:ds:end)./1e6, cp(:, 1:ds:end), 'ko', f./1e6, cp_kk, 'k-');
set(gca, 'XLim', [0 f_max]);
box on;
xlabel('Frequency [MHz]');
ylabel('C_p [m/s]');

% label the plots
text(40, 1517, 'y = 1.1');
text(40, 1509, 'y = 1.5');
text(40, 1500, 'y = 1.9');

repo: Aero-matlab


text:
function img = colorSegments(cube, method)
% Transforms the segmentation mask into a color image. This  function is
% just intended as a visualization tool.
% It might introduce artifacts if you try to use this function for segmenting objects.

[nrows, ncols, ncolor] = size(cube);
cube = double(reshape(cube, [nrows*ncols ncolor]));

if nargin == 1
    % Sort segments from smaller to larger and visualize with occlusions
    area = squeeze(sum(cube, 1));
    [foo, k] = sort(-area);
    cube = cube(:,k);
end

mask = zeros([nrows*ncols 1]);
for i = ncolor:-1:1
    cube(:,i) = cube(:,i).*(1-mask);
    mask = mask+cube(:,i);
end

map = hsv(ncolor);
map = map(randperm(ncolor),:);

img = cube * map;
img = reshape(img, [nrows ncols 3]);
img(img>1) = 1;
img = uint8(255*img);

author: CSAILVision

repo: LabelMeToolbox


text:
function string = utl_printapproach(app,strip_direct,indent,indent_incr,subset)
% Convert an approach to a string representation
% String = utl_printapproach(Approach)
%
% In:
%   Approach : a BCI approach, either designed in the GUI or constructed in a script
%
%   StripDirect : strip arg_direct flags (default: true)
%
%   Indent : initial indent (default: 0)
%
%   IndentIncrement : indentation increment (default: 4)
%
%   Subset : subset of parameters to display; if 'all' then all parameters are printed,
%            and if 'diff' then only those parameters are printed that differ from the defaults
%            of the respective paradigm (default: 'diff')
%
% Out:
%   String : a string representation of the approach, for use in scripts
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-10-22

% check inputs
if nargin < 2 || isempty(strip_direct)
    strip_direct = true; end
if nargin < 3 || isempty(indent)
    indent = 0; end
if nargin < 4 || isempty(indent_incr)
    indent_incr = 4; end
if nargin < 5 || isempty(subset)
    subset = 'diff'; end
if ~isa(strip_direct,'logical')
    error('The StripDirect argument must be a logical/boolean value.'); end

% get required approach properties
if ischar(app)
    paradigm = ['Paradigm' app];
    parameters = {};
elseif iscell(app)
    paradigm = ['Paradigm' app{1}];
    parameters = app(2:end);
elseif all(isfield(app,{'paradigm','parameters'}))
    paradigm = app.paradigm;
    parameters = app.parameters;
else
    error('The given data structure is not an approach.');
end

% get a handle to the paradigm's calibrate() function
try
    instance = eval(paradigm);
catch e
    if ~exist(char(paradigm),'file')
        error('A BCI paradigm with name %s does not exist.',char(paradigm));
    else
        error('Failed to instantiate the paradigm named %s with error: %s; this is likely an error in the Paradigm''s code.',char(paradigm),e.message);
    end
end
func = @instance.calibrate;

% report the defaults of the paradigm
defaults = clean_fields(arg_report('rich',func));

% report the current settings of the paradigm in form of an argument specification
settings = clean_fields(arg_report('lean',func,parameters));

% get the difference between the defaults and settings
if strcmp(subset,'diff')
    specdiff = arg_diff(defaults,settings);
elseif strcmp(subset,'all')
    specdiff = settings;
else
    error('Unsupported subset: %s',hlp_tostring(subset,100));
end

% convert to nested cell arrays of human-readable name-value pairs
difference = arg_tovals(specdiff,[],'HumanReadableCell',false);

% pre-pend the paradigm choice
paradigm_name = char(paradigm);
difference = [{'arg_selection',paradigm_name(9:end)} difference];

% and convert to string
string = arg_tostring(difference,strip_direct,indent,indent_incr);

% clean fields of x, by removing all arg_direct instances and
% all skippable fields
function x = clean_fields(x)
x(strcmp({x.first_name},'arg_direct') | [x.skippable]) = [];
try
    children = {x.children};
    empty_children = cellfun('isempty',children);
    [x(~empty_children).children] = celldeal(cellfun(@clean_fields,children(~empty_children),'UniformOutput',false));
catch %#ok<CTCH>
end

repo: REST


text:
function g = tensorKernGradient(kern, x, varargin)

% TENSORKERNGRADIENT Gradient of TENSOR kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% tensor product
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO tensorKernParamInit, kernGradient, tensorKernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

% Last of varargin is covGrad.
  g = zeros(1, kern.nParams);
  startVal = 1;
  endVal = 0;
  twoXin = 0;
  if length(varargin) > 1
    twoXin = 1;
    x2 = varargin{1};
    covGrad = varargin{end};
  else
    covGrad = varargin{end};
  end

  for i = 1:length(kern.comp)
    tempKern = tensorKernSlash(kern, i);
    if twoXin
      tempCovGrad = covGrad.*kernCompute(tempKern, x, x2);
    else
      tempCovGrad = covGrad.*kernCompute(tempKern, x);
  end
  endVal = endVal + kern.comp{i}.nParams;
  if ~isempty(kern.comp{i}.index)
    % only part of the data is involved in the kernel.
    if ~twoXin
      g(1, startVal:endVal)  = kernGradient(kern.comp{i}, ...
                                            x(:, kern.comp{i}.index), ...
                                            tempCovGrad);
    else
      g(1, startVal:endVal) = kernGradient(kern.comp{i}, ...
                                           x(:, kern.comp{i}.index), ...
                                           x2(:, kern.comp{i}.index), ...
                                           tempCovGrad);
    end
  else
    if ~twoXin
      % all the data is involved with the kernel.
      g(1, startVal:endVal)  = kernGradient(kern.comp{i}, x, ...
                                            tempCovGrad);
    else
      g(1, startVal:endVal) = kernGradient(kern.comp{i}, x, x2, ...
                                           tempCovGrad);
    end
  end
  startVal = endVal + 1;
end
g = g*kern.paramGroups;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
clear all, close all, clc
load ../../CH01_SVD/DATA/allFaces.mat
mustache = double(rgb2gray(imread('mustache.jpg'))/255);

%% Build Training and Test sets
Training = zeros(size(faces,1),30*20);
Test = zeros(size(faces,1),20*20);
for k=1:20
    baseind = 0;
    if(k>1) baseind = sum(nfaces(1:k-1));
    end
%     inds = baseind + randperm(nfaces(k));
    inds = baseind + (1:nfaces(k));
    Training(:,(k-1)*30+1:k*30) = faces(:,inds(1:30));
    Test(:,(k-1)*20+1:k*20) = faces(:,inds(31:50));
end

%% Oclude test image (image 6 of test set for person 7)
% 126 = 6 + 20*(7-1);
% b1 - clean test image
b1 = Test(:,126);
% b2 - mustache disguise
b2 = Test(:,126).*reshape(mustache,n*m,1);
randvec = randperm(n*m);
first40 = randvec(1:floor(.3*length(randvec)));
vals40 = uint8(255*rand(size(first40)));
b3 = b1;
b3(first40) = vals40;
% imagesc(reshape(b3,n,m)), colormap gray
b4 = b1 + 50*randn(size(b1));
imagesc(reshape(b4,n,m)), colormap gray

B = [b1 b2 b3 b4];
for k=1:4
    axes('position',[0  0  1  1])
    imagesc(reshape(B(:,k),n,m)), colormap gray
    axis off
    set(gcf,'Position',[100 100 168*2 192*2])
    set(gcf,'PaperPositionMode','auto')
%     print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b',num2str(k)]);
end
%% DOWNSAMPLE ALL IMAGES
M = size(Training,2);

TrainingSmall = zeros(120,M);
for k=1:M
    temp = reshape(Training(:,k),n,m);
%     tempSmall = temp(16:16:end,16:16:end);
    tempSmall = imresize(temp,[12 10],'lanczos3');
    TrainingSmall(:,k) = reshape(tempSmall,120,1);
end

Bsmall = zeros(120,4);
for k=1:4
    temp = reshape(B(:,k),n,m);
%     tempSmall = temp(16:16:end,16:16:end);
    tempSmall = imresize(temp,[12 10],'lanczos3');
    Bsmall(:,k) = reshape(tempSmall,120,1);
end

% RENORMALIZE COLUMNS OF TRAINING SMALL

for k=1:M
    scalingFactor(k) = norm(TrainingSmall(:,k));
    TrainingSmall(:,k) = TrainingSmall(:,k)/scalingFactor(k);
end

for k=1:4
    axes('position',[0  0  1  1])
    imagesc(reshape(Bsmall(:,k),12,10)), colormap gray
    axis off
    set(gcf,'Position',[100 100 168 192])
    set(gcf,'PaperPositionMode','auto')
%     print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_bsmall',num2str(k)]);
end

%% L1 SEARCH, TESTCLEAN
b1s = Bsmall(:,1);
eps = .01;
cvx_begin;
    variable a1(M);                 % alph is slope of fit to be optimized
    minimize( norm(a1,1) );
    subject to
        norm(TrainingSmall*a1 - b1s,2) < eps;
cvx_end;

figure
plot(a1,'k')
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRCa1');

figure
axes('position',[0  0  1  1])
imagesc(reshape(Training*(a1./scalingFactor'),n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b1L']);

imagesc(reshape(b1-(Training*(a1./scalingFactor')),n,m)), colormap gray
caxis([0 255])
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b1S']);

figure
binnedError = zeros(20,1);
for k=1:20
    L = (k-1)*30+1:k*30;
    binnedError(k) = norm(b1-(Training(:,L)*(a1(L)./scalingFactor(L)')))/norm(b1)
end
bar(binnedError,'FaceColor',[.7 .7 .7])
hold on
bar(7,binnedError(7),'FaceColor',[.7 0 0])
box on
xlim([0 21])
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRChistogram_b1');

%% L1 SEARCH, MUSTACHE
b2s = Bsmall(:,2);
eps = 500;
cvx_begin;
    variable a2(M);                 % alph is slope of fit to be optimized
    minimize( norm(a2,1) );
    subject to
        norm(TrainingSmall*a2 - b2s,2) < eps;
cvx_end;

figure
plot(a2,'k')
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRCa2');

figure
axes('position',[0  0  1  1])
imagesc(reshape(Training*(a2./scalingFactor'),n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b2L']);

imagesc(reshape(b2-(Training*(a2./scalingFactor')),n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b2S']);

binnedError = zeros(20,1);
for k=1:20
    L = (k-1)*30+1:k*30;
    binnedError(k) = norm(b2-(Training(:,L)*(a2(L)./scalingFactor(L)')))/norm(b2)
end
bar(binnedError,'FaceColor',[.7 .7 .7])
hold on
bar(7,binnedError(7),'FaceColor',[.7 0 0])
box on
xlim([0 21])
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRChistogram_b2');

%% L1 SEARCH, OCCLUSION
b3s = Bsmall(:,3);
eps = 1000;
cvx_begin;
    variable a3(M);                 % alph is slope of fit to be optimized
    minimize( norm(a3,1) );
    subject to
        norm(TrainingSmall*a3 - b3s,2) < eps;
cvx_end;

figure
plot(a3,'k')
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRCa3');

figure
axes('position',[0  0  1  1])
imagesc(reshape(Training*(a3./scalingFactor'),n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b3L']);

imagesc(reshape(b3-(Training*(a3./scalingFactor')),n,m)), colormap gray
caxis([0 255])
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b3S']);

binnedError = zeros(20,1);
for k=1:20
    L = (k-1)*30+1:k*30;
    binnedError(k) = norm(b3-(Training(:,L)*(a3(L)./scalingFactor(L)')))/norm(b3)
end
bar(binnedError,'FaceColor',[.7 .7 .7])
hold on
bar(7,binnedError(7),'FaceColor',[.7 0 0])
box on
xlim([0 21])
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRChistogram_b3');

%% L1 SEARCH, NOISE
b4s = Bsmall(:,4);
eps = 10;
cvx_begin;
    variable a4(M);                 % alph is slope of fit to be optimized
    minimize( norm(a4,1) );
    subject to
        norm(TrainingSmall*a4 - b4s,2) < eps;
cvx_end;

figure
plot(a4,'k')
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRCa4');

figure
axes('position',[0  0  1  1])
imagesc(reshape(Training*(a4./scalingFactor'),n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b4L']);

imagesc(reshape(b4-(Training*(a4./scalingFactor')),n,m)), colormap gray
caxis([0 255])
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b4S']);

binnedError = zeros(20,1);
for k=1:20
    L = (k-1)*30+1:k*30;
    binnedError(k) = norm(b4-(Training(:,L)*(a4(L)./scalingFactor(L)')))/norm(b4)
end
bar(binnedError,'FaceColor',[.7 .7 .7])
hold on
bar(7,binnedError(7),'FaceColor',[.7 0 0])
box on
xlim([0 21])
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRChistogram_b4');

%% LEAST SQUARES IS NO GOOD

a4L2 = pinv(Training)*b4;
figure
plot(a4L2,'k')
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRCa4L2');

figure
axes('position',[0  0  1  1])
imagesc(reshape(Training*a4L2,n,m)), colormap gray
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b4L_L2']);

imagesc(reshape(b4-(Training*a4),n,m)), colormap gray
% caxis([0 255])
axis off
set(gcf,'Position',[100 100 168*2 192*2])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_b4S_L2']);

binnedError = zeros(20,1);
for k=1:20
    L = (k-1)*30+1:k*30;
    binnedError(k) = norm(b4-(Training(:,L)*(a4L2(L))))/norm(b4)
end
bar(binnedError,'FaceColor',[.7 .7 .7])
xlim([0 21])
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/f_chCS_ex06_SRChistogram_b4L2');

%%  PLOT LIBRARY AND TEST SAMPLE
figure
axes('position',[0  0  1  1])
imagesc(TrainingSmall), colormap gray
axis off
set(gcf,'Position',[100 100 1200 240])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_LIBRARYSMALL']);

imagesc((b1s)), colormap gray, axis off
caxis([0 255])
set(gcf,'Position',[100 100 50 900])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_TESTIMAGE']);

axes('position',[0  0  1  1])

imagesc((450-a1')), colormap gray, axis off
caxis([0 550])
set(gcf,'Position',[1 1 1440 50])
set(gcf,'PaperPositionMode','auto')
% print('-dpng','-r0','-loose', ['../figures/f_chCS_EX06SRC_TESTSPARSE']);

author: dynamicslab


text:
function [nstate] = tapas_mh_mc3_hier_node_sample(data, model, ...
    inference, state, node)
%% Samples for a particular node assuming there is a function for it.
%
%
% This works as the conditional independency of the level allows one to sample
% without worrying about the values in other nodes.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

nstate = state;

% Theta is a data object
theta = inference.mh_sampler{node}.propose_sample(data, model, inference, ...
    state, node);

y = state.graph{node - 1};

nllh = model.graph{node - 1}.llh(y, theta, model.graph{node-  1}.htheta);
nlpp = model.graph{node}.llh(theta, state.graph{node + 1}, ...
    model.graph{node}.htheta);

[v] = inference.mh_sampler{node}.ar_rule(state.llh{node - 1}, ...
    state.llh{node}, nllh, nlpp, 0, model.graph{node - 1}.htheta.T);

nstate.graph{node}.y(:, v) = theta.y(:, v);
nstate.llh{node - 1}(:, v) = nllh(:, v);
nstate.llh{node}(:, v) = nlpp(:, v);

nstate.nsample = state.nsample + 1;

end

author: translationalneuromodeling

repo: tapas


text:
function mask = triangle_mask ( dim_num, triangle_order, nodes, coord )

%*****************************************************************************80
%
%% TRIANGLE_MASK is a user routine which masks triangles.
%
%  Discussion:
%
%    The region to be considered is the [0,4]x[0,4] square.
%
%    We want to remove the lower left triangular corner,
%    and part of the upper right triangular corner.
%
%    The following diagram of the 25 nodes indicates by "O" the
%    nodes that should end up being deleted, although the deletion
%    is actually done by triangles.
%
%    Before masking:
%
%      X - X - X - X - X
%      | \ | \ | \ | \ |
%      X - X - X - X - X
%      | \ | \ | \ | \ |
%      X - X - X - X - X
%      | \ | \ | \ | \ |
%      X - X - X - X - X
%      | \ | \ | \ | \ |
%      X - X - X - X - X
%
%    After masking:
%
%      X - X   O   O   O
%      | \ | \
%      X - X - X   O   O
%      | \ | \ | \
%      X - X - X - X - X
%        \ | \ | \ | \ |
%      O   X - X - X - X
%            \ | \ | \ |
%      O   O   X - X - X
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer TRIANGLE_ORDER, the number of nodes in the triangle.
%
%    Input, integer NODES(TRIANGLE_ORDER), the indices of the nodes.
%
%    Input, real COORD(DIM_NUM,TRIANGLE_ORDER), the coordinates
%    of the nodes.
%
%    Output, logical MASK, is TRUE if the triangle should be discarded,
%    and FALSE if the triangle should be retained.
%

%
%  Compute the centroid.
%
  for dim = 1 : dim_num
    centroid(dim) = sum ( coord(dim,1:triangle_order) ) / triangle_order;
  end
%
%  Remove the lower left corner
%
  if ( centroid(1) + centroid(2) < 2.0 )

    mask = 1;
%
%  Remove the upper right section.
%
  elseif ( 5.0 < centroid(1) + centroid(2) & 2.0 < centroid(2) )

    mask = 1;
%
%  Keep everything else.
%
  else

    mask = 0;

  end

  return
end


text:
%LOADPGM	Load a PGM image
%
%	I = loadpgm(filename)
%
%	Returns a matrix containing the image loaded from the PGM format
%	file filename.  Handles ASCII (P2) and binary (P5) PGM file formats.
%
%	If the filename has no extension, and open fails, a '.pgm' will
%	be appended.  If the file cannot be opened it returns [].
%
%	Wildcards are allowed in file names.  If multiple files match
%	a 3D image is returned where the last dimension is the number
%	of images contained.
%
%	I = loadpgm
%
%	Presents a file selection GUI from which the user can pick a file.
%	Uses the same path as previous call.
%
%	Second return argument is the image creation time that comes
% from a TIMESPEC header comment (local convention).

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [I,t] = loadpgm(file)
	t = [];
	fid = fopen(filename, 'r');
	if fid < 0,
		fid = fopen([filename '.pgm'], 'r');
	end
	if fid < 0,
		I = [];
		error(['Couldn''t open file ' filename]);
	end

	white = [' ' 9 10 13];	% space, tab, lf, cr
	white = setstr(white);

	% read the PX header
	magic = fread(fid, 2, 'char');

	%  read cols and process any comment field before hand
	while 1
		c = fread(fid,1,'*char');
		if c == '#',
			comment = fgetl(fid);
			time = sscanf(comment, ' TIMESPEC %d %d');
			if length(time) == 2,
				t = time(1) + time(2)*1e-9;
			end
		elseif ~any(c == white)
			fseek(fid, -1, 'cof');	% unputc()
			break;
		end
	end
	cols = fscanf(fid, '%d', 1);

	%  read rows and process any comment field before hand
	while 1
		c = fread(fid,1,'*char');
		if c == '#',
			fgetl(fid);
		elseif ~any(c == white)
			fseek(fid, -1, 'cof');	% unputc()
			break;
		end
	end
	rows = fscanf(fid, '%d', 1);

	%  read maxval and process any comment field before hand
	while 1
		c = fread(fid,1,'*char');
		if c == '#',
			fgetl(fid);
		elseif ~any(c == white)
			fseek(fid, -1, 'cof');	% unputc()
			break;
		end
	end
	maxval = fscanf(fid, '%d', 1);

	% read the newline
	c = fread(fid,1,'*char');

	% Process the header info and read the image
	if magic(1) == 'P',
		if magic(2) == '2',
			fprintf('%s: ASCII PGM file (%d x %d)\n', filename, rows, cols)
			I = fscanf(fid, '%d', [cols rows])';
		elseif magic(2) == '5',
			fprintf('%s: binary PGM file (%d x %d)\n', filename, rows, cols)
			if maxval == 1,
				fmt = 'unint1';
			elseif maxval == 15,
				fmt = 'uint4';
			elseif maxval == 255,
				fmt = 'uint8';
			elseif maxval == 2^16-1,
				fmt = 'uint16';
			elseif maxval == 2^32-1,
				fmt = 'uint32';
			end
			I = fread(fid, [cols rows], fmt)';
		else
			disp('Not a PGM file');
		end
	end
	fclose(fid);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function plotIMPB(IMCalc)
% JC
% Plot the pb information for IMCalc. IM struct.
% To make sure the coordinates are alright. and pbs are adjacent to each
% other.
%

figure;
xPosV = IMCalc.beams.xPBPosV;
yPosV = IMCalc.beams.yPBPosV;
beamlet_delta_x = IMCalc.beams.beamletDelta_x;
beamlet_delta_y = IMCalc.beams.beamletDelta_y;
w_field = ones(size(xPosV));
%for i=1:length(xPosV)
for i = 250: 251
   patch([xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) + beamlet_delta_x(i)/2 xPosV(i) - beamlet_delta_x(i)/2], [yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) + beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2 yPosV(i) - beamlet_delta_y(i)/2], w_field(i));
   hold on; plot(xPosV(i), yPosV(i), 'r.');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ HessGrad ] = lbfgs_two_loop_recursion( grad, s_array, y_array )
% Two loop recursion algorithm for L-BFGS.
%
% Reference:
%       Jorge Nocedal and Stephen Wright,
%       "Numerical optimization,"
%       Springer Science & Business Media, 2006.
%
%       Algorithm 7.4 in Section 7.2.
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created H.Kasai on Oct. 17, 2016

    if(size(s_array,2)==0)
        HessGrad = -grad;
    else
        q = grad;

        for i = size(s_array,2):-1:1
            rk(i) = 1/(y_array(:,i)'*s_array(:,i));
            a(i) = rk(i)*s_array(:,i)'*q;
            q = q - a(i)*y_array(:,i);
        end

        Hk0 = (s_array(:,end)'*y_array(:,end))/(y_array(:,end)'*y_array(:,end));
        R = Hk0.*q;

        for jj = 1:size(s_array,2)
            beta = rk(jj)*y_array(:,jj)'*R;
            R = R + s_array(:,jj)*(a(jj) - beta);
        end

        HessGrad = -R;
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [ n_data, a, b, fx ] = agm_values ( n_data )

%*****************************************************************************80
%
%% AGM_VALUES returns some values of the AGM.
%
%  Discussion:
%
%    The AGM is defined for nonnegative A and B.
%
%    The AGM of numbers A and B is defined by setting
%
%      A(0) = A,
%      B(0) = B
%
%      A(N+1) = ( A(N) + B(N) ) / 2
%      B(N+1) = sqrt ( A(N) * B(N) )
%
%    The two sequences both converge to AGM(A,B).
%
%    In Mathematica, the AGM can be evaluated by
%
%      ArithmeticGeometricMean [ a, b ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input, 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 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, real A, B, the argument ofs the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 14;

  a_vec = [ ...
     22.0, ...
     83.0, ...
     42.0, ...
     26.0, ...
      4.0, ...
      6.0, ...
     40.0, ...
     80.0, ...
     90.0, ...
      9.0, ...
     53.0, ...
      1.0, ...
      1.0, ...
      1.0, ...
      1.5 ];
  b_vec = [ ...
     96.0, ...
     56.0, ...
      7.0, ...
     11.0, ...
     63.0, ...
     45.0, ...
     75.0, ...
      0.0, ...
     35.0, ...
      1.0, ...
     53.0, ...
      2.0, ...
      4.0, ...
      8.0, ...
      8.0 ];
  fx_vec = [ ...
     52.274641198704240049, ...
     68.836530059858524345, ...
     20.659301196734009322, ...
     17.696854873743648823, ...
     23.867049721753300163, ...
     20.717015982805991662, ...
     56.127842255616681863, ...
      0.000000000000000000, ...
     59.269565081229636528, ...
     3.9362355036495554780, ...
     53.000000000000000000, ...
     1.4567910310469068692, ...
     2.2430285802876025701, ...
     3.6157561775973627487, ...
     4.0816924080221632670 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s=v_lpcconv(from,to,x,y,np)
%V_LPCCONV(from,to,x,y)->s convert between LPC parameter sets
%
% The output is a string that may be passed to eval(s)
% x and y are optionally the input and output matrices
% and np the new value of the order p.
% with one frame stored per row. from and to are taken
% from the following list which also gives the column dimension:
%
%  1 ar p+1  Autoregressive coevfficients: ar(1)=1 always.
%  2 cc  p   Complex cepstral coefficients
%  3 ls  p   Line spectrum pair frequencies (normalized Hz)
%  4 zz  p   Z-plane roots
%  5 ss  p   S-plane roots (normalized Hz)
%  6 rf  p   Reflection coefficients (= -PARCOR coefs)
%  7 ao  p   Area ratios
%  8 aa p+2  Vocal tract areas: aa(p+2)=1 always
%  9 rr p+1  Autocorrelation coefficients
% 10 dl  p   DCT of log area function
% 11 lo  p   Log area ratios
% 12 la p+1  Log areas: la(1)=0 always
% 13 ra p+1  Autocorrelation coefs of inverse filter
% 14 ff p+2  Fourier transform of forward filter (all-pole)
% 15 pf p+2  Power spectrum of forward filter (all-pole)
% 16 gc  p   Gain and cos(w) of each formant
% 17 im p+1  Impulse response of forward filter

%	   Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lpcconv.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nm=['aa';'am';'ao';'ar';'cc';'db';'dl';'ff';'fq';'im';'is';'la';'lo';'ls';'pf';'ra';'rf';'rr';'ss';'zz';];
nx=[...
      0 17  3 17 17 17  7 17  0 17 17 17 17 17 17 17 17 17 17 17;...
      0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;...
     17 17  0 17 17 17 17 17  0 17 17 17 17 17 17 17 17 17 17 17;...
     17 18 17  0  5  6 17  8  0 10 17 17 17 14 15 16 17 18 20 20;...
      4  4  4  4  5  6  4  4  0  4  4  4  4  4 15  4  4  4  4  4;...
     15 15 15 15 15  0 15 15  0 15 15 15 15 15 15 15 15 15 15 15;...
      1  1  1  1  1  1  0  1  0  1  1  1  1  1  1  1  1  1  1  1;...
     15 15 15 15 15 15 15  0  0 15 15 15 15 15 15 15 15 15 15 15;...
     20 20 20 20 20 20 20 20  0 20 20 20 20 20 20 20 20 20 20 20;...
      4  4  4  4  4  4  4  4  0  0  4  4  4  4  4  4  4  4  4  4;...
     17 17 17 17 17 17 17 17  0 17  0 17 17 17 17 17 17 17 17 17;...
     17 17 17 17 17 17 17 17  0 17 17  0 17 17 17 17 17 17 17 17;...
     17 17 17 17 17 17 17 17  0 17 17 17  0 17 17 17 17 17 17 17;...
      4  4  4  4  4  4  4  4  0  4  4  4  4  0  4  4  4  4  4  4;...
      5 18  5  5  5  5  5  5  0  5  5  5  5  5  0  5  5 18  5  5;...
     15 15 15 15 15 15 15 15  0 15 15 15 15 15 15  0 15 15 15 15;...
      1 18  3  4  4  4  1  4  0  4 11 12 13  4  4  4  0 18  4  4;...
      4  2  4  4  4  4  4  4  0  4  4  4  4  4  4  4  4  0  4  4;...
     20 20 20 20 20 20 20 20  0 20 20 20 20 20 20 20 20 20  0 20;...
      4  4  4  4  5  4  4  4  0  4  4  4  4  4  4  4  4  4 19  0;...
];
na=size(nm,1);
b=256*nm(:,1)+nm(:,2);
jf=find(b==256*from(1)+from(2));
jt=find(b==256*to(1)+to(2));
if length([jf jt])~=2
   [x,idx]=sort(b);
  error(sprintf('lpcxx2yy types are: %s',[nm(idx,:)';' '*ones(1,na)]));
end
if nargin<3 s=nm(jf,:); else s=x; end
while jf ~= jt
  jn=nx(jf,jt);
  if jn==0
     error(sprintf('cannot convert between %s and %s',nm(jf,:),nm(jt,:)));
 end
  s=sprintf('lpc%s2%s(%s)',nm(jf,:),nm(jn,:),s);
  jf=jn;
  end
if nargin<4 sn=nm(jt,:); else sn=y; end
s=sprintf('%s=%s;',sn,s);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function x = rstepHess(obj,x0,dx,grad0)
alpha = sqrt(eps);
x = x0 + alpha*dx;
[y,deriv] = obj(x);
g = deriv();
x = (g-grad0)/alpha;

author: nesl