text:
function test_example_bids_meg

% MEM 4gb
% WALLTIME 00:10:00

%
%% Converting an example MEG dataset for sharing in BIDS
%
%
% The general description of how MEG data can be shared in the [BIDS format](https://bids.neuroimaging.io) is described in the paper [MEG-BIDS, the brain imaging data structure extended to magnetoencephalography](https://doi.org/10.1038/sdata.2018.110) by Guiomar Niso, et al. Scientific Data volume 5, Article number: 180110 (2018). The BIDS standard itself with all details can be found on <https://bids-standard.org>.
%
%% # CTF
%
% The CTF MEG file format consists of recordings that are represented in a directory with the extension `*.ds`. The directory contains multiple files, where the `*.res4` contains most of the header information, and the `*.meg4` files (there can be multiple) contain the data. However, also the other files contain relevant data and/or metadata, such as markers and information on how the data has been processed.
%
% In general you **should not** simply rename the directory, or mess with the files in the directory by hand, since there are multiple files that need to have a consistent file naming structure. The CTF Linux command-line software includes tools for anonimizing and renaming datasets in a consistent fashion and [this example page](/faq/how_can_i_anonymize_a_ctf_dataset/) also presents a way of doing it directly in MATLAB.
%
%% ## Example
%
% This dataset corresponds to the paper by Erik te Woerd et al. [Impaired auditory-to-motor entrainment in Parkinson’s disease](https://doi.org/10.1152/jn.00547.2016). J Neurophysiol. 2017 May 1;117(5):1853-1864. doi: 10.1152/jn.00547.2016.
%
% The complete dataset consists of MEG data recorded from 15 Parkinson patients and 15 control participants. In three of the sessions the recording was interrupted, causing those sessions to be split over 2 runs. Consequently, there are 33 CTF `*.ds` dataset directories to start with. The data from two participants was excluded from further analysis.
%
% Control subjects are coded from S101 to S115, patients are coded from S201 to S215. Besides the pseudonimized subject identifier, there is also a run number in the original file name (which is relevant for three of the recordings), the project number, the experimenter, and the date at which the recording was done. Since the date might be recognizeable for the participants and/or their family members, we will treat that as potentially identifying information and scrub it from the dataset.
%
% The data is anonimized using the `go_anonymiseDs.m` function from [this example page](/faq/how_can_i_anonymize_a_ctf_dataset/#using-matlab). This reads in the data with a low-level CTF function, scrubs potentially identifying  information, and writes it with a low-level CTF function. Some of the auxiliary files that can optionally be present (for example when the data was opened with the CTF DataEditor application) and the headlocalizer data will be ignored upon reading, and hence also not appear in the copy that is written.
%
% The anonimized version of the CTF data has a short descriptive name. Using the **[data2bids](https://github.com/fieldtrip/fieldtrip/blob/release/data2bids.m)** we copy and rename the CTF data so that it has the right directory structure and file names, and we add the metadata.
%
%
subj = {
  'S101'
  'S102'
  'S103'
  'S104'
  'S105'
  'S106'
  'S107'
  'S108'
  'S109'
  'S110'
  'S111'
  'S112'
  'S113'
  'S114'
  'S115'
  'S201'
  'S202'
  'S203'
  'S204'
  'S205'
  'S206'
  'S207'
  'S208'
  'S209'
  'S210'
  'S211'
  'S212'
  'S213'
  'S214'
  'S215'
  };

for subindx=1:numel(subj)
  for runindx=1:2

    % the number 3018009.04 in the dataset name refers to the DCCN project identifier
    % all data acquisition was done using study-specific participant identifiers rather than names
    % the CTF dataset name contains the recording date, but here we use a wildcard instead

    d = dir(sprintf('../original/301800904eritwoe%s00%d_1200hz_*_01.ds', subj{subindx}, runindx));
    if isempty(d)
      % for most subjects the data was recorded in a single run
      % in that case run 2 does not exist
      continue
    else
      origname = fullfile(d.folder, d.name);
      anonname = fullfile(d.folder, sprintf('%s_%d.ds', subj{subindx}, runindx));
      disp(anonname); % this is just an intermediate name, the final name will be assigned by data2bids
    end

    if ~isfolder(anonname)
      go_anonymiseDs(origname, anonname);
    end

    cfg = [];

    cfg.bidsroot = '../bids';  % write to the present working directory
    cfg.sub = subj{subindx};
    cfg.run = runindx;
    cfg.dataset = anonname; % this is the intermediate name

    cfg.datatype = 'meg';
    cfg.method = 'copy'; % the original data is in a BIDS-compliant format and can simply be copied

    cfg.InstitutionName             = 'Radboud University';
    cfg.InstitutionalDepartmentName = 'Donders Institute for Brain, Cognition and Behaviour';
    cfg.InstitutionAddress          = 'Kapittelweg 29, 6525 EN, Nijmegen, The Netherlands';

    % required for dataset_description.json
    cfg.dataset_description.Name                = 'Impaired auditory-to-motor entrainment in Parkinson''s disease';
    cfg.dataset_description.BIDSVersion         = 'v1.5.0';

    % optional for dataset_description.json
    cfg.dataset_description.License             = 'RU-DI-HD-1.0';
    cfg.dataset_description.Authors             = 'Te Woerd ES, Oostenveld R, de Lange FP, Praamstra P';
    cfg.dataset_description.ReferencesAndLinks  = {'https://doi.org/10.1152/jn.00547.2016'};
    cfg.dataset_description.EthicsApprovals     = 'DCCN 3018009.04';

    cfg.TaskDescription = 'The experiment consisted of an auditory target detection task. Three cue tones of different frequencies (900, 1,100 and 1,300 Hz; 50-ms duration) predicted the probability (10, 30, 50%) of the next stimulus being the target tone (2,000 Hz; 50-ms duration) (Fig. 1A). The pitch of standard tones was chosen to be separated by 200 Hz and the target differed 700 Hz from the highest standard tone, which enabled subjects to detect the target tone easily. Participants were informed about the meaning of the standard tones and instructed to press the response button as swift as possible with the index finger of their dominant hand.';
    cfg.task = 'auditory';

    try
      data2bids(cfg);
    catch
      % this is probably because the output dataset already exists
      % this is due to running the script multiple times
      disp(lasterr)
    end

  end % for each run
end % for each subject
%
% Following this conversion, the `participant.tsv` file only contains the subject identifiers. There is more information in an Excel file, such as the age, sex, clinical state, questionaire data, lab notes. In principle it would be possible to copy that into the MATLAB script and use data2bids to export that, but in this case it was easier to use "save-as" and then select "tab-delimited-text" as the output format.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function estimation_results = FEKF_SLAM(data)
% FEJ-EKF SLAM
% load pre-given data: odometry and observations
if nargin < 1
    load('./data.mat');
end

data_matrix = data.state;
odo_cov = data.odom_cov;   % constant variable
obs_cov = data.obse_cov;   % constant variable

odom_sigma = data.odom_sigma;
obsv_sigma = data.obsv_sigma;

%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% In every step, all elements of Estimation_X will be changed %%%%%%%%%%%%
estimation_x.orientation = data.poses.orientation(1:3,1:3);
estimation_x.position    = data.poses.position(:,1);
estimation_x.cov         = sparse(6,6);
estimation_x.landmarks   = [];       % the landmarks observed until this step (included), 4*N format, the 4-th row is the index
estimation_x0.IndexObservedNew=[];
estimation_x0.IndexObservedAlreadyThis=[];
%Estimation_X.IndexOfFeature=[];     % the names(indexes) of the landmarks observed until this step (included)
%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%
FirstPosition=[]; % store  p_{n|n-1}
FirstLandmarks=[]; % store f_{k0}

% Initialize
n_steps = max(data_matrix(:,4));  % step instead of pose,  hence, it does not include pose 0
estimation_results = cell(1, n_steps+1);
estimation_results{1} = estimation_x;

for i = 0:n_steps
    IndexOfCurrentStepInDataMatrix = find(data_matrix(:,4) == i);
    m = size(IndexOfCurrentStepInDataMatrix, 1);
    if ( mod(i, 50) == 0 )
        disp(['Processing pose ', int2str(i)]);
    end
    % det(Estimation_X.cov)
    if i ~= n_steps
        OdometryFromThis2Next = data_matrix(IndexOfCurrentStepInDataMatrix(m-5):IndexOfCurrentStepInDataMatrix(m),1);
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 1 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 3 )];
            [estimation_x,FirstLandmarks] = FEKF_update(estimation_x, CameraMeasurementThis, obsv_sigma, FirstLandmarks );
        end

        estimation_results{i+1} = estimation_x;

        [estimation_x, FirstPosition] = FEKF_propagate(estimation_x, OdometryFromThis2Next, odom_sigma, FirstPosition );

    else
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 1 ) , data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 3 )];
           [estimation_x,FirstLandmarks] = FEKF_update(estimation_x, CameraMeasurementThis, obsv_sigma, FirstLandmarks);
        end
        estimation_results{i+1} = estimation_x;
    end
end

clearvars -except estimation_results;

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
% The COBRAToolbox: testMoomin.m
%
% Purpose:
%     - Tests the moomin function
%
% Author:
%     - Original file: Taneli Pusa 01/2020

global CBTDIR

% initialize the test
currentDir = pwd;
fileDir = fileparts(which('testMoomin.m'));
cd(fileDir);
testPath = pwd;

% glpk appears to be extremely slow with stoichiometric constraints
% gurobi needs further testing to ensure it works
solvers = prepareTest('needsMILP', true, 'requiredSolvers', {'ibm_cplex'}, 'excludeSolvers',...
	 {'glpk', 'gurobi'});

% load model and data for the test
inputModel = getDistributedModel('ecoli_core_model.mat');
load([testPath filesep 'testData_moomin.mat']);

nMets = size(inputModel.mets, 1);
nRxns = size(inputModel.rxns, 1);

% tolerance
tol = 1e-6;

% run tests with different solvers
for i=1:length(solvers.MILP)
	fprintf(' -- Running testMoomin using the solver interface: %s ... ', solvers.MILP{i});

	solverOK = changeCobraSolver(solvers.MILP{i}, 'MILP', 0);

	if solverOK
		% test stoichiometric version
		[outputModel, MILPsolutions, MILPproblem] = moomin(inputModel, expression, ...
			'enumerate', 5);

		% check that reading GPRs works (solver independent, only needs to be done once)
		if i==1
			assert(any(outputModel.inputColours));
		end

		% there should be only one solution
		assert(size(MILPsolutions, 1) == 2 && MILPsolutions{2, 1}.stat == 0);

		% check the MILP
		A = MILPproblem.A;
		x = MILPsolutions{1,1}.full;
		% stoichiometric balance
		assert(all(abs(A(1:nMets, 1:nRxns) * x(1:nRxns)) < tol));
		% binary variables worked as intended
		assert(all( abs((x(1:nRxns) > 0) - x(nRxns + 1:2 * nRxns)) < tol ));
		assert(all( abs((x(1:nRxns) < 0) - x(2 * nRxns + 1:end)) < tol ));

		% test the topological version
		[outputModel, MILPsolutions, MILPproblem] = moomin(inputModel, expression, ...
			'enumerate', 5, 'stoichiometry', 0);

		% there should be only one solution
		assert(size(MILPsolutions, 1) == 2 && MILPsolutions{2, 1}.stat == 0);

		% check the MILP
		A = MILPproblem.A;
		b = MILPproblem.b;
		x = MILPsolutions{1,1}.full;
		assert(all( A(1:nRxns + nMets, :) * x - b(1:nRxns + nMets) < tol ));
		assert(all( A(nRxns + nMets + 1:nRxns + 3 * nMets, :) * x -...
			b(nRxns + nMets + 1:nRxns + 3 * nMets) > -tol ));
	end
	fprintf('Done.\n');
end

% delete solver log files
delete('*.log');

cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function qqplot(o)
% quantile-quantile of misorientation angle against random angular
% misorientation
%
% Example
%   cs = crystalSymmetry('-43m');
%   odf1 = unimodalODF(orientation.id(cs),'halfwidth',20*degree);
%   odf2 = unimodalODF(orientation.id(cs),'halfwidth',50*degree);
%
%   qqplot(odf1.discreteSample(1000))
%   qqplot(odf2.discreteSample(1000))
%

angles = o.angle;
[pdf,omegas] = calcAngleDistribution(o.CS,o.SS);

pdf = cumsum(pdf);
pdf = pdf./pdf(end);

omegas(pdf == 0) = [];
pdf(pdf == 0) = [];

mof = cumsum(hist(angles,omegas));
mof = mof./mof(end);

figure
plot(pdf,mof,'.')
line([0 1],[0 1])

axis equal tight

author: mtex-toolbox

repo: mtex

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


text:
function pagechan(page_num,mode)
% PAGECHAN Change pages within environment.
%          PAGECHAN handles the transition between pages within the
%          Page Setup... environment.

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

% obtain handle information
f2 = gcf;
f = get(f2,'userdata');
ui_data = get(f,'userdata');
ui_han = ui_data{1};
ui_han2 = get(ui_han(49),'userdata');

% variables internal to Page Setup... interface
page_setup = get(ui_han2(1),'userdata');
ptch = get(ui_han2(2),'userdata');
txt = get(ui_han2(3),'userdata');

set(txt,'string',['Page ',int2str(page_num),' Layout']);

% first, how many plots on this page
plots = page_setup(:,page_num);
plots(find(plots==0)) = [];
set(ui_han2(4),'userdata',length(plots));

color_mat = [1,0,0;0,1,0;0,0,1;1,1,0;0,1,1;0.5,0,0;0,0.5,0;0,0,0.5];
if ~mode,
 set(ui_han2(1),'value',page_num);
 set(ui_han2(2:9),'back','w','value',0);
 for k = 1:length(plots),
  set(ui_han2(plots(k)+1),'back',color_mat(plots(k),:),'value',1);
 end
end

% locations of patchs
fig_wdt = 780;
frm_wdt = 190;

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

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

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

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

set(ptch,'vis','off');
if length(plots)==1,
 xpos1 = [axs1_l, axs1_l+axs1_w, axs1_l+axs1_w, axs1_l];
 ypos1 = [axs1_b, axs1_b, axs1_b+axs1_h, axs1_b+axs1_h];
 set(ptch(1),'vis','on','facecolor',get(ui_han2(plots(1)+1),'back'),...
             'xdata',xpos1,'ydata',ypos1);

elseif length(plots)==2,
 xpos1 = [axs1_l, axs1_l+axs1_w, axs1_l+axs1_w, axs1_l];
 ypos1 = [axs2_b, axs2_b, axs2_b+axs2_h, axs2_b+axs2_h];
 ypos2 = [axs1_b, axs1_b, axs1_b+axs2_h, axs1_b+axs2_h];
 set(ptch(1),'vis','on','facecolor',get(ui_han2(plots(1)+1),'back'),...
             'xdata',xpos1,'userdata',plots(1),'ydata',ypos1);
 set(ptch(2),'vis','on','facecolor',get(ui_han2(plots(2)+1),'back'),...
             'xdata',xpos1,'userdata',plots(2),'ydata',ypos2);

elseif length(plots)==3,
 xpos1 = [axs1_l, axs1_l+axs1_w, axs1_l+axs1_w, axs1_l];
 ypos1 = [axs2_b, axs2_b, axs2_b+axs2_h, axs2_b+axs2_h];
 xpos2 = [axs1_l, axs1_l+axs2_w, axs1_l+axs2_w, axs1_l];
 ypos2 = [axs1_b, axs1_b, axs1_b+axs2_h, axs1_b+axs2_h];
 xpos3 = [axs2_l, axs2_l+axs2_w, axs2_l+axs2_w, axs2_l];
 ypos3 = [axs1_b, axs1_b, axs1_b+axs2_h, axs1_b+axs2_h];
 set(ptch(1),'vis','on','facecolor',get(ui_han2(plots(1)+1),'back'),...
             'xdata',xpos1,'userdata',plots(1),'ydata',ypos1);
 set(ptch(2),'vis','on','facecolor',get(ui_han2(plots(2)+1),'back'),...
             'xdata',xpos2,'userdata',plots(2),'ydata',ypos2);
 set(ptch(3),'vis','on','facecolor',get(ui_han2(plots(3)+1),'back'),...
             'xdata',xpos3,'userdata',plots(3),'ydata',ypos3);

elseif length(plots)==4,
 xpos1 = [axs1_l, axs1_l+axs2_w, axs1_l+axs2_w, axs1_l];
 ypos1 = [axs2_b, axs2_b, axs2_b+axs2_h, axs2_b+axs2_h];
 xpos2 = [axs2_l, axs2_l+axs2_w, axs2_l+axs2_w, axs2_l];
 ypos2 = [axs2_b, axs2_b, axs2_b+axs2_h, axs2_b+axs2_h];
 xpos3 = [axs1_l, axs1_l+axs2_w, axs1_l+axs2_w, axs1_l];
 ypos3 = [axs1_b, axs1_b, axs1_b+axs2_h, axs1_b+axs2_h];
 xpos4 = [axs2_l, axs2_l+axs2_w, axs2_l+axs2_w, axs2_l];
 ypos4 = [axs1_b, axs1_b, axs1_b+axs2_h, axs1_b+axs2_h];
 set(ptch(1),'vis','on','facecolor',get(ui_han2(plots(1)+1),'back'),...
             'xdata',xpos1,'userdata',plots(1),'ydata',ypos1);
 set(ptch(2),'vis','on','facecolor',get(ui_han2(plots(2)+1),'back'),...
             'xdata',xpos2,'userdata',plots(2),'ydata',ypos2);
 set(ptch(3),'vis','on','facecolor',get(ui_han2(plots(3)+1),'back'),...
             'xdata',xpos3,'userdata',plots(3),'ydata',ypos3);
 set(ptch(4),'vis','on','facecolor',get(ui_han2(plots(4)+1),'back'),...
             'xdata',xpos4,'userdata',plots(4),'ydata',ypos4);
end

author: Sable

repo: mcbench-benchmarks


text:
% This demo compares exponential vs hyperbolic temporal discounting.
% This script simulates an agent making choices between two options, which
% differ in terms of the expected reward and temporal horizons, under both
% an hyperbolic and an exponantial temporal discounting model.
% The choice data are then inverted using these two models. NB: in both
% models, the discount factor and the behavioural temperature are unknown.
% The effect of simulating the data under each model onto the log-Bayes
% factor (comparison between hyperbolic vs exponential models) is assessed
% using Monte-Carlo simulations.

clear all
close all

models = {'hyperbolic','exponential'};

% initialize simulations
ntrials = 128;
g_fname = @g_discounting;
phi = [1;-1];
in.ind.logk = 1; % index of discount factor (phi)
in.ind.logb = 2; % index of behavioural temperature (phi)
in.ind.t = [1;3]; % index of temporal horizon (u)
in.ind.R = [2;4]; % index of expected reward (u)

% Build options and dim structures for model inversion
dim.n_theta = 0;
dim.n_phi = 2;
dim.n = 0;
dim.p = 1;
dim.n_t = ntrials;
options.inG = in;
options.dim = dim;
options.sources = struct('type',1 ,'out', 1); % one binomial observation;
options.DisplayWin = 0;

% Build time series of hidden states and observations
Nmcmc = 16;
nm = length(models);
p = cell(nm,nm,Nmcmc);
o = cell(nm,nm,Nmcmc);
F = zeros(nm,nm,Nmcmc);
for ii=1:Nmcmc
    R = exp(randn(2,ntrials));
    t = exp(randn(2,ntrials));
    u = zeros(4,ntrials);
    u(in.ind.t,:) = t;
    u(in.ind.R,:) = R;
    for i=1:nm
        options.inG.model = models{i};
        [y,x,x0,eta,e] = VBA_simulate (ntrials,[],g_fname,[],phi,u,[],[],options);
        % Call inversion routine
        for j=1:nm
            options.inG.model = models{j};
            [p{i,j,ii},o{i,j,ii}] = VBA_NLStateSpaceModel(y,u,[],g_fname,dim,options);
            F(i,j,ii) = o{i,j,ii}.F;
        end
    end
end

hf = figure('color',[1 1 1]);
ha = axes('parent',hf);
dF = F(:,1,:) - F(:,2,:);
mdF = mean(dF,3);
vdF = var(dF,[],3)./Nmcmc;
plotUncertainTimeSeries(mdF,vdF,[],ha);
set(ha,'xlim',[0,3],'xtick',[1,2],'xticklabels',models)
xlabel(ha,'type of simulated data')
ylabel(ha,['log p(y|',models{1},') - log p(y|',models{2},')'])
box(ha,'off')

author: MBB-team

repo: VBA-toolbox


text:
% fmri_data.subs

% This is different from fmri_data.apply_parcellation, which returns
% averages or weighted averages per image for a set of images in a data
% object.  apply_parcellation can be used to return local pattern responses
% for any parcellation, but it does not return vectors with the local
% weights themselves.

% expand into atlas subregions
% Take an image object (one image) and an integer-valued parcellation/atlas whose values
% indicate k regions, and expand the image object into a set of k images with
% the original values only for voxels in the kth region, and zero
% elsewhere. Useful for breaking a pattern mask into local atlas-defined subregions.

% define target image obj to sample to : nps

nps = replace_empty(nps); % just in case
target_obj = nps;  % this will be expanded to set of images

orig_image_vector = target_obj.dat(:, 1);

parcellation_file = which('BN_Atlas_274_noCb_uint16.nii');
parcel_obj = fmri_data(parcellation_file);

parcel_obj = resample_space(parcel_obj, nps);

names_file = fullfile(fileparts(parcellation_file), 'cluster_names.mat');
clear names
load(names_file, 'names');
names;      % returns names var. make sure: break if not

% Find which names match
%wh = ~cellfun(@isempty, strfind(names, string_to_find));

% Index values to look for
% indx = find(wh)';

%% Do expansion
% -----------------------------------------------------------------------
for indx = 1:length(names)

    output_obj = select_voxels_by_value(parcel_obj, indx);

    target_subregion_vector = orig_image_vector; % will become image values only in region, 0 elsewhere
    target_subregion_vector(output_obj.removed_voxels) = 0;

    target_obj.dat(:, indx) = target_subregion_vector;

end

% Add image for everything not already covered
% *******

% Remove empty images, attach names to object

empty_imgs = sum(abs(target_obj.dat)) < 10*eps;
target_obj.dat(:, empty_imgs) = [];

names(empty_imgs) = [];
target_obj.additional_info = names;

%% Get region object, r
% May have different length (n regions) because of contiguity requirement
% -----------------------------------------------------------------------

dat_wh = get_wh_image(target_obj, 1);
r = region(dat_wh);

region_names = repmat(names(1), 1, length(r));

for i = 2:size(target_obj.dat, 2)

    dat_wh = get_wh_image(target_obj, i);

    % threshold at 20% probability, get rid of stray voxels
    %dat_wh = threshold(dat_wh, [.2 Inf], 'raw-between', 'k', 3);

    r_wh = region(dat_wh);

    for j = 1:length(r_wh)
        [r_wh(j).title, r_wh(j).shorttitle] = deal(names{i});
    end

    r = [r r_wh];

    region_names = [region_names repmat(names(i), 1, length(r_wh))];

end
%%

% similarity_output = canlab_pattern_similarity(dat, pattern_weights, varargin)

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ctr = ctrCreate(user_defaults)
% Create ConTrac (ctr) structure
%
%   ctr = ctrCreate(user_defaults)
%
% Creates a structure of parameters needed to run ConTrac/MetroTrack.
%
% Originally by BW?
%
% Modified by DY 7/25/2008 to include calls to ctrSet to install default
% values and allows for the user to pass in a new set of default values
% (user_defaults) if they want, but more needed to make this work. Right
% now, they pass in default values to ctrSet according to what I happen to
% want to be defaults at the moment.

% Initialize the CTR struct
ctr=struct;

% Set the defaults
% If user_defaults passed in, use them. Otherwise, use defaults below.
if ~notDefined('user_defaults')
    d=user_defaults;
else
    d.version=1; % We are at version 1 on July 18, 2008
    d.image_directory='none';
    d.tensors_filename='none';
    d.fa_filename='none';
    d.pdf_filename='none';
    d.mask_filename='none';
    d.Xmask_filename='none';
    d.script_filename='none';
    d.compute_fa='false';
    d.require_way='false';
    d.sampler_type     = 'SISR';
    d.desired_samples  = 50000;
    d.burnin           = 1;
    d.max_nodes    = 240;
    d.min_nodes    = 3;
    d.step_size    = 1;
    d.skip_samples = 1;
    d.roi{1}.coords = [1 1 1; 2 2 2];
    d.roi{2}.coords = [1 1 1; 3 3 3];
    d.roi{1}.valid_cortex = 'true';
    d.roi{2}.valid_cortex = 'true';
    d.roi{1}.seed_region  = 'true';
    d.roi{2}.seed_region  = 'true';
    d.translate = 0;
    d.isr       = 0;
    d.esr       = 1;
    d.temp_swap = 0.1;
    d.inv_temp  = 1;
    d.start_path_tries = 1;
    d.save_out_spacing = 50;
    d.fa_absorb   = 0.01;
    d.abs_normal  = 0;
    d.abs_penalty = 0;
    d.smooth_std  = 14;
    d.angle_cutoff = 130;
    d.shape_params_vec = [0.175 0.15 100];
end

ctr=ctrSet(ctr,'version',d.version);
ctr=ctrSet(ctr,'tensors_filename',d.tensors_filename);
ctr=ctrSet(ctr,'image_directory',d.image_directory);
ctr=ctrSet(ctr,'fa_filename',d.fa_filename);
ctr=ctrSet(ctr,'pdf_filename',d.pdf_filename);
ctr=ctrSet(ctr,'mask_filename',d.mask_filename);
ctr=ctrSet(ctr,'Xmask_filename',d.Xmask_filename);
ctr=ctrSet(ctr,'script_filename',d.script_filename);
ctr=ctrSet(ctr,'compute_fa',d.compute_fa);
ctr=ctrSet(ctr,'require_way',d.require_way);
ctr=ctrSet(ctr,'sampler_type',d.sampler_type);
ctr=ctrSet(ctr,'desired_samples',d.desired_samples);
ctr=ctrSet(ctr,'burnin',d.burnin);
ctr=ctrSet(ctr,'max_nodes',d.max_nodes);
ctr=ctrSet(ctr,'min_nodes',d.min_nodes );
ctr=ctrSet(ctr,'step_size',d.step_size);
ctr=ctrSet(ctr,'skip_samples',d.skip_samples);
ctr=ctrSet(ctr,'roi',d.roi{1}.coords,1,'coords');
ctr=ctrSet(ctr,'roi',d.roi{2}.coords,2,'coords');
ctr=ctrSet(ctr,'roi',d.roi{1}.valid_cortex,1,'valid_cortex');
ctr=ctrSet(ctr,'roi',d.roi{2}.valid_cortex,2,'valid_cortex');
ctr=ctrSet(ctr,'roi',d.roi{1}.seed_region,1,'seed_region');
ctr=ctrSet(ctr,'roi',d.roi{2}.seed_region,2,'seed_region');
ctr=ctrSet(ctr,'translate',d.translate);
ctr=ctrSet(ctr,'isr',d.isr);
ctr=ctrSet(ctr,'esr',d.esr);
ctr=ctrSet(ctr,'temp_swap',d.temp_swap);
ctr=ctrSet(ctr,'inv_temp',d.inv_temp);
ctr=ctrSet(ctr,'start_path_tries',d.start_path_tries);
ctr=ctrSet(ctr,'save_out_spacing',d.save_out_spacing);
ctr=ctrSet(ctr,'fa_absorb',d.fa_absorb);
ctr=ctrSet(ctr,'abs_normal',d.abs_normal);
ctr=ctrSet(ctr,'abs_penalty',d.abs_penalty);
ctr=ctrSet(ctr,'smooth_std',d.smooth_std);
ctr=ctrSet(ctr,'angle_cutoff',d.angle_cutoff);
ctr=ctrSet(ctr,'shape_params_vec',d.shape_params_vec);
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function varargout = MU_Plot(varargin)
% MU_PLOT MATLAB code for MU_Plot.fig
%      MU_PLOT, by itself, creates a new MU_PLOT or raises the existing
%      singleton*.
%
%      H = MU_PLOT returns the handle to a new MU_PLOT or the handle to
%      the existing singleton*.
%
%      MU_PLOT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MU_PLOT.M with the given input arguments.
%
%      MU_PLOT('Property','Value',...) creates a new MU_PLOT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MU_Plot_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MU_Plot_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help MU_Plot

% Last Modified by GUIDE v2.5 24-Sep-2013 12:14:42

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MU_Plot_OpeningFcn, ...
                   'gui_OutputFcn',  @MU_Plot_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before MU_Plot is made visible.
function MU_Plot_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to MU_Plot (see VARARGIN)

% Choose default command line output for MU_Plot
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes MU_Plot wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = MU_Plot_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in Keep_pushbutton.
function Keep_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Keep_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function FileMenu_Callback(hObject, eventdata, handles)
% hObject    handle to FileMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function OpenMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to OpenMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
file = uigetfile('*.fig');
if ~isequal(file, 0)
    open(file);
end

% --------------------------------------------------------------------
function PrintMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to PrintMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
printdlg(handles.figure1)

% --------------------------------------------------------------------
function CloseMenuItem_Callback(hObject, eventdata, handles)
% hObject    handle to CloseMenuItem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
selection = questdlg(['Close ' get(handles.figure1,'Name') '?'],...
                     ['Close ' get(handles.figure1,'Name') '...'],...
                     'Yes','No','Yes');
if strcmp(selection,'No')
    return;
end

delete(handles.figure1)

% --- Executes on button press in Copy_pushbutton.
function Copy_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Copy_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

dataObjs = findobj(handles.Plot_axes,'Type','line');
xdata = get(dataObjs, 'XData'); %data from low-level grahics objects
ydata = get(dataObjs, 'YData');

clipboard('copy', [xdata(:), ydata(:)]);

% --- Executes on button press in List_pushbutton.
function List_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to List_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

h=figure('Resize','off','position',[100 120 200 550],'Name','Voxel List','MenuBar','none');

dataObjs = findobj(handles.Plot_axes,'Type','line');
xdata = get(dataObjs, 'XData'); %data from low-level grahics objects
ydata = get(dataObjs, 'YData');

uitable('Parent',h,'Data',[xdata(:), ydata(:)],'ColumnName',{'X-Data','Y-Data'},'RowName',[],...
        'Position',[5  10  195  530]);

% --- Executes on button press in Save_pushbutton.
function Save_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Save_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

dataObjs = findobj(handles.Plot_axes,'Type','line');
xdata = get(dataObjs, 'XData'); %data from low-level grahics objects
ydata = get(dataObjs, 'YData');

try
    Flag = MU_load_matrix('data_plt', [xdata(:), ydata(:)], 0);
catch me
    error_msg{1,1}='ERROR!!! Saving plot data aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    return;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function hStructMenu = putStructMenu(hParent)
%"putStructMenu"
%   Add structure menu to CERR slice viewer, or update the existing structure menu.
%   Uses a self-updating callback.
%
%LM:  5 Apr 02, JOD. Creation.
%    26 Mar 04, JRA, self-updating callback, added struct toggling.
%
%Usage:
%   function hStructMenu = putStructMenu()
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
global stateS

%Necessary for compiled version:
if ~exist('planC')
    planC = [];
end

structuresC = {};

IMRTPdir = fileparts(which('IMRTP'));
if ~isempty(IMRTPdir)
    pos = 7;  %position of menu
else
    pos = 6;  %position of menu
end

%Create new menu if necessary.
if isfield(stateS, 'handle') && isfield(stateS.handle, 'CERRStructMenu') && ishandle(stateS.handle.CERRStructMenu)
    hStructMenu = stateS.handle.CERRStructMenu;
else
    hStructMenu = uimenu(stateS.handle.CERRSliceViewer, 'label', '&Structures', 'callback', 'putStructMenu;', 'Interruptible', 'off');
    stateS.handle.CERRStructMenu = hStructMenu;

    %First time for menu, add permenant ui submenus.
    %Call up Mesh-representation GUI
    %uimenu(hStructMenu, 'label', 'Mesh Representation', 'callback','selectStructsToMeshGUI(''init'')','interruptible','on');

    %Call up editStructures
    uimenu(hStructMenu, 'label', 'Manual contouring', 'callback',['sliceCallBack(''contourMode'')'],'interruptible','on', 'separator', 'on');

    %Call up deepLearnSegGui
    uimenu(hStructMenu, 'label', 'Deep-learning auto-contouring', 'callback',['deepLearnSegGui'],'interruptible','on');

    %One element for Segment Labeler
    uimenu(hStructMenu, 'label', 'Segment Labeler', 'callback','segmentLabelerControl(''segmentLabeler'', ''init'');','interruptible','on');

    %Call up editStructures
    %     uimenu(hStructMenu, 'label', 'Add/Edit structures', 'callback',['editStructFields'],'interruptible','on');

    %One element for structure fusion.
    uimenu(hStructMenu, 'label', 'Derive new structure', 'callback',['sliceCallBack(''structurefusion'')'],'interruptible','on');

    %One element for copying structure from one scan to another
    uimenu(hStructMenu, 'label', 'Copy from one scan to other', 'callback',['sliceCallBack(''copyStr'')'],'interruptible','on');

    %One element for structure consensus
    uimenu(hStructMenu, 'label', 'Consensus', 'callback',['sliceCallBack(''structConsensus'')'],'interruptible','on');

    %Turn all structs on
    uimenu(hStructMenu, 'label', 'View all structures', 'callback', ['sliceCallBack(''ViewAllStructures'')'], 'interruptible', 'on', 'separator', 'on');

    uimenu(hStructMenu, 'label', 'View no structures', 'callback', ['sliceCallBack(''ViewNoStructures'')'], 'interruptible', 'on');
end

if isempty(planC)
    set(hStructMenu, 'visible', 'off');
    return;
else
    set(hStructMenu, 'visible', 'on');
end
indexS = planC{end};

numStructs = length(planC{indexS.structures});

%Init structure visible list.
if ~isfield(planC{indexS.structures}, 'visible') & length(planC{indexS.structures}) ~= 0
    [planC{indexS.structures}.visible] = deal(1);
end

%Build list of structure names.
if numStructs > 0
    [structuresC{1:numStructs}] = deal(planC{indexS.structures}.structureName);
end

%Find and remove old structure listings.
kids = get(hStructMenu, 'children');
numOldMenus = length(kids);
delete(kids(1:numOldMenus-8));

[assocScansV, relStructNum] = getStructureAssociatedScan(1:numStructs, planC);
allScans = unique(assocScansV);
for i = 1:length(allScans)
    strSep = 'off';
    if i == 1
        strSep = 'on';
    end

    scanType = planC{indexS.scan}(allScans(i)).scanType;
    if isempty(scanType)
        scanType = ['Scan ', num2str(allScans(i))];
    end
    hStrSetPannel(i)= uimenu(hStructMenu, 'label', [num2str(allScans(i)), '. ', scanType], 'callback','','interruptible','on', 'Separator', strSep);

    scanIndxV = find(assocScansV == allScans(i));

    uimenu(hStrSetPannel(i), 'label', ['Associated Scan Set: ' num2str(allScans(i))],'interruptible', 'off','Enable','off');

    %Reassign color
    uimenu(hStrSetPannel(i), 'label', 'Reassign Color', 'callback',['structColorGUI(''init'',''' num2str(allScans(i)) ''')'], 'interruptible', 'on','Separator','on');

    %Turn all structs on
    uimenu(hStrSetPannel(i), 'label', ['Turn ON Set ' num2str(allScans(i))], 'callback',['sliceCallBack(''ViewAllStructures'',''' num2str(allScans(i)) ''')'], 'interruptible', 'on','Separator','on');

    uimenu(hStrSetPannel(i), 'label', ['Turn OFF Set ' num2str(allScans(i))], 'callback',['sliceCallBack(''ViewNoStructures'',''' num2str(allScans(i)) ''')'], 'interruptible', 'on');

    maxStructToShow = 25;
    numStructs = length(scanIndxV);
    for j = 1:min(numStructs,maxStructToShow)
        %Populate with new structure menu items.

        sep = 'off';
        checked = 'off';

        str = [num2str(scanIndxV(j)) '.  ' structuresC{scanIndxV(j)}];

        %Seperator if 1st structure.
        if j == 1
            sep = 'on';
        end

        %Checkmark if visible.
        if planC{indexS.structures}(scanIndxV(j)).visible
            checked = 'on';
        end

        %Create actual menu element(s).
        uimenu(hStrSetPannel(i), 'label', str, 'callback',['sliceCallBack(''toggleSingleStruct'',''', num2str(scanIndxV(j)),''')'],...
            'interruptible','on','tag','struct menu', 'Separator', sep, 'Checked', checked);
    end

    if numStructs > maxStructToShow
        uimenu(hStrSetPannel(i), 'label', 'More Structures...', 'callback',['sliceCallBack(''selectStructMore'',',num2str(allScans(i)),')'],...
            'interruptible','on','tag','struct menu', 'Separator', sep, 'Checked', checked);
    end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
clear; clc;
tic;
% [dn,fs1] = audioread('C:\Users\Xu\Desktop\��׿��\echo cancellation\subband\book\After reading Chapter 4 -- Subband Adaptive Filtering\room1_2m_pb\001_2mplayback_1646-802.wav');
% [un,fs1] = audioread('x45.wav');
[x,fs1] = audioread('C:\Users\Xu\Desktop\��׿��\echo cancellation\subband\book\After reading Chapter 4 -- Subband Adaptive Filtering\room1_2m_pb\015_2mplayback_1646-802.wav');
dn = x(:,4);
un = x(:,5);
% Arguments:
% w0               Coefficients of FIR filter at start
% mu               Step size
% N                Number of subbands
% L                length of analysis filter
% M                length of adaptive weight vector
% H                Analysis filter bank (optional), each column represents a filter
% F                Synthesis filter bank (optional), each column represents a filter
% alpha            Adjust scaling of tap weights
% delta            Small constant
M = 512;
w0 = zeros(M,1);
mu = 0.5;
N = 8;
L = 512;
alpha = 1;
delta = 0.0001;
S = PMSAFinit(w0,mu,N,L,alpha,delta)
[en,S] = PMSAFadapt(un,dn,S);
% % result = un'-1*en;
% % plot(dn'); hold on; plot(1*en);
% % legend('dn','1 en')
toc;
soundsc(en,fs1);
pause;
clear sound
% audiowrite('4.wav',en,fs1);

author: CharlesThaCat

repo: acoustic-interference-cancellation

path: github-repos/MATLAB/CharlesThaCat-acoustic-interference-cancellation/acoustic-interference-cancellation-edb394499ea6f9c96445a3e9613bd64a854c289e/Subband processing/Book/After reading Chapter 5 -- Critically Sampled and Oversampled Subband Structures/AIC.m


text:
function results = PTKComputeVolumeFromSegmentation(segmentation_mask, reporting)
    % PTKComputeVolumeFromSegmentation. Computes volume and surface area from masks
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.

    voxel_size = segmentation_mask.VoxelSize;

    voxel_volume_mm3 = voxel_size(1) * voxel_size(2) * voxel_size(3);
    lung_volume_mm3 = sum(segmentation_mask.RawImage(:) > 0)*voxel_volume_mm3;

    surface = PTKGetSurfaceFromSegmentation(segmentation_mask.RawImage);
    surface_volume_mm3 = sum(surface(:))*voxel_volume_mm3;

    results = PTKMetrics;
    results.AddMetric('VolumeCm3', lung_volume_mm3/1000, 'Volume (cm^3)');
    results.AddMetric('SurfaceVolumeCm3', surface_volume_mm3/1000, 'Volume of surface voxels (cm^3)');
end

author: tomdoel

repo: pulmonarytoolkit


text:
function r8cc_print_some ( m, n, nz_num, colptr, rowind, a, ilo, jlo, ihi, ...
  jhi, title )

%*****************************************************************************80
%
%% R8CC_PRINT_SOME prints some of a R8CC matrix.
%
%  Discussion:
%
%    The R8CC format is the double precision sparse compressed column
%    format.  Associated with this format, we have an M by N matrix
%    with NZ_NUM nonzero entries.  We construct the column pointer
%    vector COL of length N+1, such that entries of column J will be
%    stored in positions COL(J) through COL(J+1)-1.  This indexing
%    refers to both the ROW and A vectors, which store the row indices
%    and the values of the nonzero entries.  The entries of the
%    ROW vector corresponding to each column are assumed to be
%    ascending sorted.
%
%    The R8CC format is equivalent to the MATLAB "sparse" format,
%    and the Harwell Boeing "real unsymmetric assembled" (RUA) format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in A.
%
%    Input, integer COLPTR(N+1), points to the first element of each column.
%
%    Input, integer ROWIND(NZ_NUM), contains the row indices of the elements.
%
%    Input, real A(NZ_NUM), the matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  incx = 5;
%
%  Print the columns of the matrix, in strips of 5.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%4d  ', i );
%
%  Print out (up to) 5 entries in row I, that lie in the current strip.
%
      for j = j2lo : j2hi

        aij = 0.0;

        for k = colptr(j) : colptr(j+1)-1
          if ( rowind(k) == i )
            aij = a(k);
            break;
          end
        end

        fprintf ( 1, '%12g  ', aij );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [] = assertVectorsAlmostEqual(a, b, comparetype, tolerance)
  if (max(abs(reshape(a-b,[],1))) > tolerance)
    testFailed(a,b);
  elseif (min(size(a) == size(b)) < 1)
    testFailed(a,b);
  end

function [] = testFailed(a, b)
  [ST, I] = dbstack(2);
  disp(strcat("FAILED: ",  ST(1).name));
  disp(a)
  disp("--")
  disp(b)

author: ricedsp


text:
function [w, infos] = subsamp_svrg(problem, options)
% Sabsampled SVRG algorithm.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       Subsampled SVRG:
%       R. Kolte, M. Erdogdu and A. Ozgur,
%       "Accelerating SVRG via second-order information,"
%       OPT2015, 2015.
%
%
% Created by H.Kasai on Oct. 28, 2016
% Modified by H.Kasai on Mar. 25, 2018

    % set dimensions and samples
    d = problem.dim();
    n = problem.samples();

    % extract options
    if ~isfield(options, 'stepsizefun')
        options.stepsizefun = @stepsize_alg;
    else
    end

    if ~isfield(options, 'tol_optgap')
        tol_optgap = 1.0e-12;
    else
        tol_optgap = options.tol_optgap;
    end

    if ~isfield(options, 'batch_size')
        batch_size = 10;
    else
        batch_size = options.batch_size;
    end

    if batch_size > n
        batch_size = n;
    end
    num_of_bachces = floor(n / batch_size);

    if ~isfield(options, 'max_epoch')
        max_epoch = 100;
    else
        max_epoch = options.max_epoch;
    end

    if ~isfield(options, 'r')
        r = inf;
    else
        r = options.r;
    end

    if r > d + 1
        r = d - 1;
    end

    if ~isfield(options, 'w_init')
        w = randn(d,1);
    else
        w = options.w_init;
    end

    if ~isfield(options, 'f_opt')
        options.f_opt = -Inf;
    end

    if ~isfield(options, 'permute_on')
        permute_on = 1;
    else
        permute_on = options.permute_on;
    end

    if ~isfield(options, 'verbose')
        verbose = false;
    else
        verbose = options.verbose;
    end

    if ~isfield(options, 'store_w')
        options.store_w = false;
    end

    % initialize
    total_iter = 0;
    epoch = 0;
    grad_calc_count = 0;
    w = options.w_init;
    num_of_bachces = floor(n / options.batch_size);

    % store first infos
    clear infos;
    [infos, f_val, optgap] = store_infos(problem, w, options, [], epoch, grad_calc_count, 0);

    %
    sample_size = round(10*r*log(d));

    % set start time
    start_time = tic();

    % display infos
    if verbose > 0
        fprintf('Sub-sampled SVRG: Epoch = %03d, cost = %.16e, optgap = %.4e\n', epoch, f_val, optgap);
    end

    % main loop
    while (optgap > tol_optgap) && (epoch < max_epoch)

        % permute samples
        if permute_on
            perm_idx = randperm(n);
        else
            perm_idx = 1:n;
        end

       % compute full gradient
        full_grad = problem.grad(w,1:n);
        % store w
        w0 = w;
        grad_calc_count = grad_calc_count + n;

        % calculated Hessian using subsamples every outer loop
        sub_indices = datasample((1:n), sample_size);
        H = problem.hess(w, sub_indices);
        [u, sigma,~] = svd(H);
        Q = u(:,1:r);
        gamma = sigma(r+1,r+1);
        sigma = sigma(1:r,1:r);
        Sing_inv_gamma = diag(1./diag(sigma)) - diag(ones(r,1)/gamma);

        for j = 1 : num_of_bachces

            % update step-size
            step = options.stepsizefun(total_iter, options);

            % calculate variance reduced gradient
            start_index = (j-1) * batch_size + 1;
            indice_j = perm_idx(start_index:start_index+batch_size-1);
            grad = problem.grad(w, indice_j);
            grad_0 = problem.grad(w0, indice_j);
            grad_est = full_grad + grad - grad_0;

            % update w
            v = -Q*(Sing_inv_gamma)*(Q' * grad_est) - (1/gamma)*grad_est;
            w = w + step * v;

            % proximal operator
            if ismethod(problem, 'prox')
                w = problem.prox(w, step);
            end

            total_iter = total_iter + 1;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + j * batch_size + sample_size;
        epoch = epoch + 1;

        % store infos
        [infos, f_val, optgap] = store_infos(problem, w, options, infos, epoch, grad_calc_count, elapsed_time);

        % display infos
        if verbose > 0
            fprintf('Sub-sampled SVRG: Epoch = %03d, cost = %.16e, optgap = %.4e\n', epoch, f_val, optgap);
        end
    end

    if optgap < tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', tol_optgap);
    elseif epoch == max_epoch
        fprintf('Max epoch reached: max_epochr = %g\n', max_epoch);
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function test_bug1637

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY megplanar_sincos channelconnectivity ft_prepare_neighbours ft_channelselection

% this function checks whether megplanar_sincos relies on a fixed channel
% order or whether this can be totally mixed up (it should be able to deal
% with that!)
%
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1637

% load neighbours
cfg = [];
cfg.method = 'template';
cfg.template = 'CTF275_neighb';
neighbours = ft_prepare_neighbours(cfg);

cd(dccnpath('/home/common/matlab/fieldtrip/data/test'))
load bug1637_hdr.mat
load bug1637_grad.mat

d = pwd;
ftPath = fileparts(mfilename('fullpath')); % get path, strip away 'ft_defaults'
ftPath = strrep(ftPath, '\', '\\');

cd(fullfile(ftPath, '..', 'private'));
errored = false;
try
  cfg = [];
  cfg.neighbours = neighbours;
  cfg.channel    = ft_channelselection('MEG', hdr.label);
  [neighbsel] = match_str({cfg.neighbours.label}, cfg.channel);
  cfg.neighbours = cfg.neighbours(neighbsel);
  cfg.neighbsel = channelconnectivity(cfg);

  montage_other = megplanar_sincos(cfg, hdr.grad);

  montage = [];
  montage{1} = megplanar_sincos(cfg, grad);

  [sel12, sel22] = match_str(montage{1}.labelold, montage_other.labelold);
  [sel11, sel21] = match_str(montage{1}.labelnew, montage_other.labelnew);
  if ~isequal(montage{1}.tra(sel11, sel12), montage_other.tra(sel21, sel22))
    warning('tra matrix differs - but that''s because the montage_other.tra consists of nans - checking whether nonzero elements are the same')
    idx1 = montage{1}.tra(sel11, sel12)~=0;
    idx2 = montage_other.tra(sel21, sel22)~=0;
    if ~isequal(idx1, idx2)
      errored = true;
      error('tra matrix qualitatively differ!');
    else
      disp('...passed!')
    end
  elseif ~(all(cellfun(@isequal, montage{1}.labelold(sel12), montage_other.labelold(sel22))))
    errored = true;
    error('labelold mismatch');
  elseif ~(all(cellfun(@isequal, montage{1}.labelnew(sel11), montage_other.labelnew(sel21))))
    errored = true;
    error('labelnew mismatch');
  end

  for i=1:10 % 10 is an arbitrary number here, just do it for some time
    labperm = randperm(numel(grad.label));
    % permute everything
    grad.chanori    = grad.chanori(labperm, :);
    grad.chanpos    = grad.chanpos(labperm, :);
    grad.chantype   = grad.chantype(labperm, :);
    grad.chanunit   = grad.chanunit(labperm, :);
    grad.coiloriri  = grad.coilori(labperm, :);
    grad.coilpos    = grad.coilpos(labperm, :);
    grad.label      = grad.label(labperm, :);
    grad.tra        = grad.tra(labperm, :);

    montage{2} = megplanar_sincos(cfg, grad);

    % since we now reordered everything according to cfg.channel, the
    % following should not be necessary - but hey!
    [sel12, sel22] = match_str(montage{1}.labelold, montage{2}.labelold);
    [sel11, sel21] = match_str(montage{1}.labelnew, montage{2}.labelnew);
    if ~(isequal(montage{1}.tra(sel11, sel12), montage{2}.tra(sel21, sel22)))
      errored = true;
      error('tra matrix differs')
    elseif ~(all(cellfun(@isequal, montage{1}.labelold(sel12), montage{2}.labelold(sel22))))
      errored = true;
      error('labelold mismatch');
    elseif ~(all(cellfun(@isequal, montage{1}.labelnew(sel11), montage{2}.labelnew(sel21))))
      errored = true;
      error('labelnew mismatch');
    end
  end

  for i=1:10 % 10 is an arbitrary number here, just do it for some time
    labperm = randperm(numel(grad.label));
    % permute everything
    grad.chanori    = hdr.grad.chanori(labperm, :);
    grad.chanpos    = hdr.grad.chanpos(labperm, :);
    grad.chantype   = hdr.grad.chantype(labperm, :);
    grad.chanunit   = hdr.grad.chanunit(labperm, :);
    grad.coiloriri  = hdr.grad.coilori(labperm, :);
    grad.coilpos    = hdr.grad.coilpos(labperm, :);
    grad.label      = hdr.grad.label(labperm, :);
    grad.tra        = hdr.grad.tra(labperm, :);

    montage{2} = megplanar_sincos(cfg, hdr.grad);

    % since we now reordered everything according to cfg.channel, the
    % following should not be necessary - but hey!
    [sel12, sel22] = match_str(montage{1}.labelold, montage{2}.labelold);
    [sel11, sel21] = match_str(montage{1}.labelnew, montage{2}.labelnew);
    idx1 = montage{1}.tra(sel11, sel12)~=0;
    idx2 = montage_other.tra(sel21, sel22)~=0;
    if ~isequal(idx1, idx2)
      errored = true;
      error('tra matrix qualitatively differ!');
    elseif ~(all(cellfun(@isequal, montage{1}.labelold(sel12), montage_other.labelold(sel22))))
      errored = true;
      error('labelold mismatch');
    elseif ~(all(cellfun(@isequal, montage{1}.labelnew(sel11), montage_other.labelnew(sel21))))
      errored = true;
      error('labelnew mismatch');
    end
  end

end
cd(d)

if errored
  error('channelorder matters whereas it should not!');
end


text:
classdef LogarithmNode < GraphNode
    properties
        const = 1e-2;
    end

    methods
        function obj = LogarithmNode(dimOut, const)
            obj = obj@GraphNode('Logarithm',dimOut);
            if nargin>=2
                obj.const = const;
            end
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);
            input = prev_layers{1}.a;
            [D,T,N] = size(input);

            if strcmpi(class(gather(input(1))), 'single')
                obj.const = single(obj.const);
            end

            if N==1
                obj.a = log(input+obj.const);
            else
                if obj.variableLength
                    input = obj.PadShortTrajectory(input, 0);
                    if sum(input(:)+obj.const<=0)
                        fprintf('error: input to log is negative');
                    end
                    obj.a = log(input+obj.const);
                else
                    obj.a = log(input+obj.const);
                end
            end
            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,prev_layers, future_layers)
            future_grad = obj.GetFutureGrad(future_layers);
            input = prev_layers{1}.a;
            obj.grad{1} = 1./(input+obj.const).*future_grad;
        end
    end
end

author: singaxiong

repo: SignalGraph


text:
function [M] = ROIcmap(nc,opt)
% [M] = ROIcmap([nc],[opt]);
% creates a colormap with
%           no color too close to grayscale,
%           no two consecutive colors too close
%           no colors exeedingly close to another in the map
%           no colors too close to a background color (optional)
%   nc: number of colors in map. Default is 64
%   opt: optional options structure
%      .show: Figure handle in which to show the map
%             Default is 1. Set to 0 for no show.
%             In show mode, you get to pick colors with mouse
%             and read in their values. Hit enter to exit from
%             this mode.
%      .state: State of the random number generator.
%              Default is 0.
%      .write: Name of file to write colormap into
%              Default is '', no writing. Use something like
%              ROI64s0.1D.cmap, for a 64 cols, seed 0 colormap.
%      .avoid: Color to avoid getting close to.
%      .verb: verbosioty. 1 is default. 0 is for quiet
% returns
%   M: The colormap.
%
%see also readXcol, rgbdectohex, and ScaleToMap
% Ziad S. Saad SSCC/NIMH/NIH, saadz@mail.nih.gov

if (nargin == 0),
   nc = 64;
   opt.show = 1;
elseif (nargin == 1),
   opt.show = 1;
end

if (isempty(nc)), nc = 64; end

if (~isfield(opt,'show') | isempty(opt.show)), opt.show = 1; end
if (~isfield(opt,'state') | isempty(opt.state)), opt.state = 0; end
if (~isfield(opt,'write') | isempty(opt.write)), opt.write = ''; end
if (~isfield(opt,'verb') | isempty(opt.verb)), opt.verb = 1; end
if (~isfield(opt,'avoid') | isempty(opt.avoid)), opt.avoid = []; end

%initialize rng
rand('state',opt.state);

M = zeros(nc,3);
alldiff_lim = 0.5; %between 0 and 1, controls how different all colors in map are.
               %The first few colors can be quite different, high alldiff_lim
               %The difference is adjusted as more colors are demanded.
g_lim = 0.2; %limit for too gray (0-1)
d_lim = 0.40; %limit for too dim (0-3)
b_lim = 2.2; %limit for too bright  (0-3)
for (i=1:1:nc),
   M(i,:) = rand(1,3);
   cnt = 0;
   %reject if too gray or too close to previous color
   while (  toogray(M(i,:), g_lim, d_lim, b_lim) |  tooclose(M,i, 0.6, alldiff_lim) | (~isempty(opt.avoid) & (sum(abs(M(i,:)-opt.avoid)) < 0.6))),
      M(i,:) = rand(1,3);
      cnt = cnt + 1;
      if (cnt > 2000), % too tight, relax
         alldiff_lim = max([0.95.*alldiff_lim 0.02]) ;
         d_lim = max([0.95.*d_lim 0.01]);
         b_lim = min([b_lim*1.05, 8.0]);
         if (opt.verb) fprintf(1,'Reduced alldiff_lim to %g, d_lim to %g, b_lim to %g\n', alldiff_lim, d_lim, b_lim); end
         cnt = 0;
      end
   end
   if (opt.verb) fprintf(1,'Color %d OK\n', i); end
end
if (opt.verb) fprintf(1,'alldiff_lim final was %g, d_lim final was %g, b_lim final was %g\n', alldiff_lim, d_lim, b_lim); end

if (~isempty(opt.write)),
   optw.OverWrite = 'p';
   wryte3(M, opt.write, optw);
end

if (opt.show),
   ShowCmap(M, opt.show);
end

return;

function [a] = toogray(c, g_lim, d_lim, b_lim)

   a = 0;
   dc = abs(c - mean(c));
   cs = sum(c);
   if (dc(1) < g_lim & dc(2) < g_lim & dc(3) < g_lim), a = 1; return; end
   if (cs < d_lim | cs > b_lim), a = 1; return; end
   return;

function [a] = tooclose(M,i,prev_lim, alldiff_lim)

   if (i==1), a = 0; return; end

   a = 1;

   %too close to previous ?
   dc = abs(M(i,:)-M(i-1,:));
   if (sum(dc) < prev_lim), return; end

   %too close to one before?
   if (i > 2),
      for (j=1:1:i-2),
         dc = abs(M(i,:)-M(j,:));
         if (dc(1) < alldiff_lim & dc(2) < alldiff_lim & dc(3) < alldiff_lim), return; end
      end
   end
   %OK if you get here
   a = 0;
   return;


text:
function model = ivmAddPoint(model, i)

% IVMADDPOINT Add a point into the IVM representation.
% FORMAT
% DESC incorporates the ith point from the data set into the IVM
% model.
% ARG model : the model to which the point is to be added.
% ARG index : the index of the point in the training data which is
% to be added.
% RETURN model : the returned model with the point added in.
%
% SEEALSO : ivmUpdateSites, ivmUpdateM, ivmUpdateNuG, ivmSelectPoint,
% ivmRemovePoint, ivmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% IVM

index = find(model.J == i);
if isempty(index)
  error(['Point ' num2str(i) ' is not in inactive set'])
end

%/~model = ivmUpdateNuG(model, i);
%~/
model = ivmUpdateSites(model, i);
model = ivmUpdateM(model, i);

% Remove point from the non-active set and place in the active.
model.J(index) = [];
model.I = [model.I; i];

model = ivmUpdateNuG(model, model.J);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [param]=zn2pd(input)
% [param]=zn2pd(input)
% Ziegler-Nichols PD controller for processes of 2nd order.
% Controller uses reference variable (w) just in proportional component.
% This function computes parameters of the controller (r0, q0, q1, q2, p1, p2).
% Output of the controller is calculated follows:
%
%                  r0                            q0 + q1*z^-1 + q2*z^-2
% U(z^-1) = ----------------------- * W(z^-1) - ------------------------ * Y(z^-1)
%            1 + p1*z^-1 + p2*z^-2                1 + p1*z^-1 + p2*z^-2
%
% where p1=0, p2=0, q2=0
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2
% Gs(z^-1) = -----------------------
%             1 + a1*z^-1 + a2*z^-2
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... sample time T0
% Output: param ... controller parameters
%                   param(1) ... r0
%                   param(2) ... q0
%                   param(3) ... q1
%                   param(4) ... q2 (0)
%                   param(5) ... p1 (0)
%                   param(6) ... p2 (0)

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
T0 = input(5);

% compute ultimate gain and frequency
[Kpu, Tu] =  ultim([b1 b2],[a1 a2],T0);

Kp = 0.4*Kpu;
Td = Tu/20;

r0 = Kp;
q0 = Kp*(1+Td/T0);
q1 = -Kp*(Td/T0);
q2 = 0;
p1 = 0;
p2 = 0;

param=[r0; q0; q1; q2; p1; p2];


text:
function sequences = sequence_load(directory, varargin)
% sequence_load Load a set of sequences
%
% Create a cell array of new sequence structures for sequences, specified in a listing file.
%
% Input:
% - directory: Path to the directory with sequences.
% - varargin[List]: Name of the file that lists all the sequences. By default `list.txt` is used.
%
% Output:
% - sequences: A cell array of new sequence structures.

list = 'list.txt';

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'list'
            list = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i},'!']) ;
    end
end

list_file = fullfile(directory, list);

sequences = cell(0);

mkpath(directory);

bundle_url = get_global_variable('bundle');

if ~exist(list_file, 'file') && ~isempty(bundle_url)

    if strsuffix(bundle_url, '.zip')

        print_text('Downloading sequence bundle from "%s". This may take a while ...', bundle_url);
        bundle = [tempname, '.zip'];
        try
            urlwrite(bundle_url, bundle);
            unzip(bundle, directory);
		    delete(bundle);
            list_file = fullfile(directory, 'list.txt');
        catch
            print_text('Unable to retrieve sequence bundle from the server. This is either a connection problem or the server is temporary offline.');
            print_text('Please try to download the bundle manually from %s and uncompress it to %s', bundle_url, directory);
            return;
        end;

    elseif strsuffix(bundle_url, '.json')
        listing = tempname;
        urlwrite(bundle_url, listing);
        meta = json_decode(fileread(listing));

        print_text('Downloading sequence dataset "%s" with %d sequences.', meta.name, numel(meta.sequences));

        slashes = strfind(bundle_url, '/');
        base_url = bundle_url(1:slashes(end));

        print_indent(1);

        for i = 1:numel(meta.sequences)
            sequence = meta.sequences{i};
            print_text('Downloading sequence "%s" ...', sequence.name);
            sequence_directory = fullfile(directory, sequence.name);
            %if exist(sequence_directory, 'dir')
            %    continue;
            %end;
            mkpath(sequence_directory);
            data = struct('name', sequence.name, 'fps', sequence.fps, 'format', 'default');
            data.channels = struct();

            if strncmp(sequence.annotations.url, 'http://', 7) || strncmp(sequence.annotations.url, 'https://', 8)
                annotations_url = sequence.annotations.url;
            else
                annotations_url = [base_url, sequence.annotations.url];
            end

            try
                bundle = [tempname, '.zip'];
                urlwrite(annotations_url, bundle);
                unzip(bundle, sequence_directory);
                delete(bundle);
            catch
                print_text('Unable to retrieve sequence bundle from the server. This is either a connection problem or the server is temporary offline.');
                return;
            end;

            channels = fieldnames(sequence.channels);

            for c = 1:numel(channels)
                channel = sequence.channels.(channels{c});
                channel_directory = fullfile(sequence_directory, channels{c});
                mkpath(channel_directory);
                if strncmp(channel.url, 'http://', 7) || strncmp(channel.url, 'https://', 8)
                    channel_url = channel.url;
                else
                    channel_url = [base_url, channel.url];
                end
                try
                    bundle = [tempname, '.zip'];
                    urlwrite(channel_url, bundle);
                    unzip(bundle, channel_directory);
                    delete(bundle);
                catch
                    print_text('Unable to retrieve sequence bundle from the server. This is either a connection problem or the server is temporary offline.');
                    return;
                end;

                if isfield(channel, 'pattern')
                    data.channels.(channels{c}) = [channels{c}, filesep, channel.pattern];
                else
                    data.channels.(channels{c}) = [channels{c}, filesep];
                end
            end

            writestruct(fullfile(sequence_directory, 'sequence'), data);

        end

        list_file = fullfile(directory, 'list.txt');
        % Save list file
        fid = fopen(list_file, 'w');
        for i = 1:numel(meta.sequences)
            fprintf(fid, '%s\n', meta.sequences{i}.name);
        end
        fclose(fid);

        print_indent(-1);

    else
        error('Unknown dataset type!');
    end;

end;

fid = fopen(list_file, 'r');

while true
    sequence_name = fgetl(fid);
    if sequence_name == -1
        break;
    end

    if exist(fullfile(directory, sequence_name, 'sequence'), 'file')
        sequence_path = fullfile(directory, sequence_name, 'sequence');
    elseif exist(fullfile(directory, sequence_name), 'dir')
        sequence_path = fullfile(directory, sequence_name);
    else
        continue;
    end;

    print_debug('Loading sequence %s', sequence_name);

    sequences{end+1} = sequence_create(sequence_path); %#ok<AGROW>

end;

fclose(fid);

author: votchallenge

repo: toolkit-legacy


text:
function [aVal,aJacob,aHess,papt]=aPolarCoordTurn2DOmega(x,tauInv,tauAccelInv)
%%APOLARCOORDTURN2DOMEGA The continuous-time drift function for a 2D
%           coordinated turn model with the velocity expressed in terms of
%           a heading angle and a speed rather than in terms of a Cartesian
%           velocity vector. Additionally a turn rate is part of the state
%           (unlike in aPolarLin2D) and a linear acceleration term can be
%           given, which acts in the direction of motion. The turn rate and
%           linear acceleration can optionally have time constants
%           associated with them, like in the Singer model, modelling a
%           tendancy to eventually return to non-accelerating, straight-
%           line motion.
%
%INPUTS: xState The 5X1 or 6X1 target state for 2D motion. If there is no
%               linear acceleration (acceleration along the direction of
%               motion), then the state is [xPos;yPos;theta;v;omega],
%               where (xPos,yPos) are the Cartesian position, theta is the
%               angle of the velocity vector in radians counterclockwise
%               from the x axis, v is the speed, and omega is the turn
%               rate, typically in radians per second. If a linear
%               acceleration component is provided, then the state is
%               [xPos;yPos;theta;v;omega,al].
%        tauInv The inverse of the correlation time constant tau for the
%               turn rate in seconds. tauInv must be positive. The default
%               if omitted or an empty matrix is passed is 0.
%   tauAccelInv The inverse of the correlation time constant for the linear
%               acceleration in seconds. This parameter is not used if
%               there is no linear acceleration. the default if omitted or
%               an empty matrix is passed is 0.
%
%OUTPUTS: aVal The 5X1 (or 6X1 with linear acceleration) time-derivative of
%              the state.
%       aJacob The 5X5 (or 6X6) matrix of partial derivatives of aVals
%              such that aJacob(:,k) is the partial derivative of
%              aVals(:,k) with respect to x(k).
%        aHess The 5X5X5 (or 6X6X6) matrix of second derivatives of aVals
%              such that aHess(:,k1,k2) is the second partial derivative of
%              aVals with respect to x(k1) and x(k2).
%         papt The 5X1 or 6X1 derivative with resect to time of aVals.
%              This is all zeros, because the model is time invariant.
%
%The basic 2D coordinated turn model in Cartesian coordinates is described
%in Section VA of [1]. When the turn rate is something that must be
%estimated, it is assumed that the continuous-time turn rate model is
%omegaDot=-(1/tauTurn)*Omega+noise
%Note that the ordering of the state elements assumed by this function
%differs from the ordering of the state elements assumed in [1].
%
%The 2D coordinates turn model in Cartesian coordinates is also described
%in Chapter 4.2.3 of [2].
%
%The corresponding diffusion matrix is given by the function DCoordTurn2D.
%The corresponding discrete-time functions are FCoordTurn2D and
%QCoordTurn. However, note that the discrete-time functions are
%direct-discrete models and not discretizations of the continuous-time
%models as the propagated PDF does not remain Gaussian over time.
%
%REFERENCES:
%[1] X. R. Li and V. P. Jilkov, "Survey of maneuvering target tracking.
%    Part I: Dynamic models," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 39, no. 4, pp. 1333-1364, Oct. 2003.
%[2] S. Blackman and R. Popoli, Design and Analysis of Modern Tracking
%    Systems. Norwood, MA: Artech House, 1999.
%
%January 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(tauInv))
    tauInv=0;
end

if(nargin<3||isempty(tauAccelInv))
    tauAccelInv=0;
end

%Extract the heading and speed (polar velocity components).
theta=x(3);
v=x(4);
omega=x(5);%The turn rate

xDim=size(x,1);

sinTheta=sin(theta);
cosTheta=cos(theta);
switch(xDim)
    case 5%There is no linear acceleration
        %From Equations 61 and 67 in Li's paper.
        aVal=[v*cosTheta;%Position dervative
              v*sinTheta;%Position derivative
              omega;%Heading derivative
              0;%Speed derivative
              -tauInv*omega];%Turn rate derivative

        if(nargout>1)
            dXdTheta=-v*sinTheta;
            dYdTheta=v*cosTheta;
            dXdv=cosTheta;
            dYdv=sinTheta;

            aJacob=[0,0,dXdTheta,dXdv,0;
                    0,0,dYdTheta,dYdv,0;
                    0,0,0,          0,1;
                    0,0,0,          0,0;
                    0,0,0,          0,-tauInv];

            if(nargout>2)
                aHess=zeros(5,5,5);

                dXdThetadTheta=-v*cosTheta;
                dXdvdTheta=-sinTheta;

                dXdThetadv=dXdvdTheta;
                dXdvdv=0;

                %%%
                dYdThetadTheta=-v*sinTheta;
                dYdvdTheta=cosTheta;

                dYdThetadv=dYdvdTheta;
                dYdvdv=0;

                aHess(:,:,3)=[0,0,dXdThetadTheta,dXdvdTheta,0;
                              0,0,dYdThetadTheta,dYdvdTheta,0;
                              zeros(3,5)];
                aHess(:,:,4)=[0,0,dXdThetadv,dXdvdv,0;
                              0,0,dYdThetadv,dYdvdv,0;
                              zeros(3,5)];

                if(nargout>3)
                    papt=zeros(5,1);
                end
            end
        end
    case 6%There is a linear acceleration component.
        al=x(6);%The linear acceleration (vDot)

        aVal=[v*cosTheta;%Position dervative
              v*sinTheta;%Position derivative
              omega;%Heading derivative
              al;%Speed derivative
              -tauInv*omega;%Turn rate derivative
              -tauAccelInv*al];%Linear acceleration derivative

        if(nargout>1)
            dXdTheta=-v*sinTheta;
            dYdTheta=v*cosTheta;
            dXdv=cosTheta;
            dYdv=sinTheta;

            aJacob=[0,0,dXdTheta,dXdv,        0,0;
                    0,0,dYdTheta,dYdv,        0,0;
                    0,0,       0,   0,        1,0;
                    0,0,       0,   0,        0,1;
                    0,0,       0,   0,  -tauInv,0;
                    0,0,       0,   0,        0,-tauAccelInv];

            if(nargout>2)
                aHess=zeros(6,6,6);

                dXdThetadTheta=-v*cosTheta;
                dXdvdTheta=-sinTheta;

                dXdThetadv=dXdvdTheta;
                dXdvdv=0;

                %%%
                dYdThetadTheta=-v*sinTheta;
                dYdvdTheta=cosTheta;

                dYdThetadv=dYdvdTheta;
                dYdvdv=0;

                aHess(:,:,3)=[0,0,dXdThetadTheta,dXdvdTheta,0,0;
                              0,0,dYdThetadTheta,dYdvdTheta,0,0;
                              zeros(4,6)];

                aHess(:,:,4)=[0,0,dXdThetadv,dXdvdv,0,0;
                              0,0,dYdThetadv,dYdvdv,0,0;
                              zeros(4,6)];

                if(nargout>3)
                    papt=zeros(6,1);
                end
            end
        end
    otherwise
        error('The dimensionality of the state is neither 5 nor 6.');
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function options = ppcaOptions(mappingType, latentPoints, mappingOptions)

% PPCAOPTIONS Options for probabilistic PCA.
% FORMAT
% DESC returns the default options for probabilistic PCA.
% RETURN options : default options structure for probabilistic PCA.
%
%
% SEEALSO : ppcaCreate
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

options = [];


text:
function checkEnvironment(param,num_arg)
% checkEnvironment checks if the required tools are present in the host system.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Check for input arguments
if(num_arg==0)
    error('Type help startSimulation')
end

if(num_arg == 3 && param{1}.OperatingMode == 1)
    error('For constant input currents, please provide the current density as a parameter to startSimulation')
end

% Check if SUNDIALS is installed correctly
disp('Checking for SUNDIALS availability...');
result = exist('IDAFree');
if(result~=2)
    error('The SUNDIALS software does not seem to be present in the search path or installed on your system.\nPlease go to:\nhttps://computation.llnl.gov/casc/sundials/main.html\nand follow the instruction for installing SUNDIALS in MATLAB.');
else
    disp('SUNDIALS Package FOUND!');
end

%Check for CasADi
disp('Checking for CasADi availability...');

try
    import casadi.*
catch e
    error('The CasADi package does not seem to be present in the search path or installed on your system.\nPlease go to:\nhttps://github.com/casadi/casadi/wiki/InstallationInstructions\nand follow the instruction for installing CasADi in MATLAB.');
end

disp('CasADi Package FOUND!');

% For each parameter structure, check if all the fields are set.
for i=1:length(param)
    %Check if the parameters have been set correctly
    [result, missing] = checkBatteryParameters(param{i});
    if(result==0)
        disp('Warning, there are missing parameters in the param array.')
        disp('Please find below, the list of missing parameters:')
        for jj=1:length(missing)
            disp(missing{jj});
        end
        error('Please fix the problem and restart the script')
    end
end

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
%
% Usage:  V=mexSparseProject(U,param);
%
% Name: mexSparseProject
%
% Description: mexSparseProject solves various optimization
%     problems, including projections on a few convex sets.
%     It aims at addressing the following problems
%     for all columns u of U in parallel
%       1) when param.mode=1 (projection on the l1-ball)
%           min_v ||u-v||_2^2  s.t.  ||v||_1 <= thrs
%       2) when param.mode=2
%           min_v ||u-v||_2^2  s.t. ||v||_2^2 + lamuda1||v||_1 <= thrs
%       3) when param.mode=3
%           min_v ||u-v||_2^2  s.t  ||v||_1 + 0.5lamuda1||v||_2^2 <= thrs
%       4) when param.mode=4
%           min_v 0.5||u-v||_2^2 + lamuda1||v||_1  s.t  ||v||_2^2 <= thrs
%       5) when param.mode=5
%           min_v 0.5||u-v||_2^2 + lamuda1||v||_1 +lamuda2 FL(v) + ...
%                                                   0.5lamuda_3 ||v||_2^2
%          where FL denotes a "fused lasso" regularization term.
%       6) when param.mode=6
%          min_v ||u-v||_2^2 s.t lamuda1||v||_1 +lamuda2 FL(v) + ...
%                                             0.5lamuda3||v||_2^2 <= thrs
%
%        When param.pos=true and param.mode <= 4,
%        it solves the previous problems with positivity constraints
%
% Inputs: U:  double m x n matrix   (input signals)
%               m is the signal size
%               n is the number of signals to project
%         param: struct
%           param.thrs (parameter)
%           param.lambda1 (parameter)
%           param.lambda2 (parameter)
%           param.lambda3 (parameter)
%           param.mode (see above)
%           param.pos (optional, false by default)
%           param.numThreads (optional, number of threads for exploiting
%             multi-core / multi-cpus. By default, it takes the value -1,
%             which automatically selects all the available CPUs/cores).
%
% Output: V: double m x n matrix (output matrix)
%
% Note: this function admits a few experimental usages, which have not
%     been extensively tested:
%         - single precision setting
%
% Author: Julien Mairal, 2009

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
classdef PTKDependency < handle
    % PTKDependency. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     PTKDependency sores a dependency object, uniquely tagging a particular
    %     result generated from running a plugin. This dependency is stored
    %     alongside the plugin result and stored in the
    %     dependency list of every plugin which uses this result. This allows us
    %     to determine that a result is still valid, because its dependencies
    %     still match the dependency objects in the cache.
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        PluginName
        Context
        Uid
        DatasetUid
        Attributes
    end

    methods
        function obj = PTKDependency(name, context, uid, dataset_uid, attributes)
            obj.PluginName = name;
            obj.Context = context;
            obj.Uid = uid;
            obj.DatasetUid = dataset_uid;
            obj.Attributes = attributes;
        end
    end
end


text:
function varargout = isoutside(varargin)
%ISOUTSIDE Define avoidance constraint on polytope
%
% F = ISOUTSIDE(P)
%
% Input
%  P : Linear inequalities defining polytope
%
% Since the modelling is based on big-M formulations, all
% involved variables should be explicitly bounded.

% Here is the real overloaded ismember
switch class(varargin{1})
    case {'lmi','constraint'}

        P = varargin{1};
        if ~(all(is(P,'linear')) & all(is(P,'elementwise')))
            error('ISOUTSIDE only applicable to linear inequalities')
        end
        varargout{1} = (yalmip('define',mfilename,varargin{:}) == 1);
        varargout{1} = setupMeta(lmi([]), mfilename,varargin{:});

    case 'char'
        varargout{1} = isoutside_internal(varargin{3});
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function geometry_test20321 ( )

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

  tetra = [ ...
     0.577350269189626,  0.0, 0.0; ...
    -0.288675134594813,  0.5, 0.0; ...
    -0.288675134594813, -0.5, 0.0; ...
     0.0,                0.0, 0.816496580927726 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST20321\n' );
  fprintf ( 1, '  For a tetrahedron in 3D,\n' );
  fprintf ( 1, '  TETRAHEDRON_EDGE_LENGTH_3D computes the edge lengths;\n' );

  r8mat_transpose_print ( dim_num, 4, tetra, '  Tetrahedron vertices:' );

  edge_length = tetrahedron_edge_length_3d ( tetra );

  r8vec_print ( 6, edge_length, '  Edge lengths:' );

  return
end


text:
function [engine, loglik] = initialize_engine(engine)
%initialize
bnet = bnet_from_engine(engine);
ns = bnet.node_sizes(:);
N = length(bnet.dag);

pot_type = 'scg'
check_for_cd_arcs([], bnet.cnodes, bnet.dag);

% Evaluate CPDs with evidence, and convert to potentials
pot = cell(1, N);
C = length(engine.cliques);
inited = zeros(1, C);
clpot = cell(1, C);
evidence = cell(1, N);
for n=1:N
  fam = family(bnet.dag, n);
  e = bnet.equiv_class(n);
  pot{n} = CPD_to_scgpot(bnet.CPD{e}, fam, ns, bnet.cnodes, evidence);
  cindex = engine.clq_ass_to_node(n);
  if inited(cindex)
      %clpot{cindex} = direct_combine_pots(clpot{cindex}, pot{n});
      clpot{cindex} = direct_combine_pots(pot{n}, clpot{cindex});
  else
      clpot{cindex} = pot{n};
      inited(cindex) = 1;
  end
end

for i=1:C
    if inited(i) == 0
        clpot{i} = scgpot([], [], [], []);
    end
end

seppot = cell(C, C);
% separators are is not need to initialize

% collect to root (node to parents)
for n=engine.postorder(1:end-1)
  for p=parents(engine.jtree, n)
      [margpot, comppot] = complement_pot(clpot{n}, engine.separator{p,n});
      margpot = marginalize_pot(clpot{n}, engine.separator{p,n});
      clpot{n} = comppot;
      %seppot{p, n} = margpot;
      clpot{p} = combine_pots(clpot{p}, margpot);
      %clpot{p} = combine_pots(margpot, clpot{p});
  end
end

temppot = clpot;
%temppot = clpot{engine.root};
for n=engine.preorder
  for c=children(engine.jtree, n)
    seppot{n,c} = marginalize_pot(temppot{n}, engine.separator{n,c});
    %seppot{n,c} = marginalize_pot(clpot{n}, engine.separator{n,c});
    %clpot{c} = direct_combine_pots(clpot{c}, seppot{n,c});
    temppot{c} = direct_combine_pots(temppot{c}, seppot{n,c});
  end
end

engine.clpot = clpot;
engine.seppot = seppot;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [affineMat,pixDim,planC] = getScanAffineMat(planC,scanNum)

if ischar(planC)
    planC = loadPlanC(planC);
end

indexS = planC{end};

if ~exist('scanNum','var') || isempty(scanNum)
    scanNum = 1;
end

iHat = [1; 0; 0; 0];
jHat = [0; 1; 0; 0];
kHat = [0; 0; 1; 0];

N = numel(planC{indexS.scan}(scanNum).scanInfo);

try
    iop = planC{indexS.scan}(scanNum).scanInfo(1).imageOrientationPatient;
    if isempty(iop)
        disp('defaulting to HFS orientation');
        iop = [1 0 0 0 1 0]';
    end
    ipp = (planC{indexS.scan}(scanNum).scanInfo(end).imagePositionPatient - planC{indexS.scan}(scanNum).scanInfo(1).imagePositionPatient)/(N-1);
    if isempty(ipp)
            ipp = [0 0 -planC{indexS.scan}(scanNum).scanInfo(1).sliceThickness*10]';
    end
catch err
    disp(err);
    disp('defaulting to HFS orientation');
    iop = [1 0 0 0 1 0]';
    ipp = [0 0 -planC{indexS.scan}(scanNum).scanInfo(1).sliceThickness*10]';
end
iop = iop(:);
ipp = ipp(:);

pixsp = 10*[planC{indexS.scan}(scanNum).scanInfo(1).grid1Units planC{indexS.scan}(scanNum).scanInfo(1).grid2Units];
sliceThickness = planC{indexS.scan}(scanNum).scanInfo(1).sliceThickness * 10;
% voxel_size = [pixsp sliceThickness];
planeMat = [pixsp(2)*iop(4:end) pixsp(1)*iop(1:3)]; %.*[-1 -1;-1 -1; 1 1];

%planeMat = [pixsp(2)*iop(1:3) pixsp(1)*iop(4:end)].*[-1 -1;-1 -1; 1 1];
% [~,orientationStr,~] = returnViewerAxisLabels(planC,scanNum);
% if strcmpi('FFP',orientationStr) || strcmpi('FFS',orientationStr)
%     originLPS = planC{indexS.scan}(scanNum).scanInfo(1).imagePositionPatient;
% else
    originLPS = planC{indexS.scan}(scanNum).scanInfo(end).imagePositionPatient;
    if isempty(originLPS)
        originLPS = [0; 0; 0];
    end
% end
originLPS = originLPS(:);
rawAffineMat = [planeMat ipp originLPS; 0 0 0 1];

rawPixDim = [pixsp(2) pixsp(1) sliceThickness];
[~,xCol] = max(abs(rawAffineMat * iHat)); %[1; 0; 0; 0]))
[~,yCol] = max(abs(rawAffineMat * jHat)); %[0; 1; 0; 0]))
[~,zCol] = max(abs(rawAffineMat * kHat)); %[0; 0; 1; 0]))

pixDim = [rawPixDim(xCol) rawPixDim(yCol) rawPixDim(zCol)];

% LIA -> RAS
affIdent = eye(4);
affIdent(xCol,xCol) = -1;
affIdent(yCol,yCol) = -1;
% ##affIdent

affineMat = rawAffineMat * affIdent;

zCorrect = [9 10 3 7];
affineMat(zCorrect) = - affineMat(zCorrect);


text:
function [K, sK] = indexardKernCompute(kern, x, x2)

% INDEXARDKERNCOMPUTE Compute the INDEXARD kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the index ard based covariance function
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the index ard based covariance function
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : indexardKernParamInit, kernCompute, kernCreate, indexardKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN
  if size(x, 2)>1
    error('Index kernel requires 1-dimensional input.')
  end

  if nargin<3
    x2 = x;
  end
  K = zeros(size(x, 1), size(x2, 1));
  for i = 1:size(x, 1)
    for j = 1:size(x2, 1)
      if round(x(i)) == round(x2(j))
        ind = find(round(x(i))==kern.indices);
        if isempty(ind)
          error('Unknown index in input');
        end
        K(i, j) = kern.indexScales(ind);
      end
    end
  end
end


text:
function IntegrateEuler(j)
% following Featherstone "Robot Dynamics Algorithm" p.103
% '04 May 8 s.k AIST
global uLINK Dtime

if j == 0 return; end
if j == 1
    [uLINK(j).p, uLINK(j).R] = SE3exp(j, Dtime);
    uLINK(j).vo = uLINK(j).vo + Dtime * uLINK(j).dvo;
    uLINK(j).w  = uLINK(j).w  + Dtime * uLINK(j).dw;
    %IntegrateSE3(j, Dtime);
else
    uLINK(j).q  = uLINK(j).q  + Dtime * uLINK(j).dq;
    uLINK(j).dq = uLINK(j).dq + Dtime * uLINK(j).ddq;
end

IntegrateEuler(uLINK(j).sister);
IntegrateEuler(uLINK(j).child);

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
%                                                     [lab,q,f] = eigK(x,K)
%
% EIGK  Computes the spectral values ("eigenvalues") or even the complete
%       spectral decomposition  of a vector x with respect to a self-dual
%       homogeneous cone K.
%
% > LAB = EIGK(x,K) This yield the spectral values of x with respect to
%       the self-dual homogeneous cone that you describe in the structure
%       K. Up to 3 fields can be used, called K.l, K.q and K.s, for
%       Linear, Quadratic and Semi-definite. Type `help sedumi' for more
%       details on this structure.
%
%       The length of the vector LAB is the order of the cone. Remark that
%       x in K if and only if LAB>=0, and x in int(K) if and only if LAB>0.
%
% > [LAB,Q,F] = EIGK(x,K) Also produces the eigenvectors for the symmetric and
%       Hermitian parts of x (corresponding to K.s), and the Lorentz frame F
%       for the Lorentz blocks in x (corresponding to K.q). This extended
%       version of EIGK is intended for INTERNAL USE BY SEDUMI.
%
% See also sedumi, mat, vec, eyeK.

function [lab,q,f] = eigK(x,K) %#ok

% THE M-FILE VERSION OF THIS FUNCTION IS HERE ONLY AS ILLUSTRATION.
% SEE THE C-SOURCE FOR THE MEX-VERSION.

%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

lab = zeros(K.l + 2*length(K.q) + sum(K.s),1); %#ok
% ----------------------------------------
% LP: lab = x
% ----------------------------------------
lab(1:K.l) = x(1:K.l);
firstk = K.l+1;
nextlab = K.l+1;
% ----------------------------------------
% LORENTZ: (lab, f) = qeig(x)
% ----------------------------------------
if nargout < 3
    for k=1:length(K.q)
        nk = K.q(k); lastk = firstk + nk - 1;
        lab(nextlab:nextlab+1) = qeig(x(firstk:lastk));
        firstk = lastk + 1; nextlab = nextlab + 2;
    end
else
    nextf = 1;
    for k=1:length(K.q)
        nk = K.q(k); lastk = firstk + nk - 1;
        [labk, fk] = qeig(x(firstk:lastk));
        lab(nextlab:nextlab+1) = labk;
        f(nextf:nextf+nk-2) = fk;
        firstk = lastk + 1; nextlab = nextlab + 2; nextf = nextf + nk-1;
    end
end
% ----------------------------------------
% SDP: [q,lab] = eig(x)
% ----------------------------------------
offq = firstk - 1;
q = zeros(length(x)-offq,1);
for k = 1:K.rsdpN
    nk = K.s(k); lastk = firstk + nk*nk - 1;
    Xk = mat(x(firstk:lastk),nk);
    Xk = (Xk + Xk') / 2;
    [Qk, Labk] = eig(Xk);
    lab(nextlab:nextlab+nk-1) = diag(Labk);
    q(firstk-offq:lastk-offq) = Qk;
    firstk = lastk + 1; nextlab = nextlab + nk;
end
for k = (1+K.rsdpN):length(K.s)
    nk = K.s(k); ifirstk = firstk + nk*nk; lastk = ifirstk + nk*nk - 1;
    Xk = mat(x(firstk:ifirstk-1)+sqrt(-1)*x(ifirstk:lastk),nk);
    [Qk, Labk] = eig(Xk);
    lab(nextlab:nextlab+nk-1) = real(diag(Labk));
    q(firstk-offq:ifirstk-1-offq) = real(Qk);
    q(ifirstk-offq:lastk-offq) = imag(Qk);
    firstk = lastk + 1; nextlab = nextlab + nk;
end

repo: REST


text:
function newImg = dtiResliceVolume(img, xform, boundingBox, mmPerVoxOut)
%
% newImg = dtiResliceVolume(img, xform, boundingBox, mmPerVoxOut)
%
% HISTORY:
%   2004.01.13 RFD (bob@white.stanford.edu) wrote it.

if(~exist('mmPerVoxOut') | isempty(mmPerVoxOut))
    mmPerVoxOut = [1 1 1];
end

if(~exist('boundingBox') | isempty(boundingBox))
    % create a default bounding box for the resliced data.
    % The following is the spm default bounding box. Note that the box is
    % defined in Talairach space (units = mm).
    boundingBox = [-78 -120 -60;
                    78  80   85];
end

% myCinterp3 likes [rows,columns,slices], so we need a permute here
img = double(permute(img,[2,1,3]));

x = (boundingBox(1,1):mmPerVoxOut(1):boundingBox(2,1));
y = (boundingBox(1,2):mmPerVoxOut(2):boundingBox(2,2));
z = (boundingBox(1,3):mmPerVoxOut(3):boundingBox(2,3));
[X,Y,Z] = meshgrid(x, y, z);
newSize = size(X);
clear x y z;
talCoords = [X(:) Y(:) Z(:)];
newSize = size(X);
clear X Y Z;

imgCoords = mrAnatXformCoords(xform, talCoords);
newImg = myCinterp3(img, [size(img,1) size(img,2)], size(img,3), imgCoords, 0.0);
newImg = squeeze(reshape(newImg,newSize));
return;

author: vistalab

repo: vistasoft