text:
% This function recursively cuts a sentence into shorter segments no longer
% than a predefined maximum length, based purely on the information about
% short pause durations and locations.
% Inputs:
%   sent_len: length of the sentence in terms of frames
%   sp_len: length of the short pauses in the sentence
%   middle_frame: middle of the short pauses
%   max_seg_len: maximum length of every resulting segments.
%   weight_location: the weight of short pause location. This is relative to the
%       short pause duration
% Outputs:
%   cutting_point: an array of suggested cutting points
%
% Authors: Xiong Xiao
% Date Created: 2009
% Last Modified: 16 April 2014
%
function cutting_point = choose_cutting_points_by_sp(sent_len, sp_len, middle_frame, max_seg_len, loc_weight)
if nargin<5
    loc_weight = 1;
end
if sent_len < max_seg_len
    cutting_point = [];
elseif length(sp_len)<2
    cutting_point = [];
else
    % Choose the sp for splitting
    % Two rules to select a sp: 1) longer sp segments should be
    % given higher priority; 2) sp near to center of the sentence
    % should be given higher priority. The weight of these two
    % criteria will be adjsuted by a weighting loc_weight
    sent_center = round(sent_len/2);
    score = [];
    for k = 1:length(sp_len)
        if middle_frame(k) < max_seg_len/5 || sent_len-middle_frame(k)<max_seg_len/5    % we don't cut by short pause that is too close to the both ends
            score(k) = 0;
        else
            dist_to_center = abs(middle_frame(k)-sent_center)/sent_center;      % how far is current short pause to the sentence center relatively, which is in [0 0.5] interval.
            % we change the distance to log scale, multiple it by -3 to convert it to positive number. As log(0) is infinity, we need to set a upper limit to the log distance.
            location_score = min(20, -3*log(dist_to_center) );
            score(k) = sp_len(k) + loc_weight * location_score;     % we still have loc_weight to tune the relative weight of the two factors
        end
    end
    [best_sp, idx] = max(score);    % everytime we only cut the current sentence once.
    cutting_point = middle_frame(idx);

    sub_points = choose_cutting_points_by_sp(cutting_point, sp_len(1:idx-1), middle_frame(1:idx-1), max_seg_len, loc_weight);     % See whether left segment needs to be cut further
    cutting_point = ([sub_points cutting_point]);
    sub_points = choose_cutting_points_by_sp(sent_len - cutting_point(end), sp_len(idx+1:end), middle_frame(idx+1:end)-cutting_point(end), max_seg_len, loc_weight);     % See whether left segment needs to be cut further
    cutting_point = ([cutting_point sub_points+cutting_point(end)]);
end

author: singaxiong

repo: SignalGraph


text:
function M = oned_bilinear ( kernel, phi, test, w_g )

%*****************************************************************************80
%
%% ONE_BILINEAR evaluates a bilinear form.
%
%  Discussion:
%
%    This function computes integral ( kernel * phi * test }
%
%    The calculation that is carried out is equivalent to:
%
%      [ n_gauss, n_test ] = size ( test );
%      [ n_gauss, n_dof ] = size ( phi );
%
%      M = zeros ( n_test, n_dof );
%      for i = 1 : n_test
%        for j = 1 : n_dof
%          M(i,j) = ( kernel' .* test(:,i)' ) * ( phi(:,j) .* w_g );
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Parameters:
%
%    Input, real KERNEL(N_GAUSS), the kernel function in the integral,
%    evaluated at the Gauss points.
%
%    Input, real PHI(N_GAUSS,N_DOF), the element test functions evaluated
%    at the Gauss points.
%
%    Input, real TEST(N_GAUSS,N_TEST), the test functions evaluated at the
%    Gauss points.
%
%    Input, real W_G(N_GAUSS,1), a column vector of Gauss weights.
%
%    Output, real M(N_TEST,N_DOF), the bilinear form.
%
  M = test' * diag ( kernel .* w_g ) * phi;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function models = sample_models(n_models, n_vars, values)
% SAMPLE_MODELS: Sample candidate model vectors with categorical values
%
% Inputs: n_models - number of vectors to return
%         n_vars - number of variables for each vector
%         values - {n_vars x 1} cell - numeric variables for each dimension

% check values and default to binary if not specified
if nargin < 3
	values = repmat({[0,1]},n_vars,1);
end

% check if n_vars and values have the same dimension
if length(values) ~= n_vars
    error('Number of variables does not match size of values cell')
end

% Generate matrix of zeros to store models
models = zeros(n_models, n_vars);

% Construct model vectors by sampling each dimension independently
% and with replacement
for j=1:n_vars
	models(:,j) = randsample(values{j}, n_models, true);
end

% -- END OF FILE --

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function varargout = plot_cubic(C,pe,p)
  % PLOT_CUBIC  Plot a cubic Bezier curve
  %
  % [pe,p] = plot_cubic(C)
  %
  % Inputs:
  %   C  4 by dim list of control points
  % Outputs:
  %   pe  plot handles for UI
  %   p  plot handle for curve
  %
  % Example:
  %   clf;hold on;arrayfun(@(c) set(plot_cubic(P(C(c,:),:)),'Color','b'),1:size(C,1));hold off;
  t = linspace(0,1)';
  P = cubic_eval(C,t);
  if ~exist('p','var') || isempty(p)
    p = plt(P,'-k','LineWidth',2);
  else
    set(p,'XData',P(:,1),'YData',P(:,2));
  end
  ish = ishold;
  hold on;
  aiblue = hex2dec(['4D';'80';'FF'])'/255;
  if ~exist('pe','var') || isempty(pe)
    pe = plot_edges(C,[1 2;3 4],'-o','Color',aiblue,'LineWidth',1);
  else
    set(pe(1),'XData',C(1:2,1),'YData',C(1:2,2));
    set(pe(2),'XData',C(3:4,1),'YData',C(3:4,2));
  end
  hold off;
  if ish
    hold on
  end
  if nargout>=1
    varargout{1} = pe;
    if nargout >= 2
      varargout{2} = p;
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
function newval = i2e_data(mpc, val, oldval, ordering, dim)
%I2E_DATA   Converts data from internal to external indexing.
%
%   VAL = I2E_DATA(MPC, VAL, OLDVAL, ORDERING)
%   VAL = I2E_DATA(MPC, VAL, OLDVAL, ORDERING, DIM)
%
%   For a case struct using internal indexing, this function can be
%   used to convert other data structures as well by passing in 3 or 4
%   extra parameters in addition to the case struct. If the value passed
%   in the 2nd argument (VAL) is a column vector or cell array, it will
%   be converted according to the ordering specified by the 4th argument
%   (ORDERING, described below). If VAL is an n-dimensional matrix or
%   cell array, then the optional 5th argument (DIM, default = 1) can be
%   used to specify which dimension to reorder. The 3rd argument (OLDVAL)
%   is used to initialize the return value before converting VAL to
%   external indexing. In particular, any data corresponding to off-line
%   gens or branches or isolated buses or any connected gens or branches
%   will be taken from OLDVAL, with VAL supplying the rest of the
%   returned data.
%
%   The ORDERING argument is used to indicate whether the data
%   corresponds to bus-, gen- or branch-ordered data. It can be one
%   of the following three strings: 'bus', 'gen' or 'branch'. For
%   data structures with multiple blocks of data, ordered by bus,
%   gen or branch, they can be converted with a single call by
%   specifying ORDERING as a cell array of strings.
%
%   Any extra elements, rows, columns, etc. beyond those indicated
%   in ORDERING, are not disturbed.
%
%   Examples:
%       A_ext = i2e_data(mpc, A_int, A_orig, {'bus','bus','gen','gen'}, 2);
%
%       Converts an A matrix for user-supplied OPF constraints from
%       internal to external ordering, where the columns of the A
%       matrix correspond to bus voltage angles, then voltage
%       magnitudes, then generator real power injections and finally
%       generator reactive power injections.
%
%       gencost_ext = i2e_data(mpc, gencost_int, gencost_orig, {'gen','gen'}, 1);
%
%       Converts a GENCOST matrix that has both real and reactive power
%       costs (in rows 1--ng and ng+1--2*ng, respectively).
%
%   See also E2I_DATA, I2E_FIELD, INT2EXT.

%   MATPOWER
%   Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if ~isfield(mpc, 'order')
    error('i2e_data: mpc does not have the ''order'' field required for conversion back to external numbering.');
end
o = mpc.order;
if o.state ~= 'i'
    error('i2e_data: mpc does not appear to be in internal order');
end
if nargin < 5
    dim = 1;
end
if ischar(ordering)         %% single set
    if strcmp(ordering, 'gen')
        v = get_reorder(val, o.(ordering).e2i, dim);
    else
        v = val;
    end
    newval = set_reorder(oldval, v, o.(ordering).status.on, dim);
else                        %% multiple sets
    be = 0;  %% base, external indexing
    bi = 0;  %% base, internal indexing
    for k = 1:length(ordering)
        ne = size(o.ext.(ordering{k}), 1);
        ni = size(mpc.(ordering{k}), 1);
        v = get_reorder(val, bi+(1:ni), dim);
        oldv = get_reorder(oldval, be+(1:ne), dim);
        new_v{k} = i2e_data(mpc, v, oldv, ordering{k}, dim);
        be = be + ne;
        bi = bi + ni;
    end
    ni = size(val, dim);
    if ni > bi              %% the rest
        v = get_reorder(val, bi+1:ni, dim);
        new_v{length(new_v)+1} = v;
    end
    newval = cat(dim, new_v{:});
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function scnlab_pca_check1(imgs, realign_files, X, spersess)
% :Usage:
% ::
%
%     function scnlab_pca_check1(imgs, realign_files or params (t x 6) across all runs, X, spersess)
%
%
% :Inputs:
%
%   **imgs:**
%        list of all image names in order
%
%   **realign_files:**
%        movement param file for each session, names in a cell array, OR
%
%        a t x 6 matrix of realignment parameters across all sessions
%
%   **X:**
%        design matrix; no intercept is needed

% :Examples:
% ::
%
%    % setup code for auditory oddball data
%    cd('/Users/tor/Documents/Tor_Documents/Coursework_and_Teaching/Mind_Res_Net_fMRI_Course_2008/data/auditory_oddball/2subjects-processed/s01/')
%
%    imgs = filenames('*/sw*img','absolute','char')
%    realign_files = filenames('*/rp*txt')
%
%    % LOAD TASK ONSETS and CREATE DESIGN MATRIX
%    onsets{1} = load('novel_stimuli_run1.asc');
%    onsets{2} = load('target_stimuli_run1.asc');
%    onsets{3} = load('standard_stimuli_run1.asc');
%    onsets{4} = load('novel_stimuli_run2.asc');
%    onsets{5} = load('target_stimuli_run2.asc');
%    onsets{6} = load('standard_stimuli_run2.asc');
%
%    regs_per_sess = 3;
%    nsess = 2;
%    for i = 1:length(onsets), onsets{i} = onsets{i}'; end
%    X = cell(1, nsess);
%    X{1} = onsets2delta(onsets(1:3), 1, 249);
%    X{1} = X{1}(:, 1:end-1);
%    X{2} = onsets2delta(onsets(4:6), 1, 249);
%    X{2} = X{2}(:, 1:end-1);
%    X = blkdiag(X{:});

% ..
%    PCA ANALYSIS USING FMRISTAT
%
%    intercept, for later removal
% ..
Xint = intercept_model(spersess, 1:10);

scum = cumsum(spersess);

if ~exist('qc_images', 'dir'), mkdir('qc_images'); end

figure; set(gcf,'Position',[307         310        1109         470]);
mask_thresh = fmri_mask_thresh(imgs(1,:));
% [V, D] = pca_image(imgs, [], 5, imgs(1,:), mask_thresh, 'pca_raw', 0);

%print('-dpsc2', '-append', 'qc_report');
%scn_export_papersetup(800); saveas(gcf,['qc_images' filesep 'pca_image_raw'],'png');

[V, D] = pca_image(imgs, [], 5, imgs(1,:), mask_thresh, 'pca_noint', 0, Xint);
subplot(1, 2, 1);
for i = 1:length(spersess), plot_vertical_line(scum(i)); end
scn_export_papersetup(800); saveas(gcf,['qc_images' filesep 'pca_image_no_intcpt'],'png');

% -------------------------------------------------------------------------
% LOAD REALIGNMENT PARAMS
% -------------------------------------------------------------------------
if iscell(realign_files)
    disp('LOADING REALIGNMENT PARAMS')

    for i = 1:length(realign_files), rparams{i} = load(realign_files{i}); end;
    rp = cat(1, rparams{:});

    % get RP by session -- for relation to PCs
    for i = 1:length(rparams), rparams{i} = zscore(rparams{i}); end
    rpz = blkdiag(rparams{:});

else
    disp('I THINK MOVEMENT PARAMS ARE ENTERED IN A SINGLE MATRIX (T X 6)')
    rp = realign_files;
    rpz = rp;

end

create_figure('Realignment params');
plot(rp);
for i = 1:length(spersess), plot_vertical_line(scum(i)); end

%OPT = scnlab_outlier_id('setup', 'tr', 2, 'spersess', spersess, 'dummy', 1:3, 'hp', 100, 'mad', 5, 'niter', 3, 'mvmt', rpzcat);
OPT = scnlab_outlier_id('setup', 'tr', 2, 'spersess', spersess, 'hp', 100, 'mad', 5, 'niter', 3, 'mvmt', rp);

m = size(V, 2); n = size(rpz, 2); p = size(X, 2);
t = size(V, 1);

% ------------------------------------------------------------------------
% Movement plot
% -------------------------------------------------------------------------
fprintf('%% Var Explained by mvmt:\n');
% for j = 1:4
%
%     [b, bint, r, rint, stats] = regress(V(:, j), rpz);
%     fprintf('\tComp %3.0f : %3.0f%%\n', j, stats(1)*100);
% end

create_figure('Scatterplots');

[H, Ax,BigAx] = plotmatrix(rp, V);
xlabel('Movement param') ; ylabel('Component')

print('-dpsc2', '-append', 'qc_report');
scn_export_papersetup(800); saveas(gcf,['qc_images' filesep 'motion_component_scatter'],'png');

Xm = intercept(rpz, 'end');
rsquare = zeros(1, 4);

disp('Variance in each component explained by realignment parameters:')

for j = 1:m
    % Center component to avoid counting intercept in r-square?
    % Don't have to, doesn't matter because var operator is 2nd moment
    y = V(:, j);
    b =  Xm \ y;

    fits = Xm * b;

    rsquare(j) = var(fits) / var(y);

    fprintf('\tComp %3.0f : %3.0f%%\n', j, rsquare(j)*100);
end

%%
% -------------------------------------------------------------------------
% Get VIFs
% -------------------------------------------------------------------------

% Task: Variance inflation factors
vifs = getvif(X);

% Variance inflation factors, considering movement params too
vifsm = getvif([X rpz]);
vifsm = vifsm(1:size(X,2));

% -------------------------------------------------------------------------
% PLOT OF RELATIONSHIPS
% -------------------------------------------------------------------------

cc = corrcoef([rpz V]);
mbyc = cc(1:n, n+1:end);

cc = corrcoef([rpz X]);
mbyt = cc(1:n, n+1:end);

cc = corrcoef([X V]);
tbyc = cc(1:p, p+1:end);

create_figure('Relationships', 3, 3);
for i = 1:9, axh(i) = subplot(3, 3, i); end

% % for i = 1:9
% %     pp = get(axh(i), 'Position'); pp(2) = pp(2) - .05; pp(4) = pp(4) + .05; pp(3) = pp(3) + .05; set(axh(i), 'Position', pp);
% % end
axes(axh(1))
axis off

% marginals
axes(axh(2)); cla
plot_matrix_cols(V, 'vertical')
set(gca,'YLim', [0 t+1], 'XLim', [0 m+1], 'XAxisLocation', 'top');
title('Component'), drawnow

axes(axh(3)); cla
plot_matrix_cols(X, 'vertical')
set(gca,'YLim', [0 t+1], 'XLim', [0 p+1], 'XAxisLocation', 'top');
title('Task Design'), drawnow

axes(axh(4)); cla
plot_matrix_cols(rpz)
set(gca,'YLim', [0 n+1], 'XLim', [0 t+1]);
ylabel('Movement'), drawnow

axes(axh(7)); cla
plot_matrix_cols(X)
set(gca,'YLim', [0 p+1], 'XLim', [0 t+1]);
ylabel('Task Design'), drawnow

% bivariate relationships
axes(axh(5));
imagesc(mbyc, [-1 1]); %title('Component'); ylabel('Movement');
set(gca,'YLim', [1 n], 'XLim', [.5 m+.5], 'YDir', 'Reverse');
axis off
axes(axh(6));
imagesc(mbyt, [-1 1]); %title('Task');
set(gca,'YLim', [1 n], 'XLim', [.5 p+.5], 'YDir', 'Reverse');
axis off
axes(axh(8));
imagesc(tbyc, [-1 1]); xlabel('Component'); %ylabel('Task');
set(gca,'YLim', [1 p], 'XLim', [.5 m+.5], 'YDir', 'Reverse');
set(gca,'YColor',[1 1 1])

axes(axh(9));
axis off

% multivariate relationships
h3 = axes('position',[0.4108    0.7093-.08    0.2134    0.01]);
rsq_mvmt = rsquare_calc(rpz, V);
imagesc(rsq_mvmt, [-1 1]);
title('R^2 by Mvmt', 'FontSize', 14);
axis off

h4 = axes('position',[0.4108    0.4096-.08    0.2134    0.01]);
rsq_task = rsquare_calc(X, V);
imagesc(rsq_task, [-1 1]);
title('R^2 by Task', 'FontSize', 14);
axis off

h4 = axes('position',[0.6916    0.7093-.08    0.2134    0.01]);
rsq_mvmt = rsquare_calc(rpz, X);
imagesc(rsq_mvmt, [-1 1]);
title('R^2 by Mvmt', 'FontSize', 14);
axis off

% axis colorbar
h = axes('Position', [0.11    0.7093 .2 .01]);
imagesc([0], [-1 1]); colorbar('horiz');
title('Correlation', 'FontSize', 14);
axis off

% Var Inflation
h2 = axes('Position', [.69 .27 .2134 .04], 'FontSize', 14, 'YColor', [1 1 1]);
imagesc(vifs, [0 10]);
%colorbar('horiz')
title('Var. Inflation (VIF)');
axis off

% Var Inflation
h2 = axes('Position', [.69 .19 .2134 .04], 'FontSize', 14, 'YColor', [1 1 1]);
imagesc(vifsm, [0 10]);
colorbar('horiz')
title('VIF w/Mvmt');

%cmap = colormap_tor([0 0 1], [1 0 0], [1 1 1]);
cmap = colormap_tor([0 0 1], [1 1 0], [.5 .2 1], [1 1 1], [1 0 0]);
colormap(cmap)

print('-dpsc2', '-append', 'qc_report');
scn_export_papersetup(800); saveas(gcf,['qc_images' filesep 'task_mvmt_component_summary'],'png');

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef PathType
    % This class (enumerator) is meant to manage the  path types:
    % fillet, dubins, straight_rrt
   enumeration
      fillet, dubins, straight_rrt
   end
end

author: lis-epfl

repo: swarmlab

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


text:
function rgb = ImGray2Pseudocolor(gim, map, n)
% IMGRAY2PSEUDOCOLOR transform a gray image to pseudocolor image
%   GIM is the input gray image data
%   MAP is the colormap already defined in MATLAB, for example:
%      'Jet','HSV','Hot','Cool','Spring','Summer','Autumn','Winter','Gray',
%      'Bone','Copper','Pink','Lines'
%   N specifies the size of the colormap
%   rgb is the output COLOR image data
%
% Main codes stolen from:
%       http://www.alecjacobson.com/weblog/?p=1655
%       %% rgb = ind2rgb(gray2ind(im,255),jet(255));                      %
%

[nr,nc,nz] = size(gim);
rgb = zeros(nr,nc,3);

if ( ~IsValidColormap(map) )
    disp('Error in ImGray2Pseudocolor: unknown colormap!');
elseif (~(round(n) == n) || (n < 0))
    disp('Error in ImGray2Pseudocolor: non-integer or non-positive colormap size');
else
    fh = str2func(ExactMapName(map));
    rgb = ind2rgb(gray2ind(gim,n),fh(n));
    rgb = uint8(rgb*255);
end

if (nz == 3)
    rgb = gim;
    disp('Input image has 3 color channel, the original data returns');
end

function y = IsValidColormap(map)

y = strncmpi(map,'Jet',length(map)) | strncmpi(map,'HSV',length(map)) |...
    strncmpi(map,'Hot',length(map)) | strncmpi(map,'Cool',length(map)) |...
    strncmpi(map,'Spring',length(map)) | strncmpi(map,'Summer',length(map)) |...
    strncmpi(map,'Autumn',length(map)) | strncmpi(map,'Winter',length(map)) |...
    strncmpi(map,'Gray',length(map)) | strncmpi(map,'Bone',length(map)) |...
    strncmpi(map,'Copper',length(map)) | strncmpi(map,'Pink',length(map)) |...
    strncmpi(map,'Lines',length(map));

function emapname = ExactMapName(map)

if strncmpi(map,'Jet',length(map))
    emapname = 'Jet';
elseif strncmpi(map,'HSV',length(map))
    emapname = 'HSV';
elseif strncmpi(map,'Hot',length(map))
    emapname = 'Hot';
elseif strncmpi(map,'Cool',length(map))
    emapname = 'Cool';
elseif strncmpi(map,'Spring',length(map))
    emapname = 'Spring';
elseif strncmpi(map,'Summer',length(map))
    emapname = 'Summer';
elseif strncmpi(map,'Autumn',length(map))
    emapname = 'Autumn';
elseif strncmpi(map,'Winter',length(map))
    emapname = 'Winter';
elseif strncmpi(map,'Gray',length(map))
    emapname = 'Gray';
elseif strncmpi(map,'Bone',length(map))
    emapname = 'Bone';
elseif strncmpi(map,'Copper',length(map))
    emapname = 'Copper';
elseif strncmpi(map,'Pink',length(map))
    emapname = 'Pink';
elseif strncmpi(map,'Lines',length(map))
    emapname = 'Lines';
end

author: Sable

repo: mcbench-benchmarks


text:
function f = addBreaks(f, breaks, tol)
%ADDBREAKS   Add breakpoints to the domain of a CHEBFUN.
%   ADDBREAKS(F, BREAKS) attempts to insert breakpoints in F at the points in
%   the vector BREAKS.  BREAKS need not be sorted or have only unique values,
%   but it should consist only of points in the domain of F.  Breakpoints will
%   not be inserted if they are too close together or too close to existing
%   breakpoints.
%
%   ADDBREAKS(F, BREAKS, TOL) does the same but uses the tolerance TOL as a
%   lower bound for the tolerance used in deciding if breakpoints are too close
%   to each other or to existing ones to qualify for insertion.
%
% See also ADDBREAKSATROOTS.

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

% Lower bound for tolerance:
if ( nargin == 2 )
    tol = 0;
end

% Since each column of an array-valued CHEBFUN must have the same breakpoints,
% we simply take unique(breaks(:)).
breaks = unique(breaks(:));

% Get rid of new breakpoints at infinity:
breaks( isinf(breaks) ) = [];

% Get rid of new breakpoints at infinity:
domain = f.domain( isfinite(f.domain) );

% Set a tolerance:
breakTol = max(100*eps*max(min(diff(domain)), 1), tol);

% Distance between new breakpoints and existing ones:
dist = abs(bsxfun(@minus, breaks, domain));

% Avoid introducing new breakpoints close to existing ones:
% Note that: 1. In case of [-Inf Inf], DIST is empty.
%            2. In case of [a Inf] or [-Inf b], BREAKTOL is empty.
if ( ~isempty(dist) && ~isempty(breakTol) )
    breaks( any(dist < breakTol, 2) ) = [];
end

% Add new breaks if required:
if ( ~isempty(breaks) )
    % Get the domain with the new breakpoints: (union is not required, by above)
    dom = unique([f.domain, breaks.']);

    % Handle periodic functions:
    if ( isPeriodicTech(f) )
        f = chebfun(f);
    end

    % Introduce these breakpoints into f:
    f = restrict(f, dom);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ccf,pst] = spm_csd2ccf(csd,Hz,dt)
% Converts cross spectral density to cross covariance function
% FORMAT [ccf,pst] = spm_csd2ccf(csd,Hz,dt)
%
% csd  (n,:,:)          - cross spectral density (cf, mar.P)
% Hz   (n x 1)          - vector of frequencies (Hz)
% dt                    - samping interval (default = 1/(2*Hz(end)))
%
% ccf                   - cross covariance functions
% pst  (N,1)            - vector of lags for evaluation (seconds)
%
% Note that because this scheme uses FFT one can only change dt.
%
% See also:
%  spm_ccf2csd.m, spm_ccf2mar, spm_csd2ccf.m, spm_csd2mar.m, spm_mar2csd.m,
%  spm_csd2coh.m, spm_Q.m, spm_mar.m and spm_mar_spectral.m
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_csd2ccf.m 6395 2015-03-26 15:05:04Z adeel $

% Nyquist
%--------------------------------------------------------------------------
if nargin < 3, dt  = 1/(2*Hz(end)); end

% unpack cells
%--------------------------------------------------------------------------
if iscell(csd)
    for i = 1:length(csd)
       [ccfi,pst] = spm_csd2ccf(csd{i},Hz,dt);
       ccf{i}     = ccfi;
    end
    return
end

% unpack time bins (for time-frequency responses)
%--------------------------------------------------------------------------
if ndims(csd) == 4
    for i = 1:size(csd,1)
       [ccfi,pst]   = spm_csd2ccf(squeeze(csd(i,:,:,:)),Hz,dt);
       ccf(i,:,:,:) = ccfi;
    end
    return
end

% indices for FFT
%--------------------------------------------------------------------------
dw    = Hz(2) - Hz(1);
Hz    = Hz/dw;
ns    = 1/dt;
N     = ceil(ns/2/dw);
gj    = find(Hz > 0 & Hz < (N + 1));
gi    = gj + ceil(Hz(1)) - 1;
g     = zeros(N,1);

% Fourier transform cross-spectral density
%==========================================================================
for i = 1:size(csd,2)
    if ismatrix(csd)
        g(gi)      = csd(gj,i);
        f          = ifft([0; g; flipud(conj(g))]);
        ccf(:,i)   = real(fftshift(f))*N*dw;
    else
        for j = 1:size(csd,3)
            g(gi)      = csd(gj,i,j);
            f          = ifft([0; g; flipud(conj(g))]);
            ccf(:,i,j) = real(fftshift(f))*N*dw;
        end
    end
end

% Compute time bins
%--------------------------------------------------------------------------
pst = dt*(-N:N);

author: spm


text:
function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly
g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the sigmoid of each value of z (z can be a matrix,
%               vector or scalar)

g = 1./(1 + exp(-z));

% =============================================================

end

author: UtkarshPathrabe


text:
% gabor2d() - generate a two-dimensional gabor matrice.
%
% Usage:
%   >> [ matrix ] = gabor2d(rows, columns);
%   >> [ matrix ] = gabor2d( rows, columns, freq, ...
%                             angle, sigmaR, sigmaC, meanR, meanC, dephase, cut)
% Example :
%	>> imagesc(gabor2d( 50, 50))
%
% Inputs:
%   rows        - number of rows
%   columns     - number of columns
%   freq        - frequency of the sinusoidal function in degrees (default: 360/rows)
%   angle       - angle of rotation of the resulting 2-D array in
%                 degrees of angle {default: 0}.
%   sigmaR      - standard deviation for rows {default: rows/5}
%   sigmaC      - standard deviation for columns {default: columns/5}
%   meanR       - mean for rows {default: center of the row}
%   meanC       - mean for columns {default: center of the column}
%   dephase     - phase offset in  degrees {default: 0}. A complex Gabor wavelet
%                 can be build using gabor2dd(...., 0) + i*gabor2d(...., 90),
%                 0 and 90 being the phase offset of the real and imaginary parts
%   cut	        - percentage (0->1) of maximum value below which to remove values
%                 from the matrix {default: 0}
% Ouput:
%   matrix - output gabor matrix
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function mat = gabor2d( sizeX, sizeY, freq, angle, sigmaX, sigmaY, meanX, ...
meanY, dephase, cut);

if nargin < 2
	help gabor2d
	return;
end;
if nargin < 3
	freq = 360/sizeX;
end;
if nargin < 4
	angle = 0;
end;
if nargin < 5
	sigmaX = sizeX/5;
end;
if nargin < 6
	sigmaY = sizeY/5;
end;
if nargin < 7
	meanX = (sizeX+1)/2;
end;
if nargin < 8
	meanY = (sizeY+1)/2;
end;
if nargin < 9
	dephase = 0;
end;
if nargin < 10
	cut = 0;
end;
freq = freq/180*pi;

X = linspace(1, sizeX, sizeX)'* ones(1,sizeY);
Y = ones(1,sizeX)'   		  * linspace(1, sizeY, sizeY);
%[-sizeX/2:sizeX/2]'*ones(1,sizeX+1);
%Y = ones(1,sizeY+1)'   *[-sizeY/2:sizeY/2];

rotatedmat = ((X-meanX)+i*(Y-meanY)) * exp(i*angle/180*pi);
mat = sin(real(rotatedmat)*freq + dephase/180*pi).*exp(-0.5*(  ((X-meanX)/sigmaX).*((X-meanX)/sigmaX)...
				+((Y-meanY)/sigmaY).*((Y-meanY)/sigmaY)))...
            			/((sigmaX*sigmaY)^(0.5)*pi);

if cut > 0
	maximun = max(max(mat))*cut;
	I = find(mat < maximun);
	mat(I) = 0;
end;

return;

% other solution
% --------------

for X = 1:sizeX
    for Y = 1:sizeY
        mat(X,Y) = sin(real((X-meanX+j*(Y-meanY))*exp(i*angle/180*pi))*freq + dephase/180*pi) ...
            .*exp(-0.5*(  ((X-meanX)/sigmaX).*((X-meanX)/sigmaX)...
                          +((Y-meanY)/sigmaY).*((Y-meanY)/sigmaY)))...
            			/((sigmaX*sigmaY)^(0.5)*pi);
    end;
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
function A = myACF(V)
%%
%% Provides all N-1 autocovariance estimates from a vector of observations
%%
N = length(V);
A = xcov(V(~isnan(V)), 'biased')';
NV = sum(~isnan(V));
A = [repmat(NaN, 1, N-NV+1) A(NV:(2*(NV-1)))];


text:
function computePredictionCT_batchHN(pathWORK,fSetName,outcomes,freedomMat,maxOrder,nBoot,imbalance,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% function computePredictionCT_batchHN(pathWORK,fSetName,outcomes,freedomMat,maxOrder,nBoot,imbalance,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for the 'CT'
% feature set only and for all outcomes. See ref. [1,2] for more details.
% Goal: comparison with the "Radiomics signature" of (Aerts et al., Nat Commun, 2014)
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the HN WORKSPACE directory.
% - fSetName: Cell of strings specifying the name of the type of feature
%             set analyzed (e.g., {'PET', 'SEPARATE', 'FUSED'})
% - outcomes: Structure specifying the status (1 or 0) for different
%             outcomes in HN cancer. Contains: outcomes.Failure,
%             outcomes.Locoregional, outcomes.Distant, outcomes.Death. See
%             ref.[1] for more details.
% - freedomMat:  Cell vector of matrices of row vectors of 1's and 0's to
%                specify the degree of freedom on texture extraction
%                parameters for all experiments. For example, for an nth
%                experiment of an ith scan where extraction parameters 1, 2
%                and 4 in paramCell are allowed to vary, use
%                freedomMat{i}(n,:) = [1,1,0,1].
% - maxOrder: Integer specifying the maximal multivariable model order
%             to construct.
% - nBoot: Number of bootstrap samples to use.
% - imbalance: String specifying the type of imbalance-adjustement strategy
%              employed. Either 'IABR' for imbalance-adjusted bootstrap
%              resampling (see ref.[2]), or 'IALR' for imbalance-adjusted
%              logistic regression (see ref.[1]).
% - nBatch: Number of parallel batch.
% - matlabPATH: Full path to the MATLAB excutable on the system.
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named
%          'RESULTS' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathWORK), cd('RESULTS')
nameOutcomes = fieldnames(outcomes);
nOutcomes = length(nameOutcomes);
outcomeSep = cell(1,nOutcomes);
for o = 1:nOutcomes
    temp.(nameOutcomes{o}) = outcomes.(nameOutcomes{o});
    outcomeSep{o} = temp;
    temp = [];
end
nType = length(fSetName);
[param] = batchExperiments(nType,outcomes,nBatch); nBatch = length(param);
cd('SCRIPTS'), pathScripts = pwd;
time = 60; % Number of seconds to wait before checking if parallel computations are done

% PRODUCE BATCH COMPUTATIONS
save('workspace'), pause(5);
for i = 1:nBatch
    nameScript = ['Script_PredPerfEst_CT',num2str(i),'.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:size(param{i},1)
        fprintf(fid,['computeAllPrediction_HN(pathWORK,fSetName{param{',num2str(i),'}(',num2str(j),',1)},outcomeSep{param{',num2str(i),'}(',num2str(j),',2)},freedomMat{param{',num2str(i),'}(',num2str(j),',1)},maxOrder,nBoot,imbalance,',num2str(i),')\n']);
    end
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
nameVerif = {};
for i = 1:nType
    nExp = size(freedomMat{i},1);
    nParam = size(freedomMat{i},2);
    for j = 1:nExp
        temp = ['RESULTS_',fSetName{i},'_'];
        for k = 1:nParam
            temp = [temp,num2str(freedomMat{i}(j,k))];
        end
        nameVerif = [nameVerif;temp];
    end
end
nVerif = numel(nameVerif);
while 1
    pause(time);
    check = zeros(nVerif,1);
    for o = 1:nOutcomes
        cd([pathWORK,'/RESULTS/',nameOutcomes{o}])
        for i = 1:nVerif
            check(i) = check(i) + exist([nameVerif{i},'.mat']);
        end
    end
    if sum(check) == nVerif*nOutcomes*2
        break
    end
end

cd(pathScripts)
delete('workspace.mat')
cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function pass = test_plus( )
% Test diskfun plus() command

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

f1 = @(x,y) sin(pi*x.*y);  % Strictly even/pi-periodic
f2 = @(x,y) sin(pi*x);  % Strictly odd/anti-periodic
g1 = diskfun(f1);
g2 = diskfun(f2);
gplus = g1 + g2;
fplus = @(t, r) diskfun.pol2cartf(@(x,y) f1(x, y) + ...
    f2(x, y), t, r);
t = pi*(2*rand-1);
r = rand;
pass(1) = abs(feval(gplus, t, r, 'polar') - fplus(t, r)) < tol;

f1 = @(x,y) sin(pi*x.*y) + sin(pi*x);   % Mixed symmetric terms
f2 = @(x,y) sin(pi*(x-.35).*(y-2.3)) + sin(pi*(x-.35)) ;  % Mixed symmetric terms
g1 = diskfun(f1);
g2 = diskfun(f2);
gplus = g1 + g2;
fplus = @(t, r) diskfun.pol2cartf( @(x,y) f1(x, y) + f2(x, y), t, r);
t = pi*(2*rand-1);
r = rand;
pass(2) = abs(feval(gplus, t, r, 'polar') - fplus(t, r)) < tol;

% Check that compression is working:
f = diskfun(@(x,y) x.^2 + y.^2);
r = rank(f);
g = f;
for k = 1:10
    g = g + f;
end
pass(3) = norm(g - 11*f, inf) < vscale(g)*tol;
pass(4) = (rank(g) - r ) == 0;

% Check what happens with cancellation errors:
f = diskfun(@(x,y) sin(x.*y));
g = 2*f;
pass(5) = ( norm(g - f - f, inf) < tol );

% Check for the case where f and g are non-zero at the poles bug
% f+g is zero at the poles.
f = diskfun(@(x,y) x+1);
g = diskfun(@(x,y) cos(y));
h = f - g;
pass(6) = h.nonZeroPoles == 0;

% Check for the case where f non-zero at the poles but g is zero at the
% poles.
f = diskfun(@(x,y) cos(x));
g = diskfun(@(x,y) sin(y));
h = f + g;
pass(7) = h.nonZeroPoles == 1;
h = g + f;
pass(8) = h.nonZeroPoles == 1;

end


text:
function [oceanloadcorr] = ocean_loading_correction(time, XR, XS)

% SYNTAX:
%   [oceanloadcorr] = ocean_loading_correction(time, XR, XS);
%
% INPUT:
%   time = GPS time
%   XR   = receiver position  (X,Y,Z)
%   XS   = satellite position (X,Y,Z)
%
% OUTPUT:
%   oceanloadcorr = ocean loading correction terms (along the satellite-receiver line-of-sight)
%
% DESCRIPTION:
%   Computation of the ocean loading displacement terms.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%ocean loading displacements matrix, station-dependent (see http://holt.oso.chalmers.se/loading/)
global ol_disp zero_time

oceanloadcorr = zeros(size(XS,1),1);
if (isempty(ol_disp))
    return
end

%terms depending on the longitude of the lunar node (see Kouba and Heroux, 2001)
fj = 1; %(at 1-3 mm precision)
uj = 0; %(at 1-3 mm precision)

%ref: http://202.127.29.4/cddisa/data_base/IERS/Convensions/Convension_2003/SUBROUTINES/ARG.f
tidal_waves = [1.40519E-4, 2.0,-2.0, 0.0, 0.00; ... % M2  - semidiurnal
               1.45444E-4, 0.0, 0.0, 0.0, 0.00; ... % S2  - semidiurnal
               1.37880E-4, 2.0,-3.0, 1.0, 0.00; ... % N2  - semidiurnal
               1.45842E-4, 2.0, 0.0, 0.0, 0.00; ... % K2  - semidiurnal
               0.72921E-4, 1.0, 0.0, 0.0, 0.25; ... % K1  - diurnal
               0.67598E-4, 1.0,-2.0, 0.0,-0.25; ... % O1  - diurnal
               0.72523E-4,-1.0, 0.0, 0.0,-0.25; ... % P1  - diurnal
               0.64959E-4, 1.0,-3.0, 1.0,-0.25; ... % Q1  - diurnal
               0.53234E-5, 0.0, 2.0, 0.0, 0.00; ... % Mf  - long-period
               0.26392E-5, 0.0, 1.0,-1.0, 0.00; ... % Mm  - long-period
               0.03982E-5, 2.0, 0.0, 0.0, 0.00];    % Ssa - long-period

refdate = datenum([1975 1 1 0 0 0]);

[week, sow] = time2weektow(zero_time + time);
dateUTC = datevec(gps2utc(datenum(gps2date(week, sow))));

%separate the fractional part of day in seconds
fday = dateUTC(4)*3600 + dateUTC(5)*60 + dateUTC(6);
dateUTC(4:end) = 0;

%number of days since reference date (1 Jan 1975)
days = (datenum(dateUTC) - refdate);

capt = (27392.500528 + 1.000000035*days)/36525;

%mean longitude of the Sun at the beginning of day
H0 = (279.69668 + (36000.768930485 + 3.03e-4*capt)*capt)*pi/180;

%mean longitude of the Moon at the beginning of day
S0 = (((1.9e-6*capt - 0.001133)*capt + 481267.88314137)*capt + 270.434358)*pi/180;

%mean longitude of the lunar perigee at the beginning of day
P0 = (((-1.2e-5*capt - 0.010325)*capt + 4069.0340329577)*capt + 334.329653)*pi/180;

corr = zeros(3,1);
for k = 1 : 11
    angle = tidal_waves(k,1)*fday + tidal_waves(k,2)*H0 + tidal_waves(k,3)*S0 + tidal_waves(k,4)*P0 + tidal_waves(k,5)*2*pi;
    corr  = corr + fj*ol_disp(1).matrix(1:3,k).*cos(angle + uj - ol_disp(1).matrix(4:6,k)*pi/180);
end
corrENU(1,1) = -corr(2,1); %east
corrENU(2,1) = -corr(3,1); %north
corrENU(3,1) =  corr(1,1); %up

%displacement along the receiver-satellite line-of-sight
XRcorr = local2globalPos(corrENU, XR);
corrXYZ = XRcorr - XR;
for s = 1 : size(XS,1)
    LOS  = XR - XS(s,:)';
    LOSu = LOS / norm(LOS);
    % oceanloadcorr(s,1) = dot(corrXYZ,LOSu);
    oceanloadcorr(s,1) = sum(conj(corrXYZ).*LOSu);
end

author: goGPS-Project


text:
function y = sinc(x)

y = pi * x;

sel = abs(y) > 1e-15;

y(sel) = sin(y(sel)) ./ y(sel);
y(~sel) = 1;

author: thomaskuestner


text:
function CPD = gaussian_CPD(bnet, self, varargin)
% GAUSSIAN_CPD Make a conditional linear Gaussian distrib.
%
% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,
% where node is the number of a node in this equivalence class.

% To define this CPD precisely, call the continuous (cts) parents (if any) X,
% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:
% - no parents: Y ~ N(mu, Sigma)
% - cts parents : Y|X=x ~ N(mu + W x, Sigma)
% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))
% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
%
% The list below gives optional arguments [default value in brackets].
% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)
% Parameters will be reshaped to the right size if necessary.
%
% mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]
% cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(100*eye(Y,Y), [1 1 Q]) ]
% weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]
% cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]
% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]
% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]
% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]
% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]
% cov_prior_weight - weight given to I prior for estimating Sigma [0.01]
% cov_prior_entropic - if 1, we also use an entropic prior for Sigma [0]
%
% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 1)

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = init_fields;
  clamp = 0;
  CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));
  return;
elseif isa(bnet, 'gaussian_CPD')
  % This might occur if we are copying an object.
  CPD = bnet;
  return;
end
CPD = init_fields;

CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));

args = varargin;
ns = bnet.node_sizes;
ps = parents(bnet.dag, self);
dps = myintersect(ps, bnet.dnodes);
cps = myintersect(ps, bnet.cnodes);
fam_sz = ns([ps self]);

CPD.self = self;
CPD.sizes = fam_sz;

% Figure out which (if any) of the parents are discrete, and which cts, and how big they are
% dps = discrete parents, cps = cts parents
CPD.cps = find_equiv_posns(cps, ps); % cts parent index
CPD.dps = find_equiv_posns(dps, ps);
ss = fam_sz(end);
psz = fam_sz(1:end-1);
dpsz = prod(psz(CPD.dps));
cpsz = sum(psz(CPD.cps));

% set default params
CPD.mean = randn(ss, dpsz);
CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]);
CPD.weights = randn(ss, cpsz, dpsz);
CPD.cov_type = 'full';
CPD.tied_cov = 0;
CPD.clamped_mean = 0;
CPD.clamped_cov = 0;
CPD.clamped_weights = 0;
CPD.cov_prior_weight = 0.01;
CPD.cov_prior_entropic = 0;
nargs = length(args);
if nargs > 0
  CPD = set_fields(CPD, args{:});
end

% Make sure the matrices have 1 dimension per discrete parent.
% Bug fix due to Xuejing Sun 3/6/01
CPD.mean = myreshape(CPD.mean, [ss ns(dps)]);
CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]);
CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]);

% Precompute indices into block structured  matrices
% to speed up CPD_to_lambda_msg and CPD_to_pi
cpsizes = CPD.sizes(CPD.cps);
CPD.cps_block_ndx = cell(1, length(cps));
for i=1:length(cps)
  CPD.cps_block_ndx{i} = block(i, cpsizes);
end

%%%%%%%%%%%
% Learning stuff

% expected sufficient statistics
CPD.Wsum = zeros(dpsz,1);
CPD.WYsum = zeros(ss, dpsz);
CPD.WXsum = zeros(cpsz, dpsz);
CPD.WYYsum = zeros(ss, ss, dpsz);
CPD.WXXsum = zeros(cpsz, cpsz, dpsz);
CPD.WXYsum = zeros(cpsz, ss, dpsz);

% For BIC
CPD.nsamples = 0;
switch CPD.cov_type
 case 'full',
  % since symmetric
    %ncov_params = ss*(ss-1)/2;
    ncov_params = ss*(ss+1)/2;
  case 'diag',
    ncov_params = ss;
  otherwise
    error(['unrecognized cov_type ' cov_type]);
end
% params = weights + mean + cov
if CPD.tied_cov
  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params;
else
  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params;
end

% for speeding up maximize_params
CPD.useC = exist('rep_mult');

clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights;
CPD = set_clamped(CPD, clamped);

%%%%%%%%%%%

function CPD = init_fields()
% This ensures we define the fields in the same order
% no matter whether we load an object from a file,
% or create it from scratch. (Matlab requires this.)

CPD.self = [];
CPD.sizes = [];
CPD.cps = [];
CPD.dps = [];
CPD.mean = [];
CPD.cov = [];
CPD.weights = [];
CPD.clamped_mean = [];
CPD.clamped_cov = [];
CPD.clamped_weights = [];
CPD.cov_type = [];
CPD.tied_cov = [];
CPD.Wsum = [];
CPD.WYsum = [];
CPD.WXsum = [];
CPD.WYYsum = [];
CPD.WXXsum = [];
CPD.WXYsum = [];
CPD.nsamples = [];
CPD.nparams = [];
CPD.cov_prior_weight = [];
CPD.cov_prior_entropic = [];
CPD.useC = [];
CPD.cps_block_ndx = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear all, close all, clc
A = imread('../../CH01_SVD/DATA/dog.jpg');
B = rgb2gray(A);

for j=1:size(B,1);
    Cshift(j,:) = fftshift(fft(B(j,:)));
    C(j,:) = (fft(B(j,:)));
end

for j=1:size(C,2);
    D(:,j) = fft(C(:,j));
end

subplot(1,3,1)
imagesc(B);
subplot(1,3,2)
imagesc(log(abs(Cshift)))
subplot(1,3,3)
imagesc(fftshift(log(abs(D))))
colormap gray
%%
figure
imagesc(B)
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../figures/2DFFTa');

figure
imagesc(log(abs(Cshift)))
colormap gray

axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../figures/2DFFTb');
figure
imagesc(fftshift(log(abs(D))))
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../figures/2DFFTc');

%% FFT Compression
figure
Bt=fft2(B);
F = log(abs(fftshift(Bt))+1);  % put FFT on log-scale

% Zero out all small coefficients and inverse transform
Btsort = sort(abs(Bt(:)));
keepvec = [.1 .05 .01 .005 .002 .001];
for k=4:6
    keep = keepvec(k);
    % keep = 0.05;
    thresh = Btsort(floor((1-keep)*length(Btsort)));
    ind = abs(Bt)>thresh;
    Atlow = Bt.*ind;
    Flow = log(abs(fftshift(Atlow))+1);  % put FFT on log-scale
    imshow(mat2gray(Flow),[]);

    % Plot Reconstruction
    Alow=uint8(ifft2(Atlow));
    imshow(Alow)

    axis off
    set(gcf,'PaperPositionMode','auto')
    %     print('-depsc2', '-loose', ['../figures/2DFFT_Compress',num2str(k)]);
    print('-dpng', '-loose', ['../figures/2DFFT_Compress',num2str(k)]);

end

%% Denoise
B = Bold(2:2:end,2:2:end);
ind = rand(size(B));
Bnoise = B + uint8(200*randn(size(B)));
% Bnoise = B + uint8(200*rand(size(B)));
% Bnoise = B.*uint8(ind<.2);
imagesc(Bnoise)
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
% print('-dpng', '-loose', '../figures/2DFFT_Noise1');

Bt=fft2(Bnoise);
Btshift = fftshift(Bt);
[nx,ny] = size(B);

figure
F = log(abs(Btshift)+1);  % put FFT on log-scale
imagesc(F)
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
% print('-dpng', '-loose', '../figures/2DFFT_Noise2');

[X,Y] = meshgrid(-ny/2+1:ny/2,-nx/2+1:nx/2);
R2 = X.^2+Y.^2;
ind = R2<150^2;
Btshiftfilt = Btshift.*ind;

figure
Ffilt = log(abs(Btshiftfilt)+1);  % put FFT on log-scale
imagesc(Ffilt)
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
% print('-dpng', '-loose', '../figures/2DFFT_Noise3');

figure
Btfilt = ifftshift(Btshiftfilt);
Bfilt = ifft2(Btfilt);
imagesc(uint8(Bfilt))
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
% print('-dpng', '-loose', '../figures/2DFFT_Noise4');

%% Wavelet Compression
figure
[C,S] = wavedec2(B,4,'db1');

% Zero out all small coefficients and inverse transform
Csort = sort(abs(C(:)));
keepvec = [.1 .05 .01 .005 .002 .001];
for k=1:6
    figure
    keep = keepvec(k);
    % keep = 0.05;
    thresh = Csort(floor((1-keep)*length(Csort)));
    ind = abs(C)>thresh;
    Cfilt = C.*ind;

    % Plot Reconstruction
    Alow=uint8(waverec2(Cfilt,S,'db1'));
    imagesc(uint8(Alow))
    colormap gray

    colormap gray
    axis off
    axis tight
    set(gcf,'Position',[100 100 600 800])
    set(gcf,'PaperPositionMode','auto')
    %     print('-depsc2', '-loose', ['../figures/2DFFT_Compress',num2str(k)]);
        print('-dpng', '-loose', ['../figures/2DWavelet_Compress',num2str(k)]);

end

author: dynamicslab


text:
function [ d2, p2 ] = poly_power ( d1, p1, n )

%*****************************************************************************80
%
%% POLY_POWER computes a power of a polynomial.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_integrals/poly_power.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D1, the degree of the polynomial.
%
%    Input, real P1(M1), the polynomial coefficients.
%    M1 = ((D1+1)*(D1+2))/2.
%
%    Input, integer N, the nonnegative integer power.
%
%    Output, integer D2, the degree of the power polynomial.
%    D2 = N * D1.
%
%    Output, real P2(M2), the polynomial power.
%    M2 = ((D2+1)*(D2+2))/2.
%

%
%  Create P2, a polynomial representation of 1.
%
  d2 = 0;
  m2 = ( ( d2 + 1 ) * ( d2 + 2 ) ) / 2;
  p2 = zeros ( m2, 1 );
  p2(1) = 1.0;
%
%  Iterate N times:
%    P2 <= P2 * P1
%
  for i = 1 : n
    [ d2, p2 ] = poly_product ( d2, p2, d1, p1 );
  end

  return
end


text:
%[2009]-"Text feature selection using ant colony optimization"

% (9/12/2020)

function ACO = jAntColonyOptimization(feat,label,opts)
% Parameters
tau   = 1;      % pheromone value
eta   = 1;      % heuristic desirability
alpha = 1;      % control pheromone
beta  = 0.1;    % control heuristic
rho   = 0.2;    % pheromone trail decay coefficient

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'tau'), tau = opts.tau; end
if isfield(opts,'alpha'), alpha = opts.alpha; end
if isfield(opts,'beta'), beta = opts.beta; end
if isfield(opts,'rho'), rho = opts.rho; end
if isfield(opts,'eta'), eta = opts.eta; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial Tau & Eta
tau = tau * ones(dim,dim);
eta = eta * ones(dim,dim);
% Pre
fitG = inf;
fit  = zeros(1,N);

curve = inf;
t = 1;
% Iterations
while t <= max_Iter
	% Reset ant
	X = zeros(N,dim);
	for i = 1:N
    % Random number of features
    num_feat = randi([1,dim]);
    % Ant start with random position
    X(i,1)   = randi([1,dim]);
    k        = [];
    if num_feat > 1
      for d = 2:num_feat
        % Start with previous tour
        k      = [k(1:end), X(i, d-1)];
        % Edge/Probability Selection (2)
        P      = (tau(k(end),:) .^ alpha) .* (eta(k(end),:) .^ beta);
        % Set selected position = 0 probability (2)
        P(k)   = 0;
        % Convert probability (2)
        prob   = P ./ sum(P(:));
        % Roulette Wheel selection
        route  = jRouletteWheelSelection(prob);
        % Store selected position to be next tour
        X(i,d) = route;
      end
    end
  end
  % Binary
  X_bin = zeros(N,dim);
  for i = 1:N
    % Binary form
    ind           = X(i,:);
    ind(ind == 0) = [];
    X_bin(i,ind)  = 1;
  end
  % Fitness
  for i = 1:N
    % Fitness
    fit(i) = fun(feat,label,X_bin(i,:),opts);
    % Global update
    if fit(i) < fitG
      Xgb  = X(i,:);
      fitG = fit(i);
    end
  end
%---// [Pheromone update rule on tauK] //
  tauK = zeros(dim,dim);
  for i = 1:N
    % Update Phromones
    tour = X(i,:);
    tour(tour == 0) = [];
    % Number of features
    len_x = length(tour);
    tour  = [tour(1:end), tour(1)];
    for d = 1:len_x
      % Feature selected on graph
      x = tour(d);
      y = tour(d + 1);
      % Update delta tau k on graph (3)
      tauK(x,y) = tauK(x,y) + (1 / (1 + fit(i)));
    end
  end
%---// [Pheromone update rule on tauG] //
  tauG = zeros(dim,dim);
  tour = Xgb;
  tour(tour == 0) = [];
  % Number of features
  len_g = length(tour);
  tour  = [tour(1:end), tour(1)];
  for d = 1:len_g
    % Feature selected on graph
    x = tour(d);
    y = tour(d + 1);
    % Update delta tau G on graph
    tauG(x,y) = 1 / (1 + fitG);
  end
%---// Evaporate pheromone // (4)
  tau = (1 - rho) * tau + tauK + tauG;
  % Save
  curve(t) = fitG;
  fprintf('\nIteration %d Best (ACO)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Sf = Xgb;
Sf(Sf == 0) = [];
sFeat = feat(:,Sf);
% Store results
ACO.sf = Sf;
ACO.ff = sFeat;
ACO.nf = length(Sf);
ACO.c  = curve;
ACO.f  = feat;
ACO.l  = label;
end

%// Roulette Wheel Selection //
function Index = jRouletteWheelSelection(prob)
% Cummulative summation
C = cumsum(prob);
% Random one value, most probability value [0~1]
P = rand();
% Roulette wheel
for i = 1:length(C)
	if C(i) > P
    Index = i;
    break;
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function mypatch=patchtest(some_results, col_heading)

    %{
    some_results is the results from one of the ZmapGridFunctions. It will contain:
    1.  a field called "values" which is a table.  This table will have columns for 'x', 'y' and
        col_heading (where col_heading is the name of one of the columns).
    2.  a field called Grid, which is a ZmapGrid and contains fields X, Y

    %}
    gr_s = some_results.Grid;
    results = some_results.values;

    myX=some_results.values.x;
    myX=reshape(myX,size(some_results.Grid.X));
    hold on;
    % assume: dx is constant for every latitude
    dx_at_lat = min(diff(myX,[],2),[],2);   % Nx1

    % fill all X position nans, because they will cause holes
    for i=1:size(myX,1)
        AnchorIdx = find(~ismissing(myX(i,:)),1,'first');
        AnchorVal = myX(i,AnchorIdx);
        missIdx = find(ismissing(myX(i,:)));
        myX(i,missIdx) = (missIdx - AnchorIdx) .* dx_at_lat(i) + AnchorVal;
    end

    % assume dy is constant everywhere
    dy = mean(min(diff(gr_s.Y)));
    shifted_X = myX - repmat(dx_at_lat ./2, 1, size(myX,2)) ;
    shifted_Y = gr_s.Y-dy./2;

    myresults=results.(col_heading);
    myresults=reshape(myresults,size(shifted_X));

    %% because surfaces and patches are based on the lower-left corner, add col & row.

    % add row to top with same values as existing last (top)row
    myresults(end+1,:)=myresults(end,:);
    shifted_X(end+1,:)=shifted_X(end,:);
    shifted_Y(end+1,:)=shifted_Y(end,:) + dy;

    % add column to end with same values as existing last (right) column
    myresults(:,end+1)=myresults(:,end);
    shifted_X(:,end+1)=shifted_X(:,end) + [dx_at_lat; dx_at_lat(end)];
    shifted_Y(:,end+1)=shifted_Y(:,end);

    pa=surf2patch(shifted_X,shifted_Y,zeros(size(shifted_X)),myresults);

    mypatch=patch(pa);
    mypatch.Faces(end,:)=[]; %this point was made up.
    mypatch.Tag='the_patch';
    mypatch.HitTest='off';
    shading faceted;

    set(gca,'Children',circshift(get(gca,'Children'),-1)); % put this patch at the end
end

author: CelsoReyes


text:
function linpack_z_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests ZPPCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  For a double precision complex (C)\n' );
  fprintf ( 1, '  Hermitian positive definite packed matrix (PP),\n' );
  fprintf ( 1, '  ZPPCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  a(1) = complex ( 2.5281,  0.0000 );

  a(2) = complex ( 2.1341, -0.2147 );
  a(3) = complex ( 3.0371,  0.0000 );

  a(4) = complex ( 2.4187,  0.2932 );
  a(5) = complex ( 2.0905,  1.1505 );
  a(6) = complex ( 2.7638,  0.0000 );
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition number.\n' );

  [ a, rcond, info ] = zppco ( a, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reciprocal condition number = %f\n', rcond );

  return
end

author: johannesgerer


text:
function [node,elem,face]=cgals2m(v,f,opt,maxvol,varargin)
%
% [node,elem,face]=cgals2m(v,f,opt,maxvol)
%
% wrapper for CGAL 3D mesher (CGAL 3.5 and newer)
% convert a triangular surface to tetrahedral mesh
%
% http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input:
%	 v: the node coordinate list of a surface mesh (nn x 3)
%	 f: the face element list of a surface mesh (be x 3)
%	 opt: parameters for CGAL mesher, if opt is a structure, then
%	     opt.radbound: defines the maximum surface element size
%	     opt.angbound: defines the miminum angle of a surface triangle
%	     opt.distbound: defines the maximum distance between the
%		 center of the surface bounding circle and center of the
%		 element bounding sphere
%	     opt.reratio:  maximum radius-edge ratio
%	     if opt is a scalar, it only specifies radbound.
%	 maxvol: target maximum tetrahedral elem volume
%
% output:
%	 node: output, node coordinates of the tetrahedral mesh
%	 elem: output, element list of the tetrahedral mesh, the last
%	      column is the region id
%	 face: output, mesh surface element list of the tetrahedral mesh
%	      the last column denotes the boundary ID
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1,'creating surface and tetrahedral mesh from a polyhedral surface ...\n');

exesuff=fallbackexeext(getexeext,'cgalpoly');

ang=30;
ssize=6;
approx=0.5;
reratio=3;

flags=varargin2struct(varargin{:});

if(~isstruct(opt))
	ssize=opt;
end

if(isstruct(opt) && length(opt)==1)  % does not support settings for multiple labels
	if(isfield(opt,'radbound'))   ssize=opt.radbound; end
	if(isfield(opt,'angbound'))   ang=opt.angbound; end
	if(isfield(opt,'distbound'))  approx=opt.distbound; end
	if(isfield(opt,'reratio'))    reratio=opt.reratio; end
end
if(getoptkey('DoRepair',0,flags)==1)
    [v,f]=meshcheckrepair(v,f);
end
saveoff(v,f,mwpath('pre_cgalpoly.off'));
deletemeshfile(mwpath('post_cgalpoly.mesh'));

randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"

if(~isempty(getvarfrom({'caller','base'},'ISO2MESH_RANDSEED')))
        randseed=getvarfrom({'caller','base'},'ISO2MESH_RANDSEED');
end

cmd=sprintf('"%s%s" "%s" "%s" %.16f %.16f %.16f %.16f %.16f %d',mcpath('cgalpoly'),exesuff,...
    mwpath('pre_cgalpoly.off'),mwpath('post_cgalpoly.mesh'),ang,ssize,...
    approx,reratio,maxvol,randseed);
system(cmd);
if(~exist(mwpath('post_cgalpoly.mesh'),'file'))
    error(sprintf('output file was not found, failure was encountered when running command: \n%s\n',cmd));
end
[node,elem,face]=readmedit(mwpath('post_cgalpoly.mesh'));

fprintf(1,'node number:\t%d\ntriangles:\t%d\ntetrahedra:\t%d\nregions:\t%d\n',...
    size(node,1),size(face,1),size(elem,1),length(unique(elem(:,end))));
fprintf(1,'surface and volume meshes complete\n');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dataCorners2 = atlasMergeAdjacentCorners(dataCorners,dataCorners2);
%
%   dataCorners2 = atlasMergeAdjacentCorners(dataCorners,dataCorners2);
%
% Author:  Wandell
% Purpose:
%     When adding an atlas segment, sometimes we want the corners from
%     adjacent segments to be the same.  This routine finds the two closest
%     corners and makes the corner positions in the second segment equal to
%     the closest ones in the first segment
%
%Examples
%

% First, we check the distance from each corner in corners2 to each corner
% in corners.  This is a 4x4 distance matrix.
%
for ii=1:4
    for jj=1:4
        d(ii,jj) = norm(dataCorners(ii,:) - dataCorners2(jj,:));
    end
end

% Find the two points with the closest values
[v,idx] = sort(d(:));

% Find the row and column coordinates of these two points
[r1,c1] = ind2sub(size(d),idx(1));

% Set the dataCorners2(c,:) equal to the positions of dataCorners(r,:)
dataCorners2(c1,:) = dataCorners(r1,:);

% Choose the second point that is closest with the constraint that this
% second point in dataCorners2 cannot be attached to the point we have
% already chosen.  Another way to do this is to remove the row containing
% point at r1, recalculale the distance matrix, and then find the minimum.
for ii=2:length(idx)
    [r2,c2] = ind2sub(size(d),idx(ii));
    if r1 ~= r2
        dataCorners2(c2,:) = dataCorners(r2,:);
        break;
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
clear all, close all
n1=10; n2=5;                       % number of data points from each class n1/n2
sc=1;                                  % change the scale of the initial problem

% 1) draw data from two Gaussians
S1 = eye(2); S2 = [1 0.95; 0.95 1];
m1 = [0.75;  0]; m2 = [-0.75; 0];
randn('seed',17); rand('seed',17)
% coordinates
x1 = chol(S1)'*randn(2,n1)+repmat(m1,1,n1);
x2 = chol(S2)'*randn(2,n2)+repmat(m2,1,n2);
x1(:,5) = x1(:,2); x1(:,1) = x1(:,4); x2(:,1) = x2(:,2);  % rank deficiency in K
xtr = [x1 x2]';
% class labels
ytr = [repmat(-1,[1,n1]) repmat(1,[1,n2])]';
tt = (-5:.3:5)*sc; [t1 t2] = meshgrid(tt,tt);
x1 = x1*sc; x2=x2*sc; xtr=xtr*sc; m1=m1*sc; m2=m2*sc;
t  = [t1(:) t2(:)];
z1 = exp(-sum((t-repmat(m1',length(t),1))*inv(S1).*(t-repmat(m1',length(t),1)),2)/2);
z2 = exp(-sum((t-repmat(m2',length(t),1))*inv(S2).*(t-repmat(m2',length(t),1)),2)/2);
xte = t; clear t
xtr(end,:)=[2.9,2.3]; ytr(end)=-1;                                       % modif

% 2) set GP parameters
cov = {@covSum,{@covSEiso,@covNoise}}; hyp.cov = [0; 2; -Inf];       % ell,sf,sn
lik =  {'likLogistic'};  hyp.lik  = [];                  % likLogistic or likErf
mn  = {'meanZero'};      hyp.mean = [];
inf = 'infEP';

% 3) predict using approximated inference
[nlZ,dnlZ,post] = gp(hyp, inf, mn, cov, lik, xtr, ytr);
[ymu,ys2,fmu,fs2,junk,post] = gp(hyp,inf,mn,cov,lik,xtr,post,xte);

% 4) prediction using MCMC sampling
% set MCMC parameters, see some more details in inf/infMCMC.m
% We have two samplers implemented, namely
%  hmc - Hybrid Monte Carlo, and
%  ess - Elliptical Slice Sampling.
% par.sampler = 'hmc'; par.Nsample = 20;
par.sampler = 'ess'; par.Nais = 5; par.Nsample = 20; par.Nburnin = 20; par.Nskip = 2;

hyp = gpminimize(hyp,@gp,-200,0,inf,mn,cov,lik,xtr,ytr);
tic
[posts,nlZs,dnlZs] = infMCMC(hyp,mn,cov,lik,xtr,ytr,par);
posts
toc
[ymus,ys2s,fmus,fs2s,junk,posts] = gp(hyp,@infMCMC,mn,cov,lik,xtr,posts,xte);

% 5a) echo results
fprintf('nlZ-EP=%f, nlZ-AIS=%f\n', nlZ, nlZs)
fprintf('acceptance rate (MCMC) = %1.2f%%\n',100*posts.acceptance_rate_MCMC)
for r=1:length(posts.acceptance_rate_AIS)
  fprintf('acceptance rate (AIS) = %1.2f%%\n',100*posts.acceptance_rate_AIS(r))
end
% 5b) print results
figure
subplot(221)
  plot([-12,12],[-12,12],'r'), hold on, plot(fmus,fmu,'k.'), title('\mu_f')
  xlabel('MCMC'), ylabel(inf)
subplot(222)
  plot([0,10],[0,10],'r'), hold on, plot(sqrt(fs2s),sqrt(fs2),'k.')
  title('\sigma_f')
  xlabel('MCMC'), ylabel(inf)
subplot(223)
  plot([-1,1],[-1,1],'r'), hold on, plot(ymus,ymu,'k.'), title('\mu_y')
  xlabel('MCMC'), ylabel(inf)
subplot(224)
  plot([0,1],[0,1],'r'), hold on, plot(sqrt(ys2s),sqrt(ys2),'k.')
  title('\sigma_y')
  xlabel('MCMC'), ylabel(inf)

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [data,colorRange,minData,maxData] = scaleData(data,varargin)

data = real(data);

%if check_option(varargin,{'log','logarithmic'}), data(data==0)=NaN; end

% min and max
if check_option(varargin,{'log','logarithmic'})
  minData = min(data(~isnan(data) & ~isinf(log(data))));
else
  minData = min(data(~isnan(data) & ~isinf(data)));
end
maxData = max(data(~isnan(data) & ~isinf(data)));

% get colorrange from data
colorRange = [minData,maxData];
minData = nanmin(data(:));
maxData = nanmax(data(:));
if minData == maxData
  if minData == 0
    maxData = 1;
  else
    minData = 0;
  end
end

% from options
if check_option(varargin,{'contourf'},'double')

  contours = get_option(varargin,{'contourf','contour'},[],'double');
  colorRange = [contours(1),contours(end)];

elseif check_option(varargin,'colorRange','double')

  colorRange = get_option(varargin,'colorrange',[],'double');

  if isinf(colorRange(1))
    if isfinite(minData)
      colorRange(1) = minData;
    else
      colorRange(1) = colorRange(2)-2*eps;
    end
  end

end

% correct for allmost constant data
if isempty(colorRange)
  colorRange = [0,1];
elseif colorRange(1)>0 && ...
    ((colorRange(2)-colorRange(1))/colorRange(1) < 1e-15)
  colorRange = [min(colorRange(1),0),max(colorRange(2),1)];
end

author: mtex-toolbox

repo: mtex

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


text: %MIN Datafile overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function quadrule_test25 ( )

%*****************************************************************************80
%
%% TEST25 tests LEGENDRE_SET_SQRTX2_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2009
%
%  Author:
%
%    John Burkardt
%
  order_max = 20;

  a = 0.0;
  b = 1.0;

  nsub = 1;

  xlo = 0.0;
  xhi = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST25\n' );
  fprintf ( 1, '  LEGENDRE_SET_SQRTX2_01 sets up Gauss-Legendre\n' );
  fprintf ( 1, '  quadrature over [0,1] with weight function 1/SQRT(X);\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for iexp = 0 : 3

      norder = 2^iexp;

      fprintf ( 1, '  %2d', norder );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = legendre_set_sqrtx2_01 ( norder );

        result = 0.0;
        for j = 1 : norder
          result = result + weight(j) * func ( xtab(j) );
        end

        fprintf ( 1, '  %12f', result );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end

repo: jburkardt-m


text:
function a = bab ( n, alpha, beta )

%*****************************************************************************80
%
%% BAB returns the BAB matrix.
%
%  Example:
%
%    N = 5
%    ALPHA = 5, BETA = 2
%
%    5  2  .  .  .
%    2  5  2  .  .
%    .  2  5  2  .
%    .  .  2  5  2
%    .  .  .  2  5
%
%  Properties:
%
%    A is banded, with bandwidth 3.
%
%    A is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is Toeplitz: constant along diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    CM da Fonseca, J Petronilho,
%    Explicit Inverses of Some Tridiagonal Matrices,
%    Linear Algebra and Its Applications,
%    Volume 325, 2001, pages 7-21.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA, BETA, the parameters.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    a(i,i) = alpha;
  end

  for i = 1 : n - 1
    a(i,i+1) = beta;
    a(i+1,i) = beta;
  end

  return
end


text:
%
% FILE NAME:    mergeGlutamateSubsystem.m
%
% PURPOSE: Currently, there is only one reaction in the "Glutamate
%          metabolism" subsystem. The subsystem for this reaction will be
%          changed to "Alanine, aspartate and glutamate metabolism", which
%          contains 39 reactions.
%

% load model
load('humanGEM.mat');

% convert subsystems to cell array
subSystems = cellfun(@(x) x, ihuman.subSystems);

% find reaction in Glutamate metabolism subsystem
rxnInd = ismember(subSystems,'Glutamate metabolism');

% change subsystem to Alanine, aspartate and glutamate metabolism
ihuman.subSystems{rxnInd} = {'Alanine, aspartate and glutamate metabolism'};

% export model
exportHumanGEM(ihuman,'humanGEM','../../',{'mat','yml'},false,false);

author: SysBioChalmers

repo: Human-GEM


text:
function [hmm,fehist,feterms,rho] = hmmstrain(Xin,T,hmm,info,options)
% Stochastic HMM variational inference
%
% INPUTS
% Xin: cell with strings referring to the files containing each subject's data,
%       or cell with with matrices (time points x channels) with each
%       subject's data
% T: cell of vectors, where each element has the length of each trial per
%       subject. Dimension of T{n} has to be (1 x nTrials)
% options: HMM options for both the subject and the group runs
%
% Diego Vidaurre, OHBA, University of Oxford (2016)

N = length(Xin); K = length(hmm.state);
X = loadfile(Xin{1},T{1},options); ndim = size(X,2); XW = [];
orders = formorders(hmm.train.order,hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag);
pcapred = hmm.train.pcapred>0;
p = hmm.train.lowrank; do_HMM_pca = (p > 0); rangeK = 1:K;
if pcapred
    npred = hmm.train.pcapred;
else
    if isfield(hmm.train,'B') && ~isempty(hmm.train.B), Q = size(hmm.train.B,2);
    else Q = ndim;
    end
    npred = length(orders)*Q;
end
subjfe = info.subjfe;
loglik = info.loglik;
statekl = info.statekl;
Dir2d_alpha = info.Dir2d_alpha;
Dir_alpha = info.Dir_alpha;
fehist = info.fehist;

hmm_best = hmm;
% Dir2d_alpha_best = Dir2d_alpha;
% Dir_alpha_best = Dir_alpha;
cyc_best = 1;
tp_less = max(hmm.train.embeddedlags) + max(-hmm.train.embeddedlags);

clear info;

% init stochastic learning stuff
nUsed = zeros(1,N);
sampling_weights = options.BIGbase_weights;
undertol = 0; count = 0;
Tfactor = N/options.BIGNbatch;
if options.downsample > 0
    downs_ratio = (options.downsample/options.Fs);
else
    downs_ratio = 1;
end

% Stochastic learning
for cycle = 2:options.BIGcyc

    % sampling batch
    I = datasample(1:N,options.BIGNbatch,'Replace',false,'Weights',sampling_weights);
    nUsed(I) = nUsed(I) + 1;
    nUsed = nUsed - min(nUsed) + 1;
    sampling_weights = options.BIGbase_weights.^nUsed;

    % read data for this batch
    Tbatch = []; Tbatch_list = cell(length(I),1);
    for ii = 1:length(I)
        i = I(ii);
        Tbatch = [Tbatch; ceil(downs_ratio*(T{i}-tp_less)) ];
        Tbatch_list{ii} = ceil(downs_ratio*(T{i}-tp_less));
    end

    [X,XX,Y] = loadfile(Xin(I),T(I),options);
    if hmm.train.lowrank > 0,  XX = X; end

    % local parameters (Gamma, Xi, P, Pi, Dir2d_alpha and Dir_alpha),
    % and free energy relative these local parameters
    tacc = 0; t2acc = 0;
    Gamma = cell(options.BIGNbatch,1); Xi = cell(options.BIGNbatch,1);
    XXGXX = cell(K,1);
    for k = 1:K, XXGXX{k} = zeros(size(XX,2)); end
    for ii = 1:length(I)
        i = I(ii);
        t = (1:sum(Tbatch_list{ii})) + tacc;
        t2 = (1:(sum(Tbatch_list{ii})-length(Tbatch_list{ii})*options.order)) + t2acc;
        tacc = tacc + sum(Tbatch_list{ii});
        t2acc = t2acc + sum(Tbatch_list{ii}) - length(Tbatch_list{ii})*options.order;
        data = struct('X',X(t,:),'C',NaN(sum(Tbatch_list{ii})-length(Tbatch_list{ii})*options.order,K));
        XX_i = XX(t2,:); Y_i = Y(t2,:);
        hmm_i = hmm;
        [Gamma{ii},~,Xi{ii}] = hsinference(data,Tbatch_list{ii},hmm_i,Y_i,[],XX_i); % state time courses
        checkGamma(Gamma{ii},Tbatch_list{ii},hmm_i.train,i);
        if isempty(Xi{ii}) % id_mixture
            Xi{ii} = approximateXi(Gamma{ii},Tbatch_list{ii},hmm_i);
        end
        for k = 1:K
            XXGXX{k} = XXGXX{k} + (XX_i' .* repmat(Gamma{ii}(:,k)',size(XX_i,2),1)) * XX_i;
        end
        % update transition probabilities
        Dir_alpha(:,i) = 0;
        for trial = 1:length(Tbatch_list{ii})
            t3 = sum(Tbatch_list{ii}(1:trial-1)) - options.order*(trial-1) + 1;
            Dir_alpha(:,i) = Dir_alpha(:,i) + Gamma{ii}(t3,:)';
        end
        Dir2d_alpha(:,:,i) = squeeze(sum(Xi{ii},1));
    end
    %if length(unique(hmm.train.grouping))==1
    %    hmm.Dir_alpha = sum(Dir_alpha,2)' + hmm.prior.Dir_alpha;
    %    hmm.Dir2d_alpha = sum(Dir2d_alpha,3) + hmm.prior.Dir2d_alpha;
    %else
    %    g = unique(hmm.train.grouping);
    %    for i = 1:length(g)
    %        hmm.Dir_alpha(:,i) = sum(Dir_alpha(:,hmm.train.grouping==g(i)),2) + hmm.prior.Dir_alpha';
    %        hmm.Dir2d_alpha(:,:,i) = sum(Dir2d_alpha(:,:,hmm.train.grouping==g(i)),3) + hmm.prior.Dir2d_alpha;
    %    end
    %end
    hmm.Dir_alpha = sum(Dir_alpha,2)' + hmm.prior.Dir_alpha;
    hmm.Dir2d_alpha = sum(Dir2d_alpha,3) + hmm.prior.Dir2d_alpha;
    [hmm.P,hmm.Pi] = computePandPi(hmm.Dir_alpha,hmm.Dir2d_alpha);
    subjfe(:,3,cycle) = evalfreeenergy([],[],[],[],hmm,[],[],[0 0 0 1 0]) / N; % "shared" P/Pi KL
    for ii = 1:length(I)
        i = I(ii); % Gamma entropy&LL
        subjfe(i,1:2,cycle) = evalfreeenergy([],Tbatch_list{ii},Gamma{ii},Xi{ii},hmm,[],[],[1 0 1 0 0]);
    end

    % global parameters (hmm), and collect state free energy
    rho(cycle) = (cycle + options.BIGdelay)^(-options.BIGforgetrate);
    Gamma = cell2mat(Gamma);
    if do_HMM_pca
        hmm_noisy = updatePCAparam (hmm,sum(Gamma),XXGXX,Tfactor,rangeK);
        hmm = states_supdate(hmm,hmm_noisy,rho(cycle),[1 2]);
    else
        % W
        if isfield(hmm.state(1),'W') && ~isempty(hmm.state(1).W.Mu_W)
            [hmm_noisy,XW] = updateW(hmm,Gamma,Y,XX,XXGXX,Tfactor);
            hmm = states_supdate(hmm,hmm_noisy,rho(cycle),1);
        end
        % Omega
        if isfield(hmm.state(1),'Omega') || isfield(hmm,'Omega')
            hmm_noisy = updateOmega(hmm,Gamma,Y,XX,XXGXX,XW,Tfactor);
            hmm = states_supdate(hmm,hmm_noisy,rho(cycle),2);
        end
        % Priors
        if pcapred
            % beta
            hmm_noisy = updateBeta(hmm);
            hmm = states_supdate(hmm,hmm_noisy,rho(cycle),5);
        else
            % sigma
            if ~isempty(orders)
                hmm_noisy = updateSigma(hmm);
                hmm = states_supdate(hmm,hmm_noisy,rho(cycle),3);
            end
            % alpha
            if ~isempty(orders)
                hmm_noisy = updateAlpha(hmm);
                hmm = states_supdate(hmm,hmm_noisy,rho(cycle),4);
            end
        end
    end

    % rest of the free energy (states' KL and data loglikelihood)
    [fe,ll] = evalfreeenergy(X,Tbatch,Gamma,cell2mat(Xi),hmm,Y,XX,[0 1 0 0 1]); % state KL & loglik
    statekl(cycle) = sum(fe(2:end));
    tacc = 0;
    for ii = 1:length(I)
        i = I(ii);
        t = (1:sum(Tbatch_list{ii})-length(Tbatch_list{ii})*options.order) + tacc; tacc = tacc + length(t);
        loglik(i,cycle) = sum(ll(t));
    end

    % bring from last iteration whatever was not updated
    for i = setdiff(1:N,I)
        loglik(i,cycle) = loglik(i,cycle-1);
    end

    fehist(cycle) = (-sum(loglik(:,cycle)) + statekl(cycle) + sum(sum(subjfe(:,:,cycle))));
    ch = (fehist(end)-fehist(end-1)) / abs(fehist(end)-fehist(1));
    if min(fehist)==fehist(cycle)
        hmm_best = hmm;
        cyc_best = cycle;
        %Dir2d_alpha_best = Dir2d_alpha; Dir_alpha_best = Dir_alpha;
        count = 0;
    else
        count = count + 1;
    end

    if options.BIGverbose
        fprintf('Cycle %d, free energy: %g (relative change %g), rho: %g \n', ...
            cycle,fehist(end),ch,rho(cycle));
    end

    %if hmm.train.tudamonitoring
    %    hmm.tudamonitor.synch(cycle+1,:) = getSynchronicity(Gamma,T);
    %    hmm.tudamonitor.accuracy(cycle+1,:) = getAccuracy(XX,residuals,T,Gamma,[],0);
    %end

    if cycle>5 && abs(ch) < options.BIGtol
        undertol = undertol + 1;
    else
        undertol = 0;
    end
    if cycle>options.BIGmincyc
        if undertol > options.BIGundertol_tostop, break; end
        if count >= options.BIGcycnobetter_tostop, break; end
    end

end

hmm = hmm_best;
fehist = fehist(1:cyc_best);
loglik = loglik(:,1:cyc_best);
subjfe = subjfe(:,:,1:cyc_best);
statekl = statekl(1:cyc_best);
rho = rho(1:cyc_best);

feterms = struct();
feterms.loglik = loglik;
feterms.subjfe = subjfe;
feterms.statekl = statekl;

for fn = fieldnames(options)'
    if ~isfield(hmm.train,fn{1})
        hmm.train.(fn{1}) = options.(fn{1});
    end
end

if options.BIGverbose
    if do_HMM_pca
        fprintf('Model: %d states, %d subjects, batch size %d, covariance: %s \n', ...
            K,length(T),options.BIGNbatch,hmm.train.covtype);
    else
        fprintf('Model: %d states, %d subjects, batch size %d, covariance: %s \n', ...
            K,length(T),options.BIGNbatch,hmm.train.covtype);
    end
    if hmm.train.exptimelag>1
        fprintf('Exponential lapse: %g, order %g, offset %g \n', ...
            hmm.train.exptimelag,hmm.train.order,hmm.train.orderoffset)
    else
        fprintf('Lapse: %d, order %g, offset %g \n', ...
            hmm.train.timelag,hmm.train.order,hmm.train.orderoffset)
    end
    if hmm.train.useMEX==0
        fprintf('MEX file was not used \n')
    else
        fprintf('MEX file was used for acceleration \n')
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [wordCounts] = updateCounts(wordCounts, data)
  wordCounts.word = wordCounts.word + data.numWords;
  if isfield(data, 'numChars') && data.numChars > 0
    wordCounts.char = wordCounts.char + data.numChars;
  end
end

author: lmthang


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the following exercises:
%   A) A simulation of a DC motor.
%
% Help: Function prototype
% [y, t] = runge_kutta(f, y0, [t0 tfinal], timestep)
% where f is the function being integrated as dy/dt = f(t, y).
% y0 are the initial conditions
% t0: initial time
% tfinal: final time.
% h: time step for the calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function simulate_motor_runge_kutta()
close all;
t0 = 0;
tfinal = 15;

[t, y] = runge_kutta(@motor, [0 0 0]', [t0 tfinal], 0.01);
y = y(:, 1:length(t));
plot(t, y(1,:), 'r'), hold
plot(t, y(2,:), 'g')
plot(t, y(3,:), 'b')
legend('Position theta (rad)', 'Speed (rad/s)', 'Current (A)')

figure,
plot(y(1,:), y(2,:), 'r')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function for Runge-Kutta.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = motor(t, x)
% We must return the solution of
% [dx1/dt; dx2/dt; dx3/dt dx4/dt]
Jm=500*1e-3; %kg*m2
R=0.566;%(Ohm)
L=0.172e-3;%(H)
Kb=1/(183*(pi/30)); % V/rad/s
Kp= 52.2/1000;
Bm = 0.001;
G=1;
tau = 0; % Nm
V = 5; % Volt

thetadd = (1/Jm)*(-Bm*x(2)+Kp*x(3)-tau/G);
currentd = (1/L)*(V-Kb*x(2)-R*x(3));

xd(1) = x(2);
xd(2) = thetadd;
xd(3) = currentd;
xd = xd(:);

repo: ARTE


text:
function [cl] = oz2cl(oz)
% Convert volume from US liquid ounces to centiliters.
% Chad Greene 2012
cl = oz*2.9573529563;

author: Sable


text:
function [Btu] = Nm2Btu(Nm)
% Convert energy or work from newton-meters to British thermal units.
% Chad A. Greene 2012
Btu = Nm*0.00094781707775;


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)  %'

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

p = round(sigmoid(X * theta));

% =========================================================================

end

author: gopaczewski

repo: coursera-ml


text:
function T=Entropy_Kapur(X,Trange);
hg=hist(X,Trange);
pg=hg./sum(hg);
Ti=0;
TrangeSpan=Trange(pg>0);
JT=zeros(length(TrangeSpan),1);
for Tk=TrangeSpan
    Ti=Ti+1;
    PT=sum(pg(Trange<=Tk));
    H=-pg./PT.*log10(pg./PT);
    isnanind=~isnan(H);
    Hf=sum(H(Trange<=Tk&isnanind));
    Hb=sum(H(Trange>Tk&isnanind));
    JT(Ti)=Hf+Hb;
end
[~,Tind]=max(JT);
T=TrangeSpan(Tind);
end

author: radishgiant

repo: ThresholdAndSegment


text:
function determ = carry_determinant ( n, alpha )

%*****************************************************************************80
%
%% CARRY_DETERMINANT returns the determinant of the CARRY matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ALPHA, the numeric base being used in the addition.
%
%    Output, real DETERM, the determinant.
%
  power = ( n * ( n - 1 ) ) / 2;
  determ = 1.0 / alpha^power;

  return
end


text:
function varargout = text2( varargin )
% Wrapper for text.m that ensures displayed text fits in figure.
%
% When text is called, Matlab displays the text, but does nothing to ensure
% that all of the text fits in the figure.  This function, after calling
% text, shrinks the axes until the text is fully visible.  Note that since
% font doesn't resize as the figure does, just because the text fully fits
% after text2 is called, no guaranties are made after the figure is
% resized.  Hence it is a good idea to resize the figure appropriately
% before calling text2.
%
% Same input/output options as Matlab's text.m command.
%
% USAGE
%  varargout = text2( varargin )
%
% INPUTS
%  varargin   - input to Matlab's text
%
% OUTPUTS
%  varargout  - output to Matlab's text
%
% EXAMPLE
%  figure(1); clf; text( -1 , 2, 'hello world' )  % will not appear
%  figure(2); clf; text2( -1 , 2, 'hello world' ) % will appear
%
% See also TEXT, IMLABEL
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

%%% call text normally
hs = text( varargin{:} );
if(nargout>=1); varargout={hs}; end;

%%% now make sure text fits
%%% update the OuterPosition (not regular Position) so good resize behavior
nhs = length(hs);
for i=1:nhs;  set(hs(i),'Units','Normalized');  end;
oldunits = get(gca,'Units'); set(gca,'Units','Normalized');
opos = get(gca,'OuterPosition'); %x y w h
updated = 1; updatedpos=1;
while( updated ); updated=0;
  for p=1:4
    if( updatedpos ); updatedpos=0;
      % get overall pos of text in normalized AXIS coordinates
      pos = [inf inf -inf -inf]; % [l b r t]
      for i=1:nhs
        extent = get(hs(i),'Extent'); % slow [l b w h]
        pos(1:2) = min( pos(1:2), extent(1:2) );
        pos(3:4) = max( pos(3:4), extent(1:2) + extent(3:4) );
      end;
      pos = pos * 1.05;
      pos(3:4) = pos(3:4) - pos(1:2);  % [l b w h]

      % convert pos to normalized FIGURE coordinates
      % must move origin, and also account for scaling
      axisPos = get(gca,'Position');
      pos(1) = axisPos(1) + pos(1) * axisPos(3);
      pos(2) = axisPos(2) + pos(2) * axisPos(4);
      pos(3) = pos(3) * axisPos(3);
      pos(4) = pos(4) * axisPos(4);
    end

    % if necessary, move OuterPosition of axes accordingly
    if( p<=2 && pos(p)<0 )
      step = min(.1,max(.001,-pos(p)/50));
      opos(p) = opos(p)+step;
      opos(p+2) = opos(p+2)-step;
      set(gca,'OuterPosition',opos);
      updated=1; updatedpos=1;
    elseif( p>2 && (pos(p-2)+pos(p))>1 )
      step = min(.1,max(.001,pos(p)/50));
      opos(p)=opos(p)-step;
      set(gca,'OuterPosition',opos);
      updated=1; updatedpos=1;
    end
  end
end
set(gca,'Units',oldunits);

%   if necessary, move OuterPosition of axes accordingly
%   if(p<=2 && pos(p)<0 )
%       step = max(.01,abs(pos(p))/1.5);
%       opos, step, pos, p
%       opos(p)=opos(p)+step;
%       opos(p+2)=opos(p+2)-step;
%       updated=1; updatedpos=1;
%       set(gca,'OuterPosition',opos);
%   elseif( p>2 && pos(p)>1 )
%       step = max(.01,(pos(p)-1)/1.5);
%       opos(p)=opos(p)-step;
%       updated=1; updatedpos=1;
%       set(gca,'OuterPosition',opos);
%   end

author: zouchuhang

repo: LayoutNet


text:
function [betas,stats] = weighted_glmfit(Y,varargin)
% Calculate weighted average using weighted linear least squares
% See examples below for usage
%
% Model:
% Y_i = 1*Ypop + noise
%
% :Inputs:
%
%   **Y:**
%        data matrix (nsub x T)
%
%   **w:**
%        weights
%
%   **varY:**
%        variance of data at each time point (nsub x T) + var between
%
% :Outputs:
%
%   **Ymean:**
%        weighted mean of each column of Y
%
%   **dfe:**
%        error degrees of freedom, adjusted for inequality of variance
%        (Sattherwaite) and pooled across data columns
%
% :Extended output in stats structure:
%
%   **stats.t:**
%        t-values for weighted t-test
%
%   **stats.p:**
%        2-tailed p-values for weighted t-test
%
%   **r:**
%        weighted correlation coeff across columns of Y
%
%   **xy:**
%        weighted covariance matrix
%
%   **v:**
%        weighted variance estimates for each column of Y
%          - sqrt(v) is the standard error of the mean (or grp difference)
%
%
%   **stats.fits:**
%        fits  for each group (Ymean by group), low contrast weight group then high
%
% Fastest if no stats are asked for.
%
% :Computation time:
% For FULL stats report
%   - Triples from 500 -> 1000 columns of Y, continues to increase
%
% For mean/dfe only, fast for full dataset (many columns of Y)
%
% :Examples:
%
% Basic multivariate stats for 1000 columns of dat, no weighting
% Multivariate covariances are meaningful if cols of Y are organized, e.g.,
% timeseries
% ::
%
%    [means,stats] = weighted_glmfit(dat(:,1:1000));
%
% The same, but return univariate stats only (good for large Y)
% ::
%
%    [means,stats] = weighted_glmfit(dat,'uni');
%
% A weighted version, where we put in the weights, and with a design matrix too:
% ::
%
%    [means,stats] = weighted_glmfit(X,dat,'uni','w',weights);
%
% A weighted version, where weights are determined from w/i subject variances:
% ::
%
%    [means,stats] = weighted_glmfit(X,dat,'uni','vary',variances);
%
% ..
%    NOTE: TOR CHANGED INPUT TO ASSUME THAT WE SHOULD ENTER VARWI + VARBETWEEN
% ..

if nargin == 0, error('Must at least enter data as 1st argument.'); end

% --------------------------------------
% * Set up arguments
% --------------------------------------

domultivariate = 0;     % multivariate covariance est for Y

zpdiff = []; w = []; varY = []; X = [];

for i = 1:length(varargin)
    arg = varargin{i};
    if ischar(arg)
        switch lower(arg)
            case 'w', w = varargin{i+1};
            case 'vary', varY = varargin{i+1};
            case 'uni', domultivariate = 0;
            case 'multi', domultivariate = 1;

            case {'X','x'}, X = varargin{i+1};
        end
    end
end

% fill in missing inputs with default values

if ~is_entered(w), w = ones(m,1);  end
if ~is_entered(varY), varY = ones(m,1); end
if ~is_entered(X), X = ones(m,1);  end

% Get rid of missing values
nancols = find(any(isnan(Y) | Y==0,1));
Y(:,nancols) = [];

[m,n] = size(Y);

% --------------------------------------
% * Means and contrast
% --------------------------------------

% pool weights across all voxels
[betas,invxwx,bform,fits] = get_betas_singleweight(X,Y,w);

if nargout == 1, return, end

% --------------------------------------
% * Residuals
% --------------------------------------

e = Y - fits;         % residuals

% --------------------------------------
% * Degrees of freedom
% --------------------------------------

[dfe,dfediff] = get_dfe(m,n,X,bform,varY,0);

if ~domultivariate
    % ======================================
    %
    %
    % Univariate stats: MSE, t, and p-values
    %
    %
    % ======================================

    % --------------------------------------
    % * Mean squared error
    % --------------------------------------

    W = diag(w);                    % Weight matrix
    % Loop version of MSE: avoids out of memory errors for large voxel sets
    MSE = zeros(1,n); for i=1:n, MSE(i) = e(:,i)'*W*e(:,i); end, MSE = MSE/dfe;

    k = size(betas,1);
    v = repmat(diag(invxwx),1,n) .* repmat(MSE,k,1);
    %v = invxwx * MSE;       % variances for mean

    % output
    stats.descrip1 = 'Univariate stats for test against zero:';
    stats.v = v;
    stats.v_descrip = 'V = ste^2; variance of mean estimate';
    stats.t = betas ./ sqrt(v);
    stats.p = 2 * ( 1 - tcdf(abs(stats.t),dfe) );
    stats.dfe = dfe;

else
    % ======================================
    %
    %
    % Multivariate stats: MSE, cov(Y), r(Y)
    % Useful for simulating t-values under dependence
    %
    % ======================================

    % --------------------------------------
    % * Mean squared error
    % --------------------------------------

    % additional output: covariance matrix for betas and zdiff across time
    % (columns)
    % and correlation matrix for betas and zdiff
    % used in Monte Carlo simulations for controlling false positives
    % across columns

    MSE = (e'*W*e)/dfe;             % Mean square error

    if dobtwn
        MSEdiff = (ediff'*W*ediff)/dfediff;
    end

    % --------------------------------------
    % * Estimated covariance and correlation
    %   Estimated between-subjects variance (v)
    % --------------------------------------

    xy = invxwx * MSE;           % Covariance matrix for betas;

    xy = 0.5*(xy+xy');              % Remove rounding error

    if dobtwn
        xydiff = inv(bcon'*W*bcon)*MSEdiff;           % Covariance matrix for betas;
        xydiff = 0.5*(xydiff+xydiff');
    end

    v = diag(xy);                   % Variance for betas

    if dobtwn
        vdiff = diag(xydiff);
    end

    r = xy./sqrt(v*v');             % Correlation matrix for betas

    if dobtwn
        rdiff = xydiff./sqrt(vdiff*vdiff');         % Correlation matrix for betas
    end

    stats.descrip1 = 'Multivariate stats for test against zero:';
    stats.r = r;
    stats.v = v;
    stats.xy = xy;

    stats.t = betas ./ sqrt(v');
    stats.p = 2 * ( 1 - tcdf(abs(stats.t),dfe) );

end

return

function [dfe,dfediff] = get_dfe(m,n,X,bform,varY,dobtwn,hatdiff,bcon)

dfediff = [];

% Set up residual-forming matrix
% --------------------------------------
dfe_v = zeros(n,1);
R = eye(m) - X*bform;    % residual inducing matrix

% contrast, if entered
if dobtwn
    dfe_vdiff = zeros(n,1);
    Rdiff = eye(m) - bcon * hatdiff;
end

% Calculate effective degrees of freedom
% --------------------------------------

have_unique_vars = size(varY,2) == n;

if ~have_unique_vars

    % Only one (pooled?) vector of variance estimates
    % --------------------------------------
    V = diag(varY(:,1));
    dfe = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation
    if dobtwn, dfediff = (trace(Rdiff*V)^2)/trace(Rdiff*V*Rdiff*V); end
else
    % Variance estimates for each data vector
    % --------------------------------------
    for i=1:n,

        % make diagonal matrix of variances
        V = diag(varY(:,i));

        dfe_v(i) = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation

        if dobtwn
            dfe_vdiff(i) = (trace(Rdiff*V)^2)/trace(Rdiff*V*Rdiff*V);
        end

    end

    dfe = mean(dfe_v);               % Calculate average df over all columns (pool over data vectors)
    if dobtwn
        dfediff = mean(dfe_vdiff);
    end

end

return

function bool = is_entered(x)

bool = exist('x','var') && ~isempty(x);

return

%
%
% Duplicated in robust_reg_pooled
%
%

function [betas,invxwx,bform,fits] = get_betas_singleweight(X,Y,w)

W = diag(w);                    % Weight matrix

%X = repmat(1,m,1);              % Design matrix - 1 column of all ones to calculate average
% and, separately, use bcon if that's entered

invxwx = inv(X'*W*X);
bform = invxwx * X'* W;         % beta-forming matrix.  hat = X * bform

% rows are columns of design (X), cols are Y variables
betas = bform*Y;

if nargout > 3
    fits = X * betas;
end

return

function w = bisquare_weight(r,radjust,xrank)
% r is residuals
% radjust is adjustment factor: DuMouchel & O'Brien
% xrank is rank of weighted X matrix (design)
% w is weights from bisquare function
% n is number of Y variables to replicate weights over
tuneconst = 4.685;

r = r .* radjust;
s = mad_sigma_pooled(r,xrank);
r = r ./ (s*tuneconst);
w = (abs(r)<1) .* (1 - r.^2).^2;
return

function s = mad_sigma_pooled(r,xrank)
%    Compute std estimate using MAD of residuals from 0
rsort = sort(abs(r));
rsort = rsort(xrank:end,:); % eliminate smallest; like reducing df
s = median(rsort(:)) / 0.6745;
return

function str = display_string(str)
str = sprintf(str); fprintf(1,'%s',str);
return

function erase_string(str)

len = length(str);
str2 = repmat('\b',1,len);

fprintf(1,str2);

return

author: canlab

repo: CanlabCore


text:
function geometry_test0202 ( )

%*****************************************************************************80
%
%% TEST0202 tests CYLINDER_POINT_INSIDE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 6;

  inside_test = [ 0, 1, 0, 0, 1, 1 ];
  p_test = [ ...
      4.0,   0.5,  0.0; ...
     -0.5,  -1.0,  0.0; ...
      4.0,   6.0,  0.0; ...
      0.75, -10.0, 0.0; ...
      0.0,   0.0,  0.0; ...
      0.25,  1.75, 0.0 ]';
  p1 = [ 0.0, -2.0, 0.0 ];
  p2 = [ 0.0,  2.0, 0.0 ];
  r = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0202\n' );
  fprintf ( 1, '  CYLINDER_POINT_INSIDE_3D determines if a point is\n' );
  fprintf ( 1, '  inside a cylinder.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radius R = %f\n', r );
  fprintf ( 1, '  Center of bottom disk = %10f  %10f  %10f\n', p1(1:dim_num) );
  fprintf ( 1, '  Center of top disk =    %10f  %10f  %10f\n', p2(1:dim_num) );

  for test = 1 : test_num

    p(1:dim_num) = p_test(1:dim_num,test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  P =    %10f  %10f  %10f\n', p(1:dim_num) );

    inside = cylinder_point_inside_3d ( p1, p2, r, p );

    fprintf ( 1, '  INSIDE (computed) = %d\n', inside );
    fprintf ( 1, '  INSIDE (exact) =    %d\n', inside_test(test) );

  end

  return
end


text:
function S = globMatrixPGIFE3DMass(fun,mesh,femI,fem1,fem2)

%% USAGE: generate stiffness global matrix on a 3D mesh
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem1 --- global DoF for test function space
% fem2 --- global DoF for trial function space
%
% OUTPUTS:
% [IN JN XN] --- triplets of the sparse matrix from regular elements.
% [II JI XI] --- triplets of the sparse matrix from interface elements.

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Initializaiton
if strcmp(fem1.type,'P1')||strcmp(fem1.type,'DGP1')||strcmp(fem1.type,'CR')
    feEvalBas1 = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem1.type,'DGP2')
    feEvalBas1 = @evalP2Bas3D;
end

if strcmp(fem2.type,'P1')||strcmp(fem2.type,'DGP1')||strcmp(fem2.type,'CR')
    feEvalBas2 = @evalP1Bas3D;
elseif strcmp(fem2.type,'P2')||strcmp(fem2.type,'DGP2')
    feEvalBas2 = @evalP2Bas3D;
end

%% 1. Matrix on noninterface elements
dof1 = fem1.ldof; dof2 = fem2.ldof; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);
AN = fem1.area(ntID);
gxN = fem1.gx(ntID,:); gyN = fem1.gy(ntID,:); gzN = fem1.gz(ntID,:); gw = fem1.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(fun,gxN,gyN,gzN);
Ibas = cell(dof1,1);
Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(fem1.bas(ntID,:,i), gxN, gyN, gzN, [0,0,0]);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(fem2.bas(ntID,:,j), gxN, gyN, gzN, [0,0,0]);
end

IN = reshape(repmat(fem1.t(ntID,:),4,1),nloc*ntN,1);
JN = repmat(reshape(fem2.t(ntID,:),dof2*ntN,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibas{i}.*(coefN.*Jbas{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
SN = sparse(IN(ID),JN(ID),XN(ID),size(fem1.p,1),size(fem2.p,1));

%% 2. Matrix on interface elements
AI = femI.area; gw = femI.gw; gxI = femI.gx; gyI = femI.gy; gzI = femI.gz;
ntI = size(femI.t,1); % not number of interface element, but quadrature element
XI = zeros(nloc*ntI, 1);
tID = femI.QelemID;

coefI = feval(fun,gxI,gyI,gzI);
Ibas = cell(dof1,1);
Jbas = cell(dof2,1);
for i = 1:dof1
    Ibas{i} = feEvalBas1(fem1.bas(tID,:,i), gxI, gyI, gzI, [0,0,0]);
end
for j = 1:dof2
    Jbas{j} = feEvalBas2(femI.bas(:,:,j), gxI, gyI, gzI, [0,0,0]);
end

II = reshape(repmat(femI.t,4,1),nloc*ntI,1);
JI = repmat(reshape(femI.t,dof2*ntI,1),4,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XI(ind+1:ind+ntI) = AI.*(sum(((Ibas{i}.*(coefI.*Jbas{j})).*gw'),2));
        ind = ind + ntI;
    end
end
SI = sparse(II,JI,XI,size(fem1.p,1),size(fem2.p,1));
S = SN + SI;

repo: ifem