text:
% ----------------------------------------------------------------------- %
% Function plot_areaerrorbar plots the mean and standard deviation of a   %
% set of data filling the space between the positive and negative mean    %
% error using a semi-transparent background, completely customizable.     %
%                                                                         %
%   Input parameters:                                                     %
%       - data:     Data matrix, with rows corresponding to observations  %
%                   and columns to samples.                               %
%       - options:  (Optional) Struct that contains the customized params.%
%           * options.handle:       Figure handle to plot the result.     %
%           * options.color_area:   RGB color of the filled area.         %
%           * options.color_line:   RGB color of the mean line.           %
%           * options.alpha:        Alpha value for transparency.         %
%           * options.line_width:   Mean line width.                      %
%           * options.x_axis:       X time vector.                        %
%           * options.error:        Type of error to plot (+/-).          %
%                   if 'std',       one standard deviation;               %
%                   if 'sem',       standard error mean;                  %
%                   if 'var',       one variance;                         %
%                   if 'c95',       95% confidence interval.              %
% ----------------------------------------------------------------------- %
%   Example of use:                                                       %
%       data = repmat(sin(1:0.01:2*pi),100,1);                            %
%       data = data + randn(size(data));                                  %
%       plot_areaerrorbar(data);                                          %
% ----------------------------------------------------------------------- %
%   Author:  Victor Martinez-Cagigal                                      %
%   Date:    30/04/2018                                                   %
%   E-mail:  vicmarcag (at) gmail (dot) com                               %
% ----------------------------------------------------------------------- %
function plot_areaerrorbar(data, options)
    % Default options
    if(nargin<2)
        options.handle     = figure(1);
        options.color_area = [128 193 219]./255;    % Blue theme
        options.color_line = [ 52 148 186]./255;
        %options.color_area = [243 169 114]./255;    % Orange theme
        %options.color_line = [236 112  22]./255;
        options.alpha      = 0.5;
        options.line_width = 2;
        options.error      = 'std';
    end
    if(isfield(options,'x_axis')==0), options.x_axis = 1:size(data,2); end
    options.x_axis = options.x_axis(:);

    % Computing the mean and standard deviation of the data matrix
    data_mean = mean(data,1);
    data_std  = std(data,0,1);

    % Type of error plot
    switch(options.error)
        case 'std', error = data_std;
        case 'sem', error = (data_std./sqrt(size(data,1)));
        case 'var', error = (data_std.^2);
        case 'c95', error = (data_std./sqrt(size(data,1))).*1.96;
    end

    % Plotting the result
    figure(options.handle);
    x_vector = [options.x_axis', fliplr(options.x_axis')];
    patch = fill(x_vector, [data_mean+error,fliplr(data_mean-error)], options.color_area);
    set(patch, 'edgecolor', 'none');
    set(patch, 'FaceAlpha', options.alpha);
    hold on;
    plot(options.x_axis, data_mean, 'color', options.color_line, ...
        'LineWidth', options.line_width);
    hold off;

end

author: kavli-ntnu


text:
function [x,D_struc,problem,r,res,solvertime,prob] = call_mosek_dual(model)

% Convert if the caller is bnb or bmibnb which might have appended bounds
% Sure, we could setup model with bounds, but...
[model.F_struc,model.K] = addStructureBounds(model.F_struc,model.K,model.ub,model.lb);

param = model.options.mosek;

prob.c = model.F_struc(1:model.K.f+model.K.l+sum(model.K.q)+3*model.K.e,1);
prob.a = -model.F_struc(1:model.K.f+model.K.l+sum(model.K.q)+3*model.K.e,2:end)';
prob.blc = -model.c;
prob.buc = -model.c;
prob.blx = -inf(size(prob.a,2),1);
prob.bux = inf(size(prob.a,2),1);
top = model.K.f+model.K.l;
prob.blx(1+model.K.f:model.K.f+model.K.l) = 0;

if model.K.q(1)>0 || model.K.e > 0
    prob.cones.type = [];
    prob.cones.subptr = [];
    prob.cones.sub = [];
end

if model.K.q(1)>0
    nq = length(model.K.q);
    prob.cones.type = zeros(nq, 1);
    prob.cones.subptr = zeros(nq, 1);
    prob.cones.sub = zeros(sum(model.K.q), 1);
    top0 = top;
    for i = 1:length(model.K.q)
        prob.cones.subptr(i) = top - top0 + 1;
        prob.cones.sub(top-top0+1:top-top0+model.K.q(i)) = top+1:top+model.K.q(i);
        top = top + model.K.q(i);
    end
end

if model.K.e>0
    for i = 1:model.K.e
        prob.cones.type = [prob.cones.type(:)' 3];
        prob.cones.subptr = [prob.cones.subptr(:)' length(prob.cones.sub)+1];
        prob.cones.sub = [prob.cones.sub(:)' top+3 top+2 top+1];
        top = top + 3;
    end
end

if model.K.s(1)>0
   prob = appendMosekSDPdata(model.F_struc,model.K,prob);
end

if model.options.savedebug
    ops = model.options;
    save mosekdebug prob param
end

if model.options.mosektaskfile
    mosekopt(sprintf('min write(%s) echo(0)', model.options.mosektaskfile), prob, param);
end

[r,res,solvertime] = doCall(prob,param,model.options);

try
    x = res.sol.itr.y;
catch
    if ~isempty(model.options.mosek) & isequal(model.options.mosek.MSK_IPAR_OPTIMIZER,'MSK_OPTIMIZER_FREE_SIMPLEX')
        x = res.sol.bas.y;
    else
        x = nan(length(model.c),1);
    end
end

if model.options.saveduals & ~isempty(x)
    try
        D_struc_SDP = zeros(sum(model.K.s.^2),1);
        top = 1;
        dtop = 1;
        for i = 1:length(model.K.s)
            n = model.K.s(i);
            I = find(tril(ones(n)));
            v = res.sol.itr.barx(top:((top+n*(n+1)/2)-1));
            D_struc_SDP(dtop + I - 1) = v;
            in = ceil(I/n);
            jn = mod(I-1,n)+1;
            D_struc_SDP(dtop + (jn-1)*n+in - 1) = v;
            top = top + n*(n+1)/2;
            dtop = dtop  + n^2;
        end
        D_struc = [res.sol.itr.xx;D_struc_SDP];
    catch
        D_struc = [];
    end
else
    D_struc = [];
end

problem = MosekYALMIPError(res);

function [res,sol,solvertime] = doCall(prob,param,options)

showprogress('Calling Mosek',options.showprogress);
if options.verbose == 0
    solvertime = tic;
    [res,sol] = mosekopt('minimize echo(0)',prob,param);
    solvertime = toc(solvertime);
else
    solvertime = tic;
    [res,sol] = mosekopt('minimize info',prob,param);
    solvertime = toc(solvertime);
end

function problem = MosekYALMIPError(res)

if res.rcode == 2001
    problem = 1;
    return
elseif res.rcode == 1305
    problem = -4;
    return
elseif res.rcode == 10007
    problem = 16;
    return
elseif res.rcode == 1400
    problem = 20;
    return
elseif res.rcode == 1001
    problem = -11;
    return;
elseif res.rcode == 1008
    problem = -12;
    return;
end

try
    solinfo = res.sol.itr;
catch
    solinfo = res.sol.bas;
end

switch solinfo.prosta
    case 'PRIMAL_AND_DUAL_FEASIBLE'
        problem = 0;
    case 'DUAL_INFEASIBLE'
        problem = 1;
    case 'PRIMAL_INFEASIBLE'
        problem = 2;
    case 'MSK_RES_TRM_USER_CALLBACK'
        problem = 16;
    case 'MSK_RES_TRM_STALL'
        problem = 4;
    case 'UNKNOWN'
        try
            if isequal(res.rcodestr,'MSK_RES_TRM_STALL')
                problem = 4;
            elseif isequal(res.rcodestr,'MSK_RES_OK')
                problem = 0;
            else
                problem = 11;
            end
        catch
            problem = 9;
        end
    otherwise
        problem = -1;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Nonnegative Matrix Factorization with Quadratic Programming
% The function iptrnr(Y,X,no_iter,lambdaA,rho,eta,epsil_A) returns
% the matrix A for the constrained system of linear equations: AX = Y

function A = iptrnr(Y,X,no_iter,lambdaA,rho,eta,epsil_A)

% INPUTS:
% > Y:        observation matrix [I by K] of mixed signals (images)
% > X:        source component matrix [J by K]
% > no_iter:  number of inner iterations
% > lambdaA:  Tikhonov regularization parameter
% > rho:      auxiliary parameter for defining the log barrier parameter
% > eta:      parameter in stopping rule for iterations
% > epsil_A:  threshold for active-set variables

% OUTPUT:
% > A:        mixing matrix [I by J]

[J,K] = size(X); [I,Ky] = size(Y); M = I*J;

% Initialization
A = ones(I,J);
Q_bar = spalloc(M,M,M*J); Q_tilde = spalloc(M,M,M*J);

B = X*X'; C_bar = X*Y';
Q_bar = kron(speye(I),B); % Hessian of D(Y||AX) with respect to A
Go = (C_bar - B*A'); % Gradient of D(Y||AX) with respect to A^T
At = A';
theta = -(rho/(M))*abs(Go(:)'*At(:)); % initialization for
                                      %log barrier parameter
l = 0;
while l < no_iter % inner iterations

    l = l + 1;
    a = At(:);
    active = find(a < epsil_A); % active set
    inactive = find(a >= epsil_A); % inactive set

    if ~isempty(active)
        a_tilde = a;             a_tilde(active) = [];
        N = length(a_tilde);
        c_tilde = C_bar(:);      c_tilde(active) = [];
        c_tilde = 2*theta*1./a_tilde - c_tilde;
        Q_tilde = Q_bar;
        Q_tilde(:,active) = [];  Q_tilde(active,:) = [];
        Q_tilde = Q_tilde - theta*spdiags(1./a_tilde.^2,0,N,N);
    else
        a_tilde = a;
        Q_tilde = Q_bar - theta*spdiags(1./a.^2,0,M,M) + lambdaA*speye(M);
        c_tilde = 2*theta*(1./a) - C_bar(:);
    end

     [Q,R] = qr(Q_tilde,-c_tilde); % Q-less QR factorization
     z_tilde = R\Q; % Gaussian elimination
     h_tilde = z_tilde - a_tilde;
     beta = min(1, 0.9995*min(a_tilde./abs(h_tilde)));
     a_tilde = a_tilde + beta*h_tilde;

     theta_tilde = theta*(2*(1./a_tilde) - (1./a_tilde.^2).*z_tilde);
     theta = ( theta_tilde'*a_tilde);

     if ~isempty(active)
        an = zeros(M,1);   an(inactive) = a_tilde;   a = an;
     else
        a = a_tilde;
     end

     if abs(theta) < eta*norm(a) % stopping rule
        break
     end
     theta = (rho/sqrt(M))*theta;
     At = reshape(a,J,I);

end
A = At';

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
B.e0  = 12.69; % [V]
B.e1  = -3.14; % [V]
B.e2  = 1.58;  % [V]
B.A   = 1.53;  % [V]
B.B   = 29.89; % [Ah^(-1)]
B.V0  = 12.59; % [V]
B.Q0  = 0.034; % [Ah]
B.Qf  = 1.20;  % [Ah]
B.R   = 0.061; % [Ohm]
B.tau = 1.95;  % [s]

B.p1  = 1.31e-7;
B.p2  = 4.03e-15;
B.p3  = -1.22e-23;
B.p4  = 1.65e-31;

Pboard = 0;

% Pboard = ;

nominal_V = @(Q)   B.e0 + B.e1*(Q/B.Qf) + B.e2*(Q/B.Qf)^2; %[V]
power   = @(rpm)   0.73*(B.p1*rpm^2 + B.p2*rpm^4 + B.p3*rpm^6 + B.p4*rpm^8); %[W]

Q       = B.Q0;
V       = B.V0;
i       = zeros(5,1);

w1 = 800*30/pi;  %rad/s to rpm
w2 = 800*30/pi;
w3 = 800*30/pi;
w4 = 800*30/pi;

window_size = 5;
b = (1/window_size)*ones(1,window_size);
a = 1;

dt = 0.1;
t = 0;
figure;

while Q < B.Qf
    % Sum board and motors power
    t           = t + dt;
    Pow         = Pboard + power(w1) + power(w2) + power(w3) + power(w4);
    i           = circshift(i, -1);
    i(end)      = Pow/V;
    signal_filt = filter(b,a,i);
    i_filt      = signal_filt(end);

    Q       = Q + i_filt*dt/3600; %[Ah]
    Vf      = B.A*exp(-B.B*(B.Qf-B.Q0-Q));
    V_nom   = nominal_V(Q);
    V       = V_nom - B.R*i_filt - Vf;

    plot(t,V,'bo', t,Q,'ro');
    hold on;

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function c = model_center_of_mass(model)

c = [mean(model.x),mean(model.y),mean(model.z)];

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function dummyObject = fixDummyObjectSize(dummyObject,originalObject);

frameUpdateInterval = dummyObject.frameUpdateInterval;
endMargin = size(originalObject.spectrogram,2)*frameUpdateInterval-max(originalObject.anchorTimeLocation);
if size(dummyObject.spectrogram,2)*frameUpdateInterval < max(dummyObject.anchorTimeLocation)+endMargin
    dummyFrameSize = max(dummyObject.anchorTimeLocation)+endMargin;
    dimmyFrequencySize = size(originalObject.spectrogram,1);
    dummyObject.spectrogram = ones(dimmyFrequencySize,dummyFrameSize);
    dummyObject.aperiodicityIndex = ones(dimmyFrequencySize,dummyFrameSize);
    dummyObject.F0 = ones(1,dummyFrameSize);
    dummyObject.vuv = ones(1,dummyFrameSize);
end;

author: HidekiKawahara


text:
% The COBRAToolbox: testConnectedComponents.m
%
% Purpose:
%     - tests the Connected Components functionality using a toy model
%
%
% Authors:
%     - Thomas Pfau

%% The model created looks as follows:
%
%
%   <-> A -> B -> C <->
%
%
%              -> F -
%             /      \
%   <-> D -> E -> G ---> H <->
%            ^
%            |
%            v
%
% And should yield 2 connected components, along with a bunch of exchangers
% which are ignored.

currentDir = pwd;
testdir = fileparts(which('testConnectedComponents.m'));
cd(testdir)

requiredToolboxes = {'image_toolbox'};

prepareTest('toolboxes',requiredToolboxes);

model = createToyModelForConnectedComponentAnalysis();
[groups,orphans,R,C] = connectedComponents(model);

%The orphan reactions, are the exchangers, which by definition are dropped
%from the connected components.Those are the last 5 reactions)
assert(isequal(orphans,numel(model.rxns)-4:numel(model.rxns)));

assert(isequal(model.rxns(1:2),model.rxns(groups(1).elements)));
assert(isequal(model.rxns(3:6),model.rxns(groups(2).elements)));

[groups,orphans] = connectedComponents(model,'largestComponent',0,1);

assert(isequal(numel(groups),1));
assert(isequal(model.rxns(3:6),model.rxns(groups.elements)));

delete('reactionsNotConnectedByAnything.txt');
delete('reactionAdjacencyOtherThanCofactors.txt')

cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [wp_lat_deg, wp_lon_deg, wp_alt_ft]=get_wp(wp, host, port)
% [WP_LAT_DEG, WP_LON_DEG, WP_ALT_FT]=GET_WP(WP, HOST, PORT)
%
% Read waypoint WP from waypoint stack in Flightgear over socket
% connection (HOST, PORT).
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".
%
% Function returns:
% - Lattitude degree
% - Longitude degree
% - Altitude feet

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);
in = BufferedReader(InputStreamReader(fg_socket.getInputStream()));

% Read waypoint wp
% Latitude degree
if wp == 1
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp/latitude-deg\r\n'));
else
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp[%d]/latitude-deg\r\n'...
    ,wp-1));
end
wp_lat_deg_str = char(in.readLine());
j = findstr(wp_lat_deg_str,'''');
wp_lat_deg = str2double(wp_lat_deg_str(j(1)+1:j(2)-1));

% Longitude degree
if wp == 1
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp/longitude-deg\r\n'));
else
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp[%d]/longitude-deg\r\n'...
    ,wp-1));
end
wp_lon_deg_str = char(in.readLine());
j = findstr(wp_lon_deg_str,'''');
wp_lon_deg = str2double(wp_lon_deg_str(j(1)+1:j(2)-1));

% Altitude feet
if wp == 1
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp/altitude-ft\r\n'));
else
    out.println(sprintf(...
    'get /autopilot/route-manager/route/wp[%d]/altitude-ft\r\n'...
    ,wp-1));
end
wp_alt_ft_str = char(in.readLine());
j = findstr(wp_alt_ft_str,'''');
wp_alt_ft = str2double(wp_alt_ft_str(j(1)+1:j(2)-1));

out.close();
in.close();
fg_socket.close();

end

author: Sable

repo: mcbench-benchmarks


text:
function hash = calculate_results_fingerprint(tracker, experiment, sequences)
% calculate_results_fingerprint Calculate results hash
%
% Calculates a hash fingerprint based on timestamps of result files for a tracker
%
% Input:
% - tracker: Tracker structure.
% - experiment: Experiment structure.
% - sequences: Cell array of sequence structures.
%
% Output:
% - hash: A string containing hash fingerprint.

time_string = iterate(experiment, tracker, sequences, 'iterator', @fingerprint_iterator, 'context', []);

hash = md5hash(time_string);

end

function context = fingerprint_iterator(event, context)
% fingerprint_iterator Iterator function
%
% Iterator function that calculates parts of the fingerprint.
%
% Input:
% - event (structure): Event structure.
% - context (structure): Context structure.
%
% Output:
% - context (structure): Modified context structure.
%

    switch (event.type)
        case 'sequence_enter'

            files = tracker_evaluate(event.tracker, event.sequence, event.experiment, ...
                'scan', true);

            dates = zeros(1, numel(files));

            for j = 1:numel(files)
                stat = dir(files{j});
                dates(j) = stat.datenum;
            end;

            context = [context, dates];

    end;

end

author: votchallenge

repo: toolkit-legacy


text:
function i4vec_mask_print ( n, a, mask_num, mask, title )

%*****************************************************************************80
%
%% I4VEC_MASK_PRINT prints a masked I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of components of the vector.
%
%    Input, integer A(N), the vector to be printed.
%
%    Input, integer MASK_NUM, the number of masked elements.
%
%    Input, integer MASK(MASK_NUM), the indices of the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Masked vector printout:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = 1 : mask_num
    fprintf ( 1, '  %6d: %6d  %10d\n', i, mask(i), a(mask(i)) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [img_files, pos, target_sz, ground_truth, video_path] = load_video_info(base_path, video)
%LOAD_VIDEO_INFO
%   Loads all the relevant information for the video in the given path:
%   the list of image files (cell array of strings), initial position
%   (1x2), target size (1x2), the ground truth information for precision
%   calculations (Nx2, for N frames), and the path where the images are
%   located. The ordering of coordinates and sizes is always [y, x].
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	%see if there's a suffix, specifying one of multiple targets, for
	%example the dot and number in 'Jogging.1' or 'Jogging.2'.
	if numel(video) >= 2 && video(end-1) == '.' && ~isnan(str2double(video(end))),
		suffix = video(end-1:end);  %remember the suffix
		video = video(1:end-2);  %remove it from the video name
	else
		suffix = '';
	end

	%full path to the video's files
	if base_path(end) ~= '/' && base_path(end) ~= '\',
		base_path(end+1) = '/';
	end
	video_path = [base_path video '/'];

	%try to load ground truth from text file (Benchmark's format)
	filename = [video_path 'groundtruth_rect' suffix '.txt'];
	f = fopen(filename);
	assert(f ~= -1, ['No initial position or ground truth to load ("' filename '").'])

	%the format is [x, y, width, height]
	try
		ground_truth = textscan(f, '%f,%f,%f,%f', 'ReturnOnError',false);
	catch  %#ok, try different format (no commas)
		frewind(f);
		ground_truth = textscan(f, '%f %f %f %f');
	end
	ground_truth = cat(2, ground_truth{:});
	fclose(f);

	%set initial position and size
	target_sz = [ground_truth(1,4), ground_truth(1,3)];
	pos = [ground_truth(1,2), ground_truth(1,1)] + floor(target_sz/2);

	%from now on, work in the subfolder where all the images are
	video_path = [video_path 'img/'];

	%for these sequences, we must limit ourselves to a range of frames.
	%for all others, we just load all png/jpg files in the folder.
	frames = {'David', 300, 465; %770
			  'Football1', 1, 74;
			  'Freeman3', 1, 460;
			  'Freeman4', 1, 283};

	idx = find(strcmpi(video, frames(:,1)));

	if isempty(idx),
		%general case, just list all images
		img_files = dir([video_path '*.png']);
		if isempty(img_files),
			img_files = dir([video_path '*.jpg']);
			assert(~isempty(img_files), 'No image files to load.')
		end
		img_files = sort({img_files.name});
        target_sz = [ground_truth(1,4), ground_truth(1,3)];
	    pos = [ground_truth(1,2), ground_truth(1,1)] + floor(target_sz/2);
    else

        target_sz = [ground_truth(frames{idx, 2},4), ground_truth(frames{idx, 2},3)];
        pos = [ground_truth(frames{idx, 2},2), ground_truth(frames{idx, 2},1)] + floor(target_sz/2);
		%list specified frames. try png first, then jpg.
		if exist(sprintf('%s%04i.png', video_path, frames{idx,2}), 'file'),
			img_files = num2str((frames{idx,2} : frames{idx,3})', '%04i.png');

		elseif exist(sprintf('%s%04i.jpg', video_path, frames{idx,2}), 'file'),
			img_files = num2str((frames{idx,2} : frames{idx,3})', '%04i.jpg');

		else
			error('No image files to load.')
		end

		img_files = cellstr(img_files);
    end

    	if size(ground_truth,1) == 1,
		%we have ground truth for the first frame only (initial position)
		ground_truth = [];
	else
		%store positions instead of boxes
		ground_truth = ground_truth(:,[2,1]) + ground_truth(:,[4,3]) / 2;
	end

end

repo: KCF


text:
function [M,fig] = plot_marginals(sampleCell,T,truespikes,int_show)

% Plots marginal posterior empirical pdfs for # of spikes for each timebin
% similar to figure 1B in Pnevmatikakis et al., Neuron 2016

% Inputs:
% sampleCell:   Cell array with spike times in continuous time (SAMPLES.ss)
% T:            Number of timebins
% truespikes:   Number of true spikes per timebin (vector of size T x 1)
% int_show:     Show only a specified interval (default: [1,T])

% Output:
% M:            matrix of empirical posterior pdfs for each timebin

% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation

if nargin < 4
    int_show = 1:T;
end
nT = length(int_show);

if nargin == 2
    truespikes = -0.5*ones(1,nT);
end

if length(truespikes) == T
    truespikes = truespikes(int_show);
end

Mat = samples_cell2mat(sampleCell,T);
Mat = Mat(:,int_show);

mS = min(max([Mat(:);truespikes(:)])+1,6);
M = zeros(mS,nT);
for i = 1:nT
    M(:,i) = hist(Mat(:,i),0:mS-1)/size(Mat,1);
end

cmap = bone(100);
cmap(2:26,:) = [];
fig = figure;
imagesc(M); axis xy;
%set(gca,'Ytick',[0.5:(mS+.5)],'Yticklabel',[-1:(mS)]);  %hold all; plot(traceData.spikeFrames+1); set(gca,'YLim',[1,5])
%set(gca,'YLim',[1.25,mS+.25]);
hold all; scatter(1:nT,truespikes+1,[],'m');
set(gca,'YLim',[1.5,mS+.125]);
pos = get(gca,'Position');
set(gca,'Ytick',0.5+[-0.5:mS-.5],'Yticklabel',[-1+(0:mS)]);
colormap(cmap);
ylabel('# of Spikes ','fontweight','bold','fontsize',14);
xlabel('Timestep ','fontweight','bold','fontsize',14);
title('Posterior Spike Histogram (MCMC) ','fontweight','bold','fontsize',14);
cbar = colorbar('Location','East');
cpos = get(cbar,'Position');
set(cbar,'Position',[pos(1)+pos(4),cpos(2:4)]);
%set(gca,'Xtick',[])
%set(cbar,'Color',[1,1,1]);
set(cbar,'Fontsize',12);

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [ fea, out ] = ex_poisson1( varargin )
%EX_POISSON1  1D Poisson equation example.
%
%   [ FEA, OUT ] = EX_POISSON1( VARARGIN ) Poisson equation on a line with a
%   constant source term equal to 1, homogenous boundary conditions, and
%   exact solution (-x^2+x)/2. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       hmax        scalar {1/10}          Grid cell size
%       sfun        string {sflag1}        Finite element shape function
%       iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
%                                          or directly define fea.eqn/bdr fields (=0)
%                                          or use core assembly functions        (<0)
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'hmax',   1/10;
            'sfun',   'sflag1';
            'refsol', '(-x^2+x)/2';
            'fsrc',   '1';
            'iphys',  1;
            'icub',   2;
            'iplot',  1;
            'tol',    2e-2;
            'fid',    1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Grid generation.
if( opt.hmax>0 )
  nx = round( 1/opt.hmax );
  fea.grid = linegrid( nx, 0, 1 );
else   % Scrambled testing grid.
  fea.grid.p = [ 0 1/10 4/10 1/3 1 1-1/3 ];
  fea.grid.c = [ 1 4 2 6 3 ;
                 2 3 4 5 6 ];
  fea.grid.a = [ 0 2 1 4 3 ;
                 2 4 3 0 5 ];
  fea.grid.b = [ 1 1 1 -1 ;
                 4 2 2  1 ]';
  fea.grid.s = ones(1,5);
end
n_bdr = 2;   % Number of boundaries.

% Problem definition.
fea.sdim  = { 'x' };                    % Coordinate name.
switch opt.iphys

 case 0   % Directly define fea.eqn/bdr fields.

  fea.dvar  = { 'u' };                  % Dependent variable name.
  fea.sfun  = { opt.sfun  };            % Shape function.

  % Define equation system.
  fea.eqn.a.form = { [2;2] };           % First row indicates test function space   (2=x-derivative),
                                        % second row indicates trial function space (2=x-derivative).
  fea.eqn.a.coef = { 1 };               % Coefficient used in assembling stiffness matrix.

  fea.eqn.f.form = { 1 };               % Test function space to evaluate in right hand side (1=function values).
  fea.eqn.f.coef = { opt.fsrc };        % Coefficient used in right hand side.

  % Define boundary conditions.
  if( strcmp(opt.sfun(end-1:end),'H3') )   % Prescribed derivatives at end points for Hermite elements.
    fea.bdr.d = {{ 0 0 ; 1/2 -1/2 }};
  else
    fea.bdr.d     = cell(1,n_bdr);
    [fea.bdr.d{:}] = deal(0);              % Assign zero to all boundaries (homogenous Dirichlet conditions).
  end

  fea.bdr.n     = cell(1,n_bdr);        % No Neumann boundaries ('fea.bdr.n' empty).

  % Parse and solve problem.
  fea       = parseprob(fea);           % Check and parse problem struct.
  fea.sol.u = solvestat(fea,'fid',fid,'icub',opt.icub); % Call to stationary solver.

 case 1   % Use physics mode.

  fea = addphys(fea,@poisson);          % Add Poisson equation physics mode.
  fea.phys.poi.sfun = { opt.sfun };     % Set shape function.
  fea.phys.poi.eqn.coef{3,4} = { opt.fsrc };   % Set source term coefficient.
  fea.phys.poi.bdr.coef{1,end} = repmat({0},1,n_bdr);   % Set Dirichlet boundary coefficient to zero.
  fea = parsephys(fea);                 % Check and parse physics modes.
  if( strcmp(opt.sfun(end-1:end),'H3') )   % Prescribed derivatives at end points for Hermite elements.
    fea.bdr.d = {{ 0 0 ; 1/2 -1/2 }};
  end

  % Parse and solve problem.
  fea       = parseprob(fea);           % Check and parse problem struct.
  fea.sol.u = solvestat(fea,'fid',fid,'icub',opt.icub); % Call to stationary solver.

 otherwise   % Use core assembly functions.

  fea.dvar  = { 'u' };                  % Dependent variable name.
  fea.sfun  = { opt.sfun  };            % Shape function.
  fea       = parseprob(fea);           % Check and parse problem struct.

  % Assemble stiffness matrix.
  form  = [2;2];
  sfun  = {opt.sfun;opt.sfun};
  coefa = 1;
  sind  = 1;
  i_cub = opt.icub;

  [vRowInds,vColInds,vAvals,n_rows,n_cols] = ...
    assemblea(form,sfun,coefa,i_cub,fea.grid.p,fea.grid.c,fea.grid.a,fea.grid.s,[]);
  A = sparse(vRowInds,vColInds,vAvals,n_rows,n_cols);

  % Check and compare with finite difference stencil.
  if (strcmp(opt.sfun,'sflag1'))
    h = 1/nx;
    n = nx+1;
    e = ones(n,1);
    A_ref = 1/h*spdiags([-e 2*e -e], -1:1, n, n);
    A_ref(1) = A_ref(1)/2;
    A_ref(end) = A_ref(end)/2;

    err = norm(A(:)-A_ref(:));
    if err>opt.tol
      out.err  = err;
      out.pass = -1;
      return
    end
  end

  form  = 1;
  sfun  = sfun{1};
  coeff = 1;

  f = assemblef(form,sfun,coeff,i_cub,fea.grid.p,fea.grid.c,fea.grid.a,fea.grid.s,[]);

  % Check and compare with finite difference stencil.
  if (strcmp(opt.sfun,'sflag1'))
    f_ref          = coeff*h*ones(n,1);
    f_ref([1 end]) = coeff*1/2*h;

    err = norm(f-f_ref);
    if err>1e-6
      out.err  = err;
      out.pass = -2;
      return
    end
  end

  % Set homogenous Dirichlet boundary conditions on first and last dof/node.
  bind = [1 nx+1];
  A = A';   %'
  A(:,bind) = 0;         % Zero out Dirichlet BC rows.
  for i=1:length(bind)   % Loop to set diagonal entry to 1.
    i_a = bind(i);
    A(i_a,i_a) = 1;
  end
  A = A';   %'
  f(bind) = 0;           % Set corresponding source term entries to Dirichlet BC values.

  % Solve problem.
  fea.sol.u = A\f;

end

% Postprocessing.
if ( opt.iplot>0 )
  x = linspace( 0, 1, 41 );
  u = evalexpr( 'u', x, fea )';
  figure
  subplot(3,1,1)
  plot( x, u )
  axis( [0 1 0 0.2])
  grid on
  title('Solution u')
  subplot(3,1,2)
  ux = (-x.^2+x)/2;
  plot( x, ux )
  axis( [0 1 0 0.2])
  grid on
  title('Exact solution')
  subplot(3,1,3)
  plot( x, abs(ux-u) )
  title('Error')
end

% Error checking.
xi = [1/2; 1/2];
s_err = ['abs(',opt.refsol,'-u)'];
err = evalexpr0(s_err,xi,1,1:size(fea.grid.c,2),[],fea);
ref = evalexpr0('u',xi,1,1:size(fea.grid.c,2),[],fea);
err = sqrt(sum(err.^2)/sum(ref.^2));

if( ~isempty(fid) )
  fprintf(fid,'\nL2 Error: %e\n',err)
  fprintf(fid,'\n\n')
end

out.err  = err;
out.tol  = opt.tol;
out.pass = out.err<out.tol;
if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function days = days_before_month_common ( y, m )

%*****************************************************************************80
%
%% DAYS_BEFORE_MONTH_COMMON returns the number of days before a Common month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of
%    days in the year before the first day of the given month.
%
  mdays = [ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 ];
%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_common ( y2, m2 );

  if ( ierror ~= 0 )
    days = 0;
    return
  end

  days = mdays ( m2 );

  if ( 2 < m2 && year_is_leap_common ( y2 ) )
    days = days + 1;
  end

  return
end


text:
function f = f_equilibrium_1d(r,u,x,t,R)
% Compute equilibrum in 1d cases for:
%
% MB:  Maxwell-Boltzmann, theta =  0
% FD:  Fermi-Diract,      theta = +1
% BE:  Bose-Einstein,     theta = -1
%
% inputs:
% u: macroscopic velocity
% x: microscopic velocity
% t: temperature

% f = 1./(exp( (x-u).^2 ./ t)./r + theta);

f = r*(1/2*pi*R*t)./(exp(- (x-u).^2 ./ (2*R*t)))

repo: Aero-matlab


text:
function [duplicated_elems] = duplicateElems(elems, num_duplicates)
%DUPLICATEELEMS takes an elems vector of length N and num_duplicates vector
% of length N, and produces a duplicated_elems vector of length
% sum(num_duplicates). Each element in num_duplicates specifies how many
% times an element in elems needs to be replicated.
%
% E.g. on how to do it in a loopy way:
%
%   elems = [2,4,6,7,10,9,20,15];
%   num_duplicates = [0,2,4,0,3,2,0,0];
%   duplicated_elems = [];
%   % loopy way to do it
%   for idx = 1:length(elems)
%       duplicated_elems = [duplicated_elems, ...
%                           repmat(elems(idx), 1, num_duplicates(idx))];
%   end
%   disp(duplicated_elems)
%
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    assert(isvector(elems), 'elems should be a vector');
    assert(isvector(num_duplicates), 'num_duplicates should be a vector');
    assert(numel(elems) == numel(num_duplicates), ...
        'elems and num_duplicates should be of the same size');

    num_duplicates = num_duplicates(:);

    % remove elements which weren't replicated
    elems(num_duplicates <= 0) = [];
    num_duplicates(num_duplicates <= 0) = [];

    % if no duplication element is present in num_duplicates
    if isempty(num_duplicates)
        duplicated_elems = [];
        return;
    end

    duplicate_idxs = zeros(sum(num_duplicates), 1);
    % mark the starting locations
    duplicate_idxs([1; cumsum(num_duplicates(1:end-1))+1]) = 1;
    duplicate_idxs = cumsum(duplicate_idxs);
    duplicated_elems = elems(duplicate_idxs);
end

author: Cloud-CV

repo: object-proposals


text:
function code = tree_rooted_code ( nnode, parent )

%*****************************************************************************80
%
%% TREE_ROOTED_CODE returns the code of a rooted tree.
%
%  Discussion:
%
%    This code for a rooted tree depends on the node ordering, so it's actually
%    the code for a labeled rooted tree.  To eliminate the effects of node
%    labeling, one could choose as the code for a tree the maximum of all
%    the codes associated with the different possible labelings of the tree.
%    There are more effective ways of arriving at this code than simply
%    generating all possible codes and comparing them.
%
%    For a tree with NNODES, the code is a list of 2*NNODE 0's and 1's,
%    describing a traversal of the tree starting at an imaginary node 0,
%    moving "down" to the root (a code entry of 1), and then moving
%    "down" (1) or "up" (0) as the tree is traversed in a depth first
%    manner.  The final move must be from the root up to the imaginary
%    node 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NNODE, the number of nodes.
%
%    Input, integer PARENT(NNODE), is the parent node of each node.
%    The node with parent 0 is the root.
%
%    Output, integer CODE(2*NNODE), the code for the tree.
%

%
%  Find the root.
%
  father = 0;
  for i = 1 : nnode
    if ( parent(i) == 0 )
      k = 1;
      code(1) = 1;
      father = i;
      break
    end
  end

  if ( father == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TREE_ROOTED_CODE - Fatal error!\n' );
    fprintf ( 1, '  Could not find the root.\n' );
    error ( 'TREE_ROOTED_CODE - Fatal error!' );
  end

  while ( father ~= 0 )

    k = k + 1;
    code(k) = 0;

    for son = 1 : nnode
      if ( parent(son) == father )
        code(k) = 1;
        father = son;
        break
      end
    end

    if ( code(k) == 0 )
      parent(father) = - parent(father);
      father = - parent(father);
    end

  end

  parent(1:nnode) = - parent(1:nnode);

  return
end


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] = ubah_run(fid, data, tc, opts)
% This file is the run core for the market strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] = market_run(fid, data, tc, opts)
%
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
% daily_portfolio: daily portfolio, achieved by the strategy
%
% data: market sequence vectors
% fid: handle for write log file
% tc: transaction fee rate
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret] ...
%          = market_run(fid, data, 0, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m]=size(data);

% Variables for return, start with uniform weight
cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);

% portfolio at the beinning (end) of a period
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% print log file head
fprintf(fid, '-------------------------------------\n');
fprintf(fid, 'Parameters [tc: %f].\n', tc);
fprintf(fid, 'day\t Daily Return\t Total return\n');

fprintf(1, '-------------------------------------\n');
if(~opts.quiet_mode)
    fprintf(1, 'Parameters [tc: %f].\n', tc);
    fprintf(1, 'day\t Daily Return\t Total return\n');
end

% Backtests
for t = 1:1:n,

    % Calculate t's portfolio at the beginning of t-th trading day
    if (t >= 2)
        [day_weight] = ubah_kernel(data(1:t-1, :), day_weight_o);
    end

    % Normalize the constraint, always useless
    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    % Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % Adjust weight(t, :) for the transaction cost issue
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Log information
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
        end
    end
end

% Output the cumulative return and log it.
fprintf(fid, 'Market(tc=%.4f), Cumulative return: %.2f\n',tc, cum_ret);
fprintf(fid, '-------------------------------------\n');
fprintf(1, 'Market(tc=%.4f), Cumulative return: %.2f\n',tc, cum_ret);
fprintf(1, '-------------------------------------\n');

end
%%%%%%%%%%%%%%End%%%%%%%%%%%%%%%%%%%%%%

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function out=asin(x)

out=-i.*log(i.*x+sqrt(1-x.^2));


text:
function test_suite = test_binomial1

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_BINOMIAL1
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictiveMeanVarianceAndDensity(testCase)
  % Test predictive mean, variance and density
  verifyVarsEqual(testCase, getName(), ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

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


text:
function [h,g,a,info]=wfilt_lemarie(N)
%WFILT_LEMARIE  Battle and Lemarie filters
%   Usage: [h,g,a]=wfilt_lemarie(N)
%
%   Input parameters:
%         N     : Filter length, must be even.
%
%   `[h,g,a]=wfilt_lemarie(N)` calculates $N$ (even) truncated coeficients
%   of orthonormal Battle-Lemarie wavelets. Filter coefficients are obtained
%   by frequency domain sampling and trunctating the impulse response.
%   Due to the truncation, the filterbank might not achieve a perfect
%   reconstruction. The filetrs are included nevertheless since they were
%   the original ones used in the first MRA paper.
%
%   Examples:
%   ---------
%   :::
%     wfiltinfo('lemarie50');
%
%   References: mallat89atheory

% Original copyright goes to:
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
% Author: Jose Martin Garcia
% e-mail: Uvi_Wave@tsc.uvigo.es

if rem(N,2)~=0
    error('%s: Filter length must be even.',upper(mfilename));
end

num_coefs = N;
L = 1024;
H = wfreq_lemarie(L);
hh=real(ifft(H{1},L));
hh=[ hh(L-floor(num_coefs/2)+1:L) hh(1:ceil(num_coefs/2))];
hh=hh/norm(hh);

g{1} = (hh);
g{2} = -(-1).^(1:length(hh)).*g{1}(end:-1:1);

g = cellfun(@(gEl) struct('h',gEl,'offset',-floor(numel(gEl)/2)),g,'UniformOutput',0);

h = g;
a= [2;2];
info.istight = 1;

function [H,G] = wfreq_lemarie(L)
%WFREQ_LEMARIE  Battle and Lemarie filters frequency resp. sampling
%   Usage: [H,G]=wfreq_lemarie(L)
%
%   Input parameters:
%         N     : Number of samples of the frequency response.
%
%   `[H,G]=wfreq_lemaire(L)` calculates $L$ samples of the Battle and
%   Lemarie filters frequency responses.
%
%   References: mallat89atheory
%
%

% Original copyright goes to:
% Copyright (C) 1994, 1995, 1996, by Universidad de Vigo
% Author: Jose Martin Garcia
% e-mail: Uvi_Wave@tsc.uvigo.es

% frequency axis
w=[0:2*pi/L:2*pi*(1-1/L)];
w(1)=eps;
w(L/2+1)=w(L/2+1)+1e-15;

% calculation of frequency response of analysis lowpass filter
num=0;den=0;
K=36;
for k=-K:K,
	num=1./((w+2*pi*k).^8)+num;
	den=1./((2*w+2*pi*k).^8)+den;
end
H = cell(2,1);
H{1}=sqrt(num./(2.^8*den));
H{1}(1)=1;

H{2} = fftshift(H{1});
G = cell(2,1);
G{1} = fliplr(H{1});
G{2} = fliplr(H{2});

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function y = tifs2cv(f,m,d,q)
%TIFS2CV Compresses a multi-frame TIFF image sequence.
%   Y = TIFS2CV(F,M,D,Q) compresses multiframe TIFF F using motion
%   compensated frames, 8 x 8 DCT transforms, and Huffman coding. If
%   parameter Q is omitted or is 0, only Huffman encoding is used and
%   the compression is lossless; for Q > 0, lossy JPEG encoding is
%   performed. The inputs are:
%
%      F     A multi-frame TIFF file        (e.g., 'file.tif')
%      M     Macroblock size                (e.g., 8)
%      D     Search displacement            (e.g., [16 8])
%      Q     JPEG quality for IM2JPEG       (e.g., 1)
%
%   Output Y is an encoding structure with fields:
%
%      Y.blksz      Size of motion compensation blocks
%      Y.frames     The number of frames in the image sequence
%      Y.quality    The reconstruction quality
%      Y.motion     Huffman encoded motion vectors
%      Y.video      An array of MAT2HUFF or IM2JPEG coding structures
%
%   See also CV2TIFS.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% The default reconstruction quality is lossless.
if nargin < 4
    q = 0;
end

% Compress frame 1 and reconstruct for the initial reference frame.
if q == 0
    cv(1) = mat2huff(imread(f,1));
    r = double(huff2mat(cv(1)));
else
    cv(1) = im2jpeg(imread(f,1),q);
    r = double(jpeg2im(cv(1)));
end
fsz = size(r);

% Verify that image dimensions are multiples of the macroblock size.
if ((mod(fsz(1),m) ~= 0) || (mod(fsz(2),m) ~= 0))
    error('Image dimensions must be multiples of the block size.');
end

% Get the number of frames and preallocate a motion vector array.
fcnt = size(imfinfo(f),1);
mvsz = [fsz/m 2 fcnt];
mv = zeros(mvsz);

% For all frames except the first, compute motion conpensated
% prediction residuals and compress with motion vectors.
for i = 2:fcnt
    frm = double(imread(f,i));
    frmC = im2col(frm,[m m],'distinct');
    eC = zeros(size(frmC));

    for col = 1:size(frmC,2)
        lookfor = col2im(frmC(:,col),[m m],[m m],'distinct');

        x = 1 + mod(m * (col - 1),fsz(1));
        y = 1 + m * floor((col - 1) * m / fsz(1));
        x1 = max(1,x - d(1));
        x2 = min(fsz(1),x + m + d(1) - 1);
        y1 = max(1,y - d(2));
        y2 = min(fsz(2),y + m + d(2) - 1);

        here = r(x1:x2,y1:y2);
        hereC = im2col(here,[m m],'sliding');
        for j = 1:size(hereC,2)
            hereC(:,j) = hereC(:,j) - lookfor(:);
        end
        sC = sum(abs(hereC));
        s = col2im(sC,[m m],size(here),'sliding');
        mins = min(min(s));
        [sx,sy] = find(s == mins);

        ns = abs(sx) + abs(sy);         % Get the closest vector
        si = find(ns == min(ns));
        n = si(1);

        mv(1 + floor((x - 1)/m), 1 + floor((y - 1)/m), 1:2, i) = ...
            [x - (x1 + sx(n) - 1) y - (y1 + sy(n) - 1)];
        eC(:,col) = hereC(:, sx(n) + (1 + size(here,1) - m) ...
            * (sy(n) - 1));
    end

    % Code the prediction residual and reconstruct it for use in
    % forming the next reference frame.
    e = col2im(eC,[m m],fsz,'distinct');
    if q == 0
        cv(i) = mat2huff(int16(e));
        e = double(huff2mat(cv(i)));
    else
        cv(i) = im2jpeg(uint16(e + 255),q,9);
        e = double(jpeg2im(cv(i)) - 255);
    end

    % Decode the next reference frame. Use the motion vectors to get
    % the subimages needed to subtract from the prediction residual.
    rC = im2col(e,[m m],'distinct');
    for col = 1:size(rC,2)
        u = 1 + mod(m * (col - 1),fsz(1));
        v = 1 + m * floor((col - 1) * m / fsz(1));
        rx = u - mv(1 + floor((u - 1)/m), 1 + floor((v - 1)/m), 1, i);
        ry = v - mv(1 + floor((u - 1)/m), 1 + floor((v - 1)/m), 2, i);
        temp = r(rx:rx + m - 1,ry:ry + m - 1);
        rC(:,col) = temp(:) - rC(:,col);
    end
    r = col2im(double(uint16(rC)),[m m],fsz,'distinct');
end

y = struct;
y.blksz = uint16(m);
y.frames = uint16(fcnt);
y.quality = uint16(q);
y.motion = mat2huff(mv(:));
y.video = cv;

author: dipum

repo: dipum-toolbox

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


text:
function y = R3R(x,batch)
% R3R is a repeated running medians of 3.
%
% Usage: y = R3R(x)
%
%   where  x is the input vector with length at least 3
%          y is the return vector of the same length as x
%

% See John Wilder Tukey, "EXPLORATORY DATA ANALYSIS".
% Addison-Wesley Publishing Co. 1977. Chpt 7A, 7D

if min(size(x)) ~= 1, error('Input must be a vector'); end
x = x(:); xlen = length(x); cont = 1; y = x;
if nargin<2, batch = 0; end
if xlen < 3,
   if ~batch, warning('Input vector must be at least 3 elements long'); end
   cont = 0;
end
if cont,
   % running median
   for ii = 2:xlen-1,
      y(ii) = median(x(ii-1:ii+1));
   end
   % compute the low-end point
   delta = abs( y(3)-y(2) );
   slope = sign( y(3)-y(2) );
   switch slope
   case 1
      if y(2)-y(1)<=0 | y(2)-y(1)>=2*delta, y(1) = y(2) - 2*delta; end
   case -1
      if y(1)-y(2)<=0 | y(1)-y(2)>=2*delta, y(1) = y(2) + 2*delta; end
   case 0
      y(1) = y(2);
   end
   % compute the high-end point
   delta = abs( y(xlen-1)-y(xlen-2) );
   slope = sign( y(xlen-1)-y(xlen-2) );
   switch slope
   case 1
      if y(xlen)-y(xlen-1)<=0 | y(xlen)-y(xlen-1)>=2*delta,y(xlen)=y(xlen-1)+2*delta;end
   case -1
      if y(xlen-1)-y(xlen)<=0 | y(xlen-1)-y(xlen)>=2*delta,y(xlen)=y(xlen-1)-2*delta;end
   case 0
      y(xlen) = y(xlen-1);
   end
end
% 3-ing to death?
if ~isequal(y(2:xlen-1),x(2:xlen-1)) & xlen > 3, y = R3R(y); end

return


text:
function [w,mu,P]=RunnalsGaussMixRed(w,mu,P,K,gammaBound,KMax)
%%RUNNALSGAUSSMIXRED Perform Gaussian mixture reduction using the greedy
%                    merging algorithm by Runnals in [1].
%
%INPUTS: w An NX1 or 1XN vector of weights of the components of the
%          original Gaussian mixture.
%       mu An xDimXN matrix of the means of the vector components of the
%          original Gaussian mixture.
%        P An xDim XxDim XN hypermatrix of the covariance matrices for the
%          components of the original Gaussian mixture.
%        K The minimum desired number of components desired in the mixture
%          after reduction. This will be always achieved unless gammaBound
%          is set. Often, if gammaBound is set, one will just make K=1.
% gammaBound If the number of components is less than or equal to KMax,
%          then reduction will continue until the lowest cost between two
%          components is greater than gammaBound or there are K components
%          left. The default if this parameter is omitted or an empty
%          matrix is passed is Inf, which means that reduction will
%          continue until K components is reached.
%     KMax Regardless of the gammaBound setting, it won't be allowed that
%          more than KMax components are present. The default if omitted or
%          an empty matrix is passed is Inf.
%
%OUTPUTS: w The length K vector of weights of the mixture after reduction.
%           (the length can be less than K if the input mixture had fewer
%           than K components).
%        mu The xDimXK means of the mixture after reduction.
%         P The xDimXxDimXK covariance matrices of the mixture after
%           reduction.
%
%This implements the suboptimal Gaussian mixture reduction algorithm of
%[1]. Care is taken in the implementation to avoid copies and reallocations
%associated with actually resizing the cost array every time something is
%eliminated.
%
%EXAMPLE 1:
%This scalar example plots a full 9-component PDF and then the PDF reduced
%to 5 components using sequential brute-force reduction and Runnals'
%algorithm. Increasing the number of components maintained to 6 makes
%the results of Runnals' algorithm much closer to that of the brute-force
%solution.
% w=[0.03, 0.18, 0.12, 0.19, 0.02, 0.16, 0.06, 0.1, 0.08, 0.06];
% mu=[1.45, 2.20, 0.67, 0.48, 1.49, 0.91, 1.01, 1.42, 2.77, 0.89];
% P=reshape([0.0487, 0.0305, 0.1171, 0.0174, 0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679],[1,1,10]);
%
% K=6;
% %Runnals' cost function.
% [wRed,muRed,PRed]=RunnalsGaussMixRed(w,mu,P,K,[],[],0);
% %Sequential brute-force reduction.
% [wRedBF,muRedBF,PRedBF]=bruteForceGaussMixRed(w,mu,P,K,true);
% numPoints=500;
% xVals=linspace(0,3,numPoints);
% PDFVals1=GaussianMixtureD.PDF(xVals,w,mu,P);
% PDFVals2=GaussianMixtureD.PDF(xVals,wRedBF,muRedBF,PRedBF);
% PDFVals3=GaussianMixtureD.PDF(xVals,wRed,muRed,PRed);
% figure(1)
% clf
% hold on
% plot(xVals,PDFVals1,'-k','linewidth',4)
% plot(xVals,PDFVals2,'--r','linewidth',2)
% plot(xVals,PDFVals3,'-m','linewidth',1)
% legend('Full Mixture','Brute-Force','Cost Function of Runnalls')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%This is similar to example 1, but we only consider Runnal's cost function
%and we set K=1 and use gammaVal to control how many elements it ultimately
%has in the end. In this example, it reduced to from 10 to KNew=6
%components.
% w=[0.03, 0.18, 0.12, 0.19, 0.02, 0.16, 0.06, 0.1, 0.08, 0.06];
% mu=[1.45, 2.20, 0.67, 0.48, 1.49, 0.91, 1.01, 1.42, 2.77, 0.89];
% P=reshape([0.0487, 0.0305, 0.1171, 0.0174, 0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679],[1,1,10]);
%
% K=1;
% gammaVal=0.05;
% %Runnals' cost function.
% [wRed,muRed,PRed]=RunnalsGaussMixRed(w,mu,P,K,gammaVal);
% KNew=length(wRed)%The number of reduced components.
% numPoints=500;
% xVals=linspace(0,3,numPoints);
% PDFVals1=GaussianMixtureD.PDF(xVals,w,mu,P);
% PDFVals2=GaussianMixtureD.PDF(xVals,wRed,muRed,PRed);
% figure(1)
% clf
% hold on
% plot(xVals,PDFVals1,'-k','linewidth',4)
% plot(xVals,PDFVals2,'-m','linewidth',1)
% legend('Full Mixture','Reduced Mixture')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] A. R. Runnalls, "Kullback-Leibler approach to Gaussian mixture
%    reduction," IEEE Transactions on Aerospace and Electronic Systems,
%    vol. 43, no. 3, pp. 989-999, Jul. 2007.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<5||isempty(gammaBound))
        gammaBound=Inf;
    end

    if(nargin<6||isempty(KMax))
        KMax=Inf;
    end

    N=length(w);

    %If no reduction is necessary.
    if(N<=K)
        return;
    end

    %We will only be using one triangular of this matrix.
    M=Inf*ones(N,N);%This is the cost matrix.
    wLogDetPVals=zeros(N,1);
    for k=1:N
        wLogDetPVals(k)=w(k)*log(det(P(:,:,k)));
    end

    %We shall fill the cost matrix with the cost of all pairs.
    for cur1=1:(N-1)
        for cur2=(cur1+1):N
            M(cur1,cur2)=BDist(w(cur1),w(cur2),mu(:,cur1),mu(:,cur2),P(:,:,cur1),P(:,:,cur2),wLogDetPVals(cur1),wLogDetPVals(cur2));
        end
    end

    KCur=N;
    selIdxPresent=true(N,1);
    for mergeRound=1:(N-K)
        [Mm,minRows]=min(M);%Minimize over the rows
        [minVal,minCol]=min(Mm);%Minimize over the columns.
        minRow=minRows(minCol);

        if(minVal>gammaBound&&KCur<=KMax)
            %If the distribution has been sufficiently reduced in terms of
            %cost.
            break;
        end

        %Now we know which two hypotheses to merge: The ones with indices
        %minRow and minCol. We will merge those hypotheses and put the
        %results in minRow.
        wSum=w(minRow)+w(minCol);
        w1=w(minRow)/wSum;
        w2=w(minCol)/wSum;
        muMerged=w1*mu(:,minRow)+w2*mu(:,minCol);
        diff1=mu(:,minRow)-muMerged;
        diff2=mu(:,minCol)-muMerged;
        PMerged=w1*(P(:,:,minRow)+diff1*diff1')+w2*(P(:,:,minCol)+diff2*diff2');
        w(minRow)=wSum;
        mu(:,minRow)=muMerged;
        P(:,:,minRow)=PMerged;
        wLogDetPVals(minRow)=wSum*log(det(PMerged));

        %The column is removed.
        selIdxPresent(minCol)=false;

        %Make all the costs Inf so that nothing will be assigned to the
        %removed column.
        M(minCol,:)=Inf;
        M(:,minCol)=Inf;

        %We must now fill in the costs for the merged estimate, which is in
        %minRow. We shall fill the cost matrix with the cost of all pairs.
        for cur1=1:(minRow-1)
            if(selIdxPresent(cur1))
                M(cur1,minRow)=BDist(w(cur1),w(minRow),mu(:,cur1),mu(:,minRow),P(:,:,cur1),P(:,:,minRow),wLogDetPVals(cur1),wLogDetPVals(minRow));
            end
        end

        for cur2=(minRow+1):N
            if(selIdxPresent(cur2))
                M(minRow,cur2)=BDist(w(minRow),w(cur2),mu(:,minRow),mu(:,cur2),P(:,:,minRow),P(:,:,cur2),wLogDetPVals(minRow),wLogDetPVals(cur2));
            end
        end

        KCur=KCur-1;
    end

    %The merged items are marked with selIdxPresent and can now be all
    %grouped together to be returned.
    w=w(selIdxPresent);
    w=w/sum(w);%Guarantee continued normalization.
    mu=mu(:,selIdxPresent);
    P=P(:,:,selIdxPresent);
end

function val=BDist(w1,w2,mu1,mu2,P1,P2,wLogDetP1,wLogDetP2)
%%BDIST This is the distance measure given in Equation 21 in Section VI-B
%       in [1]. This is a bound related the Kullback-Leibler divergence.
%
%REFERENCES:
%[1] A. R. Runnalls, "Kullback-Leibler approach to Gaussian mixture
%    reduction," IEEE Transactions on Aerospce and Electronic Systems,
%    vol. 43, no. 3, pp. 989-999, Jul. 2007.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    wSum=w1+w2;
    w1m=w1/wSum;
    w2m=w2/wSum;

    %In [1], there is a constant factor of (1/2) out front, but since all
    %components have it, we omit it here. Also, the dterminant is a slow
    %step, so we put in explicit determinants for things less than 4D.
    switch(size(mu1,1))
        case 1
            diff=mu1-mu2;
            P=w1m*P1+w2m*P2+w1m*w2m*(diff*diff);
            val=wSum*log(P)-wLogDetP1-wLogDetP2;
        case 2
            diff1=mu1(1)-mu2(1);
            diff2=mu1(2)-mu2(2);
            wmProd=w1m*w2m;
            val=wSum*log((w1m*P1(1,1)+w2m*P2(1,1)+wmProd*(diff1*diff1))*(w1m*P1(2,2)+w2m*P2(2,2)+wmProd*(diff2*diff2))-(w1m*P1(1,2)+w2m*P2(1,2)+wmProd*(diff1*diff2)).^2)-wLogDetP1-wLogDetP2;
        case 3
            diff=mu1-mu2;
            P=w1m*P1+w2m*P2+w1m*w2m*(diff*diff');
            detVal=-P(1,3)*P(2,2)*P(3,1)+P(1,2)*P(2,3)*P(3,1)+P(1,3)*P(2,1)*P(3,2)-P(1,1)*P(2,3)*P(3,2)-P(1,2)*P(2,1)*P(3,3)+P(1,1)*P(2,2)*P(3,3);
            val=wSum*log(detVal)-wLogDetP1-wLogDetP2;
        otherwise
            diff=mu1-mu2;
            P12=w1m*P1+w2m*P2+w1m*w2m*(diff*diff');
            val=wSum*log(det(P12))-wLogDetP1-wLogDetP2;
    end

    %Deal with the case where w1 and w2 are both essentially zero.
    if(~isfinite(val))
        val=0;
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [fa,md,rd,ad] = dtiComputeFA(eigVal)
% Computes fractional anisotropy (FA) from tensor eigenvalues
%
%   [fa,md,rd,ad] = dtiComputeFA(eigVal)
%
% eigVal  XxYxZx3xN (or nx3xN) array of tensor eigenvalues (N is No. of subjects)
% fa      XxYxZxN (or nxN) array of FA values
% md      XxYxZxN (or nxN) array of Mean Diffusivity values (10^-^6 mm^2/sec)
% rd      XxYxZxN (or nxN) array of Radial Diffusivity values
%
% HISTORY:
% 2003.12.04 ASH (armins@stanford.edu) wrote it.
% 2003.12.04 RFD (bob@white.stanford.edu) added check to avoid divide-by-zero.
% 2004.02.19 ASH added extra dimension for subjects
% 2005.05.02 RFD: added option to pass in a dt6 array.
% 2006.07.12 ASH added indexed format
% 2006.07.13 RFD: fixed the features that ASH so rudely clobbered yesterday.
% 2007.02.28 RFD: added rd return arg.

% Programming:  Since this returns MD as an option, how about renaming the
% function.  Or giving it a flag about which to return.  Or something.

% Check inputs
if (ndims(eigVal)==2)
    Ind = 1;    % eigVal in indexed nx3 format
    eigVal = shiftdim(eigVal, -2);
elseif(ndims(eigVal)==3),
    if(size(eigVal,2)==3 && size(eigVal,3)~=3)
        Ind = 1;    % Data in indexed nx3xN format
        eigVal = shiftdim(eigVal, -2);
    else
        Ind = 2;    % Data in XxYx3 format
        eigVal = shiftdim(eigVal, -1);
    end
else
    Ind = 0;    % eigVal in XxYxZx3xN format
end
if (ndims(eigVal)<4 || ndims(eigVal)>5),
    error('Wrong input format');
end

if(size(eigVal,4) ~= 3),
    if(size(eigVal,4)==6)
        [eigVec, eigVal] = dtiEig(eigVal);
        clear eigVec;
    else
    error('Wrong input format');
    end
end

% Main
epsilon = 1e-10;  % Avoid division by zero

% Mean diffusivity
md = sum(eigVal,4)/3;

%
stdevDiffusivity = sqrt(sum((eigVal-repmat(md,[1,1,1,3,1])).^2,4));
normDiffusivity = sqrt(sum(eigVal.^2,4));
nz = normDiffusivity > epsilon;

% Formula for fractional anisotropy
fa = repmat(NaN, size(md));
fa(nz) = sqrt(3/2).*(stdevDiffusivity(nz)./normDiffusivity(nz));
fa = squeeze(fa);

% Mean diffusivity
md = squeeze(md);

% If requested, return rd and ad
if(nargout>2)
    % Radial diffusivity is the average of the two smaller eigenvalues
    rd = squeeze((eigVal(:,:,:,2)+eigVal(:,:,:,3))/2);
end
if(nargout>3)
    % Axial diffusivity - the largest eigenvalue
    ad = squeeze(eigVal(:,:,:,1));
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x,iflaw]=fmpar(N,p1,p2,p3)
%FMPAR	Parabolic frequency modulated signal.
%	[X,IFLAW]=FMPAR(N,P1,P2,P3) generates a signal with
%	parabolic frequency modulation law.
%	X(T) = exp(j*2*pi(A0.T + A1/2.T^2 +A2/3.T^3))
%
%	N  : the number of points in time
%	P1 : if NARGIN=2, P1 is a vector containing the three
%	    coefficients [A0 A1 A2] of the polynomial instantaneous phase.
%	    If NARGIN=4, P1 (as P2 and P3) is a time-frequency point of
%	    the form [Ti Fi].
%	    The coefficients (A0,A1,A2) are then deduced such that
%	    the frequency modulation law fits these three points.
%	P2,P3 : same as P1 if NARGIN=4.       (optional)
%	X     : time row vector containing the modulated signal samples
%	IFLAW : instantaneous frequency law
%
%	Examples :
%	 [X,IFLAW]=fmpar(128,[1 0.4],[64 0.05],[128 0.4]);
%	 subplot(211);plot(real(X));subplot(212);plot(IFLAW);
%	 [X,IFLAW]=fmpar(128,[0.4 -0.0112 8.6806e-05]);
%	 subplot(211);plot(real(X));subplot(212);plot(IFLAW);
%
%	See also FMCONST, FMHYP, FMLIN, FMSIN, FMODANY, FMPOWER.

%	P. Goncalves - October 1995, O. Lemoine - November 1995
%	Copyright (c) 1995 Rice University
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin <= 1),
 error ( 'The number of parameters must be at least 2.' );
elseif (N <= 0),
 error ('The signal length N must be strictly positive' );
elseif nargin == 2 ;
  if length(p1) ~= 3
    error('Bad number of coefficients for P1');
  end
  a0 = p1(1) ; a1 = p1(2) ; a2 = p1(3) ;
elseif nargin == 4 ;
  if (length(p1) ~= 2) |(length(p2) ~= 2) |(length(p3) ~= 2),
    error('Bad number of coefficients for P1, P2, P3');
  end
  if p1(1)>N | p1(1)<1,
   error ('P1(1) must be between 1 and N');
  elseif p2(1)>N | p2(1)<1,
   error ('P2(1) must be between 1 and N');
  elseif p3(1)>N | p3(1)<1,
   error ('P3(1) must be between 1 and N');
  elseif p1(2)<0,
   error ('P1(2) must be > 0');
  elseif p2(2)<0,
   error ('P2(2) must be > 0');
  elseif p3(2)<0,
   error ('P3(2) must be > 0');
  end
  Y = [p1(2) p2(2) p3(2)] ;
  X = [1 1 1;p1(1) p2(1) p3(1);p1(1)^2 p2(1)^2 p3(1)^2] ;
  coef = Y*inv(X) ;
  a0 = coef(1) ;
  a1 = coef(2) ;
  a2 = coef(3) ;
end

t=1:N;

phi = 2*pi*(a0*t + a1/2*t.^2 + a2/3*t.^3) ;
iflaw = (a0 + a1*t + a2*t.^2).' ;

aliasing = find(iflaw<0 | iflaw>0.5) ;
if isempty(aliasing) == 0
  disp(['!!! WARNING: signal is undersampled or has negative frequencies']) ;
end

x = exp(i*phi).';

author: HeLiangHIT


text:
% Script: call_Omorifit.m
% -----------------------
% Script to input parameters for a modified Omori law fit and calculate fit to data
% Calculation is done in function plot_llkstest
% The function works on the catalog newt2!!!
% J.Woessner
% last update: 19.08.03

report_this_filefun(mfilename('fullpath'));
% Get input parameters
prompt  = {,'Enter length of learning period (days)','Enter number of bootstraps:'};
title   = 'Parameters ';
lines= 1;
if exist('time') &&  exist('bootloops')
    time = num2str(time);
    bootloops = num2str(bootloops);
    def     = {time,bootloops};
    answer  = inputdlg(prompt,title,lines,def);
    time = str2double(answer{1});
    timef = 1;
    bootloops = str2double(answer{2});

else
    def     = {'50','50'};
    answer  = inputdlg(prompt,title,lines,def);
    time = str2double(answer{1});
    timef = 1;
    bootloops = str2double(answer{2});
end

% maepi is the mainshock of sequence
plot_llkstest(newt2,time,timef,bootloops,maepi);

author: CelsoReyes


text:
function d = getdict
% Dictionary of NIFTI stuff
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

persistent dict;
if ~isempty(dict),
    d = dict;
    return;
end;

% Datatype
t = true;
f = false;
table = {...
    0   ,'UNKNOWN'   ,'uint8'   ,@uint8  ,1,1  ,t,t,f
    1   ,'BINARY'    ,'uint1'   ,@logical,1,1/8,t,t,f
    256 ,'INT8'      ,'int8'    ,@int8   ,1,1  ,t,f,t
    2   ,'UINT8'     ,'uint8'   ,@uint8  ,1,1  ,t,t,t
    4   ,'INT16'     ,'int16'   ,@int16  ,1,2  ,t,f,t
    512 ,'UINT16'    ,'uint16'  ,@uint16 ,1,2  ,t,t,t
    8   ,'INT32'     ,'int32'   ,@int32  ,1,4  ,t,f,t
    768 ,'UINT32'    ,'uint32'  ,@uint32 ,1,4  ,t,t,t
    1024,'INT64'     ,'int64'   ,@int64  ,1,8  ,t,f,f
    1280,'UINT64'    ,'uint64'  ,@uint64 ,1,8  ,t,t,f
    16  ,'FLOAT32'   ,'float32' ,@single ,1,4  ,f,f,t
    64  ,'FLOAT64'   ,'double'  ,@double ,1,8  ,f,f,t
    1536,'FLOAT128'  ,'float128',@crash  ,1,16 ,f,f,f
    32  ,'COMPLEX64' ,'float32' ,@single ,2,4  ,f,f,f
    1792,'COMPLEX128','double'  ,@double ,2,8  ,f,f,f
    2048,'COMPLEX256','float128',@crash  ,2,16 ,f,f,f
    128 ,'RGB24'     ,'uint8'   ,@uint8  ,3,1  ,t,t,f};

dtype = struct(...
    'code'     ,table(:,1),...
    'label'    ,table(:,2),...
    'prec'     ,table(:,3),...
    'conv'     ,table(:,4),...
    'nelem'    ,table(:,5),...
    'size'     ,table(:,6),...
    'isint'    ,table(:,7),...
    'unsigned' ,table(:,8),...
    'min',-Inf,'max',Inf',...
    'supported',table(:,9));
for i=1:length(dtype),
    if dtype(i).isint
        if dtype(i).unsigned
            dtype(i).min =  0;
            dtype(i).max =  2^(8*dtype(i).size)-1;
        else
            dtype(i).min = -2^(8*dtype(i).size-1);
            dtype(i).max =  2^(8*dtype(i).size-1)-1;
        end;
    end;
end;
% Intent
table = {...
    0   ,'NONE'         ,'None',{}
    2   ,'CORREL'       ,'Correlation statistic',{'DOF'}
    3   ,'TTEST'        ,'T-statistic',{'DOF'}
    4   ,'FTEST'        ,'F-statistic',{'numerator DOF','denominator DOF'}
    5   ,'ZSCORE'       ,'Z-score',{}
    6   ,'CHISQ'        ,'Chi-squared distribution',{'DOF'}
    7   ,'BETA'         ,'Beta distribution',{'a','b'}
    8   ,'BINOM'        ,'Binomial distribution',...
        {'number of trials','probability per trial'}
    9   ,'GAMMA'        ,'Gamma distribution',{'shape','scale'}
    10  ,'POISSON'      ,'Poisson distribution',{'mean'}
    11  ,'NORMAL'       ,'Normal distribution',{'mean','standard deviation'}
    12  ,'FTEST_NONC'   ,'F-statistic noncentral',...
        {'numerator DOF','denominator DOF','numerator noncentrality parameter'}
    13  ,'CHISQ_NONC'   ,'Chi-squared noncentral',{'DOF','noncentrality parameter'}
    14  ,'LOGISTIC'     ,'Logistic distribution',{'location','scale'}
    15  ,'LAPLACE'      ,'Laplace distribution',{'location','scale'}
    16  ,'UNIFORM'      ,'Uniform distribition',{'lower end','upper end'}
    17  ,'TTEST_NONC'   ,'T-statistic noncentral',{'DOF','noncentrality parameter'}
    18  ,'WEIBULL'      ,'Weibull distribution',{'location','scale','power'}
    19  ,'CHI'          ,'Chi distribution',{'DOF'}
    20  ,'INVGAUSS'     ,'Inverse Gaussian distribution',{'mu','lambda'}
    21  ,'EXTVAL'       ,'Extreme Value distribution',{'location','scale'}
    22  ,'PVAL'         ,'P-value',{}
    23  ,'LOGPVAL'      ,'Log P-value',{}
    24  ,'LOG10PVAL'    ,'Log_10 P-value',{}
    1001,'ESTIMATE'     ,'Estimate',{}
    1002,'LABEL'        ,'Label index',{}
    1003,'NEURONAMES'   ,'NeuroNames index',{}
    1004,'MATRIX'       ,'General matrix',{'M','N'}
    1005,'MATRIX_SYM'   ,'Symmetric matrix',{}
    1006,'DISPLACEMENT' ,'Displacement vector',{}
    1007,'VECTOR'       ,'Vector',{}
    1008,'POINTS'       ,'Pointset',{}
    1009,'TRIANGLE'     ,'Triangle',{}
    1010,'QUATERNION'   ,'Quaternion',{}
    1011,'DIMLESS'      ,'Dimensionless',{}
};
intent = struct('code',table(:,1),'label',table(:,2),...
    'fullname',table(:,3),'param',table(:,4));

% Units
table = {...
     0,   1,'UNKNOWN'
     1,1000,'m'
     2,   1,'mm'
     3,1e-3,'um'
     8,   1,'s'
    16,1e-3,'ms'
    24,1e-6,'us'
    32,   1,'Hz'
    40,   1,'ppm'
    48,   1,'rads'};
units = struct('code',table(:,1),'label',table(:,3),'rescale',table(:,2));

% Reference space
% code  = {0,1,2,3,4};
table = {...
    0,'UNKNOWN'
    1,'Scanner Anat'
    2,'Aligned Anat'
    3,'Talairach'
    4,'MNI_152'};
anat  = struct('code',table(:,1),'label',table(:,2));

% Slice Ordering
table = {...
    0,'UNKNOWN'
    1,'sequential_increasing'
    2,'sequential_decreasing'
    3,'alternating_increasing'
    4,'alternating_decreasing'};
sliceorder = struct('code',table(:,1),'label',table(:,2));

% Q/S Form Interpretation
table = {...
    0,'UNKNOWN'
    1,'Scanner'
    2,'Aligned'
    3,'Talairach'
    4,'MNI152'};
xform = struct('code',table(:,1),'label',table(:,2));

dict = struct('dtype',dtype,'intent',intent,'units',units,...
    'space',anat,'sliceorder',sliceorder,'xform',xform);

d = dict;
return;

function varargout = crash(varargin)
error('There is a NIFTI-1 data format problem (an invalid datatype).');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function histo_gram = r8vec_histogram ( n, a, a_lo, a_hi, histo_num )

%*****************************************************************************80
%
%% R8VEC_HISTOGRAM computes a histogram of the elements of an R8VEC.
%
%  Discussion:
%
%    Values between A_LO and A_HI will be histogrammed into the bins
%    1 through HISTO_NUM.  Values below A_LO are counted in bin 0,
%    and values greater than A_HI are counted in bin HISTO_NUM+1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Input, real A(N), the array to examine.
%
%    Input, real A_LO, A_HI, the lowest and highest
%    values to be histogrammed.  These values will also define the bins.
%
%    Input, integer HISTO_NUM, the number of bins to use.
%
%    Output, integer HISTO_GRAM(0:HISTO_NUM+1), contains the number of
%    entries of A in each bin.
%
  histo_gram(1:histo_num+2) = 0;

  delta = ( a_hi - a_lo ) / ( 2 * histo_num );

  for i = 1 : n

    if ( a(i) < a_lo )

      histo_gram(1) = histo_gram(2) + 1;

    elseif ( a(i) <= a_hi )

      j = round ( ...
        ( ( a_hi -       delta - a(i)        ) * ( 1         )   ...
        + (      -       delta + a(i) - a_lo ) * ( histo_num ) ) ...
        / ( a_hi - 2.0 * delta        - a_lo ) );

      histo_gram(j+1) = histo_gram(j+1) + 1;

    elseif ( a_hi < a(i) )

      histo_gram(histo_num+2) = histo_gram(histo_num+2) + 1;

    end

  end

  return
end


text:
function flip = spm_flip_analyze_images
% Do Analyze format images need to be left-right flipped? The default
% behavior is to have the indices of the voxels stored as left-handed and
% interpret the mm coordinates within a right-handed coordinate system.
%
% Note: the behavior used to be set in spm_defaults.m, but this has now
% been changed.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_flip_analyze_images.m 1418 2008-04-15 19:15:16Z john $

flip = 1;


text:
function [L1t,L2t,L3t]=ongrid(V1,V2,V3,cm,vm)

L1t=ones(size(V1,1),size(cm,1));
L2t=ones(size(V2,1),size(cm,1));
L3t=ones(size(V3,1),size(cm,1));

for i=1:1:size(cm,1)

    c=cm(i,:);
    v=vm(i,:);

    [X1,Y1,Z1]=snap2grid(V1(:,1),V1(:,2),V1(:,3),c,v);
    X1=(X1-c(1))./v(1); Y1=(Y1-c(2))./v(2); Z1=(Z1-c(3))./v(3);
    V1s=[X1(:) Y1(:) Z1(:)];

    [X2,Y2,Z2]=snap2grid(V2(:,1),V2(:,2),V2(:,3),c,v);
    X2=(X2-c(1))./v(1); Y2=(Y2-c(2))./v(2); Z2=(Z2-c(3))./v(3);
    V2s=[X2(:) Y2(:) Z2(:)];

    [X3,Y3,Z3]=snap2grid(V3(:,1),V3(:,2),V3(:,3),c,v);
    X3=(X3-c(1))./v(1); Y3=(Y3-c(2))./v(2); Z3=(Z3-c(3))./v(3);
    V3s=[X3(:) Y3(:) Z3(:)];

    Vsnap=[V1s;V2s;V3s];
    Vsnap=1+(Vsnap-ones(size(Vsnap,1),1)*min(Vsnap,[],1));
    vsiz=max(Vsnap,[],1);
    Vsnap1=Vsnap(1:size(V1,1),:);
    Vsnap2=Vsnap(size(V1,1)+1:size(V1,1)+size(V2,1),:);
    Vsnap3=Vsnap(size(V1,1)+size(V2,1)+1:size(V1,1)+size(V2,1)+size(V3,1),:);

    INDsnap1=sub2ind(vsiz,Vsnap1(:,1),Vsnap1(:,2),Vsnap1(:,3));
    INDsnap2=sub2ind(vsiz,Vsnap2(:,1),Vsnap2(:,2),Vsnap2(:,3));
    INDsnap3=sub2ind(vsiz,Vsnap3(:,1),Vsnap3(:,2),Vsnap3(:,3));

    L1t(:,i)=ismembc(INDsnap1,sort(INDsnap2)) & ismembc(INDsnap1,sort(INDsnap3));
    L2t(:,i)=ismembc(INDsnap2,sort(INDsnap1)) & ismembc(INDsnap2,sort(INDsnap3));
    L3t(:,i)=ismembc(INDsnap3,sort(INDsnap1)) & ismembc(INDsnap3,sort(INDsnap2));

end

L1t=any(L1t,2);
L2t=any(L2t,2);
L3t=any(L3t,2);

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%SampleCnt calculates number of samples needed to be done

function SampleCnt = nsamples(ni, ptNum, pf, conf)
q  = prod ([(ni-pf+1) : ni] ./ [(ptNum-pf+1) : ptNum]);

if (1 -q) < eps
   SampleCnt = 1;
else
   SampleCnt  = log(1 - conf) / log(1 - q);
end

if SampleCnt < 1
   SampleCnt = 1;
end

author: strawlab

repo: MultiCamSelfCal


text:
function varargout = transpose(varargin)
%.'   DISKFUN transpose.
%   F.' is the non-conjugate transpose of the underlying quasimatrix for
%   the DISKFUN object F. Given F in polar coordinates, F(theta, r),
%   F.' = F(r, theta). To transpose a DISKFUN with respect to Cartesian
%   coordinates, use the command FLIPXY.
%
% See also DISKFUN/CTRANSPOSE and DISKFUN/FLIPXY.

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

[varargout{1:nargout}] = transpose@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function image = autoCropImage(image)
% input an image
% output an image with white space removed
% assuming the top left pixel is the background color that you don't want

mask = mean(image,3);
mask = mask == mask(1,1);
isGood1 = find(any(~mask,1));
isGood2 = find(any(~mask,2));
image = image(min(isGood2):max(isGood2),min(isGood1):max(isGood1), :);

author: jianxiongxiao

repo: ProfXkit


text:
function catalan_row_next_test ( )

%*****************************************************************************80
%
%% CATALAN_ROW_NEXT_TEST tests CATALAN_ROW_NEXT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CATALAN_ROW_NEXT_TEST\n' );
  fprintf ( 1, '  CATALAN_ROW_NEXT computes a row of Catalan''s triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  First, compute row 7:\n' );

  ido = 0;
  i = 7;
  c = [];
  c = catalan_row_next ( ido, i, c );
  fprintf ( 1, '  %2d  ', i );
  for j = 0 : i
    fprintf ( 1, '  %4d', c(j+1) );
  end
  fprintf ( 1, '\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now compute rows one at a time:\n' );
  fprintf ( 1, '\n' );

  ido = 0;
  c = [];

  for i = 0 : n
    c = catalan_row_next ( ido, i, c );
    ido = 1;
    fprintf ( 1, '  %2d  ', i );
    for j = 0 : i
      fprintf ( 1, '  %4d', c(j+1) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function []= sounding_powerlaw_sens(hydro,wet)
% function that computes the power law coefficents and performs a
% sensitivity analysis.
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By David Bekaert -- University of Leeds 2013
% modifications
% 05/2013   DB:     Compute the power law coefficients for dedicated days.
% 08/2013   DB:     Fix the case when a netdelay cannot be computed.
% 02/2014   DB:     Add non-stamps file support of ifgs-based estimation.
% 03/2014   DB:     Fix the computation of the sensitivity analysis in
%                   absence of soundign for specific SAR dates.
% 08/2014   DB:     Expand for PS network and SB network for SAR date estimation
% 03/2016   DB:     remove warning on ifgday for non-stamps processing

if nargin<1
    hydro=1;
    wet=1;
end

% getting the data from the parms_aps file
n_months = getparm_aps('sounding_months');
sounding_start_date = getparm_aps('sounding_start_date');
sounding_end_date = getparm_aps('sounding_end_date');
sounding_dir=getparm_aps('sounding_dir');
time_stamp = getparm_aps('sounding_time_stamp');
sounding_ifg_dates = getparm_aps('sounding_ifg_dates');
time_stamp_str = [];
for k=1:size(time_stamp,1)
    if k>1
        time_stamp_str = [time_stamp_str '_' time_stamp(k,:)];
    else
        time_stamp_str = [time_stamp(k,:)];
    end
end

% checking if the directory exists
if exist([sounding_dir filesep],'dir')~=7
        error('myApp:argChk', ['The specified filepath of the sounding data does not exist,...  \nAbort,... \n'])
end
% current directry
curdir = pwd;

% getting the interferogram data information
% loading the data
stamps_processed = getparm_aps('stamps_processed');
if strcmp(stamps_processed,'y')
   fprintf('Stamps processed structure \n')
   ll_matfile = getparm_aps('ll_matfile');
   ps = load(ll_matfile);
   dates = ps.day;
   lonlat = ps.lonlat;

    % getting the parms file list from stamps to see the final ifg list

    % constructing the matrix with master and slave dates
    if strcmp(getparm('small_baseline_flag'),'y')
        % for SB
        ifgs_ix = ps.ifgday_ix;

    else
        n_ifg = ps.n_ifg;
        % slightly different for PS.
        date_slave_ix = [1:n_ifg]';
        % the master dates
        date_master_ix = repmat(ps.master_ix,size(date_slave_ix,1),1);
        % ix interferograms
        ifgs_ix = [date_master_ix date_slave_ix];
    end
else
    % getting the dates in jullian format
    ifgday_matfile = getparm_aps('ifgday_matfile');
    ifgs_dates = load(ifgday_matfile);
    ifgs_dates = ifgs_dates.ifgday;
    dates = reshape(ifgs_dates,[],1);
    dates = unique(dates);
    dates = datenum(num2str(dates),'yyyymmdd');
    dates = sort(dates);        % dates increasing with time

    % getting the ix position for the master and slave dates with respect
    % to the times
    date_master = datenum(num2str(ifgs_dates(:,1)),'yyyymmdd');
    date_slave = datenum(num2str(ifgs_dates(:,2)),'yyyymmdd');

    for k=1:size(date_master,1)
        [date_master_ix(k,1)] = find(date_master(k,1)==dates);
        [date_slave_ix(k,1)] = find(date_slave(k,1)==dates);
    end

    % ix interferograms
    ifgs_ix = [date_master_ix date_slave_ix];
    clear date_master_ix date_slave_ix
end

if strcmp(sounding_ifg_dates,'y')
    % estimate for interferogram dates
    stamps_processed = getparm_aps('stamps_processed');
    if strcmp(stamps_processed,'y')
        ps = load(ll_matfile);
        date_start_vector = datestr(ps.day-15,'yyyymmdd');
        date_end_vector = datestr(ps.day+15,'yyyymmdd');

    else

        ifgday_matfile = getparm_aps('ifgday_matfile');
        ifgs_dates = load(ifgday_matfile);
        ifgs_dates = ifgs_dates.ifgday;
        dates = reshape(ifgs_dates,[],1);
        dates = unique(dates);
        dates = datenum(num2str(dates),'yyyymmdd');
        dates = sort(dates);        % dates increasing with time

        date_start_vector =  datestr(dates-15,'yyyymmdd');
        date_end_vector =  datestr(dates+15,'yyyymmdd');
    end

    % save name of the output data
    if hydro==1 && wet==0
        save_name = ['Powerlaw_sensitivity_hydro_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
    elseif hydro==0 && wet==1
        save_name = ['Powerlaw_sensitivity_wet_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
    else
        save_name = ['Powerlaw_sensitivity_SAR_dates_1month_' time_stamp_str 'Hr.mat'  ];
    end
    if exist([sounding_dir filesep 'Powerlaw'],'dir')~=7
        mkdir([sounding_dir filesep 'Powerlaw']);
    end

else
    % estimate is on fixed intervals
    % in case no start or end time is given, get it from the data files
    if isempty(sounding_start_date) || isempty(sounding_end_date)

        if exist('sounding.list','file')~=2
            % making a list of all the sounding files
            [dummy dummy2] = system('echo sounding_list > sounding.list');
            clear dummy dummy2
            for k=1:size(time_stamp,1)
                command_str = ['ls [0-9]???????_' time_stamp(k,:) '.mat >> sounding.list'];
                [dummy dummy2] = system(command_str);
                clear dummy dummy2
            end
        end
        temp = tdfread('sounding.list');
        [dummy dummy2] = system('rm sounding.list');
        clear dummy dummy2
        date_list_temp = temp.sounding_list(:,[1:8]);

        % selecting a date range when requested
        clear ix

        if isempty(start_date)
            sounding_start_date = date_list_temp(1,:);
        end
        if isempty(end_date)
            sounding_end_date = date_list_temp(end,:);
        end
        clear date_list_temp
    end

    % putting the variables in the right set-up
    start_year = str2num(sounding_start_date(1:4));
    end_year = str2num(sounding_end_date(1:4));
    start_str = sounding_start_date(5:6);
    end_str = sounding_end_date(5:6);
    start_month = str2num(start_str);
    end_month = str2num(end_str);

    % save name of the output data
   if hydro==1 && wet==0
        save_name = ['Powerlaw_sensitivity_hydro_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
    elseif hydro==0 && wet==1
        save_name = ['Powerlaw_sensitivity_wet_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
   else
        save_name = ['Powerlaw_sensitivity_' num2str(n_months) 'month_' time_stamp_str 'Hr_' num2str(start_year) start_str '_' num2str(end_year) end_str '.mat'  ];
    end
    if exist([sounding_dir filesep 'Powerlaw'],'dir')~=7
        mkdir([sounding_dir filesep 'Powerlaw']);
    end

    if exist([sounding_dir filesep 'Powerlaw'],'dir')~=7
        mkdir([sounding_dir filesep 'Powerlaw']);
    end

    % generating the periods
    month_str = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12'];
    counter = 1;
    % runnign the computation of the powerlaw in batches
    for k=1:end_year-start_year+1
        if k==end_year-start_year+1 & k>1
            for l=1:end_month
                date_start_temp = datenum([num2str(start_year+k-1)  month_str(l,:) '01'],'yyyymmdd');
                date_start_month_vector(counter,:) = datestr(date_start_temp,'yyyymmdd');
                counter = counter+1;
            end
        elseif k==end_year-start_year+1 & k==1
            for l=start_month:end_month
                date_start_temp = datenum([num2str(start_year+k-1)  month_str(l,:) '01'],'yyyymmdd');
                date_start_month_vector(counter,:) = datestr(date_start_temp,'yyyymmdd');
                counter = counter+1;
            end
        else
            for l=1:12
                date_start_temp = datenum([num2str(start_year+k-1)  month_str(l,:) '01'],'yyyymmdd');
                date_start_month_vector(counter,:) = datestr(date_start_temp,'yyyymmdd');
                counter = counter+1;
            end
        end
    end

    n_months_total = size(date_start_month_vector,1);
    for k=1:ceil(n_months_total/n_months)
        k_lower = n_months*(k-1)+1;
        k_upper = n_months*(k)+1;

        if k_upper>n_months_total
            k_upper = n_months_total;
        end
        date_start_vector(k,:)=date_start_month_vector(k_lower,:);
        date_end_vector(k,:) = datestr(datenum(date_start_month_vector(k_upper,:),'yyyymmdd')-1,'yyyymmdd');
    end

end

% Remove those months outside the users request
if ~isempty(sounding_start_date)
    ix_drop = find(datenum(date_start_vector,'yyyymmdd')<datenum(sounding_start_date,'yyyymmdd'));
else
    ix_drop = [];
end
if ~isempty(sounding_end_date)
    ix_drop = [ix_drop ; find(datenum(date_end_vector,'yyyymmdd')>datenum(sounding_end_date,'yyyymmdd'))];
end
ix_drop = unique(ix_drop);
date_end_vector(ix_drop,:)=[];
date_start_vector(ix_drop,:)=[];

% Computing the power law coefficients
for k=1:size(date_start_vector,1)
    fprintf(['\n' num2str(k) '/' num2str(size(date_start_vector,1)) ' completed \n']);
	% When having sounding data estimate the powerlaw and b coefficients

    [alpha_all,alpha_hc,h_0_threshold,n_soundings] = sounding(date_start_vector(k,:),date_end_vector(k,:),[],hydro,wet);
	alpha_vector_all(k,1) = alpha_all;
    alpha_vector(k,1) = alpha_hc;
	n_soundings_vector(k,1)=n_soundings;
	h0_vector(k,1)=h_0_threshold;

end
save([sounding_dir filesep 'Powerlaw' filesep save_name],'alpha_vector_all','alpha_vector','n_soundings_vector','h0_vector','date_start_vector','date_end_vector')

cd(curdir)
if strcmp(sounding_ifg_dates,'y')
    % estimate for interferogram dates
    stamps_processed = getparm_aps('stamps_processed');
    if strcmp(stamps_processed,'y')
        fprintf('Updating powerlaw parameters with new values. \n')

        % removing NaN by replacing them with the other SAR date estimates
        ix_alpha_fix = find(isnan(alpha_vector(:,1)));
        if ~isempty(ix_alpha_fix)
           if length(ix_alpha_fix)~=length(alpha_vector)
                fprintf([num2str(length(ix_alpha_fix)) ' out of ' num2str(length(alpha_vector)) ' SAR dates where set to the mean alpha due to lack of sounding data. \n'])
                alpha_vector(isnan(alpha_vector))=nanmean(alpha_vector);
                alpha_vector_fix = alpha_vector;
                save([sounding_dir filesep 'Powerlaw' filesep save_name],'-append','alpha_vector_fix','ix_alpha_fix')

           else
              fprintf('None of the SAR dates where estimate using sounding data \n')
           end

        end

        ix_height_fix = find(isnan(h0_vector(:,1)));
        if ~isempty(ix_height_fix)
           if length(ix_height_fix)~=length(h0_vector)
                fprintf([num2str(length(ix_height_fix)) ' out of ' num2str(length(alpha_vector)) ' SAR dates where set to the mean h0 due to lack of sounding data. \n'])
                h0_vector(isnan(h0_vector))=nanmean(h0_vector);
                h0_vector_fix = h0_vector;
                save([sounding_dir filesep 'Powerlaw' filesep save_name],'-append','h0_vector_fix','ix_height_fix')
           else
              fprintf('None of the SAR dates where estimate using sounding data \n')
           end
        end

        alpha_SAR = [alpha_vector(ifgs_ix(:,1)) alpha_vector(ifgs_ix(:,2))];
        h0_SAR = [h0_vector(ifgs_ix(:,1)) h0_vector(ifgs_ix(:,2))];

        % computing the mean between two SAR dates
        h0_InSAR = mean(h0_SAR,2);
        alpha_InSAR = mean(alpha_SAR,2);

        setparm_aps('powerlaw_h0',h0_InSAR');
        setparm_aps('powerlaw_alpha',alpha_InSAR');
    end
else

    % computing the mean between two SAR dates
    ix = isnan(h0_vector);
    powerlaw_h0 = mean(h0_vector(~ix));
    powerlaw_alpha = mean(alpha_vector(~ix));

    setparm_aps('powerlaw_h0',powerlaw_h0');
    setparm_aps('powerlaw_alpha',powerlaw_alpha');
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


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

%%
clear; close all; clc;

%% Syntax
% |writeMtrFile_tetGen(inputStruct,bOpt);|

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


text:
function [fAphi] = calc_FaultStyle(rStress)
    % Determine the faulting style based on Quantifying Anderson's fault type
% [fAphi] = calc_FaultStyle(rStress)
% ------------------------------------------
% Determining the faulting style based on Quantifying Anderson's fault type
% by B. Simpson, 1997, JGR
%
% Incoming:
% rStress : Structure from script gui_NodeCalcStressInv giving stress
%           tensor inversion by Michael
% Output:
% fAphi   : 0 <= fAphi < 1 Normal faulting
%           1 <= fAphi < 2 Strike slip
%           2 <= fAphi < 3 Thrust faulting
%
% jowoe@gps.caltech.edu
% 19.05.2006

% Stress tensor from Michael inversion
mS=[rStress.fS11 rStress.fS12 rStress.fS13; rStress.fS12 rStress.fS22 rStress.fS23; rStress.fS13 rStress.fS23 rStress.fS33];

% Eigenvector; sorted from max. compressive to minimum (S1, S2, S3)
% Maximum compressive is the most negative value
vEig = eig(mS);

% Check for "most" vertical axis
vPlunge = [rStress.fS1Plunge; rStress.fS2Plunge; rStress.fS3Plunge];

vSel = (vPlunge == max(vPlunge));
fSigvert = vEig(vSel);

% This is because maximum compressive is negative
fSighmax = min(vEig(~vSel,:));
fSighmin = max(vEig(~vSel,:));

if fSighmax > fSigvert
    n = 0;
elseif fSigvert > fSighmin
    n = 2;
else
    n = 1;
end
% Formula 2, Simpson, JGR, 1997
fAphi = (n+0.5)+(-1)^n*(rStress.fPhi-0.5);


text:
% This function normalize a set of vectors
% Parameters:
%   v     the set of vectors to be normalized (column stored)
%   nr    the norm for which the normalization is performed (Default: Euclidean)
%   rval  replace value in case the vector is 0-norm
%
% Output:
%   vout  the normalized vector
%   vnr   the norms of the input vectors
%
% Remark: the function return Nan for vectors of null norm
function [vout, vnr] = yael_vecs_normalize (v, nr, rval)

if nargin < 2, nr = 2; end

% norm of each column
vnr = (sum (v.^nr)) .^ (1 / nr);

% sparse multiplication to apply the norm
vout = bsxfun (@times, v, 1 ./ vnr);

if exist('rval')
  [~, ko] = find (isnan (vout));
  ko = unique (ko);
  vout (:, ko) = rval;
end

author: hpatches

repo: hpatches-benchmark


text:
function V = spm_vol(P)
% Get header information for images
% FORMAT V = spm_vol(P)
% P - a char or cell array of filenames
% V - a structure array containing image volume information
%     The elements of the structures are:
%       V.fname - the filename of the image.
%       V.dim   - the x, y and z dimensions of the volume
%       V.dt    - A 1x2 array.  First element is datatype (see spm_type).
%                 The second is 1 or 0 depending on the endian-ness.
%       V.mat   - a 4x4 affine transformation matrix mapping from
%                 voxel coordinates to real world coordinates.
%       V.pinfo - plane info for each plane of the volume.
%              V.pinfo(1,:) - scale for each plane
%              V.pinfo(2,:) - offset for each plane
%                 The true voxel intensities of the jth image are given
%                 by: val*V.pinfo(1,j) + V.pinfo(2,j)
%              V.pinfo(3,:) - offset into image (in bytes).
%                 If the size of pinfo is 3x1, then the volume is assumed
%                 to be contiguous and each plane has the same scalefactor
%                 and offset.
%__________________________________________________________________________
%
% The fields listed above are essential for the mex routines, but other
% fields can also be incorporated into the structure.
%
% Note that spm_vol can also be applied to the filename(s) of 4-dim
% volumes. In that case, the elements of V will point to a series of 3-dim
% images.
%__________________________________________________________________________
% Copyright (C) 1999-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_vol.m 5958 2014-04-16 17:13:41Z guillaume $

if ~nargin
    V = struct('fname',   {},...
               'dim',     {},...
               'dt',      {},...
               'pinfo',   {},...
               'mat',     {},...
               'n',       {},...
               'descrip', {},...
               'private', {});

elseif isempty(P)
    V = spm_vol;
    if iscell(P), V = {V}; end

elseif isstruct(P)
    V = P;

elseif iscell(P)
    V = cellfun(@spm_vol,P, 'UniformOutput',false);

else
    V = spm_vol;
    cnt = 0;
    for i=1:size(P,1)
        v = spm_vol_hdr(deblank(P(i,:)));

        f = fieldnames(v);
        for j=1:numel(f)
            [V(cnt+1:cnt+size(v,2),1).(f{j})] = deal(v.(f{j}));
        end
        cnt = cnt + size(v,2);
    end
end

%==========================================================================
% function V = spm_vol_hdr(p)
%==========================================================================
function V = spm_vol_hdr(p)
[pth,nam,ext,n] = spm_fileparts(p);
p = fullfile(pth,[nam ext]);
n = str2num(n);
if ~spm_existfile(p)
    error('File "%s" does not exist.', p);
end
switch ext
    case {'.nii','.NII'}
        % Do nothing

    case {'.img','.IMG'}
        if ~spm_existfile(fullfile(pth,[nam '.hdr'])) && ...
           ~spm_existfile(fullfile(pth,[nam '.HDR']))
            error('File "%s" does not exist.', fullfile(pth,[nam '.hdr']));
        end

    case {'.hdr','.HDR'}
        ext = '.img';
        p   = fullfile(pth,[nam ext]);
        if ~spm_existfile(p)
            error('File "%s" does not exist.', p);
        end

    case {'.gz','.GZ'}
        fprintf('Compressed NIfTI files are not supported.\n');
        tmpname = tempname;
        try
            tmpname = char(gunzip(p,tmpname));
        catch
            error('Cannot uncompress "%s".',p);
        end
        try
            if isempty(n), n = ''; else n = [',' num2str(n)]; end
            V = spm_vol_hdr([tmpname n]);
            for i=1:numel(V)
                V(i).dat   = spm_read_vols(V(i));
                V(i).private.dat.fname = spm_file(p,'ext','');
                V(i).fname = p;
                V(i).dt(1) = 64;
                V(i).pinfo = [1 0 0]';
            end
        catch
            warning('Cannot read uncompressed file "%s".',p);
        end
        spm_unlink(tmpname);
        rmdir(fileparts(tmpname));
        return

    otherwise
        error('File "%s" is not of a recognised type.', p);
end

V = spm_vol_nifti(p,n);

if isempty(n) && length(V.private.dat.dim) > 3
    V0(1) = V;
    for i = 2:V.private.dat.dim(4)
        V0(i) = spm_vol_nifti(V.private, i);
    end
    V = V0;
end


text:
function A5=inflowclassstat1(file)
%INFLOWCLASSSTAT1 berfungsi untuk mendistribusikan matriks inflow dalam
%berbagai kelas
%INFLOWCLASSSTAT1 dipergunakan ketika memilih mendistribusikan kelas dalam
%format statistikal
%dalam hal ini kelas dipastikan hanya 3 buah
A1=load(file); %untuk membuka file txt disimpan di A1
[m,n]=size(A1);
%penentuan rerata
A2=mean(A1);
A3=std(A1);
%penentuan nilai maksimum per kelas
A4(1,:)=A2-0.5*A3;
A4(2,:)=A2+0.5*A3;
A4(3,:)=max(A1);
for a=1:m
    for b=1:n
        c=1;
        for d=1:3+1
        if A1(a,b)<=A4(c,b)
           A5(a,b)=c;
        else
           c=d;
        end
        end
    end
end
A4;
A5

author: Sable

repo: mcbench-benchmarks


text:
function [gx] = g_BSL(x,phi,u,in)
% observation function for a Bayesian sequence-learner (BSL)
% function [gx] = g_BSL(x,theta,u,in)
% BSL is guessing the next outcome y_t based upon P(y_t=1|y_{t-1}), ie her
% bet actually depends upon the past sequence of outcomes.
% IN:
%   - x: sufficient statistics of log-odds of P(o=1):
%       x(1:2^K)= E[log-odds]
%       x((2^K)+1:2^(K+1))= log V[log-odds]
%   - phi: phi(1) = log-temperature and phi(2) = bias
%   - u: u(1:K)= sequence of K past outcomes
%   - in: depth of sequence learning
% OUT:
%   - gx: P(y_t=1|y_{t-1})

if VBA_isWeird (u) % e.g., 1st trial
    gx = 0.5;
    return
end

a = 0.36; % for E[s(x)] when x~n(mu,Sig)
K = in.K; % sequence depth
% yb = u(2:K+1); % previous outcomes
yb = u(1:K); % previous outcomes
if K >0
    indSeq = bin2dec(num2str(yb'))+1; % index of sequence of previous outcomes
else
    indSeq = 1;
end
m = x(indSeq);
v = exp(x((2^K)+indSeq));
gx = VBA_sigmoid(phi(2)+exp(phi(1)).*m./sqrt(1+a*v)); % E[sigm(log-odds of P(y))]

author: MBB-team

repo: VBA-toolbox