text:
function [ly] = pm2ly(pm)
% Convert length from picometers to light years.
% Never know when you might need this.
% Chad A. Greene 2012
ly = pm*1.05702341e-28;


text:
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return
  end

  column_num = s_word_count ( line );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function v = legendre_van ( m, a, b, n, x )

%*****************************************************************************80
%
%% LEGENDRE_VAN returns the LEGENDRE_VAN matrix for [A,B].
%
%  Discussion:
%
%    The LEGENDRE_VAN matrix is the Legendre Vandermonde-like matrix.
%
%    Normally, the Legendre polynomials are defined on -1 <= XI <= +1.
%    Here, we assume the Legendre polynomials have been defined on the
%    interval A <= X <= B, using the mapping
%      XI = ( - ( B - X ) + ( X - A ) ) / ( B - A )
%    so that
%      Lab(A,B;X) = L(XI).
%
%    if ( I = 1 ) then
%      V(1,1:N) = 1
%    else if ( I = 2 ) then
%      V(2,1:N) = XI(1:N)
%    else
%      V(I,1:N) = ( (2*I-1) * XI(1:N) * V(I-1,1:N) - (I-1)*V(I-2,1:N) ) / I
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%
%    Input, real A, B, the limits of the interval.
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, real X(N), the abscissas.
%
%    Output, real V(M,N), the matrix.
%

%
%  Force X to be a row vector.
%
  x = ( x(:) .' );
%
%  Compute the normalized abscissas in [-1,+1].
%
  xi(1,1:n) = ( - 1.0 * ( b - x(1,1:n)     )   ...
                + 1.0 * (     x(1,1:n) - a ) ) ...
              /         ( b            - a );
%
%  Set up the matrix.
%
  v = zeros ( m, n );

  for i = 1 : m

    if ( i == 1 )
      v(i,1:n) = 1.0;
    elseif ( i == 2 )
      v(i,1:n) = xi(1,1:n);
    else
      v(i,1:n) = ( ( 2 * i - 1 ) * xi(1,1:n) .* v(i-1,1:n)   ...
                 + (   - i + 1 ) *              v(i-2,1:n) ) ...
                 / (     i );
    end

  end

  return
end


text:
function y=fdb(x,a,b)
if x>b
    y=1;
elseif x<a
    y=-1;
elseif x=b
    y=1/2;
elseif x=a
    y=-1/2;
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
%
%  Context-Aware Correlation Filters
%
%  Written by Joao F. Henriques, 2014
%  Adapted by Matthias Mueller, 2016
%
%  This function takes care of setting up parameters, loading video
%  information and computing precisions. For the actual tracking code,
%  check out the TRACKER function.
%
%  RUN_TRACKER
%    Without any parameters, will ask you to choose a video,
%    and show the results in an interactive figure.
%    Press 'Esc' to stop the tracker early. You can navigate the
%    video using the scrollbar at the bottom.
%
%  RUN_TRACKER VIDEO
%    Allows you to select a VIDEO by its name. 'all' will run all videos
%    and show average statistics. 'choose' will select one interactively.
%
%  RUN_TRACKER(VIDEO, SHOW_VISUALIZATION, SHOW_PLOTS)
%    Decide whether to show the scrollable figure, and the precision plot.

function [precision, fps] = run_tracker(video,show_visualization, show_plots)

	%path to the videos (you'll be able to choose one with the GUI).
	base_path = 'sequences/';

	%default settings
	if nargin < 1, video = 'Human7'; end
	if nargin < 2, show_visualization = ~strcmp(video, 'all'); end
	if nargin < 3, show_plots = ~strcmp(video, 'all'); end

    %default settings
    kernel.type =  'linear';

    padding = 2;  %extra area surrounding the target
    lambda1 = 1e-4;     %regularization
    lambda2 = 25;
    interp_factor = 0.015;  %linear interpolation factor for adaptation
    output_sigma_factor = 0.1; %spatial bandwidth (proportional to target)

    features.gray = false;
    features.hog = true;
    features.hog_orientations = 9;

    cell_size = 4;

	switch video
	case 'choose',
		%ask the user for the video, then call self with that video name.
		video = choose_video(base_path);
		if ~isempty(video),
			[precision, fps] = run_tracker(video,show_visualization, show_plots);

			if nargout == 0,  %don't output precision as an argument
				clear precision
			end
		end

	case 'all',
		%all videos, call self with each video name.

		%only keep valid directory names
		dirs = dir(base_path);
		videos = {dirs.name};
		videos(strcmp('.', videos) | strcmp('..', videos) | ...
			strcmp('anno', videos) | ~[dirs.isdir]) = [];

		%the 'Jogging' sequence has 2 targets, create one entry for each.
		%we could make this more general if multiple targets per video
		%becomes a common occurence.
		videos(strcmpi('Jogging', videos)) = [];
		videos(end+1:end+2) = {'Jogging.1', 'Jogging.2'};

		all_precisions = zeros(numel(videos),1);  %to compute averages
		all_fps = zeros(numel(videos),1);

		if ~exist('matlabpool', 'file'),
			%no parallel toolbox, use a simple 'for' to iterate
			for k = 1:numel(videos),
				[all_precisions(k), all_fps(k)] = run_tracker(videos{k}, show_visualization, show_plots);
			end
		else
			%evaluate trackers for all videos in parallel
			if matlabpool('size') == 0,
				matlabpool open;
			end
			parfor k = 1:numel(videos),
				[all_precisions(k), all_fps(k)] = run_tracker(videos{k}, show_visualization, show_plots);
			end
		end

		%compute average precision at 20px, and FPS
		mean_precision = mean(all_precisions);
		fps = mean(all_fps);
		fprintf('\nAverage precision (20px):% 1.3f, Average FPS:% 4.2f\n\n', mean_precision, fps)
		if nargout > 0,
			precision = mean_precision;
		end

	case 'benchmark',
		%running in benchmark mode - this is meant to interface easily
		%with the benchmark's code.

		%get information (image file names, initial position, etc) from
		%the benchmark's workspace variables
		seq = evalin('base', 'subS');
		target_sz = seq.init_rect(1,[4,3]);
		pos = seq.init_rect(1,[2,1]) + floor(target_sz/2);
		img_files = seq.s_frames;
		video_path = [];

		%call tracker function with all the relevant parameters
		positions = tracker(video_path, img_files, pos, target_sz, ...
			padding, kernel, lambda, output_sigma_factor, interp_factor, ...
			cell_size, features, false);

		%return results to benchmark, in a workspace variable
		rects = [positions(:,2) - target_sz(2)/2, positions(:,1) - target_sz(1)/2];
		rects(:,3) = target_sz(2);
		rects(:,4) = target_sz(1);
		res.type = 'rect';
		res.res = rects;
		assignin('base', 'res', res);

	otherwise
		%we were given the name of a single video to process.

		%get image file names, initial state, and ground truth for evaluation
		[img_files, pos, target_sz, ground_truth, video_path] = load_video_info(base_path, video);

		%call tracker function with all the relevant parameters
		[positions, time] = tracker(video_path, img_files, pos, target_sz, ...
			padding, kernel, lambda1, lambda2, output_sigma_factor, interp_factor, ...
			cell_size, features, show_visualization);

		%calculate and show precision plot, as well as frames-per-second
		precisions = precision_plot(positions, ground_truth, video, show_plots);
		fps = numel(img_files) / time;

		fprintf('%12s - Precision (20px):% 1.3f, FPS:% 4.2f\n', video, precisions(20), fps)

		if nargout > 0,
			%return precisions at a 20 pixels threshold
			precision = precisions(20);
		end

	end
end

repo: Context-Aware-CF-Tracking


text:
% KALMAN_FORWARD_BACKWARD		Forward Backward Propogation in Information Form
%
%
% Note	:
%
%  M file accompanying my technical note
%
%    A Technique for Painless Derivation of Kalman Filtering Recursions
%
% available from http://www.mbfys.kun.nl/~cemgil/papers/painless-kalman.ps
%

% Uses :

% Change History :
% Date		Time		Prog	Note
% 07-Jun-2001	 2:24 PM	ATC	Created under MATLAB 5.3.1.29215a (R11.1)

% ATC = Ali Taylan Cemgil,
% SNN - University of Nijmegen, Department of Medical Physics and Biophysics
% e-mail : cemgil@mbfys.kun.nl

A = [1 1;0 1];
C = [1 0];
Q = eye(2)*0.01^2;
R = 0.001^2;
mu1 = [0;1];
P1 = 3*Q;

inv_Q = inv(Q);
inv_R = inv(R);

y = [0 1.1 2 2.95 3.78];

T = length(y);
L = size(Q,1);

%%%%% Forward message Passing
h_f = zeros(L, T);
K_f = zeros(L, L, T);
g_f = zeros(1, T);
h_f_pre = zeros(L, T);
K_f_pre = zeros(L, L, T);
g_f_pre = zeros(1, T);

K_f_pre(:, :, 1) = inv(P1);
h_f_pre(:,1) = K_f_pre(:, :, 1)*mu1;
g_f_pre(1) = -0.5*log(det(2*pi*P1)) - 0.5*mu1'*inv(P1)*mu1;

for i=1:T,
  h_f(:,i) = h_f_pre(:,i) + C'*inv_R*y(:,i);
  K_f(:,:,i) = K_f_pre(:,:,i) + C'*inv_R*C;
  g_f(i) = g_f_pre(i) -0.5*log(det(2*pi*R)) - 0.5*y(:,i)'*inv_R*y(:,i);
  if i<T,
    M = inv(A'*inv_Q*A + K_f(:,:,i));
    h_f_pre(:,i+1) = inv_Q*A*M*h_f(:,i);
    K_f_pre(:,:,i+1) = inv_Q - inv_Q*A*M*A'*inv_Q;
    g_f_pre(i+1) = g_f(i) -0.5*log(det(2*pi*Q)) + 0.5*log(det(2*pi*M)) + 0.5*h_f(:,i)'*M*h_f(:,i);
  end;
end

%%% Backward Message Passing
h_b = zeros(L, T);
K_b = zeros(L, L, T);
g_b = zeros(1, T);

h_b_post = zeros(L, T);
K_b_post = zeros(L, L, T);
g_b_post = zeros(1, T);

for i=T:-1:1,
  h_b(:,i) = h_b_post(:,i) + C'*inv_R*y(:,i);
  K_b(:,:,i) = K_b_post(:,:,i) + C'*inv_R*C;
  g_b(i) = g_b_post(i) - 0.5*log(det(2*pi*R)) - 0.5*y(:,i)'*inv_R*y(:,i);
  if i>1,
    M = inv(inv_Q + K_b(:,:,i));
    h_b_post(:,i-1) = A'*inv(Q)*M*h_b(:,i);
    K_b_post(:,:,i-1) = A'*inv_Q*(Q - M)*inv_Q*A;
    g_b_post(i-1) = g_b(i) -0.5*log(det(2*pi*Q)) + 0.5*log(det(2*pi*M)) + 0.5*h_b(:,i)'*M*h_b(:,i);
  end;
end;

%%%% Smoothed Estimates

mu = zeros(size(h_f));
Sig = zeros(size(K_f));
g = zeros(size(g_f));
lalpha = zeros(size(g_f));

for i=1:T,
  Sig(:,:,i) = inv(K_b_post(:,:,i) + K_f(:,:,i));
  mu(:,i) = Sig(:,:,i)*(h_b_post(:,i) + h_f(:,i));
  g(i) = g_b_post(i) + g_f(:,i);
  lalpha(i) = g(i) + 0.5*log(det(2*pi*Sig(:,:,i))) + 0.5*mu(:,i)'*inv(Sig(:,:,i))*mu(:,i);
end;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, space, sp_mul, eigv, eigf, gnum, dofs_ornt, gnum_mul] = ...
              mp_solve_maxwell_eig_mixed1_3d (problem_data, method_data)

warning ('geopdes:obsolete','Function MP_SOLVE_MAXWELL_EIG_MIXED1_3D is obsolete. Using MP_SOLVE_MAXWELL_EIG_MIXED1 instead')

[geometry, msh, space, sp_mul, eigv, eigf, gnum, dofs_ornt, gnum_mul] = mp_solve_maxwell_eig_mixed1 (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function result = circle_annulus ( func, xc, yc, radius1, radius2, nr )

%*****************************************************************************80
%
%% CIRCLE_ANNULUS approximates an integral in a circular annulus.
%
%  Integration region:
%
%    Points (X,Y) such that
%
%      RADIUS1**2 <= ( X - XC )**2 + ( Y - YC )**2 <= RADIUS2**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Peirce,
%    Numerical Integration Over the Planar Annulus,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 5, Issue 2, June 1957, pages 66-73.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied function of two
%    variables which is to be integrated, of the form:
%      function value = func ( x, y )
%
%    Input, real XC, YC, the center of the circle.
%
%    Input, real RADIUS1, RADIUS2, the radii of the circles.
%
%    Input, integer NR, the order of the rule.  This quantity specifies
%    the number of distinct radii to use.  The number of angles used will
%    be 4*NR, for a total of 4*NR**2 points.
%
%    Output, real RESULT, the approximation to the integral.
%

%
%  Choose radial abscissas and weights.
%
  [ ra, rw ] = legendre_set ( nr );
  a = -1.0E+00;
  b = +1.0E+00;
  c = radius1^2;
  d = radius2^2;
  [ ra, rw ] = rule_adjust ( a, b, c, d, nr, ra, rw );
  ra(1:nr) = sqrt ( ra(1:nr) );
  rw(1:nr) = rw(1:nr) / ( radius2^2 - radius1^2 );
%
%  Set angular abscissas and weights.
%
  nt = 4 * nr;

  tw = 1.0E+00 / nt;
%
%  Approximate the integral.
%
  quad = 0.0E+00;
  for i = 1 : nt
    t = 2.0E+00 * pi * ( i - 1 ) / nt;
    for j = 1 : nr
      x = xc + ra(j) * cos ( t );
      y = yc + ra(j) * sin ( t );
      quad = quad + tw * rw(j) * feval ( func, x, y );
    end
  end

  area = circle_annulus_area_2d ( radius1, radius2 );
  result = quad * area;

  return
end


text:
function PDIFF = chisquare_solve(XGUESS,P,V);
%CHISQUARE_SOLVE  Internal function used by CHISQUARE_INV
%
%   PDIFF=chisquare_solve(XGUESS,P,V)  Given XGUESS, a percentile P,
%   and degrees-of-freedom V, return the difference between
%   calculated percentile and P.

% Uses GAMMAINC
%
% Written January 1998 by C. Torrence

% extra factor of V is necessary because X is Normalized
	PGUESS = gammainc(V*XGUESS/2,V/2);  % incomplete Gamma function

	PDIFF = abs(PGUESS - P);            % error in calculated P

	TOL = 1E-4;
	if (PGUESS >= 1-TOL)  % if P is very close to 1 (i.e. a bad guess)
		PDIFF = XGUESS;   % then just assign some big number like XGUESS
	end

author: grinsted

repo: wavelet-coherence


text:
% If a different type of PDM is to be trained, prepare training data here

clear
load('menpo_68_pts_flip.mat');
addpath('../PDM_helpers');

xs = all_pts(1:end/2,:);
ys = all_pts(end/2+1:end,:);
num_imgs = size(xs, 1);

pdmLoc = ['../../models/pdm/pdm_68_aligned_wild.mat'];
% pdmLoc = ['pdm_68_aligned_menpo_v1.mat'];

load(pdmLoc);

pdm = struct;
pdm.M = double(M);
pdm.E = double(E);
pdm.V = double(V);
errs_poss = [];

all_shapes = -200 * ones(num_imgs, 68 * 2);
min_x = 0;
max_x = 0;
min_y = 0;
max_y = 0;

for i=1:num_imgs

    shape2D = cat(2, all_pts(i,:)', all_pts(i+num_imgs,:)');

    M_n = M;
    ind_rem = find(shape2D(:,1) == -1);
    to_keep = setdiff(1:68, ind_rem);

    if(sum(shape2D(:)==-1) > 0)

        hidden = true;
        % which indices to remove
        inds_to_rem = shape2D(:,1) == -1 | shape2D(:,2) == -1;

        shape2D = shape2D(~inds_to_rem,:);

        inds_to_rem = repmat(inds_to_rem, 3, 1);

        M_n = M(~inds_to_rem);

    end

    % To deal with really extreme cases of roll
    M2D = cat(2, M_n(1:end/3), M_n(end/3+1:2*end/3));
    [ A, t, error, alignedShape, s ] = AlignShapesWithScale(M2D, shape2D);
    R = A/s;

    % Transform the shape
    shape2D(:,1) = shape2D(:,1) - t(1);
    shape2D(:,2) = shape2D(:,2) - t(2);

    shape2D = (R' * shape2D')/s;
    shape2D = shape2D';

%     all_pts(i,all_pts(i,:)~=-1) = shape2D(:,1);
%     all_pts(i+num_pts,all_pts(i+num_pts,:)~=-1) = shape2D(:,2);
%
    all_shapes(i,[to_keep,to_keep+68]) = shape2D(:);

    if(min(shape2D(:,1)) < min_x)
        min_x = min(shape2D(:,1));
    end
    if(min(shape2D(:,2)) < min_y)
        min_y = min(shape2D(:,2));
    end
    if(max(shape2D(:,1)) > max_x)
        max_x = max(shape2D(:,1));
    end
    if(max(shape2D(:,2)) > max_y)
        max_y = max(shape2D(:,2));
    end
end

% Scaling and shifting the model so that all points are from -1 to 1
MD = all_shapes == -200;

xs = 2 * (all_shapes(:,1:68) - min_x)/(max_x-min_x) - 1;
ys = 2 * (all_shapes(:,69:end) - min_y)/(max_y - min_y) - 1;

all_shapes = cat(2, xs, ys);
all_shapes(MD) = -2;

xs = all_shapes(:,1:68);
ys = all_shapes(:,69:end);

save('menpo_train.mat', 'all_shapes');

%%
load('menpo_68_pts_valid.mat');

xs = all_pts(1:end/2,:);
ys = all_pts(end/2+1:end,:);
num_imgs = size(xs, 1);

all_shapes = -200 * ones(num_imgs, 68 * 2);

for i=1:num_imgs

    shape2D = cat(2, all_pts(i,:)', all_pts(i+num_imgs,:)');

    M_n = M;
    ind_rem = find(shape2D(:,1) == -1);
    to_keep = setdiff(1:68, ind_rem);

    if(sum(shape2D(:)==-1) > 0)

        hidden = true;
        % which indices to remove
        inds_to_rem = shape2D(:,1) == -1 | shape2D(:,2) == -1;

        shape2D = shape2D(~inds_to_rem,:);

        inds_to_rem = repmat(inds_to_rem, 3, 1);

        M_n = M(~inds_to_rem);

    end

    % To deal with really extreme cases of roll
    M2D = cat(2, M_n(1:end/3), M_n(end/3+1:2*end/3));
    [ A, t, error, alignedShape, s ] = AlignShapesWithScale(M2D, shape2D);
    R = A/s;

    % Transform the shape
    shape2D(:,1) = shape2D(:,1) - t(1);
    shape2D(:,2) = shape2D(:,2) - t(2);

    shape2D = (R' * shape2D')/s;
    shape2D = shape2D';

%     all_pts(i,all_pts(i,:)~=-1) = shape2D(:,1);
%     all_pts(i+num_pts,all_pts(i+num_pts,:)~=-1) = shape2D(:,2);
%
    all_shapes(i,[to_keep,to_keep+68]) = shape2D(:);

end

% Scaling and shifting the model so that all points are from -1 to 1
MD = all_shapes == -200;

xs = 2 * (all_shapes(:,1:68) - min_x)/(max_x-min_x) - 1;
ys = 2 * (all_shapes(:,69:end) - min_y)/(max_y - min_y) - 1;

all_shapes = cat(2, xs, ys);
all_shapes(MD) = -2;

xs = all_shapes(:,1:68);
ys = all_shapes(:,69:end);

save('menpo_valid.mat', 'all_shapes');

author: TadasBaltrusaitis

repo: OpenFace


text:
function F = ctranspose( F )
%CTRANSPOSE   Conjugate transpose of a DISKFUNV.
%    F = CTRANSPOSE(F) is the same as F'. The orientation of the DISKFUNV
%    is transposed with complex conjugation.
%
%    Since all DISKFUNV objects are real-valued, this command is the same
%    as TRANSPOSE(F)
%
% See also DISKFUNV/CTRANSPOSE

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

F = transpose( F );
F = conj( F );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nodes, edges, vSize] = readGrayGraph(filename,hemi)
% Read a gray graph written out by mrGray
%
%    [nodes, edges, vSize] = readGrayGraph(filename,[hemisphere='l'])
%
% AUTHOR:  Teo
% DATE:    1996
% ARGUMENT:
%   name of gray graph file written by mrGray
%
% RETURNS:
% nodes:  8xN array of
%    Nx(x,y,z,num_edges,edge_offset,layer,dist,pqindex).
% edges:  1xM array of node indices.  The edge_offset of
%    each node points into the starting location of its set
%    of edges.
% where N, M are the number of nodes, edges in the graph.
% vSize = size of the original volume of data containing the anatomicals
%
% Modified 8.20.98  WP/BW
%   Added the big-endian flag and the check for existence step.
% This was done for compatibility with the PCs that normally read
% in using little-endian.
%
% DJH, 8/4/99
%   Add 1 to coordinates because MATLAB starts indexing from 1 instead
% of 0.  This fixes a mis-registration problem in mrLoadRet in the Gray
% view.  The only effect that this has on mrUnfold is that the gLocs3d
% are now incremented by one.
% TODO: Return a timestamp that is saved in the .gray file. This will be
%       used to maintain consistency between gray and flat maps.

warning('Obsolete. Vistasoft uses nifti class files, not mrGray');

if ~exist(filename,'file')
  errordlg(sprintf('Can not find file %s\n',filename));
  % error(str);
end

[p f e] =  fileparts(filename);
if(strcmpi(f, 'coords') || strcmpi(e,'.mat'))
    % Load nodes/edges from a coords file
    if(~exist('hemi','var')||isempty(hemi))
        hemi = 'l';
    end
    tmp = load(filename);
    if(lower(hemi(1))=='l')
        nodes = tmp.allLeftNodes;
        edges = tmp.allLeftEdges;
        classFileName = tmp.leftClassFile;
    elseif(lower(hemi(1))=='r')
        nodes = tmp.allRightNodes;
        edges = tmp.allRightEdges;
        classFileName = tmp.rightClassFile;
    else
        error('Unrecognized hemifield: %s. Cannot read the gray graph', hemi);
    end
    % *** TODO: confirm that we don't need to add 1 to nodes and edges!
    if(nargout>2)
        % *** TODO: check to see if vSize should be reordered to match old
        % vAnat ordering.
        ni = niftiRead(classFileName,[]);
        vSize = ni.dim;
    end
else
    % Load an old-style .Gray file

    % mrGray on all machines writes out using ieee big endian
    % format.  So, we must always open our files that way.
    %
    fid = fopen(filename,'r','b');

    % Read xsize, ysize, zsize.
    [sizes, cnt] = fread(fid, [1 3], 'int');

    % [rows, cols, planes]; [ysize, xsize, zsize].
    vSize = [sizes(2), sizes(1), sizes(3)];

    % Read in number of nodes and edges.
    [sizes, cnt] = fread(fid, [1 2], 'int');

    % Read header.
    num_nodes = sizes(1);
    num_edges = sizes(2);

    % Read nodes.
    nodes = fread(fid, [6 num_nodes], 'int');

    % Offset coordinates by 1 because MATLAB indexes arrays
    % starting from 1
    nodes(1:3,:) = nodes(1:3,:) + 1;

    % Offset edge_offset by 1 because MATLAB implements arrays
    % with base 1 *@#!#!
    nodes(5,:) = nodes(5,:) + 1;

    % Add two additional attributes for miscellaneous use.
    % These are used by the shortest path algorithm to store
    % intermediate distance calculations and to mark position
    % in the priority queue.
    nodes = [nodes ; zeros(2,num_nodes)];

    % Read edges.
    edges = fread(fid, [1 num_edges], 'int');

    % Offset all node indices by 1 because MATLAB implements arrays
    % with base 1 *@#!#!
    edges = edges + 1;

    % Insert timestamp stuff here.

    fclose(fid);

end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Fu = computeFu(m,B,I)
    Fu = [  0      0;...
            0      0;...
            0      0;...
          (1/m)  (1/m);...
            0      0;...
          (-B/I) (B/I)];
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function arctan_int_values_test ( )

%*****************************************************************************80
%
%% ARCTAN_INT_VALUES_TEST demonstrates the use of ARCTAN_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ARCTAN_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  ARCTAN_INT_VALUES stores values of \n' );
  fprintf ( 1, '  the Arctangent integral function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = arctan_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function optS = setOptsExe(opt_str,optS)
%function optS = setOptsExe(opt_str,optS)
%Parse a command that sets the field of a structure, and return the
%structure.
%LM:  J.O.Deasy, 3 Sept, 2002.
%
%copyright (c) 2001, J.O. Deasy and Washington University in St. Louis.
%Use is granted for non-commercial and non-clinical applications.
%No warranty is expressed or implied for any use whatever.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if nargin == 1
  optS = [];
end

opt_str = deblank2(opt_str);
if ~isempty(strfind(opt_str,'=')) && ~strcmp(opt_str(1),'%')
  if ~strcmp(word(opt_str,1),'function') %skip function line if present
   %Remove any comments:
   ind1 = strfind('%',opt_str);
   if ~isempty(ind1)
     if ind1 == 1
       return
     end
     opt_str = opt_str(1:ind1-1);
   end
   %remove trailing ';'
   ind1 = strfind(';',opt_str);
   if ~isempty(ind1)
     opt_str = opt_str(1:ind1-1);
   end
   %parse the command
   %num = words(opt_str);
   %Get the first wor
   firstW = word(opt_str,1);
   %Find the period:
   ind2 = strfind('.',firstW);
   leaf = firstW(ind2+1:length(firstW));
   %Get setting: everything past the '=' sign
   ind3 = strfind(opt_str,'=');
   value = opt_str(ind3+1:length(opt_str));
   value = deblank2(value);
   %Convert numerical assignments to numbers:

   %catch cell arrays
   if strcmp(value(1),'{')
     %build cell array of strings
     iV = strfind(value,'''');
     value2 = {};
     for i = 1 : length(iV)/2
       j = i * 2 - 1;
       str = value(iV(j)+1:iV(j+1)-1);
       value2 = {value2{:},str};
     end
     value = value2;
   else
     if ~isempty(str2num(value))
       value = str2num(value);
     else  %remove apostrophes
       iV = strfind(value,'''');
       if ~isempty(iV)
         value(iV) = [];
       end
     end
   end

   %Finally, make the assignment:
   optS = setfield(optS, leaf, value);
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function Script_HOG_SVM_train()

% Change to your downloaded location
addpath('C:\liblinear\matlab')
addpath('../training_code/');
addpath('../utilities/');
addpath('../../data extraction/');

%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-9:0.5:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svm_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svm_test = @svm_test_linear;

pca_loc = '../../pca_generation/generic_face_rigid.mat';

%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % load the training and testing data for the current fold
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic(train_recs, devel_recs, au, rest_aus, SEMAINE_dir, hog_data_dir);

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Cross-validate here
    [ best_params, ~ ] = validate_grid_search_no_par(svm_train, svm_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);

    model = svm_train(train_labels, train_samples, best_params);

    [prediction, a, actual_vals] = predict(valid_labels, valid_samples, model);

    % Go from raw data to the prediction
    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;

    name = sprintf('models/AU_%d_static.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    write_lin_svm(name, means, svs, b, pos_lbl, neg_lbl);

    name = sprintf('results_SEMAINE_devel/AU_%d_static.mat', au);

    tp = sum(valid_labels == 1 & prediction == 1);
    fp = sum(valid_labels == 0 & prediction == 1);
    fn = sum(valid_labels == 1 & prediction == 0);
    tn = sum(valid_labels == 0 & prediction == 0);

    precision = tp/(tp+fp);
    recall = tp/(tp+fn);

    f1 = 2 * precision * recall / (precision + recall);

    save(name, 'model', 'f1', 'precision', 'recall', 'best_params', 'valid_labels', 'prediction');

end

end


text:
function [Version, date] = CERRCurrentVersion
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

Version  = '6.0';
date =  '07-Mar-2021';


text:
function [area]=findano2() % autogenerated function wrapper
    % This script file is supposed to find an anomaly
    % and estimate the extent in time and space
    %
    % Stefan Wiemer  11/94
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));

    is3 = [];
    ni = 80;

    [i,j] = find(re3 > 5.5);
    X = reshape(loc(1,:),length(gy),length(gx));
    Y = reshape(loc(2,:),length(gy),length(gx));
    figure(map);
    hold on
    for k = 1:length(i)
        xa0 = X(i(k),j(k));
        ya0 = Y(i(k),j(k));
        plot(xa0,ya0,'xk')
        l=ZG.primeCatalog.epicentralDistanceTo(ya0,xa0);
        [s,is] = sort(l);
        is3 = [is3 ; is(1:ni)];
    end   % for k

    is3 = sort(is3);

    l = [];
    for k = 1:length(is3)-1
        if is3(k) ~= is3(k+1)
            l = [l ; is3(k)];
        end
    end

    l = sort(l);
    ZG.newt2= (a(l,:));
    figure(map);
    hold on
    plot(ZG.newt2.Longitude,ZG.newt2.Latitude,'bo');

    % estimate length of anomaly
    %
    i1 = find(ZG.newt2.Longitude == min(ZG.newt2.Longitude));i1 = max(i1);
    i2 = find(ZG.newt2.Longitude == max(ZG.newt2.Longitude));i2 = max(i2);
    di  = sqrt(((ZG.newt2.Longitude(i1)-ZG.newt2.Longitude(i2))*cosd(ya0)*111).^2 +   ((ZG.newt2.Latitude(i1)-ZG.newt2.Latitude(i2))*111).^2) ;
    li = [ZG.newt2.Longitude(i1) ZG.newt2.Latitude(i1) ; ZG.newt2.Longitude(i2) ZG.newt2.Latitude(i2)];
    plot(li(:,1),li(:,2))

    i1 = find(ZG.newt2.Latitude == min(ZG.newt2.Latitude));i1 = max(i1);
    i2 = find(ZG.newt2.Latitude == max(ZG.newt2.Latitude));i2 = max(i2)
    di2  = sqrt(((ZG.newt2.Longitude(i1)-ZG.newt2.Longitude(i2))*cosd(ya0)*111).^2 +   ((ZG.newt2.Latitude(i1)-ZG.newt2.Latitude(i2))*111).^2) ;

    li = [ZG.newt2(i1,1) ZG.newt2(i1,2) ; ZG.newt2(i2,1) ZG.newt2(i2,2)];
    plot(li(:,1),li(:,2))

    di = max([di; di2])
    area = (di/2)^2 *pi;
end

author: CelsoReyes


text:
function im = text_to_image(str,varargin)
  % TEXT_TO_IMAGE  Create an grayscale image the given text
  %
  % im = text_to_image(str)
  % im = text_to_image(str,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   str  string of text
  %   Optional:
  %     'FontName  followed by name of font. List available fonts using
  %       `/usr/local/bin/convert -list font | sed -n "s/^ *Font: //p"`
  %       {'CourierNewB'}
  %     'PointSize'  followed by height of font (in pixels) {200}
  % Outputs:
  %   im  h by w image of text
  %
  % Example:
  %   im = text_to_image('Zapfino','FontName','Zapfino');imshow(im)
  %

  % Defaults
  fontname = [];
  pointsize = 200;
  filename = 'text_to_image-tmp.png';
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'FontName','PointSize'}, ...
    {'fontname','pointsize'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  if isempty(fontname)
    fontname = 'CourierNewB';
  end

  cmd = sprintf( ...
    ['%s -size %dx%d xc:white -font "%s" -pointsize %d ' ...
     '-fill black -annotate +%d+%d "%s" -trim -bordercolor "#FFF" ' ...
     '-border 0 %s'], ...
   path_to_convert,pointsize*2*numel(str),pointsize*4,fontname,pointsize, ...
   0,2*pointsize,str, ...
   filename);
  [status,result] = system(cmd);
  if status ~= 0
    warning(cmd);
    error(result);
  end
  im = imread(filename);
  delete(filename);
end

author: alecjacobson

repo: gptoolbox


text:
function outputFileName = testSceneModel(imSet, paths, modelFileName)

	imList = getImageSet(imSet);
	dt = load(modelFileName);

	trainingParam = dt.trainingParam;

	parfor j = 1:length(imList),
		f{j} = getSceneFeatures(imList{j}, paths, trainingParam.featureParam);
	end
	F = cat(2, f{:});

	[sc2, pr2, ~] = svmMulticlassTest(dt.model, F);
	scores = pr2;
	rawScores = sc2;

	outp.imList = imList;
	outp.scores = scores;
	outp.rawScores = rawScores;

	outputFileName = fullfile(paths.sceneOutsDir, strcat(sprintf('scene-%s_%s_%s', trainingParam.featureParam.featureStr, trainingParam.fileSuffix, imSet), '.mat'));
	%Save the classifier and the predictions...

	save(outputFileName, '-STRUCT', 'outp');

end
%evalRes{i} = benchmarkScene(outputFileName, mapScene, imSet{i});

author: s-gupta

repo: rgbd

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


text:
function [prob]=mhprob3(cand,lambda,sbar,eps,Finv,n);

% compute the first part of the first exponential term
temp1=exp(-cand)-exp(-lambda);

% compute the second part of the first exponential term
% initiate the summation
temp2=0;
% loop over variables
   for jj=1:n
      % if jj=1 (first variable), the part finv*eps_i,t does not exist
      if jj==1
      temp3=(1/sbar(jj,1))*eps(1,1)^2;
      % if any other variable is considered, the term finv*eps_i,t must be taken into account
      else
      temp3=(1/sbar(jj,1))*(eps(1,1)+Finv{jj,1}'*eps(1:jj-1,1))^2;
      end
   % increment the summation
   temp2=temp2+temp3;
   end
% compute the (log of the) first exponential term
term1=-0.5*temp1*temp2;

% next compute the (log of the) second exponential term
term2=-(n/2)*(cand-lambda);

% compute the sum of the two terms to obtain the log of the acceptance prob
% and exponentiate it to obtain the actual acceptance prob
prob=min(1,exp(term1+term2));

author: european-central-bank

repo: BEAR-toolbox


text:
function [rrho, xxvec, yyvec] = density(x, y, dx, dy, xy0);
% [rho, xvec, yvec] = density(x, y, dx, dy, [x0 x1 y0 y1]);
%
%  2-D probability density plot.
%
%  Input:  Vectors x and y of equal length. If dx and/or dy are omitted, the
%          default is 75 bins in each direction.
%
%          If dx or dy is negative, then the variable is taken as the number of
%          bins rather than a grid resolution.
%
%          The vector containing the limits can be padded with NaNs if only
%          certain limits are desired, e g if x0 and y1 are wanted:
%
%            density(x, y, [.5 nan nan 45])
%
%  Output: The density matrix RHO together with vectors XVEC and YVEC.
%          If no output arguments are specified, DENSITY will plot the density
%          function with the prescribed axes using PCOLOR.
%
%  Requires bin.m. Tested under MatLab 4.2, 5.0, and 5.1.
%
%  See also bin.m for further details about dx, x0, etc, and ffgrid.m.
%

% 1.9.97 Oyvind.Breivik@gfi.uib.no.
%
% Oyvind Breivik
% Department of Geophysics
% University of Bergen
% NORWAY

DX = -75;  % Default grid size

x = x(:);

if nargin < 2
 y = x;
end

y = y(:);

xy = NaN*ones(1,4);

if (nargin < 4)
 xy0 = min(x);
end

if (nargin == 3 & length(dx) > 1)
 xy0 = dx;
 dx = DX;
end

if nargin < 3
 dx = DX;
end

if (nargin == 4 & length(dy) > 1)
 xy0 = dy;
 dy = dx;
end
if nargin < 4
 dy = dx;
end

nxy = length(xy0);
xy(1:nxy) = xy0;

if (isnan(xy(4)))
 xy(4) = max(y);
end
if (isnan(xy(3)))
 xy(3) = min(y);
end
if (isnan(xy(2)))
 xy(2) = max(x);
end
if (isnan(xy(1)))
 xy(1) = min(x);
end
x0 = xy(1); x1 = xy(2); y0 = xy(3); y1 = xy(4);

if (dx < 0)
 dx = (x1 - x0)/abs(dx);
end
if (dy < 0)
 dy = (y1 - y0)/abs(dy);
end

ix = bin(x, dx, x0, x1);
iy = bin(y, dy, y0, y1); % bin data in (x,y)-space

xvec = x0:dx:x1;
yvec = y0:dy:y1;

nx = length(xvec);
ny = length(yvec);

inx = (ix >= 1) & (ix <= nx);
iny = (iy >= 1) & (iy <= ny);
in = (inx & iny);
ix = ix(in); iy = iy(in);
N = length(ix); % how many datapoints are left now?

rho = zeros(length(xvec), length(yvec)) + eps;

for i = 1:N
 rho(ix(i), iy(i)) = rho(ix(i), iy(i)) + 1;
end

rho = rho/(N*dx*dy); % Density is n per dx per dy

rho = rho'; % Get in shape

if nargout == 0
 pcolor(xvec, yvec, rho)
 colorbar
 colormap jet
 xlabel(inputname(1))
 ylabel(inputname(2))
 dum = size(rho');
 str = sprintf('%d data points, grid: %dx%d', N, dum(1)-1, dum(2)-1);
 title(str);
end

if nargout > 0
 rrho = rho;
end

if nargout > 1
 xxvec = xvec;
 yyvec = yvec;
end

author: Sable

repo: mcbench-benchmarks


text:
function display(x)
%DISP Display TT/MPS tensor as a tensor network
%
    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    disp( x, inputname(1));
end

author: NicolasBoumal

repo: manopt


text:
% @author: Maziar Raissi

function params_list = KDV()
% quantile(params_list,[0.025 0.25 0.50 0.75 0.975])
clc; close all;

plt = 1;
save_plt = 0;

addpath ..
addpath ../Utilities
addpath ../Kernels/KDV
addpath ../Utilities/export_fig

function CleanupFun()
    rmpath ..
    rmpath ../Utilities
    rmpath ../Kernels/KDV
    rmpath ../Utilities/export_fig
end

finishup = onCleanup(@() CleanupFun());

rng('default')

set(0,'defaulttextinterpreter','latex')

%% Load Data
load('../Data/kdv.mat', 'usol', 't', 'x')
u_star = real(usol); % 512x201
t_star = t; % 201x1
x_star = x';   % 512x1
N_star = size(x_star,1);
nsteps = size(t_star,1)-1;

if plt ==1
    figure(1);
    plot_surface(t_star, x_star, u_star, '$t$', '$x$', '$u(t,x)$');
    view(3)

    drawnow()
end

%% Setup
noise = 0.00;
u_data = u_star + noise*std(u_star(:))*randn(size(u_star));

N0 = 111;
N1 = 109;
%% Optimize model
params_list = zeros(nsteps,2);
hyp = [log([1.0 1.0]) 0.0 0.0 -4.0];
idx1 = randsample(N_star, N0);
step = 1;
for i = 1:step:nsteps
    dt = t_star(i+step) - t_star(i);

    idx0 = idx1;
    x0 = x_star(idx0,:);
    u0 = u_data(idx0,i);

    idx1 = randsample(N_star,N1);
    x1 = x_star(idx1,:);
    u1 = u_data(idx1,i+step);

    model = HPM(x1, u1, x0, u0, dt, hyp);
    model = model.train(50);

    hyp = model.hyp;
    params_list(i,:) = hyp(3:4);

    [pred_n_star, var_n_star] = model.predict(x_star);
    var_n_star = abs(diag(var_n_star));

    error = norm(pred_n_star - u_star(:,i+step))/norm(u_star(:,i+step));

    fprintf(1,'=========================\n');
    fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
        t_star(i+step), model.NLML, error);

    str = sprintf('%.2f  ', params_list(i,:));
    fprintf('Parameters: %s\n\n', str)

    str = sprintf('%.2f  ', median(params_list(1:step:i,:),1));
    fprintf('Median: %s\n', str)
    fprintf(1,'=========================\n\n');

    if plt == 1
        if ~exist('fig','var')
            fig = figure(2);
        end
        set(fig,'units','normalized','outerposition',[0 0 1 1])
        clf

        subplot(3,1,1);
        tit = sprintf('Time: %.2f\n%d training points', t_star(i), N0);
        plot_data_1D(x_star, u_star(:,i), x0, u0, '$x$', '$u(t,x)$', tit);

        subplot(3,1,2);
        tit = sprintf('Time: %.2f\n%d training points', t_star(i+step), N1);
        plot_data_1D(x_star, u_star(:,i+step), x1, u1, '$x$', '$u(t,x)$', tit);

        subplot(3,1,3);
        plot_prediction_1D(x_star, u_star(:,i+step), pred_n_star, var_n_star, ...
            '$x$', '$u(t,x)$', tit);

        drawnow;
    end

end

if save_plt == 1
    export_fig ./Figures/KDV.png -r300
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
% op_ppmref.m
% Jamie Near, McGill University 2015.
%
% USAGE:
% [out,frqshift]=op_ppmref(in,ppmmin,ppmmax,ppmrefval,dimNum);
%
% DESCRIPTION:
% Search for the peak located between ppmmin and ppmmax, and then give that
% peak a new ppm reference value.
%
% INPUTS:
% in        = input data in matlab structure format.
% ppmmin    = minimum of ppm search range.
% ppmmax    = maximum of ppm search range.
% ppmrefval = new reference ppm value.
% dimNum    = which subspectrum to used for referencing (optional).
%
% OUTPUTS:
% out       = Output dataset following frequency shift.
% frqshift  = Frequency shift applied (in Hz).

function [out,frqshift]=op_ppmref(in,ppmmin,ppmmax,ppmrefval,dimNum);

if in.dims.coils>0
    error('ERROR:  Can not operate on data with multilple coils!  ABORTING!!')
end
if in.dims.averages>0
    error('ERROR:  Can not operate on data with multiple averages!  ABORTING!!');
end
if in.dims.extras>0
    error('ERROR:  Can not operate on data with extras dimension!  ABORTING!!');
end
if in.dims.subSpecs>0
    if nargin<5
        plot(in.ppm,in.specs);
        legend('subspec 1','subspec 2');
        dimNum=input('Input which subspectrum to use for referencing: ');
    end
else
    dimNum=1;
end

%Zeropad the data if it hasn't already been done
if ~in.flags.zeropadded
    in_zp=op_zeropad(in,10);
else
    in_zp=in;
end

%Find the ppm of the maximum peak magnitude within the given range:
ppmindex=find(abs(in_zp.specs(in_zp.ppm>ppmmin & in_zp.ppm<ppmmax,dimNum))==max(abs(in_zp.specs(in_zp.ppm>ppmmin & in_zp.ppm<ppmmax,dimNum))));
ppmrange=in_zp.ppm(in_zp.ppm>ppmmin & in_zp.ppm<ppmmax);
ppmmax=ppmrange(ppmindex);

%Now frequency shift the dataset so that the max peak appears at ppmrefval:
frqshift=(ppmmax-ppmrefval)*in.txfrq/1e6;
out=op_freqshift(in,(ppmmax-ppmrefval)*in.txfrq/1e6);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
% Interface lavel functions for the SLAM toolbox.
%
% Create structures
%   createControls                - Create controls structure Con.
%   createLandmarks               - Create Lmk structure array.
%   createMap                     - Create an empty Map structure.
%   createObservations            - Create Obs structure array.
%   createRobots                  - Create robots structure array.
%   createSensors                 - Create sensors structure array.
%   createSimLmk                  - Create a set of landmarks for simulation.
%   createTime                    - Create Tim structure.
%   initRobots                    - Initialize robots in Map.
%   initSensors                   - Initialize sensors in Map.
%   installSensors                - Install sensors on robots.
%
% Draw graphics
%   drawMapFig                    - Redraw the 3D map figure.
%   drawSenFig                    - Redraw one sensor figure.
%
% Synchronization with Map
%   map2rob                       - Update Rob structure from the Map information.
%   map2sen                       - Update Sen structure from the Map information.
%
% Motion and prediction
%   simMotion                     - Simulated robot motion.
%   motion                        - Robot motion.
%
% Landmark initialization
%   initNewLmk                    - Initialise one landmark.
%   initLmkParams                 - Initialize off-filter landmark parameters.
%
% Landmark correction
%   correctKnownLmks              - Correct known landmarks.
%   simObservation                - Observe simulated landmarks.
%   projectLmk                    - Project landmark estimate into sensor's measurement space.
%   selectLmksToObserve           - Select landmarks to observe.
%   matchFeature                  - Match feature.
%   observationInnovation         - Observation innovation.
%   correctLmk                    - Correct landmark.
%   ekfCorrectLmk                 - Correct landmarks.
%   reparametrizeLmk              - Reparametrize landmark.
%   updateLmkParams               - Update off-filter landmark parameters.
%
% Other
%   deleteLmk                     - Delete landmark.
%   drawObsLin                    - Draw an observed line on the pinHole sensor figure.
%   drawObsPnt                    - Redraw a landmark on the pinHole sensor figure.
%   pluckerSegment                - Segment from Plucker line and endpoint abscissas.
%   createMapFig                  - Create 3D map figure and handles.
%   createSenFig                  - Create sensors figures and handles.
%   drawLmk                       - Draw 3D landmark.
%   drawRawLines                  - Draw raw lines.
%   drawRawPnts                   - Draw raw points.
%   prepareForActiveInit          -
%
% Graph SLAM
%   addKeyFrame                   - Add key frame to trajectory.
%   addKnownLmkFactors            - CORRECTKNOWNLMKS  Correct known landmarks.
%   addNewLmkFactors              - INITNEWLMK  Initialise one landmark.
%   buildAdjacencyMatrix          -
%   checkGraphIntegrity           -
%   composeHmgPnt                 -
%   computeResidual               - Compute the residual.
%   createFactors                 - Create Fac structure array.
%   createFrames                  - Create Frm structure array.
%   errorStateJacobians           - Compute Jacobians for projection onto the manifold.
%   frm2rob                       - Frame to Robot information transfer
%   frmJacobians                  - Compute Jacobians for projection onto the manifold.
%   integrateMotion               - Integrate Robot motion, with covariance.
%   lmkJacobians                  - Compute Jacobians for projection onto the manifold.
%   makeAbsFactor                 - Make absolute factor
%   makeAbsFactorFromMotionFactor -
%   makeMeasFactor                -
%   makeMotionFactor              - Make motion factor
%   resetMotion                   - Reset motion to origin of coordinates
%   retroProjLmk                  - Retro project landmark.
%   rob2frm                       - FRM2ROB Robot to Frame information transfer
%   solveGraph                    - Solve the SLAM problem posed as a graph of states and
%   solveGraphCholesky            - Solves the SLAM graph using Cholesky decomposition.
%   updateKeyFrm                  -
%   updateLmk                     -
%   updateStates                  - Update Frm and Lmk states based on computed error.

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function y = bessels(x, type, varargin)
    lut_table = [];
    lut_signs = [];
    if numel(varargin) >= 2
        lut_table = varargin{1};
        lut_signs = varargin{2};
    end
    if strcmp(type, 'I0')

    elseif strcmp(type, 'logI0')
        idxs = x >= 120;
        idxm = x < 0;
        y = x*0;
        y(idxs) = x(idxs) - 1/2 * log(2*pi*x(idxs));

        idd = (~idxs) & (~idxm);
        if ~isempty(lut_table)
            y(idd) = lut_table(round(x(idd) * 10^lut_signs) + 1);
        else
            y(idd) = log(besseli(0, x(idd)));
        end

    elseif strcmp(type, 'I1I0')
        idxs = x >= 120;
        idxm = x < 0;
        y = x*0;
        y(idxs) = 1;

        idd = (~idxs) & (~idxm);
        if ~isempty(lut_table)
            y(idd) = lut_table(round(x(idd) * 10^lut_signs) + 1);
        else
            y(idd) = besseli(1, x(idd))./besseli(0, x(idd));
        end

    end
end

repo: pTVreg


text:
%  MULTIPLICITY --> Computing the multiplicity and multiplicity structure
%  of a system of nonlinear equations at an isolated zero.
%
%  <Synopsis>
%          m = Multiplicity(f, variables, zero)
%          m = Multiplicity(f, variables, zero, options)
%  [m, D, H] = Multiplicity(f, variables, zero)
%  [m, D, H] = Multiplicity(f, variables, zero, options)
%
%  <Input Parameters>
%    1. f         --> a cell array containing the system of nonlinear
%                     equations as strings, e.g.,
%                         >>   f = { 'x^2 + sin(y) -1',  'x-cos(y)' };
%
%    2. variables --> a cell array containing the unknown variables as
%                     strings, e.g.,
%                         >>   variables = { 'x',  'y' };
%
%    3. zero      --> a vector containing numerical zero (root) of f, e.g.,
%                         >>   zero = [1, 0];
%
%    4. options   --> an optional parameter which includes the configuration
%                     settings:
%                     Display: Set to 2 to have all output (the dual space
%                              and Hilbert function) printed to the screen,
%                              and set to 1 to have depth and Hilbert
%                              function printed to the screen. Otherwise
%                              set the default value 0;
%                         Tol: The threshold for numerical rank-revealing.
%                              Singular values above Tol are counted as
%                              nonzero. The default value is 1e-8;
%                     EqnType: The equation type for MULTIPLICITY,
%                              polynomial system ('Poly') or nonlinear
%                              functions ('Nonlinear'). The default value is
%                              'Nonlinear'. By setting the value to 'Poly',
%                              MULTIPLICITY will transfer the polynomial
%                              system to the matrix representation
%                              internally and speed up the computation.
%                      MaxInt: Maximum multiplicity allowed in the recursive
%                              computation. If a zero is not isolated, the
%                              multiplicity will be infinity.  The code can
%                              be used for identifying a nonisolated zero by
%                              setting MaxInt to a known upper bound (e.g.
%                              the Bezout number). The default value is 1000.
%
%                     All the configuration settings may be changed by using
%                     optset function, i.e.,
%                         >>   options = optset('para', value);
%                     para could be any configuration setting, value is set
%                     to para. See OPTSET for details. Any configuration
%                     setting that is not changed will be set to its default
%                     value.
%
%  <Output Parameters>
%    1. m        --> the multiplicity of f at the zero;
%
%    2. D        --> a basis for the dual space as a cell array with
%                    each component being a matrix in the Matlab form of
%                                   D{i} = [c_1, j_1;
%                                           c_2, j_2;
%                                             ...
%                                           c_n, j_n ];
%                    representing a differential functional
%                            D_i = c_1*d_{j_1} + ... + c_n*d_{j_n}
%                    where d_{j_i}'s are differential monomial functionals
%                    (e.g. For a system of equations with variables {x,y,z}
%                     at the zero x=a, y=b, z=c, the functional d_{[i,j,k]}
%                     applied to any function g is the value of the partial
%                     derivative
%
%                                                         i+j+k
%                                             1          d
%                         d_{[i,j,k]}(g) = -------- * ----------- g(a,b,c)
%                                          i! j! k!     i   j   k
%                                                     dx  dy  dz
%
%                     The dual space is the vector space consists of such
%                     differential functionals that vanish on the nonlinear
%                     system while satisfying the so-called closedness
%                     condition);
%
%    3. H        --> values of the Hilbert function in a vector.
%
%  <Examples>
%   Consider the nonlinear system
%
%               sin(x)*cos(y)-x       = 0
%               sin(y)*sin(x)^2 - y^2 = 0
%
%   at the zero (x, y) = (0, 0), the multiplicity can be computed by the
%   following statements:
%
%     >> f = {'sin(x)*cos(y) - x', 'sin(y)*sin(x)^2 - y^2'};
%     >> variables = {'x', 'y'};
%     >> zero = [0, 0];
%     >> m = Multiplicity(f, variables, zero)
%
%   To create an options structure with Tol = 1e-10, MaxInt = 100:
%     >> options = optset('Tol', 1e-10, 'MaxInt', 100);
%     >>  m = Multiplicity(f, variables, zero, options)
%
%  <Algorithm>
%   This code implements a modified closedness subspace method for
%   multiplicity identification with a newly developed equation-by-equation
%   strategy for improving efficiency.
%
%  <References>
%  [1] An algorithm and software for computing multiplicity structures at
%      zeros of nonlinear systems, W. Hao, A. J. Sommesse and Z. Zeng
%  [2] The closedness subspace method for computing the multiplicity
%      structure of a polynomial system, Z. Zeng.
%
%  See also optset, cell

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Test file for chebtech/conj.m

function pass = test_conj(pref)

if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if (n == 1)
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    % Test a scalar-valued function:
    f = testclass.make(@(x) cos(x) + 1i*sin(x), [], pref);
    g = testclass.make(@(x) cos(x) - 1i*sin(x), [], pref);
    h = conj(f);
    pass(n, 1) = norm(h.coeffs - g.coeffs, inf) < 10*vscale(h)*eps;

    % Test an array-valued function:
    f = testclass.make(@(x) [cos(x) + 1i*sin(x), -exp(1i*x)], [], pref);
    g = testclass.make(@(x) cos(x) - 1i*sin(x), [], pref);
    h = conj(f);
    pass(n, 2) = norm(h.coeffs - [g.coeffs, -(g.coeffs)], inf) < ...
        10*max(vscale(h)*eps);
end

end


text:
function obj=stat_density(obj,varargin)
% geom_point Displays an smooth density estimate of the data in x
%
% Example syntax (default arguments): gramm_object.stat_density('function','pdf','kernel','normal','npoints',100)
% the 'function','kernel', and 'bandwidth' arguments are the
% ones used by the underlying matlab function ksdensity
% 'npoints' is used to set how many x values are used to
% display the density estimates.
% 'extra_x' is used to increase the range of x values over
% which the estimated density function is displayed. Values
% will be extended to the right and to the left by extra_x
% times the range of x data.

p=inputParser;
my_addParameter(p,'bandwidth',-1);
my_addParameter(p,'function','pdf')
my_addParameter(p,'kernel','normal')
my_addParameter(p,'npoints',100)
my_addParameter(p,'extra_x',0.1)
parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_density(dobj,dd,p.Results)});
obj.results.stat_density={};
end

function hndl=my_density(obj,draw_data,params)

if obj.polar.is_polar
    %Make x data modulo 2 pi
    draw_data.x=mod(comb(draw_data.x),2*pi);
    warning('Polar density estimate is probably not proper for circular data, use custom bandwidth');
    %Let's try to make boundaries a bit more proper by
    %repeating values below 0 and above 2 pi
    draw_data.x=[draw_data.x-2*pi;draw_data.x;draw_data.x+2*pi];
    extra_x=0;
    binranges=linspace(0,2*pi,params.npoints);
else
    extra_x=(obj.var_lim.maxx-obj.var_lim.minx)*params.extra_x;
    binranges=linspace(obj.var_lim.minx-extra_x,obj.var_lim.maxx+extra_x,params.npoints);
end

if params.bandwidth>0
    [f,xi] = ksdensity(comb(draw_data.x),binranges,'function',params.function,'bandwidth',params.bandwidth,'kernel',params.kernel);
else
    [f,xi] = ksdensity(comb(draw_data.x),binranges,'function',params.function,'kernel',params.kernel);
end

obj.plot_lim.minx(obj.current_row,obj.current_column)=obj.var_lim.minx-extra_x;
obj.plot_lim.maxx(obj.current_row,obj.current_column)=obj.var_lim.maxx+extra_x;
obj.plot_lim.miny(obj.current_row,obj.current_column)=0;
if obj.firstrun(obj.current_row,obj.current_column)
    obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(f);
    %obj.firstrun(obj.current_row,obj.current_column)=0;
    obj.aes_names.y=[obj.aes_names.x ' ' params.function];
else
    if max(f)>obj.plot_lim.maxy(obj.current_row,obj.current_column)
        obj.plot_lim.maxy(obj.current_row,obj.current_column)=max(f);
    end
end

obj.results.stat_density{obj.result_ind,1}.x=xi;
obj.results.stat_density{obj.result_ind,1}.y=f;

[xi,f]=to_polar(obj,xi,f);
hndl=plot(xi,f,'LineStyle',draw_data.line_style,'Color',draw_data.color,'lineWidth',draw_data.line_size);

obj.results.stat_density{obj.result_ind,1}.handle=hndl;
end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
classdef prtClassLibSvm < prtClass
    % prtClassLibSvm  Support vector machine classifier using LibSvm
    %
    %   CLASSIFIER = prtClassLibSvm returns a SVM Classifier using the
    %   SVM toolbox "LibSvm" which provides a fast interface to training
    %   and testing support vector machines.
    %
    %   Note: requires libSvm, which should be in nfPrt\util\libsvm-mat-2.91-1
    %   On linux, you may need to re-build the LibSVM Binaries.  See the
    %   documentation for LibSvm (link below) for more information.
    %
    %    A prtClassLibSvm object inherits all properties from the abstract class
    %    prtClass. In addition is has the following properties; complete
    %    documentation for these properties can be found here:
    %
    %       http://www.csie.ntu.edu.tw/~cjlin/libsvm/
    %
    %         svmType       - Whether to use a C-SVM, nu-SVM, one-class
    %                        SVM, epsilon-SVR, or nu-SVR
    %         kernelType    - Kernel type to use - linear (0), polynomial (1),
    %                         rbf (2, default), sigmoid (3), or
    %                         user-defined (4) - see below
    %         degree        - Kernel function parameter (some kernels)
    %         gamma         - Kernel function parameter (some kernels)
    %         coef0         - Kernel function parameter (some kernels)
    %         cost          - Cost parameter
    %         nu            - nu parameter (nu-SVM's)
    %         pEpsilon      - Loss function parameter (epsilon-SVMs)
    %         cachesize     - Memory cache in MB (can affect speed,
    %                        computer dependent)
    %         eEpsilon      - Termination tolerance
    %         shrinking     - Use shrinking heuristic?
    %         probabilityEstimates - Output probability estimates?
    %         weight        - Class-specific weight parameter in C-SCM
    %
    %   Default values are:
    %     svmType = 0;
    %     kernelType = 2;
    %     degree = 3;
    %     gamma = nan;
    %     coef0 = 0;
    %     cost = 1;
    %     nu = .5;
    %     pEpsilon = .1;
    %     cachesize = 100;
    %     eEpsilon = 0.001;
    %     shrinking = 1;
    %     probabilityEstimates = 0;
    %     weight = [1 1];
    %
    %     userSpecKernel = [];  %only for kernelType = 4, see below
    %
    %   prtClassLibSvm allows the specification of user-defined kernels by
    %   setting svm.kernelType to 4.  This requires further specification
    %   of svm.userSpecKernel.  svm.userSpecKernel must be either a
    %   function handle, fn(x,y) which outputs a matrix of size
    %   size(x,1) x size(y,1), or userSpecKernel can be a prtKernel object.
    %
    %   For example:
    %     svm.kernelType = 4;
    %     svm.userSpecKernel = @(x,y) (x*y'); % correlation kernel
    %
    %     svm.kernelType = 4;
    %     svm.userSpecKernel = prtKernelHyperbolicTangent;
    %
    %
    %   Additional options can be specified by modifying the field
    %   obj.libSvmOptions using the format found here:
    %   http://www.csie.ntu.edu.tw/~cjlin/libsvm/
    %
    %   More documentation can be found here:
    %   http://www.csie.ntu.edu.tw/~cjlin/papers/libsvm.pdf
    %
    %   Note: the LibSvm will output estimated percent correct values to
    %   the screen during processing; because of the way the PRT trains and
    %   tests, these should be ignored during training and plotting. (To be
    %   fixed)
    %
    %   %Example usage:
    %     TestDataSet = prtDataGenUnimodal;       % Create some test and
    %     TrainingDataSet = prtDataGenUnimodal;   % training data
    %     classifier = prtClassLibSvm;              % Create a classifier
    %     classifier = classifier.train(TrainingDataSet);    % Train
    %     classified = run(classifier, TestDataSet);         % Test
    %     percentCorr = prtScorePercentCorrect(classified,TestDataSet);
    %     subplot(2,1,1);
    %     classifier.plot;
    %     subplot(2,1,2);
    %     [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %     h = plot(pf,pd,'linewidth',3);
    %     title('ROC'); xlabel('Pf'); ylabel('Pd');

    %   Notes for doc:
    %       1) The LIBSVM Matlab interface always outputs something like:
    %   Accuracy = 62.63% (6263/10000) (classification)
    %   when you try and run it on an unlabeled data set (in testing) or
    %   when you plot it.  This is because whoever made the interface
    %   didn't think to give us a way that I can see to turn off that
    %   output. (Actually this has now been disabled.)
    %
    %       2) All the parameters that you set below actually get set for
    %       real when we call libSvmOptionString and libSvmOptionStringTest
    %       in train.  The calling syntax is command-line-esque, but I
    %       think I check for all the right requirements.
    %
    %

    properties (SetAccess=private)
        name = 'LibSVM Support Vector Machine'  % Support Vector Machine
        nameAbbreviation = 'LibSVM'         % SVM
        isNativeMary = false;  % False
    end

    properties

        svmType = 0;
        kernelType = 2;
        degree = 3;
        gamma = nan;
        coef0 = 0;
        cost = 1;
        nu = .5;
        pEpsilon = .1;
        cachesize = 100;
        eEpsilon = 0.001;
        shrinking = 1;
        probabilityEstimates = 0;
        weight = [1 1];

        libSvmOptions = '';

        userSpecKernel = [];
    end
    properties (Hidden = true)
        %Note, the libSvm has different opinions of what classes 1 and 0
        %mean than we do; gain takes care of the difference.
        gain = 1;
        trainedSvm
        libSvmOptionsTest
    end

    methods
        function obj = set.svmType(obj,val)
            assert(isscalar(val) && ismember(val,[0 1 2 3 4]),'svmType must be one of 0,1,2,3,4; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.svmType = val;
        end
        function obj = set.kernelType(obj,val)
            assert(isscalar(val) && ismember(val,[0 1 2 3 4]),'kernelType must be one of 0,1,2,3,4; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.kernelType = val;
        end
        function obj = set.degree(obj,val)
            assert(isscalar(val) && prtUtilIsPositiveInteger(val),'degree must be a positive integer; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.degree = val;
        end
        function obj = set.gamma(obj,val)
            if isnan(val)
                return
            end
            if ischar(val)
                obj.gamma = val;
            else
                assert(isscalar(val) && val > 0,'gamma must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
                obj.gamma = val;
            end
        end
        function obj = set.coef0(obj,val)
            assert(isscalar(val),'coef0 must be a scalar; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.coef0 = val;
        end
        function obj = set.cost(obj,val)
            assert(isscalar(val) && val > 0,'cost must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.cost = val;
        end
        function obj = set.nu(obj,val)
            assert(isscalar(val) && val > 0,'nu must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.nu = val;
        end
        function obj = set.pEpsilon(obj,val)
            assert(isscalar(val) && val > 0,'pEpsilon must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.pEpsilon = val;
        end
        function obj = set.cachesize(obj,val)
            assert(isscalar(val) && val > 0,'cachesize must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.cachesize = val;
        end
        function obj = set.eEpsilon(obj,val)
            assert(isscalar(val) && val > 0,'eEpsilon must be a positive value; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.eEpsilon = val;
        end
        function obj = set.shrinking(obj,val)
            val = double(val);
            assert(isscalar(val) && ismember(val,[0 1]),'shrinking must be one of 0,1; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.shrinking = val;
        end
        function obj = set.probabilityEstimates(obj,val)
            val = double(val);
            assert(isscalar(val) && ismember(val,[0 1]),'probabilityEstimates must be one of 0,1; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.probabilityEstimates = val;
        end
        function obj = set.weight(obj,val)
            assert(numel(val)==2 && all(val > 0),'weight must be a positive 2-element vector; see the instructions at http://www.csie.ntu.edu.tw/~cjlin/libsvm/ for more information');
            obj.weight = val;
        end

        function self = prtClassLibSvm(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)

            % Its ok to have 1 class if its a 1 class classifier
            if dataSet.nClasses ~= 2 && self.svmType ~=2
                error('prt:prtClassLibSvm:UnaryData','prtClassLibSvm requires binary data for training');
            end
            if dataSet.nClasses ~= 1 && self.svmType ==2
                % You can avoid this warning using
                % prtOutlierRemovalFnTargets in your algorithm spec.
                warning('prt:prtClassLibSvm:UnaryData','prtClassLibSvm requires unary data for training when svmType = 2; using max(dataSet.uniqueClasses) as the default single class.');
                dataSet = dataSet.retainClasses(max(dataSet.uniqueClasses));
            end

            training_label_vector = dataSet.getTargetsAsBinaryMatrix;
            training_label_vector = double(training_label_vector(:,end)); %zeros and ones
            training_instance_matrix = dataSet.getObservations;
            self.libSvmOptions = self.libSvmOptionString(dataSet);
            self.libSvmOptionsTest = self.libSvmOptionStringTest(dataSet);

            if self.kernelType == 4
                %self.userSpecKernel = prtKernelRbfNdimensionScale;
                if ~isa(self.userSpecKernel,'prtKernel') && ~isa(self.userSpecKernel,'function_handle')
                    error('prtClassLibSvm:KernelType4','For kernelType = 4 (user-specified kernel), svm.userSpecKernel must be a prtKernel or function handle');
                end
                if isa(self.userSpecKernel,'prtKernel')
                    self.userSpecKernel = self.userSpecKernel.train(dataSet);
                    kernelMat = self.userSpecKernel.run(dataSet);
                    training_instance_matrix = kernelMat.data;
                else
                    training_instance_matrix = self.userSpecKernel(dataSet.X,dataSet.X);
                end
                training_instance_matrix = cat(2,(1:size(training_instance_matrix,1))',training_instance_matrix);
            end
            self.trainedSvm = prtExternal.libsvm.svmtrain(training_label_vector, training_instance_matrix, self.libSvmOptions);

            %Need to figure out whether to flip SVM outputs:
            yOut = runAction(self,dataSet);
            auc = prtScoreAuc(yOut.retainLabeled);

            %libSVM has some weird rules about target names.  the first target type it finds is called H1...
            % We can fix this either above, or here
            if auc < .5
                self.gain = -1;
            end
        end

        function DataSetOut = runAction(self,dataSet)

            testing_label_vector = double(dataSet.getTargets);
            if isempty(testing_label_vector)
                testing_label_vector = zeros(dataSet.nObservations,1);
            end

            if self.kernelType == 4
                if ~isa(self.userSpecKernel,'prtKernel') && ~isa(self.userSpecKernel,'function_handle')
                    error('prtClassLibSvm:KernelType4','For kernelType = 4 (user-specified kernel), svm.userSpecKernel must be a prtKernel or function handle');
                end
                if isa(self.userSpecKernel,'prtKernel')
                    kernelMat = self.userSpecKernel.run(dataSet);
                    testing_instance_matrix = kernelMat.data;
                else
                    %Function handle:
                    testing_instance_matrix = self.userSpecKernel(dataSet.X,self.dataSet.X);
                end
                testing_instance_matrix = cat(2,(1:size(testing_instance_matrix,1))',testing_instance_matrix);
            else
                testing_instance_matrix = dataSet.getObservations;
            end

            [dontNeed, dontNeed, decision_values] = prtExternal.libsvm.svmpredict(testing_label_vector, testing_instance_matrix, self.trainedSvm, self.libSvmOptionsTest); %#ok<ASGLU>

            DataSetOut = dataSet;
            DataSetOut = DataSetOut.setObservations(decision_values*self.gain);
        end

        function dataSetOut = runActionFast(self,x,ds)

            switch self.kernelType
                case 2
                    nFeats = size(x,2);
                    f = full(self.trainedSvm.SVs);
                    gram = prtKernelRbf.kernelFn(f,x,sqrt(nFeats));
                    yOut = ((self.trainedSvm.sv_coef'*gram)-self.trainedSvm.rho)';
                otherwise
                    error('prt:prtClassLibSvm','Unsupported kernel type for runActionFast');
            end
            dataSetOut = yOut(:)*self.gain;
        end

    end
    methods (Hidden = true)

        function str = textSummary(self)
            str = sprintf('%s: \n',class(self));
            str2 = evalc('disp(self.trainedSvm)');
            str = sprintf('%s%s',str,str2);
        end
        function optionString = libSvmOptionString(obj,dataSet)
            if isnan(obj.gamma)
                obj.gamma = 1./dataSet.nFeatures;
            end
            if ischar(obj.gamma)
                optionString = sprintf('-s %d -t %d -d %d -g %s -r %f -c %f -n %f -p %f -m %f -e %f -h %d -b %d -w0 %f -w+1 %f',...
                    obj.svmType,obj.kernelType,obj.degree,obj.gamma,obj.coef0,obj.cost,obj.nu,...
                    obj.pEpsilon,obj.cachesize,obj.eEpsilon,obj.shrinking,obj.probabilityEstimates,obj.weight(1),obj.weight(2));
            else
                optionString = sprintf('-s %d -t %d -d %d -g %f -r %f -c %f -n %f -p %f -m %f -e %f -h %d -b %d -w0 %f -w+1 %f',...
                    obj.svmType,obj.kernelType,obj.degree,obj.gamma,obj.coef0,obj.cost,obj.nu,...
                    obj.pEpsilon,obj.cachesize,obj.eEpsilon,obj.shrinking,obj.probabilityEstimates,obj.weight(1),obj.weight(2));
            end
        end
        function optionString = libSvmOptionStringTest(obj,dataSet) %#ok<INUSD>
            optionString = sprintf('-b %d',obj.probabilityEstimates);
        end
    end

    methods (Hidden)
        function str = exportSimpleText(self) %#ok<MANU>
            titleText = sprintf('%% prtClassLibSvm\n');
            svmSvText = prtUtilMatrixToText(full(self.trainedSvm.SVs),'varName','supportVectors');
            svmSvCoeffText = prtUtilMatrixToText(full(self.trainedSvm.sv_coef),'varName','svCoefficients');
            svmSigmaText = prtUtilMatrixToText(sqrt(size(self.trainedSvm.SVs,2)),'varName','sigma');
            svmRhoText = prtUtilMatrixToText(self.trainedSvm.rho,'varName','rho');
            str = sprintf('%s%s%s%s%s',titleText,svmSvText,svmSvCoeffText,svmSigmaText,svmRhoText);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [strshocks_record]=strshocks(beta_gibbs,D_record,Y,X,n,k,It,Bu,favar)

% first create the call storing the results
strshocks_record=cell(n,1);
Bgibbs=reshape(beta_gibbs,k,n,It-Bu);
Dgibbs=reshape(D_record,n,n,It-Bu);

% recall X and Y from the sampling process in this case
if favar.FAVAR==1
    if isfield(favar,'bvarXY')==1
        bvarXY=1;
        Xgibbs=reshape(favar.X_gibbs,size(X,1),size(X,2),It-Bu);
        Ygibbs=reshape(favar.Y_gibbs,size(Y,1),size(Y,2),It-Bu);
    else
        bvarXY=0;
    end
else
    bvarXY=0;
end

% then loop over iterations
for ii=1:It-Bu

    % recover the VAR coefficients, reshaped for convenience
    B=squeeze(Bgibbs(:,:,ii));

    if bvarXY==1
        X=squeeze(Xgibbs(:,:,ii));
        Y=squeeze(Ygibbs(:,:,ii));
    end

    % obtain the residuals from (XXX)
    EPS=Y-X*B;

    % then recover the structural marix D
    D=squeeze(Dgibbs(:,:,ii));

    % obtain the structural disturbances from (XXX)
    ETA=D\EPS';

    % save in struct_shocks_record
    for jj=1:n
        strshocks_record{jj,1}(ii,:)=ETA(jj,:);
    end

end


text:
function asa103_test ( )

%*****************************************************************************80
%
%% ASA103_TEST tests ASA103.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA103_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA103 library.\n' );

  asa103_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA103_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [sFile, ChannelMat] = in_fopen_eyelink(DataFile)
% IN_FOPEN_EYELINK: Open EyeLink eye tracker recordings (.edf).
%
% USAGE:  [sFile, ChannelMat] = in_fopen_eyelink(DataFile)

% @=============================================================================
% 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.
% =============================================================================@
%
% Authors: Francois Tadel & Martin Voelker, 2015

%% ===== READ HEADER =====
% Windows only
if ~ispc
    error('This format is only supported on Windows systems.');
end
% Read header and events
Trials = edfImport(DataFile, [0 1 1]);
% Get events of interest
Trials = edfExtractInterestingEvents(Trials, '^TRIALID');
Trials = edfExtractMicrosaccades(Trials);
hdr.Headers = [Trials.Header];
% Split filename
[fPath, fBase, fExt] = bst_fileparts(DataFile);
% Cannot handled multiple sampling frequencies
nEpochs = length(hdr.Headers);
allRec = [hdr.Headers.rec];
if (nEpochs > 1) && ~all([allRec.sample_rate] == hdr.Headers(1).rec.sample_rate)
    error('This function cannot handle different sampling frequencies between trials.');
end
% Sampling frequency
sfreq = double(hdr.Headers(1).rec.sample_rate);

%% ===== CREATE CHANNEL FILE =====
% List of data types
chnames = {'flags','px','py','hx','hy','pa','gx','gy','rx','ry','input','buttons','htype','hdata','errors', ...
           'gxvel','gyvel','hxvel','hyvel','rxvel','ryvel','fgxvel','fgyvel','fhxvel','fhyvel','frxvel','fryvel'};
chtypes = {'STATUS','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','STATUS','STATUS','STATUS','STATUS','STATUS', ...
           'EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE','EYE'};
% Convert to channels
hdr.chnames = {};
hdr.chtypes = {};
% hdr.chgain  = [];
for i = 1:length(chnames)
    % Channel is not available
    if ~isfield(Trials(1).Samples, chnames{i})
        continue;
    end
%     % Compute channel maximum (as a deviant from the mean)
%     chsamples = double(Trials(1).Samples.(chnames{i}));
%     champ     = max(abs(chsamples),[],2) - abs(mean(chsamples,2));
%     % Channel gain: 1/amplitude (or 0 if the value is constant)
%     chgain = zeros(size(champ));
%     chgain(champ ~= 0) = 1 ./ champ(champ ~= 0);
    % One channel of data per entry
    if (size(Trials(1).Samples.(chnames{i}),1) == 1)
        hdr.chnames{end+1} = chnames{i};
        hdr.chtypes{end+1} = chtypes{i};
%         hdr.chgain(end+1)  = chgain(1);
    % Two channels of data per entry
    else
        % Left eye only OR both eyes
        if (Trials(1).Header.rec.eye == 1) || (Trials(1).Header.rec.eye == 3)
            hdr.chnames{end+1} = [chnames{i}, '_l'];
            hdr.chtypes{end+1} = chtypes{i};
%             hdr.chgain(end+1)  = chgain(1);
        end
        % Right eye only OR both eyes
        if (Trials(1).Header.rec.eye == 2) || (Trials(1).Header.rec.eye == 3)
            hdr.chnames{end+1} = [chnames{i}, '_r'];
            hdr.chtypes{end+1} = chtypes{i};
%             hdr.chgain(end+1)  = chgain(2);
        end
    end
end
nChannels = length(hdr.chnames);
% Create structure
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'EyeLink channels';
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);
% For each channel
for i = 1:nChannels
    ChannelMat.Channel(i).Name    = hdr.chnames{i};
    ChannelMat.Channel(i).Type    = hdr.chtypes{i};
    ChannelMat.Channel(i).Loc     = [0; 0; 0];
    ChannelMat.Channel(i).Orient  = [];
    ChannelMat.Channel(i).Weight  = 1;
    ChannelMat.Channel(i).Comment = [];
end

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder = 'l';
sFile.filename  = DataFile;
sFile.format    = 'EYELINK';
sFile.device    = 'EyeLink';
sFile.header    = hdr;
sFile.comment   = fBase;
sFile.condition = [];
% Epochs
sFile.epochs = repmat(db_template('epoch'), 1, nEpochs);
for i = 1:nEpochs
    % 1 timestamp = 1 millisecond
    sFile.epochs(i).times   = round((double(hdr.Headers(i).starttime) + [0, double(hdr.Headers(i).duration)-1]) ./ 1000 .* sfreq)  ./ sfreq;
    sFile.epochs(i).label   = sprintf('Trial #%d', i);
    sFile.epochs(i).nAvg    = 1;
    sFile.epochs(i).select  = 1;
    sFile.epochs(i).bad     = 0;
end
% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq   = sfreq;
sFile.prop.times   = [min([sFile.epochs.times]), max([sFile.epochs.times])];
sFile.prop.nAvg    = 1;
% No info about bad channels
sFile.channelflag = ones(nChannels, 1);

%% ===== EVENT MARKERS =====
% Initialize structure of events
events = repmat(db_template('event'), 1, 6);
iFixL = 1;  events(iFixL).label = 'fixation L';   % eye=0
iFixR = 2;  events(iFixR).label = 'fixation R';   % eye=1
iBliL = 3;  events(iBliL).label = 'blink L';
iBliR = 4;  events(iBliR).label = 'blink R';
iSacL = 5;  events(iSacL).label = 'saccade L';
iSacR = 6;  events(iSacR).label = 'saccade R';
iUsa  = 7;  events(iUsa).label = 'microsaccade';
iBut  = 8;  events(iBut).label  = 'button';
iMessList = {}; % for storing of event numbers of message events

% Create events
for i = 1:nEpochs
    % Fixation
    if isfield(Trials, 'Fixations') && ~isempty(Trials(i).Fixations)
        Ltime = unique(double(Trials(i).StartTime + Trials(i).Fixations.sttime(Trials(i).Fixations.eye == 0)) ./ 1000);
        Rtime = unique(double(Trials(i).StartTime + Trials(i).Fixations.sttime(Trials(i).Fixations.eye == 1)) ./ 1000);
        events(iFixL).times  = [events(iFixL).times, Ltime];
        events(iFixR).times  = [events(iFixR).times, Rtime];
        events(iFixL).epochs = [events(iFixL).epochs, repmat(i, 1, length(Ltime))];
        events(iFixR).epochs = [events(iFixR).epochs, repmat(i, 1, length(Rtime))];
    end
    % Blinks
    if isfield(Trials, 'Blinks') && ~isempty(Trials(i).Blinks)
        Ltime = unique(double(Trials(i).StartTime + Trials(i).Blinks.sttime(Trials(i).Blinks.eye == 0)) ./ 1000);
        Rtime = unique(double(Trials(i).StartTime + Trials(i).Blinks.sttime(Trials(i).Blinks.eye == 1)) ./ 1000);
        events(iBliL).times  = [events(iBliL).times, Ltime];
        events(iBliR).times  = [events(iBliR).times, Rtime];
        events(iBliL).epochs = [events(iBliL).epochs, repmat(i, 1, length(Ltime))];
        events(iBliR).epochs = [events(iBliR).epochs, repmat(i, 1, length(Rtime))];
    end
    % Saccade
    if isfield(Trials, 'Saccades') && ~isempty(Trials(i).Saccades)
        Ltime = unique(double(Trials(i).StartTime + Trials(i).Saccades.sttime(Trials(i).Saccades.eye == 0)) ./ 1000);
        Rtime = unique(double(Trials(i).StartTime + Trials(i).Saccades.sttime(Trials(i).Saccades.eye == 1)) ./ 1000);
        events(iSacL).times  = [events(iSacL).times, Ltime];
        events(iSacR).times  = [events(iSacR).times, Rtime];
        events(iSacL).epochs = [events(iSacL).epochs, repmat(i, 1, length(Ltime))];
        events(iSacR).epochs = [events(iSacR).epochs, repmat(i, 1, length(Rtime))];
    end
    % Button
    if isfield(Trials, 'Buttons') && ~isempty(Trials(i).Buttons)
        Time = unique(double(Trials(i).StartTime + Trials(i).Buttons.time) ./ 1000);
        events(iBut).times  = [events(iBut).times, Time];
        events(iBut).epochs = [events(iBut).epochs, repmat(i, 1, length(Time))];
    end
    % Micro-saccade (separate detection)
    if isfield(Trials, 'Microsaccades') && ~isempty(Trials(i).Microsaccades)
        Time = unique(double(Trials(i).Microsaccades.StartTime) ./ 1000);
        events(iUsa).times  = [events(iUsa).times, Time];
        events(iUsa).epochs = [events(iUsa).epochs, repmat(i, 1, length(Time))];
    end

    % Messages send by user (e.g. as sync triggers)
    if isfield(Trials, 'Events') && ~isempty(Trials(i).Events.message)
        messageType = unique(Trials(i).Events.message);
        % Exclude "TRIALID" and empty markers
        for iMess = numel(messageType):-1:1
            if strncmpi(messageType{iMess}, 'TRIALID', 6) || isempty(messageType{iMess})
                messageType(iMess) = [];
            end
        end
        % Store messages into 'event' struct
        for iMess = 1:numel(messageType)
            if ~isempty(find(strcmp([iMessList{:}], messageType{iMess}), 1)) % if 1, the event type has already a slot in 'events'
                iEvent = find(strcmp([iMessList{:}], messageType{iMess}));
            else                                                % new event type, define index for 'events' struct and label
                iMessList{numel(iMessList)+1} = messageType{iMess};
                iEvent = numel(events)+1;
                events(iEvent).label = messageType{iMess};
            end

            Time = double(Trials(i).Events.sttime(ismember(Trials(i).Events.message, messageType{iMess})==1)) ./ 1000;
            events(iEvent).times  = [events(iEvent).times, Time];
            events(iEvent).epochs = [events(iEvent).epochs, repmat(i, 1, length(Time))];
        end
    end
end
% Additional fixes
for i = 1:length(events)
    % Rounding the times to the nearest sample
    if ~isempty(events(i).times)
        events(i).times = round(events(i).times .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    end
    % Add channels and notes fields
    events(i).channels = [];
    events(i).notes    = [];
end

% Import this list
sFile = import_events(sFile, [], events);

author: brainstorm-tools


text:
function p = imag(p);
%IMAG         Imaginary part of (interval) polynomial
%
%   r = imag(p);
%

% written  08/28/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p.c = imag(p.c);
  p = normalize(p);

repo: OpenMAS


text:
function [c, s, active_set] = oasisAR1(y, g, lam, smin, active_set)
%% Infer the most likely discretized spike train underlying an AR(1) fluorescence trace
% Solves the sparse non-negative deconvolution problem
%  min 1/2|c-y|^2 + lam |s|_1 subject to s_t = c_t-g c_{t-1} >=s_min or =0

%% inputs:
%   y:  T*1 vector, One dimensional array containing the fluorescence intensities
%withone entry per time-bin.
% OR %%
% len_active_set*4 matrix, active set

%   g:  scalar, Parameter of the AR(1) process that models the fluorescence ...
%impulse response.
%   lam:  scalar, sparsity penalty parameter lambda.
%   smin: scalar, optional, default 0
%miniumal non-zero activity within each bin (minimal 'spike size').
%   active_set: npool x 4 matrix, warm stared active sets

%% outputs
%   c: T*1 vector, the inferred denoised fluorescence signal at each time-bin.
%   s: T*1 vector, discetized deconvolved neural activity (spikes)
%   active_set: npool x 4 matrix, active sets

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% initialization
y = reshape(y, [], 1);
if isempty(y)
    T = sum(active_set(:,4));
else
T = length(y);
end
if ~exist('g', 'var') || isempty(g)
    g = estimate_time_constant(y);
elseif length(g)>1
    c = zeros(T,1);
    s = zeros(T,1);
    active_set = [];
    return;
end
if ~exist('lam', 'var') || isempty(lam);   lam = 0; end
if ~exist('smin', 'var') || isempty(smin);   smin = 0; end
if ~exist('active_set', 'var') || isempty(active_set)
    len_active_set = T;
    active_set = [y-lam*(1-g),ones(T,1),(1:T)',ones(T,1),(1:T)'-1, (1:T)'+1]; % each row is one pool: (vi, wi, t, l)
    active_set(end, :) = [y(end)-lam,1,T,1,T-1,nan] ;
    active_set(1,5) = nan;
else
    len_active_set = size(active_set,1);
    active_set(:,5) = [nan; (1:len_active_set-1)'];
    active_set(:,6) = [(2:len_active_set)';nan];
end
idx = true(len_active_set,1);

%% run OASIS
ii = 1;
ii_next = active_set(ii,6);
while ~isnan(ii_next)
    % find the active set
    while (~isnan(ii_next)) && (active_set(ii_next,1)/active_set(ii_next,2)...
            >=active_set(ii,1)/active_set(ii,2)*g^(active_set(ii,4))+smin)
        active_set(ii_next,5) = ii;
        ii = ii_next;
        ii_next = active_set(ii,6);
    end

    if isnan(ii_next); break; end

    %% merge pools
    active_set(ii,1) = active_set(ii,1) + active_set(ii_next,1)* (g^(active_set(ii,4)));
    active_set(ii,2) = active_set(ii,2) + active_set(ii_next,2)*(g^(2*active_set(ii,4)));
    active_set(ii,4) = active_set(ii,4) + active_set(ii_next,4);
    active_set(ii,6) = active_set(ii_next,6);
    idx(ii_next) = false;
    ii_next = active_set(ii,6);
    ii_prev = active_set(ii, 5);

    %% backtrack until violations fixed
    while (~isnan(ii_prev)) && (active_set(ii,1)/active_set(ii,2)<...
            max(0, active_set(ii_prev,1)/active_set(ii_prev,2)*g^(active_set(ii_prev,4)))+smin)
        ii_next = ii;
        ii = ii_prev;
        active_set(ii,1) = active_set(ii,1) + active_set(ii_next,1)* (g^(active_set(ii,4)));
        active_set(ii,2) = active_set(ii,2) + active_set(ii_next,2)*(g^(2*active_set(ii,4)));
        active_set(ii,4) = active_set(ii,4) + active_set(ii_next,4);
        active_set(ii,6) = active_set(ii_next,6);
        idx(ii_next) = false;

        ii_prev = active_set(ii, 5);
        ii_next = active_set(ii,6);
    end
end
active_set(~idx, :) = [];
len_active_set = size(active_set,1);

%% construct solution for all t
c = zeros(T, 1);
s = c;
for ii=1:len_active_set
    t0 = active_set(ii,3);
    tau = active_set(ii, 4);
    c(t0:(t0+tau-1)) = max(0,active_set(ii,1)/active_set(ii,2)) * (g.^(0:(tau-1)));
end

s(active_set(2:end,3)) = c(active_set(2:end,3)) - g*c(active_set(2:end,3)-1);

author: zhoupc


text:
function pixels = makebitmap(y, width, height, ThreeDee)
% Copyright (c) 2009, The MathWorks, Inc.
%
% Creates a bitmap image of marix y
% Size of the bitmap is determined by arguments
% width and height, and more (2D or 3D visualization)
% by argument ThreeDee.

persistent bitmapHandle;

if ishandle(bitmapHandle)
    % if figure exists, use that for plotting
    figure(bitmapHandle);
    cla;
else
    % create figure if it doesn't exist
    bitmapHandle = figure('Visible', 'off');
    set(bitmapHandle, 'Color', [1 1 1]);
    set(bitmapHandle, 'PaperPositionMode', 'auto');
end;
set(bitmapHandle, 'Visible', 'off');

% Store figure location
a = get(bitmapHandle, 'Position');
a(3) = width;
a(4) = height;

if ThreeDee
    surf(y);
else
    plot(y);
end;

axis tight;
% Resize figure
set(bitmapHandle, 'Position', a);
pixels = hardcopy(bitmapHandle,'-dOpenGL','-r0');
% Size of the bitmap produced by hardcopy may be larger than that
% specified by the figure width & height.
% Figure is transposed to get the correct bitmap orientation for .ToVector
for i = 3:-1:1,
    pixs(:,:,i) = pixels(1:height,1:width,i)'; %#ok : not growing
end;
pixels = pixs(:); % In order to use .ToVector in C#


text:
function asa007_test02 ( )

%*****************************************************************************80
%
%% TEST02 demonstrates the use of SYMINV.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2008
%
%  Author:
%
%    John Burkardt
%
  n_max = 15;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02:\n' );
  fprintf ( 1, '  SYMINV computes the inverse of a positive\n' );
  fprintf ( 1, '  definite symmetric matrix.\n' );
  fprintf ( 1, '  A compressed storage format is used.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here we look at the Hilbert matrix\n' );
  fprintf ( 1, '  A(I,J) = 1 / ( I + J - 1 )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We expect errors to grow quickly with N.\n' );

  for n = 1 : n_max
%
%  Set A to the Hilbert matrix.
%
    k = 0;
    for i = 1 : n
      for j = 1 : i
        k = k + 1;
        a(k) = 1.0 / ( i + j - 1 );
      end
    end

    [ c, nullty, ifault ] = syminv ( a, n );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Matrix order N = %d\n', n );
    fprintf ( 1, '  Maxtrix nullity NULLTY = %d\n', nullty );

    k = 0;
    for j = 1 : n
      for i = 1 : j - 1
        k = k + 1;
        cfull(i,j) = c(k);
        cfull(j,i) = c(k);
      end
      k = k + 1;
      cfull(j,j) = c(k);
    end

    k = 0;
    for j = 1 : n
      for i = 1 : j - 1
        k = k + 1;
        afull(i,j) = a(k);
        afull(j,i) = a(k);
      end
      k = k + 1;
      afull(j,j) = a(k);
    end
%
%  Compute C * A - I.
%
    b(1:n,1:n) = cfull(1:n,1:n) * afull(1:n,1:n) - eye(n);

    diff = sqrt ( sum ( sum ( b(1:n,1:n).^2 ) ) );

    fprintf ( 1,  '  RMS ( C * A - I ) = %e\n', diff );
  end

  return
end


text:
function kern = diagKernExpandParam(kern, params)

% DIAGKERNEXPANDPARAM Create kernel structure from DIAG kernel's parameters.
% FORMAT
% DESC returns a diagonal noise covariance function kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : diagKernParamInit, diagKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

kern.variance = params(1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% function x = cwt_iw(Wx, type, opt)
%
% The inverse wavelet transform of signal Wx
%
% Implements Eq. (4.67) of [1].
%
% 1. Mallat, S., Wavelet Tour of Signal Processing 3rd ed.
%
% Inputs:
%  Wx: wavelet transform of a signal, see help cwt_fw
%  type: wavelet used to take the wavelet transform,
%        see help cwt_fw and help wfiltfn
%  opt: options structure used for forward wavelet transform.
%
% Output:
%  x: the signal, as reconstructed from Wx
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function x = cwt_iw(Wx, type, opt)
    if nargin<3, opt = struct(); end

    [na, n] = size(Wx);

    [N,n1,n2] = p2up(n);
    Wxp = zeros(na, N);

    % TODO - do we want to pad the wavelet representation here?  Or not
    % cut off Wx in cwt_fw and then use that to reconstruct?
    Wxp(:, n1+1:n1+n) = Wx;
    Wx = Wxp; clear Wxp;

    % Following the same value in cwt_fw
    noct = log2(N)-1;
    nv = na/noct;
    as = 2^(1/nv) .^ (1:1:na);

    assert(mod(noct,1) == 0);
    assert(nv>0 && mod(nv,1)==0); % integer

    % Find the admissibility coefficient Cpsi
    Cpsi = cwt_adm(type, opt);

    x = zeros(1, N);
    for ai=1:na
        a = as(ai);
        Wxa = Wx(ai, :);

        psih = wfilth(type, N, a, opt);

        % Convolution theorem here
        Wxah = fft(Wxa);
        xah = Wxah .* psih;
        xa = ifftshift(ifft(xah));

        x = x + xa/a;
    end

     % Take real part and normalize by log_e(a)/Cpsi
     x = log(2^(1/nv))/Cpsi * real(x);

     % Keep the unpadded part
     x = x(n1+1: n1+n);

end % cwt_iw

author: HeLiangHIT


text:
function perm0_uniform_test ( )

%*****************************************************************************80
%
%% PERM0_UNIFORM_TEST tests PERM0_UNIFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2015
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PERM0_UNIFORM_TEST\n' );
  fprintf ( 1, '  PERM0_UNIFORM randomly selects a permutation of 0, ..., N-1.\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : 5

    [ p, seed ] = perm0_uniform ( n, seed );

    fprintf ( 1, '  ' );
    for i = 1 : n
      fprintf ( 1, '%4d', p(i) );
    end
    fprintf ( 1, '\n' );

  end

  return
end