text:
function output = callqsopt(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

n = length(c);
% Bounded variables converted to constraints
if ~isempty(ub)
    % TODO : MAKE SURE SIZES ARE CORRECT
    lb(lb==-inf) = -1e30;
    ub(ub==inf) = 1e30;
else
    lb=repmat(-1e30,n,1);    % just for sure
    ub=repmat(1e30,n,1);
end

if options.showprogress;showprogress('Calling QSOPT',options.showprogress);end

options.qsopt.verbose = options.verbose;

% Call mex-interface
if options.savedebug
    save qsoptdebug
end
solvertime = tic;
[x,lambda,STATUS] = qsopt(c,-F_struc(1+K.f:end,2:end),F_struc(1+K.f:end,1),-F_struc(1:K.f,2:end),F_struc(1:K.f,1),lb,ub,options.qsopt);
solvertime = toc(solvertime);
problem = 0;

if options.saveduals
    D_struc = -lambda;
else
    D_struc = [];
end

% Check, currently not exhaustive...
switch STATUS
    case 1
        problem = 0;
    case 2
        problem = 1;
    case 3
        problem = 2;
    case 4
        problem = 3;
    case {5,6}
        problem = 4;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,'QSOPT');

% Save all data sent to solver?
if options.savesolverinput
	solverinput = [];
else
	solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
	solveroutput = [];
else
	solveroutput = [];
end

% Standard interface
Primal      = solution.x(:);
Dual        = solution.D_struc;
problem     = solution.problem;
infostr     = yalmiperror(solution.problem,interfacedata.solver.tag);
if ~options.savesolverinput
    solverinput = [];
else
    solverinput = model;
end
if ~options.savesolveroutput
    solveroutput = [];
else
    solveroutput = solveroutput;
end
% Standard interface
Primal      = solution.x(:);
Dual        = solution.D_struc;
problem     = solution.problem;
infostr     = yalmiperror(solution.problem,interfacedata.solver.tag);
if ~options.savesolverinput
    solverinput = [];
else
    solverinput = model;
end
if ~options.savesolveroutput
    solveroutput = [];
else
    solveroutput = solveroutput;
end
% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% h = POINTLAB(P) point cloud visualizer.
%
% P: a struct where P.points contains 3D point coordinates
%
function h_out = pointlab( P )

% setup
h = gca;
hold on;

% modify the position of the axis in the current figure
hg = hgtransform;
setappdata( h, 'hg', hg );
set( h,'Position',[0 0.05 1 1-0.05 ] );
set( gcf, 'color', 'white');

% do not give output if not requested
if nargout == 0
    h_out = [];
else
    h_out = h;
end

%%%%%%%%%%%%%%%%%%%%%%%% menu setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%
menu_view = uimenu('Label','View');
    uimenu( menu_view,'Label','View Samples','Callback',@view_samples_callback);
    uimenu( menu_view,'Label','View Normals','Callback',@view_normals_callback);

menu_tools = uimenu('Label','Tools');
    uimenu( menu_tools,'Label','Pick sample', 'Callback', @point_pick_callback);
    uimenu( menu_tools,'Label','Snapshot', 'Callback', @snapshot_callback);
    uimenu( menu_tools,'Label','Record Fly Around', 'Callback', @record_callback);
% save status variables
h_samples = draw_samples();  setappdata(h,'h_samples',h_samples);

% set drawing limits and vis setup
xlim([-1,1]);
ylim([-1,1]);
axis equal
axis off

function view_samples_callback( IGNORE1, IGNORE2 ) %#ok<INUSD>
    if strcmp( get( h_samples, 'Visible' ), 'on')
        set( h_samples, 'Visible', 'off' )
    else
        set( h_samples, 'Visible', 'on' )
    end
end
% this is slighgly different...
% normals are not drawn by default (they would take too long)
% I just draw them if specifically requested.
function view_normals_callback( IGNORE1, IGNORE2 )  %#ok<INUSD>
    for i=1:size(P.points,1)
        myline( P.points(i,:), P.normals(i,:), .1, 'parent', hg );
    end
end

function point_pick_callback( IGNORE1, IGNORE2 ) %#ok<INUSD>
    % convert selection point in homogeneous and apply transform
	disp( find_closest( h, P.points ) );
end

function snapshot_callback( IGNORE1, IGNORE2 ) %#ok<INUSD>
    outfilename = sprintf('results/snapshot_%s_%s', datestr(now,'ddmmyy'), datestr(now,'HHMMss') );
    print('-dpng', outfilename);
end
function record_callback( IGNORE1, IGNORE2 ) %#ok<INUSD>
    disp('- modify the basic view direction');
    filename = input('- insert the name of the output: ', 's');

    prev_matrix = get( hg, 'matrix' );
    for alpha = linspace(0,2*pi,100);
        rot_matrix = makehgtform('yrotate',alpha);
        set(hg,'Matrix',rot_matrix*prev_matrix);
        drawnow;
        gif_add_frame(gca,filename,2);
    end
end

function h_samples = draw_samples()
    h_samples = plot3( P.points(:,1), P.points(:,2), P.points(:,3), '.b', 'parent', hg, 'MarkerSize', .25 );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ARCBALL STUFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% register callbacks ad
set( gcf, 'WindowButtonMotionFcn', @motion_callback );
set( gcf, 'WindowButtonDownFcn', @buttondown_callback );
set( gcf, 'WindowButtonUpFcn', @buttonup_callback );
mousestatus = 'buttonup';
START = [0,0,0];
M_previous = get( hg, 'Matrix' );

%%%%%%%%%%%% CALLBACKS %%%%%%%%%%%%%%
% motion callback, event "every" mouse movement
function motion_callback(src, event)
    % retrieve the current point
    currp = get(gcf,'CurrentPoint');
    % retrieve window geometry
    HGEOM = get( src, 'Position');
    % transform in sphere coordinates (3,4) = (WIDTH, HEIGHT)
    currp = point_on_sphere( currp, HGEOM(3), HGEOM(4) );

    % workaround condition (see point_on_sphere)
    if isnan(currp)
       return;
    end

    %%%%% ARCBALL COMPUTATION %%%%%
    if strcmp(mousestatus, 'buttondown')
        % compute angle and rotation axis
        rot_dir = cross( START, currp ); rot_dir = rot_dir / norm( rot_dir );
        rot_ang = acos( dot( currp, START ) );

        % convert direction in model coordinate system
        M_tr = inv( M_previous );
        rot_dir = M_tr*[rot_dir,0]';
        rot_dir = rot_dir(1:3);
        rot_dir = rot_dir / norm( rot_dir ); % renormalize
        % construct matrix
        R_matrix = makehgtform('axisrotate',rot_dir,rot_ang);
        % set hgt matrix
        set(hg,'Matrix',M_previous*R_matrix);
        % refresh drawing
        drawnow;
    end
end

% only 1 event on click
function buttondown_callback( src, evnt )
    % change status
    mousestatus = 'buttondown';
    % retrieve the current point
    currp = get(gcf,'CurrentPoint');
    % retrieve window geometry
    HGEOM = get( src, 'Position');
    % SET START POSITION
    START = point_on_sphere( currp, HGEOM(3), HGEOM(4) );
    % SET START MATRIX
    M_previous = get( hg, 'Matrix' );
end
function buttonup_callback( src, evnt )
    % change status
    mousestatus = 'buttonup';
    % reset the start position
    START = [0,0,0];
end

%%%%%%%%%%%% UTILITY FUNCTION %%%%%%%%%%%%%
function currp = point_on_sphere( currp, width, height )
    currp(3) = 0;

    % determine radius of the sphere
    R = min(width, height)/2;

    % TRANSFORM the point in window coordinate into
    % the coordinate of a sphere centered in middle window
    ORIGIN = [width/2, height/2, 0];
    currp = currp - ORIGIN;

    % normalize position to [-1:1] WRT unit sphere
    % centered at the origin of the window
    currp = currp / R;

    % if position is out of sphere, normalize it to
    % unit length
    L = sqrt( currp*currp' );
    if L > 1
       % currp = nan; % workaround to stop evaluation
       % disp('out of sphere');

       currp = currp / L;
       currp(3) = 0;
    else
       % add the Z coordinate to the scheme
       currp(3) = sqrt( 1 - currp(1)^2 - currp(2)^2 );
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END ARCBALL STUFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end % END OF CURVELAB

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function x = mysvec(X)
% This function is my implementation of the svec function of SDPT3

% SP

s = size(X,1);
ii = ones(s);
idx1 = find(triu(ii));
idx2 = find(triu(ii,1));
X(idx2) = X(idx2)*sqrt(2);
x = X(idx1);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function  res = SPIRiT(kernel,method, imSize)
%res = SPIRiT(kernel [,method, imSize])
%	Implementation of the SPIRiT kernel operator
%
%
%   Constructor inputs:
%       kernel: [kx,ky,nCoils, nCoils] the spirit 2D convolution kernel.
%               See corrMatrix.m, calibrate.m and the demos on how to
%               generate a kernel.
%       method: implementation type of the operator. There are three
%               'conv'  is a k-space convolution implementation (slowest)
%                       that operates on k-space data and produces k-space data.
%               'fft'   is an fft based implementation of the k-space
%                       convolution through image domain multiplication.
%                       It operates on k-sapce data and produces k-space data.
%               'image' The SPIRiT operator is applied to image space data.
%                       This is useful when using image-based non-cartesian
%                       reconstruction. In this case the SPIRiT operator
%                       operates on image data and produces image data.
%
%       imSize:  Size of the resulting image (only needed for 'fft' and
%                'image' modes
%
%
%   Example:
%
%   [x,y] = meshgrid(linspace(0,1,128));
%   % Generate fake Sensitivity maps
%   sMaps = cat(3,x.^2,1-x.^2,y.^2,1-y.^2);
%   % generate 4 coil phantom
%   imgs = repmat(phantom(128),[1,1,4]).*sMaps;
%   DATA = fft2c(imgs);
%   % crop 20x20 window from the center of k-space for calibration
%   kCalib = crop(DATA,[20,20,4]);
%
%   %calibrate a kernel
%   kSize = [5,5];
%   coils = 4;
%   kernel = zeros([kSize,coils,coils]);
%   [AtA,] = corrMatrix(kCalib,kSize);
%   for n=1:coils
%       kernel(:,:,:,n) = calibrate(AtA,kSize,coils,n,0.01);
%   end
%   GOP = SPIRiT(kernel, 'fft',[128,128]);
%
%   % undersample by a factor of 2
%   DATA(1:2:end,2:2:end,:) = 0;
%   DATA(2:2:end,1:2:end,:) = 0;
%
%   %reconstruct:
%   [res] = cgSPIRiT(DATA,GOP, 20, 1e-5, DATA);
%   figure, imshow(cat(2,sos(imgs), 2*sos(ifft2c(DATA)), sos(ifft2c(res))),[]);
%   title('full,  zero-fill,   result')
%
%   See Also:
%               calibrate.m, corrMatrix.m
%
% (c) Michael Lustig 2007

    if nargin < 2
	    method = 'conv';
	    KERNEL = [];
    end

    if strcmp(method,'conv')==1
	    KERNEL = [];
    end

    if strcmp(method,'fft')==1 & nargin < 3
	    error('must provide image size');
    end

    % for methods 'fft' and 'image' precompute the image domain kernel by
    % zero-padding and inverse fft

    if strcmp(method,'fft')==1 | strcmp(method,'image')==1
	    for n=1:size(kernel,4)
            KERNEL(:,:,:,n) = ifftnshift(zpad(kernel(end:-1:1,end:-1:1,:,n), imSize(1), imSize(2), size(kernel,3)),1:2);
%             KERNEL(:,:,:,n) = (ifft2c(zpad(kernel(end:-1:1,end:-1:1,:,n)*sqrt(imSize(1)*imSize(2)), imSize(1), imSize(2), size(kernel,3))));
        end
    end

    if strcmp(method,'fft')==0 & strcmp(method,'conv')==0 & strcmp(method,'image')==0
	    eror('no such method');
    end

    res.kernel = kernel;
    res.adjoint = 0;
    res.KERNEL = KERNEL;
    res.method = method;
    res.imSize = imSize;
    res = class(res,'SPIRiT');

author: thomaskuestner


text:
function bss_eval_test(T,N,W,L)
% check the behaviour of some basic functions of the toolbox BSS_EVAL
%
%
% Usage: bss_eval_test(T,N,L)
%
% Input:
%   - T: number of samples of the test signals (default 1000)
%   - N: number of sources in the test signals (default 2)
%   - W: hop size = half-size of the window (default ceil(T/20));
%   - L: number of taps of the filter (default 0)
%
% Ouput: None
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

if nargin<1
    T = 1000;
end
if nargin<2
    N=2;
end
if nargin<3
    W=ceil(T/20);
end
if nargin<4
    L=0;
end

% 1. Generate some data
S  = rand(N,T);

% 2. check that the SDR,SIR and SAR are almost infinite when performing a
% tvfilt decomposition with half overlapping rectangular windows
win = ones(1,W*2);
[starget,einterf,eartif]=bss_decomp_tvfilt(S(1,:),1,S,win,W,L);
[SDR,SIR,SAR]=bss_crit(starget,einterf,eartif);
disp(['SDR: ' num2str(SDR) ' SIR: ' num2str(SIR) 'SAR: ' num2str(SAR)]);
disp('Results for rectangular windows are expected to exceed 200 dB');

% 3. check that the SDR,SIR and SAR are almost infinite when performing a
% tvfilt decomposition with half overlapping triangular windows
win = triang(W*2)';
[starget,einterf,eartif]=bss_decomp_tvfilt(S(1,:),1,S,win,W,L);
[SDR,SIR,SAR]=bss_crit(starget,einterf,eartif);
disp(['SDR: ' num2str(SDR) ' SIR: ' num2str(SIR) 'SAR: ' num2str(SAR)]);
disp('Results for triangular windows are expected to exceed 200 dB');

%plot(starget,'r')
% hold on;
% triangsum = zeros(N*100+length(win),1);
% for i=1:100
%     idxrange = (i-1)*N+(1:length(win));
%     triangsum(idxrange) = triangsum(idxrange)+win';
% end
% plot(triangsum);

function w = triang(n)
% TRIANG Triangular window.
if rem(n,2)
% It's an odd length sequence
w = 2*(1:(n+1)/2)/(n+1);
w = [w w((n-1)/2:-1:1)]';
else
% It's even
w = (2*(1:(n+1)/2)-1)/n;
w = [w w(n/2:-1:1)]';
end

author: posenhuang

repo: deeplearningsourceseparation


text:
matlabbatch{1}.spm.tools.physio.save_dir = {'physio_out'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens_HCP';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'tfMRI_MOTOR_LR_Physio_log.txt'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'tfMRI_MOTOR_LR_Physio_log.txt'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = [];
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 91;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 0.72;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 284;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 46;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.respiratory.filter.passband = [0.01 2];
matlabbatch{1}.spm.tools.physio.preproc.respiratory.despike = false;
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

author: translationalneuromodeling

repo: tapas


text:
function mono_between_enum_test ( )

%*****************************************************************************80
%
%% MONO_BETWEEN_ENUM_TEST tests MONO_BETWEEN_ENUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_BETWEEN_ENUM_TEST\n' );
  fprintf ( 1, '  MONO_BETWEEN_ENUM can enumerate the number of monomials\n' );
  fprintf ( 1, '  in M variables, of total degree between N1 and N2.\n' );

  m = 3;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using spatial dimension M = %d\n', m );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N2:' );
  for n2 = 0 : 8
    fprintf ( 1, '  %4d', n2 );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N1 +------------------------------------------------------\n' );
  for n1 = 0 : 8
    fprintf ( 1, '  %2d |', n1 );
    for n2 = 0 : 8
      v = mono_between_enum ( m, n1, n2 );
      fprintf ( 1, '  %4d', v );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_example_bids_video

% MEM 4gb
% WALLTIME 00:10:00

%
%% Converting an example video dataset for sharing in BIDS
%
% The [BIDS standard](https://bids.neuroimaging.io) does currently not specify how to represent video recordings. This example - and the support that is implemented in the **[data2bids](https://github.com/fieldtrip/fieldtrip/blob/release/data2bids.m)** function - should be considered as a preliminary proposal to help researchers with their existing data. This example may also serve to start a discussion on whether and how this data type should be added to the [BIDS specification](http://bids-specification.readthedocs.io/).
%
% In many experiments video data can be recorded along with other measurements. This is often used to off-line screen the subjects behavior and to annotate the video, i.e. note the time segments with specific (desired or undesired) behavior. Neither the video recording  (which can be considered as raw data), nor the annotations (which can be considered as derived data) are currently formally part of BIDS, although it is not so difficult to come up with a way that these can be represented in a BIDS-like fashion.
%
%% # Example
%
% The example data is available from our [FTP server](ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/example/bids_video/).
%
% The example includes the original video recording under the “original” directory. It also has a copy of the script to do the conversion under “code”. The reorganized data is under the “bids” directory. According to [the documentation](https://bids-specification.readthedocs.io/en/stable/02-common-principles.html#source-vs-raw-vs-derived-data) the original data can be added to the BIDS dataset under the “sourcedata” directory and code can be added to the “code” directory. This way no information is lost and the conversion/reorganization is fully reproducible.
%
% In general: if your original data is in a different format than the BIDS representation (e.g., DICOM instead of NIFTI), you probably want to keep a copy of the original data, e.g., on a data acquisition collection on the Donders Repository. If it is in the same format like here - since we are not converting the video file, you could simply delete the original files after conversion. In either case - your own analyses and the shared data would be based on the BIDS representation.
%
cfg = [];

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

% this is required for dataset_description.json
cfg.dataset_description.Name                = 'Video example';
cfg.dataset_description.BIDSVersion         = 'unofficial extension';

% this is optional for dataset_description.json
cfg.dataset_description.License             = 'n/a';
cfg.dataset_description.Authors             = 'n/a';
cfg.dataset_description.Acknowledgements    = 'n/a';
cfg.dataset_description.Funding             = 'n/a';
cfg.dataset_description.ReferencesAndLinks  = 'n/a';
cfg.dataset_description.DatasetDOI          = 'n/a';

% provide some metadata about the task that was performed
cfg.TaskDescription = 'The subject was walking, stopping, turning, etcetera, while wearing a mobile EEG system on his head and a mobile NIRS system in his pocket';

cfg.method    = 'copy'; % video should simply be copied, not converted
cfg.dataset   = dccnpath('/home/common/matlab/fieldtrip/data/ftp/example/bids_video/original/Video_20190411_210224.mp4');
cfg.bidsroot  = fullfile(tempdir, 'bids');
cfg.datatype  = 'video';
cfg.sub       = '01';
cfg.task      = 'walking';

% construct an empty table that has the required columns 'onset' and 'duration' (in seconds)
% this uses the optional 'trial_type' column to indicate what happened
onset      = nan;
duration   = nan;
trial_type = {''};
cfg.events = table(onset, duration, trial_type);

% add the actual  events with onset, duration and trial_type
cfg.events(1,:) = {12, 2, 'standing up'};
cfg.events(2,:) = {15, 11, 'walking'};
cfg.events(3,:) = {26, 3, 'right turn'};
cfg.events(4,:) = {29, 10, 'walking'};

data2bids(cfg);
%
% This example also shows how to add annotations of the video to the `_events.tsv` file. If you do not have annotations (yet), you would leave the `cfg.events` field empty.
%
% In case you would have multiple subsequent videos in a single session, you would use `cfg.run=1, 2, ...` to specify the successive recordings.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function clim_positive(hax)

if nargin < 1
  hax = gca;
end

cl = get(hax, 'clim');
cl = max(cl);
set(hax, 'clim', [0 cl]);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function bisection_rc_test ( )

%*****************************************************************************80
%
%% BISECTION_RC_TEST tests the BISECTION_RC library.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/bisection_rc/bisection_rc_test.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 January 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BISECTION_RC_TEST:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the BISECTION_RC library.\n' );

  bisection_rc_test01 ( );
  bisection_rc_test02 ( );
  bisection_rc_test03 ( );
  bisection_rc_test04 ( );
  bisection_rc_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BISECTION_RC_PRB:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%%*****************************************************************
%% HSDsqlpmisc:
%% produce infeasibility certificates if appropriate
%%
%% Input: X,y,Z are the original variables, not the HSD variables.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [X,y,Z,resid,reldist,param,msg] = HSDsqlpmisc(blk,At,C,b,X,y,Z,permZ,param);

   obj         = param.obj;
   relgap      = param.relgap;
   prim_infeas = param.prim_infeas;
   dual_infeas = param.dual_infeas;
   ZpATynorm   = param.ZpATynorm;
   inftol      = param.inftol;
   m0          = param.m0;
   indeprows   = param.indeprows;
   termcode    = param.termcode;
   AX          = param.AX;
   normX0      = param.normX0;
   normZ0      = param.normZ0;
   printlevel  = param.printlevel;
%%
   resid = []; reldist = []; msg = [];
   Anorm = ops(At,'norm'); xnorm = ops(X,'norm'); ynorm = norm(y);
%%
   if (termcode <= 0)
      %%
      %% To detect near-infeasibility when the algorithm provides
      %% a "better" certificate of infeasibility than of optimality.
      %%
      err = max([prim_infeas,dual_infeas,relgap]);
      iflag = 0;
      if (obj(2) > 0)
         homRd = ZpATynorm/obj(2);
         if (homRd < 1e-2*sqrt(err*inftol))
            iflag = 1;
            termcode = 1;
            param.termcode = 1;
         end
      elseif (obj(1) < 0)
         homrp = norm(AX)/(-obj(1));
         if (homrp < 1e-2*sqrt(err*inftol))
            iflag = 1;
            termcode = 2;
            param.termcode = 2;
         end
      end
   end
   if (termcode == 1)
      rby = 1/(b'*y); y = rby*y; Z = ops(Z,'*',rby);
      resid = ZpATynorm * rby;
      reldist = ZpATynorm/(Anorm*ynorm);
      msg = 'Stop: primal problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode == 2)
      tCX = blktrace(blk,C,X);
      X = ops(X,'*',1/(-tCX));
      resid = norm(AX) /(-tCX);
      reldist = norm(AX)/(Anorm*xnorm);
      msg = 'Stop: dual problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode == 3)
      maxblowup = max(ops(X,'norm')/normX0,ops(Z,'norm')/normZ0);
      msg = sprintf('Stop: primal or dual is diverging, %3.1e',maxblowup);
      if (printlevel); fprintf('\n  %s',msg); end
   end
   [X,Z] = unperm(blk,permZ,X,Z);
   if ~isempty(indeprows)
      ytmp = zeros(m0,1);
      ytmp(indeprows) = y;
      y = ytmp;
   end
%%*****************************************************************************
%% unperm: undo the permutations applied in validate.
%%
%% [X,Z,Xiter,Ziter] = unperm(blk,permZ,X,Z,Xiter,Ziter);
%%
%% undoes the permutation introduced in validate.
%% can also be called if Xiter and Ziter have not been set as
%%
%% [X,Z] = unperm(blk,permZ,X,Z);
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

  function [X,Z] = unperm(blk,permZ,X,Z);
%%
  for p = 1:size(blk,1)
     if (strcmp(blk{p,1},'s') & ~isempty(permZ{p}))
        per = permZ{p};
        X{p} = X{p}(per,per);
        Z{p} = Z{p}(per,per);
     end
  end
%%*****************************************************************************


text:
function plotAreaUnderACurve(handles, integrationLimits, coefficients, plot1_xAxis, plot1_yAxis)

% Get the handle to the axes
axes(handles.integrationPlotAxes)

% Find the roots of the polynomial
fx = coefficients(end:-1:1);
rootsFx = roots(fx);
realRootsFx = [];
for idx = 1:length(rootsFx)
    if isreal(rootsFx(idx))
        realRootsFx(end+1,1) = rootsFx(idx);
    end
end

prevXstart = integrationLimits(1);

if ~isempty(realRootsFx)
    realRootsFx = sortrows(realRootsFx);
    % Plot the shaded area (above x-axis in blue and below x-axis in red)
    for idx = 1:length(realRootsFx)
        if ((realRootsFx(idx) > integrationLimits(1)) && (realRootsFx(idx) < integrationLimits(2)))
            % Set the x and y points
            xiPrec = (realRootsFx(idx) - prevXstart)/20;
            xi = (prevXstart:xiPrec:realRootsFx(idx));
            yi = coefficients(1)+coefficients(2)*xi+coefficients(3)*xi.^2+coefficients(4)*xi.^3+coefficients(5)*xi.^4+coefficients(6)*xi.^5;
            faceColor = [0 0.5 1];
            if sign(yi(end-1)) == -1
                faceColor = [1 0 0];
            end
            area(xi,yi,'FaceColor',faceColor);
            hold on
            prevXstart = realRootsFx(idx);
        end
    end
end

xiPrec = (integrationLimits(2) - prevXstart)/20;
xi = (prevXstart:xiPrec:integrationLimits(2));
yi = coefficients(1)+coefficients(2)*xi+coefficients(3)*xi.^2+coefficients(4)*xi.^3+coefficients(5)*xi.^4+coefficients(6)*xi.^5;
faceColor = [0 0.5 1];
if sign(yi(end-1)) == -1
    faceColor = [1 0 0];
end
area(xi,yi,'FaceColor',faceColor);
hold on

x = (integrationLimits(1):integrationLimits(2));
y = coefficients(1)+coefficients(2)*x+coefficients(3)*x.^2+coefficients(4)*x.^3+coefficients(5)*x.^4+coefficients(6)*x.^5;
plot(x,y,'LineWidth',2.5,'Color',[0 0 0]);
xlabel('x','Fontweight','b');
ylabel('f(x)','Fontweight','b');
axis([plot1_xAxis(1) plot1_xAxis(end) plot1_yAxis(1) plot1_yAxis(end)]);
title('Integration');
grid on
hold off

author: Sable

repo: mcbench-benchmarks


text:
function [interpolated,discarded] = Interpolate(samples,timestamps,varargin)

%Interpolate - Interpolate samples (positions, spikes, LFP, etc.) at given timestamps.
%
% Interpolate samples (positions, spikes, LFP, etc.) at given timestamps. As a special
% case, a point process can also be 'interpolated' by changing each original timestamp
% t to the closest target timestamp t* such that t* <= t.
%
%  USAGE
%
%    [interpolated,discarded] = Interpolate(samples,timestamps,<options>)
%
%    samples        a list of samples to interpolate
%    timestamps     a list of timestamps
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'type'        'linear' if samples are linear values (default), or
%                   'circular' otherwise
%     'trim'        discard ('on', default) or keep ('off') samples that would
%                   require extrapolation. If kept, these samples are set to
%                   NaN (this can be useful when the output is required to
%                   have the same size as the input) (default = 'on')
%     'maxGap'      time gaps between original and target times exceeding this
%                   threshold will be ignored (default = Inf)
%    =========================================================================
%
%  OUTPUT
%
%    interpolated   list of interpolated samples
%    discarded      logical vector indicating for which timestamps the samples were
%                   not interpolated (see option 'maxGap')
%
%  SEE
%
%    See also Diff.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Default values
trim = 'on';
maxGap = Inf;
type = 'linear';

if nargin < 2 | mod(length(varargin),2) ~= 0,
	error('Incorrect number of parameters (type ''help Interpolate'' for details).');
end
if ~isdvector(timestamps),
	error('Incorrect timestamps - should be a vector (type ''help Interpolate'' for details).');
end
timestamps = timestamps(:);

% Parse parameter list
for j = 1:2:length(varargin),
	if ~ischar(varargin{j}),
		error(['Parameter ' num2str(j+2) ' is not a property (type ''help Interpolate'' for details).']);
	end
	switch(lower(varargin{j})),
		case 'type',
			type = varargin{j+1};
			if ~isstring_FMAT(type,'linear','circular'),
				error('Incorrect value for property ''type'' (type ''help Interpolate'' for details).');
			end
		case 'trim',
			trim = lower(varargin{j+1});
			if ~isstring_FMAT(trim,'on','off'),
				error('Incorrect value for property ''trim'' (type ''help Interpolate'' for details).');
			end
		case 'maxgap',
			maxGap = varargin{j+1};
			if ~isdscalar(maxGap,'>0'),
				error('Incorrect value for property ''maxGap'' (type ''help Interpolate'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{j}) ''' (type ''help Interpolate'' for details).']);
	end
end

interpolated = [];
if isempty(timestamps) | isempty(samples), return; end

pointProcess = isvector(samples);

% Determine which target timestamps are too far away from sample timestamps (isolated) and should be discarded
[unused,isolated] = Match(timestamps,samples(:,1),'match','closest','error',maxGap);

if pointProcess,
	% 'Interpolate' samples at timestamps (see help above)
	[interpolated,unused,discarded] = Match(samples,timestamps,'error',maxGap);
else
	% Determine which timestamps would require extrapolation
	outside = logical(zeros(size(timestamps)));
	if strcmp(trim,'on'),
		outside = timestamps<samples(1,1)|timestamps>samples(end,1);
	end
	% Discard isolated timestamps and timestamps that would require extrapolation
	discarded = outside|isolated;
	timestamps(discarded) = [];
	% Circular data?
	if strcmp(type,'circular'),
		range = isradians(samples(:,2:end));
		samples(:,2:end) = exp(i*samples(:,2:end));
	end
	% Interpolate samples at timestamps
	interpolated = [timestamps interp1(samples(:,1),samples(:,2:end),timestamps)];
	% Circular data?
	if strcmp(type,'circular'),
		interpolated(:,2:end) = wrap(angle(interpolated(:,2:end)),range);
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% Setting of the problem
global s
pde = fracLapdata9;
% pde = checkboarddata;
pde.L = 1;
option.maxIt = 4;
option.maxN = 1e6;
option.elemType = 'P1P1';
option.solver = 'mg';
option.gNquadorder = 4;
option.tol = 1e-6;
[node,elem] = squaremesh([-1 1 -1 1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

%% s = 0.2
s = 0.2;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.4
s = 0.4;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
femfracLap(node,elem,pde,bdFlag,option);

repo: ifem


text:
function a = ne( x, y )

%   Disciplined convex/geometric programming information for NE (~=):
%      Not-equal constraints violate both the DCP and DGP rulesets. Thus
%      not-equal expressions may only appear in CVX models when both
%      sides are constant.

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '~=' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function varargout = meshReduce(nodes, varargin)
%MESHREDUCE Merge coplanar faces of a polyhedral mesh
%
%   Note: deprecated, should use "mergeCoplanarFaces" instead
%
%   [NODES FACES] = meshReduce(NODES, FACES)
%   [NODES EDGES FACES] = meshReduce(NODES, EDGES, FACES)
%   NODES is a set of 3D points (as a Nn-by-3 array),
%   and FACES is one of:
%   - a Nf-by-X array containing vertex indices of each face, with each
%   face having the same number of vertices,
%   - a Nf-by 1 cell array, each cell containing indices of a face.
%   The function groups faces which are coplanar and contiguous, resulting
%   in a "lighter" mesh. This can be useful for visualizing binary 3D
%   images for example.
%
%   FACES = meshReduce(..., PRECISION)
%   Adjust the threshold for deciding if two faces are coplanar or
%   parallel. Default value is 1e-14.
%
%   Example
%   [n e f] = createCube;
%   figure; drawMesh(n, f); view(3); axis equal;
%   f2 = meshReduce(n, f);
%   figure; drawMesh(n, f2); view(3); axis equal;
%
%   See also
%   meshes3d, mergeCoplanarFaces
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-07-05
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% 20/07/2006 add tolerance for coplanarity test
% 21/08/2006 fix small bug due to difference of methods to test
%   coplanaritity, sometimes resulting in 3 points of a face not coplanar !
%   Also add control on precision
% 14/08/2007 rename minConvexHull->meshReduce, and extend to non convex
%   shapes
% 2011-01-14 code clean up
% 2013-02-22 deprecate and rename to mergeCoplanarFaces

warning('MatGeom:meshes3d:deprecated', ...
    'Function ''meshReduce'' is deprecated, should use ''mergeCoplanarFaces'' instead');

%% Process input arguments

% set up precision
acc = 1e-14;
if ~isempty(varargin)
    var = varargin{end};
    if length(var)==1
        acc = var;
        varargin(end) = [];
    end
end

% extract faces and edges
if length(varargin)==1
    faces = varargin{1};
else
    faces = varargin{2};
end

%% Initialisations

% number of faces
Nn = size(nodes, 1);
Nf = size(faces, 1);

% compute number of vertices of each face
Fn = ones(Nf, 1)*size(faces, 2);

% compute normal of each faces
normals = faceNormal(nodes, faces);

% initialize empty faces and edges
faces2  = cell(0, 1);
edges2  = zeros(0, 2);

% Processing flag for each face
% 1: face to process, 0: already processed
% in the beginning, every triangle face need to be processed
flag = ones(Nf, 1);

%% Main iteration

% iterate on each  face
for f = 1:Nf

    % check if face was already performed
    if ~flag(f)
        continue;
    end

    % indices of faces with same normal
    ind = find(abs(vectorNorm3d(cross(repmat(normals(f, :), [Nf 1]), normals)))<acc);
    %ind = ind(ind~=i);

    % keep only coplanar faces (test coplanarity of points in both face)
    ind2 = false(size(ind));
    for j = 1:length(ind)
        ind2(j) = isCoplanar(nodes([faces(f,:) faces(ind(j),:)], :), acc);
    end
    ind2 = ind(ind2);

    % compute edges of all faces in the plane
    planeEdges  = zeros(sum(Fn(ind2)), 2);
    pos  = 1;
    for i=1:length(ind2)
        face = faces(ind2(i), :);
        faceEdges = sort([face' face([2:end 1])'], 2);
        planeEdges(pos:sum(Fn(ind2(1:i))), :) = faceEdges;
        pos = sum(Fn(ind2(1:i)))+1;
    end
    planeEdges = unique(planeEdges, 'rows');

    % relabel plane edges, and find connected components
    [planeNodes, I, J] = unique(planeEdges(:)); %#ok<ASGLU>
    planeEdges2 = reshape(J, size(planeEdges));
    component   = grLabel(nodes(planeNodes, :), planeEdges2);

    % compute degree (number of adjacent faces) of each edge.
    Npe = size(planeEdges, 1);
    edgesDegree = zeros(Npe, 1);
    for i = 1:length(ind2)
        face = faces(ind2(i), :);
        faceEdges = sort([face' face([2:end 1])'], 2);
        for j = 1:size(faceEdges, 1)
            indEdge = find(sum(ismember(planeEdges, faceEdges(j,:)),2)==2);
            edgesDegree(indEdge) = edgesDegree(indEdge)+1;
        end
    end

    % extract unique edges and nodes of the plane
    planeEdges  = planeEdges(edgesDegree==1, :);
    planeEdges2 = planeEdges2(edgesDegree==1, :);

    % find connected component of each edge
    planeEdgesComp = zeros(size(planeEdges, 1), 1);
    for e = 1:size(planeEdges, 1)
        planeEdgesComp(e) = component(planeEdges2(e, 1));
    end

    % iterate on connected faces
    for c=1:max(component)

        % convert to chains of nodes
        loops = graph2Contours(nodes, planeEdges(planeEdgesComp==c, :));

        % add a simple Polygon for each loop
        facePolygon = loops{1};
        for l = 2:length(loops)
            facePolygon = [facePolygon, NaN, loops{l}]; %#ok<AGROW>
        end
        faces2{length(faces2)+1, 1}  = facePolygon;

        % also add news edges
        edges2 = unique([edges2; planeEdges], 'rows');
    end

    % mark processed faces
    flag(ind2) = 0;
end

%% Additional processing on nodes

% select only nodes which appear in at least one edge
indNodes = unique(edges2(:));

% for each node, compute index of corresponding new node (or 0 if dropped)
refNodes = zeros(Nn, 1);
for i=1:length(indNodes)
    refNodes(indNodes(i)) = i;
end

% changes indices of nodes in edges2 array
for i=1:length(edges2(:))
    edges2(i) = refNodes(edges2(i));
end

% changes indices of nodes in faces2 array
for f=1:length(faces2)
    face = faces2{f};
    for i=1:length(face)
        if ~isnan(face(i))
            face(i) = refNodes(face(i));
        end
    end
    faces2{f} = face;
end

% keep only boundary nodes
nodes2 = nodes(indNodes, :);

%% Process output arguments

if nargout == 1
    varargout{1} = faces2;
elseif nargout == 2
    varargout{1} = nodes2;
    varargout{2} = faces2;
elseif nargout==3
    varargout{1} = nodes2;
    varargout{2} = edges2;
    varargout{3} = faces2;
end

function labels = grLabel(nodes, edges)
%GRLABEL associate a label to each connected component of the graph
%   LABELS = grLabel(NODES, EDGES)
%   Returns an array with as many rows as the array NODES, containing index
%   number of each connected component of the graph. If the graph is
%   totally connected, returns an array of 1.
%
%   Example
%       nodes = rand(6, 2);
%       edges = [1 2;1 3;4 6];
%       labels = grLabel(nodes, edges);
%   labels =
%       1
%       1
%       1
%       2
%       3
%       2
%
%   See also
%   getNeighbourNodes
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-08-14,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% init
Nn = size(nodes, 1);
labels = (1:Nn)';

% iteration
modif = true;
while modif
    modif = false;

    for i=1:Nn
        neigh = getNeighbourNodes(i, edges);
        neighLabels = labels([i;neigh]);

        % check for a modification
        if length(unique(neighLabels))>1
            modif = true;
        end

        % put new labels
        labels(ismember(labels, neighLabels)) = min(neighLabels);
    end
end

% change to have fewer labels
labels2 = unique(labels);
for i=1:length(labels2)
    labels(labels==labels2(i)) = i;
end

function nodes2 = getNeighbourNodes(node, edges)
%GETNEIGHBOURNODES find nodes adjacent to a given node
%
%   NEIGHS = getNeighbourNodes(NODE, EDGES)
%   NODE: index of the node
%   EDGES: the complete edges list
%   NEIGHS: the nodes adjacent to the given node.
%
%   NODE can also be a vector of node indices, in this case the result is
%   the set of neighbors of any input node.
%
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 16/08/2004.
%

%   HISTORY
%   10/02/2004 documentation
%   13/07/2004 faster algorithm
%   03/10/2007 can specify several input nodes

[i, j] = find(ismember(edges, node));
nodes2 = edges(i,1:2);
nodes2 = unique(nodes2(:));
nodes2 = sort(nodes2(~ismember(nodes2, node)));

function curves = graph2Contours(nodes, edges)
%GRAPH2CONTOURS convert a graph to a set of contour curves
%
%   CONTOURS = GRAPH2CONTOURS(NODES, EDGES)
%   NODES, EDGES is a graph representation (type "help graph" for details)
%   The algorithm assume every node has degree 2, and the set of edges
%   forms only closed loops. The result is a list of indices arrays, each
%   array containing consecutive point indices of a contour.
%
%   To transform contours into drawable curves, please use :
%   CURVES{i} = NODES(CONTOURS{i}, :);
%
%
%   NOTE : contours are not oriented. To manage contour orientation, edges
%   also need to be oriented. So we must precise generation of edges.
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 05/08/2004.
%

curves = {};
c = 0;

while size(edges,1)>0
	% find first point of the curve
	n0 = edges(1,1);
    curve = n0;

    % second point of the curve
	n = edges(1,2);
	e = 1;

	while true
        % add current point to the curve
		curve = [curve n];

        % remove current edge from the list
        edges = edges((1:size(edges,1))~=e,:);

		% find index of edge containing reference to current node
		e = find(edges(:,1)==n | edges(:,2)==n);
        e = e(1);

		% get index of next current node
        % (this is the other node of the current edge)
		if edges(e,1)==n
            n = edges(e,2);
		else
            n = edges(e,1);
		end

        % if node is same as start node, loop is closed, and we stop
        % node iteration.
        if n==n0
            break;
        end
	end

    % remove the last edge of the curve from edge list.
    edges = edges((1:size(edges,1))~=e,:);

    % add the current curve to the list, and start a new curve
    c = c+1;
    curves{c} = curve;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function perm=topological_order(A,varargin)
% TOPOLOGICAL_ORDER Returns the topological ordering of vertices in a dag
%
% perm=topological_order(A) generates a permutation of the vertices that
% indices a topological order.  A topological order satisfies the property
% that edge (i,j) in A implies i comes before j in the order.  A
% topological order only exists for directed acyclic graphs.  If A is not a
% directed acyclic graph, then the returned perm is empty.
%
% This method works on directed graphs.
% The runtime is O(V+E).
%
% ... = topological_order(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   There are no additional options for this function.
%
% Note: this function does not depend upon the non-zero values of A, but
% only uses the non-zero structure of A.
%
% Example:
%   n = 10; A = sparse(1:n-1, 2:n, 1, n, n); % construct a simple dag
%   p = topological_order(A);
%
% See also TEST_DAG

% David Gleich
% Copyright, Stanford University, 2007-2008

%% History
%  2008-09-23: Reformatted common section.
%  2008-10-07: Fixed example.
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end
if check, check_matlab_bgl(A,struct()); end
if trans, A = A'; end

[perm dag] = topological_order_mex(A);

% in this case it wasn't a dag
if dag == 0
    perm = [];
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function T = myrepmat(T, sizes)
% MYREPMAT Like the built-in repmat, except myrepmat(T,n) == repmat(T,[n 1])
% T = myrepmat(T, sizes)

if length(sizes)==1
  T = repmat(T, [sizes 1]);
else
  T = repmat(T, sizes(:)');
end


text:
function [ gradientParams, SigmaInvs, CholDecomps, Sigmas, bs, allXresp] = gradientCCNF( params, num_alpha, numBeta, sizeTheta, lambda_a, lambda_b, lambda_th, Precalc_Bs, x, y, Precalc_yBys, Precalc_Bs_flat, constant, num_seqs)
%gradientCCNF Summary of this function goes here
%   Detailed explanation goes here

    % pick out the relevant terms (unpack)
    alphas_init = params(1:num_alpha);
    betasInit = params(num_alpha+1:num_alpha+numBeta);
    thetasInit = reshape(params(num_alpha+numBeta+1:end), sizeTheta);

    % Compute the response from the neural layers
    allXresp = 1./(1 + exp(-thetasInit * x));
    Xt = x;

    bs = 2*alphas_init' * allXresp;

    % This is precalculated for the next step and is basically the
    % feedforward step of the neural net
    Z_precalc = 2 * (allXresp .* (1-allXresp));

    % These are the outputs weighted by the alphas (see eq TODO
    db2_precalc =  bsxfun(@times, Z_precalc, alphas_init);

    num_feats = sizeTheta(2);

    if(constant)

        seq_length = size(x,2)/num_seqs;

        % As the similarities are the same across all series we can reuse our
        % Sigma and SigmaInv calculations

        I = eye(seq_length);
        [SigmaInv] = CalcSigmaCCNFflat(alphas_init, betasInit, seq_length, Precalc_Bs_flat{1}, I, zeros(seq_length));
        CholDecomp=chol(SigmaInv);

        % This is a faster way of inverting a symmetric matrix
        Sigma=CholDecomp\(CholDecomp'\I);
        Sigma_trace = trace(Sigma);

        % mu values associated with each time step
        mus = Sigma * reshape(bs, seq_length, num_seqs);

        % difference between actual and prediction (error)
        diff = (y - mus);

        db_precalc_mult = bsxfun(@times, db2_precalc, diff(:)');

        % Equation 46 from the appendix
        gradientThetasT = Xt * db_precalc_mult';

        % Reshape into the correct format
        gradientThetasT = gradientThetasT(:)';

        gradientThetasT = reshape(gradientThetasT, sizeTheta(2), sizeTheta(1))';
        gradientThetasT = gradientThetasT(:);

        % Some useful precalculations

        % for every sequence get a dot product with itself
        yy = dot(y,y);

        % same goes for the mu
        mumu = dot(mus,mus);

        % calculating the derivative of L with respect to alpha_k (Equation 27)
        % gradientAlphas =  (-yq'*yq +(2*yq'*D')' -2 * D * mu + sum(mu.^2) + trace(Sigma));
        % allXresp is D
        gradient_alphas_add = -sum(yy) + sum(mumu) + num_seqs * Sigma_trace;
        gradient_alphas = 2 * allXresp * (y(:) - mus(:)) + gradient_alphas_add;

        gradient_betas = zeros(numBeta, 1);

        % calculating the derivative of log(L) with respect to the betas
        for k=1:numBeta

            % From Equation 38 (and 39 for gamma)
            % gradient = -yq'*B^(k)*yq + mu'*B^(k)*mu + Vec(Sigma)'*vec(B^(k)

            % We precalculate B^(k) (equation 30), as it does not change
            % over the course of optimisation
            B_k = Precalc_Bs{1}{k};

            % precalculated -yq'*B_k*yq can be used as well as it does not
            % change (stored in Precalc_yBys)
            yq_B_k_yq = sum(Precalc_yBys(:,k));

            % A vectorised version of mu'*B^(k)*mu
            B_k_mu = B_k*mus;
            mu_B_k_mu = mus(:)' * B_k_mu(:);

            % Vec(Sigma)*Vec(B^(k)) can be computed as follows:
            partition_term = num_seqs * Sigma(:)'*B_k(:);

            % Equation 38 and 39 basically
            dLdb = yq_B_k_yq + mu_B_k_mu + partition_term;

            gradient_betas(k) = dLdb;
        end

        gradientParams = [gradient_alphas;gradient_betas;gradientThetasT];

        SigmaInvs = SigmaInv;
        CholDecomps = CholDecomp;
        Sigmas = Sigma;

    else

        SigmaInvs = cell(num_seqs, 1);
        CholDecomps = cell(num_seqs, 1);
        Sigmas = cell(num_seqs, 1);
        gradients = zeros(num_seqs, numel(params));

        a_precalc = zeros(sizeTheta(2)*sizeTheta(1), size(allXresp,2));

        for i=1:size(db2_precalc,1)
            a_precalc((i-1)*num_feats+1:i*num_feats,:) = bsxfun(@times, Xt, db2_precalc(i,:));
        end

        % y can either be in cell format (diff length seqs.) or in matrix
        %, same length seqs
        beg_ind = 1;

        if(iscell(y))
            end_ind = numel(y{1});
            y_cell = true;
        else
            end_ind = size(y,1);
            y_cell = false;
        end

        % Go through every sequence summing the gradients
        for q = 1 : num_seqs

            currResp = allXresp(:,beg_ind:end_ind);
            currB = bs(beg_ind:end_ind)';

            PrecalcB = Precalc_Bs{q};
            PrecalcBFlat = Precalc_Bs_flat{q};

            if(y_cell)
                yq = y{q};
            else
                yq = y(:,q);
            end

            xq = x(2:end, beg_ind:end_ind);

            % Used for equation 46 computation
            a_precalc_curr = a_precalc(:,beg_ind:end_ind);

            precalc_eye = eye(numel(yq));
            precalc_zeros = zeros(numel(yq));

            [ gradientsAlphas, gradientsBetas, gradientsThetas, SigmaInv, CholDecomp, Sigma ] = gradientCCNF_per_seq(alphas_init, betasInit, thetasInit, PrecalcB, xq, yq, currResp, currB, Precalc_yBys(q, :), PrecalcBFlat, a_precalc_curr, precalc_eye, precalc_zeros);

            gradients(q,:) = [gradientsAlphas; gradientsBetas; gradientsThetas(:)];
            SigmaInvs{q} = SigmaInv;
            CholDecomps{q} = CholDecomp;
            Sigmas{q} = Sigma;

            % Update the references to sequence start/end
            if(q ~= num_seqs)
                beg_ind = end_ind + 1;
                if(iscell(y))
                    end_ind = end_ind + numel(y{q+1});
                else
                    end_ind = end_ind + size(y,1);
                end
            end

        end

        gradientParams = sum(gradients,1)';

    end

    % Add the regularisation term
    regAlpha = alphas_init * lambda_a;
    regBeta = betasInit * lambda_b;
    regTheta = thetasInit * lambda_th;

    gradientParams = gradientParams - [regAlpha; regBeta; regTheta(:)];
end

author: TadasBaltrusaitis

repo: OpenFace


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Goal: Convert NIfTI image to bim format
% Input:
%       inObjectName - File name of an input object with .nii or .hdr
%       extension
%
%       resampleFactor - Resampling factor to down- or up-sample a binary
%       object
%
% Output:
%       bim - a 3D matrix of binary voxels.
%       mins - lowest coordinate values in x, y, and z axis (1 x 3)
%
% 03/31/2001 - created by Li Shen
% Modified by Sungeun Kim (10-09-08)

function [bim, mins] = NIfTI2bim(inObjectName, resampleFactor)

ni = load_nii(inObjectName);
roi = ni.img;
clear ni;

DIM = size(roi);

if resampleFactor ~= 1   % Upsampling or downsampling
    roi = imresize(roi, resampleFactor, 'nearest');

    space = 1/resampleFactor;  % Interpolate binary volume along z axis
    if space < 1        % Up-sampling
        [xInt, yInt, zInt] = meshgrid(1:size(roi,1), 1:size(roi,2), 1:space:size(roi,3)+space);
    elseif space > 1    % Down-sampling
        [xInt, yInt, zInt] = meshgrid(1:size(roi,1), 1:size(roi,2), 1:space:size(roi,3));
    end
    bim = interp3(roi, xInt, yInt, zInt,'nearest');
else
    bim = roi;
end

% Make binary image
ind = find(bim>0); bim(ind) = 1;
ind = find(bim<1); bim(ind) = 0;

mins = [0 0 0];

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function out=atan2(y,x)

if isa(x,'mp')
 precision=x(1).precision;
 if ~isa(y,'mp')
  y=mp(y,precision);
 end
else
 precision=y(1).precision;
 if ~isa(x,'mp')
  x=mp(x,precision);
 end
end

if any(~isreal(x)) | any(~isreal(y))
 warning('atan2 for mp objects currently ignores imaginary parts')
 x=real(x);
 y=real(y);
end

ex=numel(x); ey=numel(y);
if ex==1
 x=ones(size(y))*x;
elseif ey==1
 y=ones(size(x))*y;
end

out=atan(y./x);
mpPi=mppi(precision);
for ii=1:numel(x)
 if x(ii)<0
  if y(ii)>0
   out(ii)=out(ii)+mpPi;
  else
   out(ii)=out(ii)-mpPi;
  end
 end
end

out=real(out);

author: opencobra

repo: cobratoolbox


text:
function tests = test_spm_dcm_fit
% Unit Tests for spm_dcm_fit
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_fit.m 7365 2018-07-03 14:10:43Z peter $

tests = functiontests(localfunctions);

function test_fmri_noparfor(testCase)

data_path = get_data_path();

% Load DCMs
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM;
GCM = GCM(1:1,1);

% Remove old posteriors
for i = 1:length(GCM)
    GCM{i} = rmfield(GCM{i},'Ep');
    GCM{i} = rmfield(GCM{i},'F');
    GCM{i}.options.nograph = true;
end

% Estimate
GCM = spm_dcm_fit(GCM);

% Check DCM
testCase.assertTrue(length(GCM) == 1);
for i = 1:length(GCM)
    testCase.assertTrue(isfield(GCM{i},'Ep'));
    testCase.assertTrue(isfield(GCM{i},'F'));
end

% -------------------------------------------------------------------------
function test_fmri_parfor(testCase)

data_path = get_data_path();

% Load DCM (one will be enough)
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
GCM = GCM.GCM;
GCM = GCM(1:1,1);

% Remove old posteriors
for i = 1:length(GCM)
    GCM{i} = rmfield(GCM{i},'Ep');
    GCM{i} = rmfield(GCM{i},'F');
    GCM{i}.options.nograph = true;
end

% Estimate
GCM = spm_dcm_fit(GCM, true);

% Check DCM
testCase.assertTrue(length(GCM) == 1);
for i = 1:length(GCM)
    testCase.assertTrue(isfield(GCM{i},'Ep'));
    testCase.assertTrue(isfield(GCM{i},'F'));
end

% -------------------------------------------------------------------------
function data_path = get_data_path()

data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region');

author: spm


text:
% .............................................................
    function [shape,naturalDerivatives]=shapeFunctionQ12(xi,eta)

    % shape function and derivatives for Q12 elements
    % shape : Shape functions
    % naturalDerivatives: derivatives w.r.t. xi and eta
    % xi, eta: natural coordinates (-1 ... +1)

    shape=[1/32*(1-xi)*(1-eta)*(-10+9*(xi^2+eta^2));
           1/32*(1+xi)*(1-eta)*(-10+9*(xi^2+eta^2));
           1/32*(1+xi)*(1+eta)*(-10+9*(xi^2+eta^2));
           1/32*(1-xi)*(1+eta)*(-10+9*(xi^2+eta^2));
           9/32*(1-xi^2)*(1-1/3*xi)*(1-eta);
           9/32*(1-xi^2)*(1+1/3*xi)*(1-eta);
           9/32*(1+xi)*(1-eta^2)*(1-1/3*eta);
           9/32*(1+xi)*(1-eta^2)*(1+1/3*eta);
           9/32*(1-xi^2)*(1+1/3*xi)*(1+eta);
           9/32*(1-xi^2)*(1-1/3*xi)*(1+eta);
           9/32*(1-xi)*(1-eta^2)*(1+1/3*eta);
           9/32*(1-xi)*(1-eta^2)*(1-1/3*eta)];

 naturalDerivatives=...
       [9/16*(1-eta)*(1-xi)*xi-(1/32)*(1-eta)*(-10+9*(eta^2+xi^2)),9/16*(1-eta)*eta*(1-xi)-1/32*(1-xi)*(-10+9*(eta^2+xi^2));
        9/16*(1-eta)*xi*(1+xi)+(1/32)*(1-eta)*(-10+9*(eta^2+xi^2)),9/16*(1-eta)*eta*(1+xi)-(1/32)*(1+xi)*(-10+9*(eta^2+xi^2)) ;
        9/16*(1+eta)*xi*(1+xi)+(1/32)*(1+eta)*(-10+9*(eta^2+xi^2)),9/16*eta*(1+eta)*(1+xi)+(1/32)*(1+xi)*(-10+9*(eta^2+xi^2)) ;
        9/16*(1+eta)*(1-xi)*xi-(1/32)*(1+eta)*(-10+9*(eta^2+xi^2)),9/16*eta*(1+eta)*(1-xi)+(1/32)*(1-xi)*(-10+9*(eta^2+xi^2)) ;
        -(9/16)*(1-eta)*(1-3*xi)*xi-(27/32)*(1-eta)*(1-xi^2),-(9/32)*(1-3*xi)*(1-xi^2);
        -(9/16)*(1-eta)*xi*(1+3*xi)+(27/32)*(1-eta)*(1-xi^2),-(9/32)*(1+3*xi)*(1-xi^2);
        (9/32)*(1-3*eta)*(1-eta^2),-(9/16)*(1-3*eta)*eta*(1+xi)-(27/32)*(1-eta^2)*(1+xi);
        (9/32)*(1+3*eta)*(1-eta^2),-(9/16)*eta*(1+3*eta)*(1+xi)+(27/32)*(1-eta^2)*(1+xi);
        -(9/16)*(1+eta)*xi*(1+3*xi)+(27/32)*(1+eta)*(1-xi^2),(9/32)*(1+3*xi)*(1-xi^2);
        -(9/16)*(1+eta)*(1-3*xi)*xi-(27/32)*(1+eta)*(1-xi^2),(9/32)*(1-3*xi)*(1-xi^2);
        -(9/32)*(1+3*eta)*(1-eta^2),-(9/16)*eta*(1+3*eta)*(1-xi)+(27/32)*(1-eta^2)*(1-xi);
        -(9/32)*(1-3*eta)*(1-eta^2),-(9/16)*(1-3*eta)*eta*(1-xi)-(27/32)*(1-eta^2)*(1-xi)];

    end % end function shapeFunctionQ12

repo: Aero-matlab


text:
function [llcfeat, x, y, wid, hgt] = llc_hog2x2(img, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~exist('c', 'var'))
	c = conf();
end

featname = 'hog2x2';
p = c.feature_config.(featname);
if(isfield(p, 'dictionary'))
  dictionary = p.dictionary;
elseif(isfield(p, 'dictionary_file'))
  try
    tmp = load(p.dictionary_file);
    dictionary = tmp.dictionary;
  catch
    error('No dictionary found!');
  end
else
  error('Must specify dictionary!');
end

[feat, x, y, wid, hgt] = extract_hog2x2(img, c);
llcfeat = sparse(LLC_coding_appr(dictionary, feat, p.llcknn));

author: adikhosla

repo: feature-extraction


text:
classdef lsl_outlet < handle
	% A stream outlet.
	% Outlets are used to make streaming data (and the meta-data) available on the lab network.

    properties (Hidden)
        LibHandle = 0;       % this is a handle to the liblsl library (to call its functions)
        OutletHandle = 0;    % this is a handle to an lsl_outlet object within the library.
        SampleRate = 0;      % copy of the sample rate...
    end

    methods

        function self = lsl_outlet(info, chunksize, maxbuffered)
            % Outlet = lsl_outlet(Streaminfo, ChunkSize, MaxBuffered)
            % Establish a new stream outlet. This makes the stream discoverable.
            %
            % In:
            %   Streaminfo : The lsl_streaminfo object to use for creating this stream. Stays constant over the lifetime of the outlet.
            %
            %   ChunkSize :  Optionally the desired chunk granularity (in samples) for transmission. If 0,
            %                each push operation yields one chunk. Stream recipients can have this setting bypassed.
            %                (default: 0)
            %
            %   MaxBuffered : Optionally the maximum amount of data to buffer (in seconds if there is a nominal
            %                 sampling rate, otherwise x100 in samples). Stream recipients can elect to use
            %                 a smaller (but not larger) buffer size than the one specified here.
            %                 (default: 360)
            %

            if ~exist('chunksize','var') || isempty(chunksize) chunksize = 0; end
            if ~exist('maxbuffered','var') || isempty(maxbuffered) maxbuffered = 360; end
            self.LibHandle = info.LibHandle;
            self.SampleRate = info.nominal_srate();
            self.OutletHandle = lsl_create_outlet(info.LibHandle,info.InfoHandle,chunksize,maxbuffered);
        end

        function delete(self)
            % Destroy the outlet when it is deleted.
            % The outlet will no longer be discoverable after destruction and all paired inlets will
            % stop delivering data.

            lsl_destroy_outlet(self.LibHandle,self.OutletHandle);
        end

        function push_sample(self,sampledata,timestamp,pushthrough)
            % push_sample(SampleData,Timestamp,Pushthrough)
            % Push a sample into the outlet.
            %
            % In:
            %   SampleData : either a numeric vector (as double, converted implicitly) with as many
            %                elements as the stream has channels, or a cell array of strings of the
            %                same size (if the stream is string-formatted).
            %
            %   Timestamp : Optionally the capture time of the most recent sample, in agreement with lsl_local_clock();
            %               if 0, the current time is used.
            %
            %   Pushthrough : Whether to push the chunk through to the receivers instead of buffering it with subsequent samples.
			%	 			  Note that the chunk_size, if specified at outlet construction, takes precedence over the pushthrough flag.
            %                 (default: 1)

            if ~exist('timestamp','var') || isempty(timestamp) timestamp = 0; end
            if ~exist('pushthrough','var') || isempty(pushthrough) pushthrough = 1; end
            lsl_push_sample(self.LibHandle,self.OutletHandle,sampledata,timestamp,pushthrough);
        end

        function push_chunk(self,chunkdata,timestamps,pushthrough)
            % push_sample(ChunkData,Timestamps,Pushthrough)
            % Push a chunk into the outlet.
            %
            % In:
            %   ChunkData : Numeric matrix with one column for each sample that shall be
            %               transmitted. May also be empty.
            %
            %   Timestamps : Optionally the capture time of the most recent sample, in agreement with lsl_local_clock();
            %                if 0, the current time is used. The time stamps of other samples are automatically
            %                derived based on the sampling rate of the stream, or replicated if the stream is irregular-rate
            %                If this is a vector, it is assumed that it contains one time stamp for each sample. (default: 0)
            %
            %   Pushthrough : Whether to push the chunk through to the receivers instead of buffering it with subsequent samples.
			%	 			  Note that the chunk_size, if specified at outlet construction, takes precedence over the pushthrough flag.
            %                 (default: 1)

            if ~exist('timestamps','var') || isempty(timestamps) timestamps = 0; end
            if ~exist('pushthrough','var') || isempty(pushthrough) pushthrough = 1; end

            if timestamps == 0
                timestamps = lsl_local_clock(self.LibHandle); end

            if isscalar(timestamps)
                if self.SampleRate == 0
                    timestamps = [timestamps; -ones(size(chunkdata,2)-1,1)];
                else
                    % otherwise infer the timestamp of the first sample in the chunk and let the
                    % library deduce the timestamps for the remaining samples (using -1.0 as the
                    % stamps).
                    timestamps = [timestamps - (size(chunkdata,2)-1)/self.SampleRate; -ones(size(chunkdata,2)-1,1)];
                end
            end

            lsl_push_chunk(self.LibHandle,self.OutletHandle,chunkdata,timestamps,pushthrough);
        end

        % === Miscellaneous Features ===

        function result = have_consumers(self)
            % Check whether consumers are currently registered.
            % While it does not hurt, there is technically no reason to push samples if there is no consumer.

            result = lsl_have_consumers(self.LibHandle,self.OutletHandle);
        end

        function success = wait_for_consumers(self,timeout)
            % Wait until some consumer shows up (without wasting resources).
            %
            % In:
            %   Timeout : Timeout for the operation. (default: 60)
            %
            % Out:
            %   Success : True if the wait was successful, false if the timeout expired.

            if ~exist('timeout','var') || isempty(timeout) timeout = 60; end
            success = lsl_wait_for_consumers(self.LibHandle,self.OutletHandle,timeout);
        end

        function result = info(self)
            % Retrieve the stream info provided by this outlet.
            % This is what was used to create the stream (and also has the Additional Network Information fields assigned).

            result = lsl_streaminfo(self.LibHandle,lsl_get_info(self.LibHandle,self.OutletHandle));
        end

        function h = get_libhandle(self)
            % get the library handle (e.g., to query the clock)
            h = self.LibHandle;
        end
    end
end

author: PatternRecognition

repo: OpenBMI


text:
function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileName, headerRow, formulaFlag, gprFlag)
% Prints a flux vector with reaction labels
%
% USAGE:
%    printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileName, headerRow, formulaFlag, gprFlag)
%
% INPUTS:
%    model:          COBRA model structure
%    fluxData:       Data matrix/vector (for example, solution.v)
%
% OPTIONAL INPUTS:
%    nonZeroFlag:    Only print nonzero rows (Default = false)
%    excFlag:        Only print exchange fluxes (Default = false)
%    sortCol:        Column used for sorting (-1, none; 0, labels; >0, data
%                    columns;) (Default = -1)
%    fileName:       Name of output file (Default = [])
%    headerRow:      Header (Default = [])
%    formulaFlag:    Print reaction formulas (Default = false)
%    gprFlag:        Print reaction GPR (Default = false)
%

% .. Authors:
%       - Markus Herrgard, Ronan Fleming

if isempty(fluxData)
    return
end
if nargin < 3
    nonZeroFlag = false;
end
if nargin < 4
    excFlag = false;
end
if nargin < 5
    sortCol = -1;
end
if nargin < 6
    fileName = [];
end
if nargin < 7
    headerRow = [];
end
if nargin < 8
    formulaFlag = false;
end
if nargin < 9
    gprFlag = false;
end

if excFlag
    bool = findExcRxns(model, true, false);
else
    bool = true(size(model.S,2),1);
    labels = model.rxns;
end

if length(nonZeroFlag)==size(model.S,2)
    bool = bool & nonZeroFlag;
    nonZeroFlag = 1;
else
    if nonZeroFlag
        bool = bool & fluxData(:,1)~=0;
    end
end

% Add reaction formulas
if formulaFlag
    if nonZeroFlag
        labels = [model.rxns,model.rxns];
        %only generate the formulas for the nonzero entries
        formulas = printRxnFormula(model, labels(bool), false, false);
        labels(bool,2) = formulas;
    else
        formulas = printRxnFormula(model, labels, false, false);
        labels = [labels, formulas];
    end
end

% Add GPR
if gprFlag
    if nonZeroFlag
        %only generate the gprs for the nonzero entries
        labels(bool,end+1) = model.grRules(bool);
    else
        labels = [labels, model.grRules];
    end
end

%only print the nonzeros
if nonZeroFlag
    labels = labels(bool,:);
    fluxData = fluxData(bool,:);
end

labels(:,2) = pad(labels(:,2));

%print the labeled data
printLabeledData(labels, fluxData, 0, sortCol, fileName, headerRow)


text:
function nh=cosmo_singleton_neighborhood(ds)
% return neighborhood where each feature is only neighbor of itself
%
% nh=cosmo_singleton_neighborhood(ds)
%
% Input:
%   ds                      dataset strucut
%
% Output:
%   nh                      neighborhood structure where each feature is
%                           only neighbor of itself.
%
%
% Example:
%     % This example shows how to run group analysis with multiple
%     % comparison correction using a clustering neighborhood defined
%     % by cosmo_singleton_neighborhood with cosmo_montecarlo_cluster_stat
%     %
%     nsubjects=10;
%     nrois=4;
%     ds=struct();
%     ds.a=struct();
%     ds.fa=struct();
%     % set random data for 10 participants and 4 ROIs
%     ds.samples=randn(nsubjects,nrois);
%     %
%     % set chunks and targets for one-sample t-test
%     ds.sa.chunks=(1:nsubjects)';
%     ds.sa.targets=ones(nsubjects,1);
%     %
%     nh=cosmo_singleton_neighborhood(ds);
%     cosmo_disp(nh);
%     %|| .origin
%     %||   .a
%     %||     struct (empty)
%     %||   .fa
%     %||     struct (empty)
%     %|| .fa
%     %||   .sizes
%     %||     [ 1         1         1         1 ]
%     %|| .a
%     %||   struct (empty)
%     %|| .neighbors
%     %||   { [ 1 ]
%     %||     [ 2 ]
%     %||     [ 3 ]
%     %||     [ 4 ] }
%     %||
%     %
%     opt=struct();
%     opt.progress=false;
%     %
%     % t-test against mean=0
%     opt.h0_mean=0;
%     %
%     % make this a fast example.
%     % usually one uses opt.niter=1000; even better is opt.niter=10000
%     opt.niter=10; % use 10000 for publication quality
%     %
%     tfce_z_scores=cosmo_montecarlo_cluster_stat(ds,nh,opt);
%     %
%     % output contains one z-score per ROI
%     cosmo_disp(size(tfce_z_scores.samples))
%     %|| [1, 4]
%
%
% Notes:
% - this function can be used for ROI group analysis with multiple
%   comparison correction using cosmo_montecarlo_cluster_stat.
%   To do so, each column in
%     ds.samples
%   should correspond to values in one ROI.
%

    cosmo_check_dataset(ds);

    nh=struct();
    nh.origin=struct();
    nh.fa=struct();
    nh.a=struct();

    if isfield(ds,'a')
        a=ds.a;

        nh.a=a;
        nh.origin.a=a;
    end

    if isfield(ds,'fa')
        fa=ds.fa;

        nh.fa=fa;
        nh.origin.fa=fa;
    end

    nfeatures=size(ds.samples,2);
    nh.fa.sizes=ones(1,nfeatures);

    nh.neighbors=num2cell((1:nfeatures)');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function test_bug542

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_multiplotER

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug542.mat'));

% this is confusing the plotting function
att_dep_ipsi = rmfield(att_dep_ipsi, 'time');

stopwatch = tic;
ft_multiplotER(cfg, att_dep_ipsi);
toc(stopwatch)


text:
function value = r8_binom ( n, m )

%*****************************************************************************80
%
%% R8_BINOM evaluates the binomial coefficient using R8 arithmetic.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, integer N, M, the arguments.
%
%    Output, real VALUE, the binomial coefficient.
%
  persistent bilnmx
  persistent fintmx
  persistent sq2pil

  sq2pil = 0.91893853320467274178032973640562;

  if ( isempty ( bilnmx ) )
    bilnmx = log ( r8_mach ( 2 ) ) - 0.0001;
    fintmx = 0.9 / r8_mach ( 3 );
  end

  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BINOM - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'R8_BINOM - Fatal error!' )
  end

  if ( m < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BINOM - Fatal error!\n' );
    fprintf ( 1, '  M < 0.\n' );
    error ( 'R8_BINOM - Fatal error!' )
  end

  if ( n < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BINOM - Fatal error!\n' );
    fprintf ( 1, '  N < M.\n' );
    error ( 'R8_BINOM - Fatal error!' )
  end

  k = min ( m, n - m );

  if ( k <= 20 && k * log ( max ( n, 1 ) ) <= bilnmx )

    value = 1.0;

    for i = 1 : k
      value = value * ( n - i + 1 ) / i;
    end

  else

    if ( k < 9 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8_BINOM - Fatal error!\n' );
      fprintf ( 1, '  Result overflows.\n' );
      fprintf ( 1, '  N or M is too big.\n' );
      error ( 'R8_BINOM - Fatal error!' )
    end

    xn =  n + 1;
    xk =  k + 1;
    xnk =  n - k + 1;

    corr = r8_lgmc ( xn ) - r8_lgmc ( xk ) - r8_lgmc ( xnk );

    value = xk * log ( xnk / xk ) ...
      - xn * r8_lnrel ( - ( xk - 1.0 ) / xn ) ...
      - 0.5 * log ( xn * xnk / xk ) + 1.0 - sq2pil + corr;

    if ( bilnmx < value )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8_BINOM - Fatal error!\n' );
      fprintf ( 1, '  Result overflows.\n' );
      fprintf ( 1, '  N or M is too big.\n' );
      error ( 'R8_BINOM - Fatal error!' )
    end

    value = exp ( value );

  end

  if ( value < fintmx )
    value = r8_aint ( value + 0.5 );
  end

  return
end


text:
function [err, C] = gen_pos_err(std, n, type, absolute)
%GEN_POS_ERR Generates position errors. We first generate i.i.d. position errors
%for each sensor. We use the first sensor as the reference sensor and substract
%the the position error of the first sensor from the remaining sensors.
%Therefore, the resulting covariance matrix is not diagonal.
%Inputs:
%   std - Standard deviation of the position errors.
%   n - Number of sensors.
%   type - 'Gaussian' or 'Uniform'.
%   absolute - If set to true, will not use the first sensor as the reference
%              sensor and the position errors will be i.i.d.
%Outputs:
%   err - 2xn vector of position errors.
%   C - Covariance matrix. (2n-2) x (2n-2) if absolute is false, (2n) x (2n) if
%       absolute is true.
if nargin < 4
    absolute = false;
end
switch lower(type)
    case 'gaussian'
        err = randn(2, n)*std;
        if nargout > 1
            if absolute
                C = std^2*eye(n + n);
            else
                C = std^2*(eye(n - 1) + ones(n - 1));
                C = blkdiag(C, C);
            end
        end
    case 'uniform'
        b = sqrt(3)*std;
        err = unifrnd(-b, b, 2, n);
        if nargout > 1
            if absolute
                C = std^2*eye(n + n);
            else
                C = std^2*(eye(n - 1) + ones(n - 1));
                C = blkdiag(C, C);
            end
        end
    otherwise
        error('Unknow type "%s".', type);
end
if ~absolute
    err(1,2:end) = err(1,2:end) - err(1,1);
    err(2,2:end) = err(2,2:end) - err(2,1);
    err(1,1) = 0;
    err(2,1) = 0;
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function [dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St] = dSbr_dV(branch, Yf, Yt, V, vcart)
%DSBR_DV   Computes partial derivatives of branch power flows w.r.t. voltage.
%
%   The derivatives can be taken with respect to polar or cartesian coordinates
%   of voltage, depending on the 5th argument.
%
%   [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V)
%   [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 0)
%
%   Returns four matrices containing partial derivatives of the complex
%   branch power flows at "from" and "to" ends of each branch w.r.t voltage
%   magnitude and voltage angle, respectively (for all buses).
%
%   [DSF_DVR, DSF_DVI, DST_DVR, DST_DVI, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 1)
%
%   Returns four matrices containing partial derivatives of the complex
%   branch power flows at "from" and "to" ends of each branch w.r.t real and
%   imaginary parts of voltage, respectively (for all buses).
%
%   If YF is a sparse matrix, the partial derivative matrices will be as well.
%   Optionally returns vectors containing the power flows themselves. The
%   following explains the expressions used to form the matrices:
%
%   If = Yf * V;
%   Sf = diag(Vf) * conj(If) = diag(conj(If)) * Vf
%
%   Polar coordinates:
%     Partials of V, Vf & If w.r.t. voltage angles
%       dV/dVa  = j * diag(V)
%       dVf/dVa = sparse(1:nl, f, j * V(f)) = j * sparse(1:nl, f, V(f))
%       dIf/dVa = Yf * dV/dVa = Yf * j * diag(V)
%
%     Partials of V, Vf & If w.r.t. voltage magnitudes
%       dV/dVm  = diag(V./abs(V))
%       dVf/dVm = sparse(1:nl, f, V(f)./abs(V(f))
%       dIf/dVm = Yf * dV/dVm = Yf * diag(V./abs(V))
%
%     Partials of Sf w.r.t. voltage angles
%       dSf/dVa = diag(Vf) * conj(dIf/dVa)
%                       + diag(conj(If)) * dVf/dVa
%               = diag(Vf) * conj(Yf * j * diag(V))
%                       + conj(diag(If)) * j * sparse(1:nl, f, V(f))
%               = -j * diag(Vf) * conj(Yf * diag(V))
%                       + j * conj(diag(If)) * sparse(1:nl, f, V(f))
%               = j * (conj(diag(If)) * sparse(1:nl, f, V(f))
%                       - diag(Vf) * conj(Yf * diag(V)))
%
%     Partials of Sf w.r.t. voltage magnitudes
%       dSf/dVm = diag(Vf) * conj(dIf/dVm)
%                       + diag(conj(If)) * dVf/dVm
%               = diag(Vf) * conj(Yf * diag(V./abs(V)))
%                       + conj(diag(If)) * sparse(1:nl, f, V(f)./abs(V(f)))
%
%   Cartesian coordinates:
%     Partials of V, Vf & If w.r.t. real part of complex voltage
%       dV/dVr  = diag(ones(n,1))
%       dVf/dVr = Cf
%       dIf/dVr = Yf
%     where Cf is the connection matrix for line & from buses
%
%     Partials of V, Vf & If w.r.t. imaginary part of complex voltage
%       dV/dVi  = j * diag(ones(n,1))
%       dVf/dVi = j * Cf
%       dIf/dVi = j * Yf
%
%     Partials of Sf w.r.t. real part of complex voltage
%       dSf/dVr = conj(diag(If)) * Cf + diag(Vf) * conj(Yf)
%
%     Partials of Sf w.r.t. imaginary part of complex voltage
%       dSf/dVi = j * (conj(diag(If)) * Cf - diag(Vf) * conj(Yf))
%
%   Derivations for "to" bus are similar.
%
%   Examples:
%       [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
%       [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = ...
%           dSbr_dV(branch, Yf, Yt, V);
%       [dSf_dVr, dSf_dVi, dSt_dVr, dSt_dVi, Sf, St] = ...
%           dSbr_dV(branch, Yf, Yt, V, 1);
%
%   For more details on the derivations behind the derivative code used
%   in MATPOWER, see:
%
%   [TN2]  R. D. Zimmerman, "AC Power Flows, Generalized OPF Costs and
%          their Derivatives using Complex Matrix Notation", MATPOWER
%          Technical Note 2, February 2010. [Online]. Available:
%          https://matpower.org/docs/TN2-OPF-Derivatives.pdf
%          doi: 10.5281/zenodo.3237866
%   [TN4]  B. Sereeter and R. D. Zimmerman, "AC Power Flows and their
%          Derivatives using Complex Matrix Notation and Cartesian
%          Coordinate Voltages," MATPOWER Technical Note 4, April 2018.
%          [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf
%          doi: 10.5281/zenodo.3237909

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

%% define named indices into bus, gen, branch matrices
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% default input args
if nargin < 5
    vcart = 0;      %% default to polar coordinates
end

%% define
f = branch(:, F_BUS);       %% list of "from" buses
t = branch(:, T_BUS);       %% list of "to" buses
nl = length(f);
nb = length(V);

%% compute intermediate values
Yfc = conj(Yf);
Ytc = conj(Yt);
Vc = conj(V);
Ifc = Yfc * Vc;     %% conjugate of "from" current
Itc = Ytc * Vc;     %% conjugate of "to" current

if issparse(Yf)             %% sparse version (if Yf is sparse)
    diagVf  = sparse(1:nl, 1:nl, V(f), nl, nl);
    diagVt  = sparse(1:nl, 1:nl, V(t), nl, nl);
    diagIfc = sparse(1:nl, 1:nl, Ifc, nl, nl);
    diagItc = sparse(1:nl, 1:nl, Itc, nl, nl);
    if ~vcart
        Vnorm       = V ./ abs(V);
        diagVc      = sparse(1:nb, 1:nb, Vc, nb, nb);
        diagVnorm   = sparse(1:nb, 1:nb, Vnorm, nb, nb);
        CVf  = sparse(1:nl, f, V(f), nl, nb);
        CVnf = sparse(1:nl, f, Vnorm(f), nl, nb);
        CVt  = sparse(1:nl, t, V(t), nl, nb);
        CVnt = sparse(1:nl, t, Vnorm(t), nl, nb);
    end
else                        %% dense version
    diagVf  = diag(V(f));
    diagVt  = diag(V(t));
    diagIfc = diag(Ifc);
    diagItc = diag(Itc);
    if ~vcart
        Vnorm       = V ./ abs(V);
        diagVc      = diag(Vc);
        diagVnorm   = diag(Vnorm);
%         CVf        = zeros(nl, nb);    CVf(sub2ind([nl,nb], (1:nl)', f)) = V(f);
%         CVnf       = zeros(nl, nb);    CVnf(sub2ind([nl,nb], (1:nl)', f)) = Vnorm(f);
%         CVt        = zeros(nl, nb);    CVt(sub2ind([nl,nb], (1:nl)', t)) = V(t);
%         CVnt       = zeros(nl, nb);    CVnt(sub2ind([nl,nb], (1:nl)', t)) = Vnorm(t);
        CVf  = full(sparse(1:nl, f, V(f), nl, nb));
        CVnf = full(sparse(1:nl, f, Vnorm(f), nl, nb));
        CVt  = full(sparse(1:nl, t, V(t), nl, nb));
        CVnt = full(sparse(1:nl, t, Vnorm(t), nl, nb));
    end
end
if vcart
    Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
    Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
    Af = diagIfc * Cf;
    Bf = diagVf * Yfc;
    At = diagItc * Ct;
    Bt = diagVt * Ytc;

    dSf_dV1 = Af + Bf;          %% dSf_dVr
    dSf_dV2 = 1j * (Af - Bf);   %% dSf_dVi
    dSt_dV1 = At + Bt;          %% dSt_dVr
    dSt_dV2 = 1j * (At - Bt);   %% dSt_dVi
else
    dSf_dV1 = 1j * (diagIfc * CVf - diagVf * Yfc * diagVc);     %% dSf_dVa
    dSf_dV2 = diagVf * conj(Yf * diagVnorm) + diagIfc * CVnf;   %% dSf_dVm
    dSt_dV1 = 1j * (diagItc * CVt - diagVt * Ytc * diagVc);     %% dSt_dVa
    dSt_dV2 = diagVt * conj(Yt * diagVnorm) + diagItc * CVnt;   %% dSt_dVm
end

if nargout > 4
    Sf = V(f) .* Ifc;
    St = V(t) .* Itc;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [marginal, msg, loglik] = filter_evidence_old(engine, evidence)
% [marginal, msg, loglik] = filter_evidence(engine, evidence) (pearl_dbn)

[ss T] = size(evidence);
bnet = bnet_from_engine(engine);
bnet2 = dbn_to_bnet(bnet, T);
ns = bnet2.node_sizes;
hnodes = mysetdiff(1:ss, engine.onodes);
hnodes = hnodes(:)';

[engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2);

msg = init_msgs(bnet2.dag, ns, evidence);
msg = init_ev_msgs(engine, evidence, msg);

verbose = 1;
if verbose, fprintf('\nold filtering\n'); end

for t=1:T
  % update pi
  for i=hnodes
    n = i + (t-1)*ss;
    ps = parents(bnet2.dag, n);
    if t==1
      e = bnet.equiv_class(i,1);
    else
      e = bnet.equiv_class(i,2);
    end
    msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg);
    %if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end
    msg{n}.pi = normalise(msg{n}.pi(:) .* msg{n}.lambda_from_self(:));
    if verbose, fprintf('%d recomputes pi\n', n); disp(msg{n}.pi); end
  end
  % send pi msg to children
  for i=hnodes
    n = i + (t-1)*ss;
    cs = children(bnet2.dag, n);
    for c=cs(:)'
      j = engine.parent_index{c}(n); % n is c's j'th parent
      pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns));
      msg{c}.pi_from_parent{j} = pi_msg;
      if verbose, fprintf('%d sends pi to %d\n', n,c); disp(pi_msg); end
    end
  end
end

marginal = cell(ss,T);
lik = zeros(1,ss*T);
for t=1:T
  for i=1:ss
    n = i + (t-1)*ss;
    %[bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda);
    [bel, lik(n)] = normalise(msg{n}.pi);
    marginal{i,t} = bel;
  end
end

loglik = sum(log(lik));

%%%%%%%

function lambda = compute_lambda(n, cs, msg, ns)
% Pearl p183 eq 4.50
lambda = prod_lambda_msgs(n, cs, msg, ns);

%%%%%%%

function pi_msg = compute_pi_msg(n, cs, msg, c, ns)
% Pearl p183 eq 4.53 and 4.51
pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c);

%%%%%%%%%

function lam = prod_lambda_msgs(n, cs, msg, ns, except)

if nargin < 5, except = -1; end

%lam = msg{n}.lambda_from_self(:);
lam = ones(ns(n), 1);
for i=1:length(cs)
  c = cs(i);
  if c ~= except
    lam = lam .* msg{n}.lambda_from_child{i};
  end
end

%%%%%%%%%%%

function msg = init_msgs(dag, ns, evidence)
% INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn)
% msg =  init_msgs(dag, ns, evidence)
%
% We assume all the hidden nodes are discrete.

N = length(dag);
msg = cell(1,N);
observed = ~isemptycell(evidence(:));

for n=1:N
  ps = parents(dag, n);
  msg{n}.pi_from_parent = cell(1, length(ps));
  for i=1:length(ps)
    p = ps(i);
    msg{n}.pi_from_parent{i} = ones(ns(p), 1);
  end

  cs = children(dag, n);
  msg{n}.lambda_from_child = cell(1, length(cs));
  for i=1:length(cs)
    c = cs(i);
    msg{n}.lambda_from_child{i} = ones(ns(n), 1);
  end

  msg{n}.lambda = ones(ns(n), 1);
  msg{n}.pi = ones(ns(n), 1);

  msg{n}.lambda_from_self = ones(ns(n), 1);
end

%%%%%%%%%

function msg = init_ev_msgs(engine, evidence, msg)
% Initialize the lambdas with any evidence

[ss T] = size(evidence);
bnet = bnet_from_engine(engine);
pot_type = 'd';
t = 1;
hnodes = mysetdiff(1:ss, engine.onodes);
for i=hnodes(:)'
  c = engine.obschild(i);
  if c > 0
    fam = family(bnet.dag, c);
    e = bnet.equiv_class(c, 1);
    CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1));
    temp = pot_to_marginal(CPDpot);
    n = i;
    msg{n}.lambda_from_self = temp.T;
  end
end
for t=2:T
  for i=hnodes(:)'
    c = engine.obschild(i);
    if c > 0
      fam = family(bnet.dag, c, 2);
      e = bnet.equiv_class(c, 2);
      CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t));
      temp = pot_to_marginal(CPDpot);
      n = i + (t-1)*ss;
      msg{n}.lambda_from_self = temp.T;
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% MAIN Illustrates how to use the EPnP algorithm described in:
%
%       Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua.
%       Accurate Non-Iterative O(n) Solution to the PnP Problem.
%       In Proceedings of ICCV, 2007.
%
% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

clear all; close all;

addpath data;
addpath error;
addpath EPnP;

fprintf('\n---------EPnP--------------\n');
%1.-Generate simulated input data------------------------------------------
load_points=0;
if ~load_points
    n=50; %number of points
    std_noise=10; %noise in the measurements (in pixels)
    [A,point,Rt]=generate_noisy_input_data(n,std_noise);
    save('data\input_data_noise.mat','A','point','Rt');
else
    load('data\input_data_noise.mat','A','point','Rt');
    n=size(point,2);
    draw_noisy_input_data(point);
end

%2.-Inputs format--------------------------------
x3d=zeros(n,4);
x2d=zeros(n,3);
A=A(:,1:3);
for i=1:n
    x3d_h(i,:)=[point(i).Xworld',1];
    x2d_h(i,:)=[point(i).Ximg(1:2)',1];

    %world and camera coordinates
    X3d_world(i,:)=point(i).Xworld';
    X3d_cam(i,:)=point(i).Xcam';
end

%3.-EPnP----------------------------------------------------
Xw=x3d_h(:,1:3);
U=x2d_h(:,1:2);

[Rp,Tp,Xc,sol]=efficient_pnp(x3d_h,x2d_h,A);

%draw Results
for i=1:n
    point(i).Xcam_est=Xc(i,:)';
end
figure; h=gcf;
plot_3d_reconstruction(point,'EPnP (Old)',h);
xlim([-2 2]); ylim([-2 2]);

%compute error
error=reprojection_error_usingRT(Xw,U,Rp,Tp,A);
fprintf('error EPnP: %.3f\n',error);

%3.-EPnP_GAUSS_NEWTON----------------------------------------------------
Xw=x3d_h(:,1:3);
U=x2d_h(:,1:2);

[Rp,Tp,Xc,sol]=efficient_pnp_gauss(x3d_h,x2d_h,A);

%draw Results
for i=1:n
    point(i).Xcam_est=Xc(i,:)';
end
figure; h=gcf;
plot_3d_reconstruction(point,'EPnP Gauss Newton',h);

%compute error
error=reprojection_error_usingRT(Xw,U,Rp,Tp,A);
fprintf('error EPnP_Gauss_Newton: %.3f\n',error);
xlim([-2 2]); ylim([-2 2]);

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function i4_log_2_test ( )

%*****************************************************************************80
%
%% I4_LOG_2_TEST tests I4_LOG_2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 17;

  x_test = [ 0, 1, 2, 3, 9, 10, 11, 99, 101, -1, -2, ...
    -3, -9, 1000, 1023, 1024, 1025 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_LOG_2_TEST\n' );
  fprintf ( 1, '  I4_LOG_2: whole part of log base 2.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       X     I4_LOG_2\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    x = x_test(test);
    fprintf ( 1, '  %6d  %12d\n', x, i4_log_2 ( x ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function U=NTds(U,X,Y,xflux,yflux,N,T,varargin)
%%
%   Performs dimensional splitting on the 2D scalar conservation law
%   U_t + xflux(U)_x + yflux(U)_y = 0 with initial data U.
%   Uses Strang splitting with timestep dt=T/N. The domain is the square
%   spanned by X and Y. varargin may contain 'wbar' if waitbar is desired,
%   'periodic' followd by [xmin xmax], [ymin ymax]. The one-dimensional
%   equations are solved using a second order central scheme.
%

%%
% Initial setup
addpath ../../AppendixA;
wb=0;
nv=length(varargin);
i=1;
boundary_cond='neumann';
while i<=nv,
	if strcmp(varargin(i),'wbar'),
		wb=1;
		i=i+1;
	elseif strcmp(varargin(i),'periodic'),
		boundary_cond='periodic';
		i=i+1;
	else
      display('Unknown option in NTds');
      i=i+1;
	end;
end;
dt=T/N;
xdir=2; ydir=1;
dx=X(2)-X(1);
dy=Y(2)-Y(1);
CFL=0.5;
limiter='superbee';
strang=dt*0.5;

%%
% The main loop
if wb
  wstr=sprintf('Computing %d steps.',N);
  wbar=waitbar(0,wstr);
end;
for k=1:N,
  U=central(xflux,U,strang,dx,xdir,boundary_cond,CFL,limiter);
  if wb,
    waitbar((k-0.5)/N,wbar);
  end;
  if (k==1),
    strang=dt;
  end;
  U=central(yflux,U,strang,dy,ydir,boundary_cond,CFL,limiter);
  if k==N,
    strang=0.5*dt;
    U=central(xflux,U,strang,dx,xdir,boundary_cond,CFL,limiter);
  end;
  if wb,
    waitbar(k/N,wbar);
  end;
end;

%%
%
if wb,
  close(wbar)
end;
rmpath ../../AppendixA;


text:
% csv2cell() - Imports the contents of a character delimited spreadsheet
%              text file as a 2D cell array, csv_matrix.
%
% Usage:
%  >> csv_matrix=csv2cell(csv_fname,delimiter);
%
% Required Input:
%   csv_fname - The name of the text file to import.  If a path is not
%               included in the name, the current working directory will
%               be used.
%
% Optional Input:
%   delimiter - The character that marks cell boundaries in the text file.
%               {default: ','}
%
% Output:
%   csv_matrix - 2D cell array of strings. csv_matrix{x,y} corresponds to
%                the cell in the xth row and yth column in the spreadsheet.
%
% Note, be sure not to use the delimiting character in the text file as
% anything but a cell boundary marker.
%
% Example:
%  >> demo_dat=round(rand(5,3)*10);
%  >> fid=fopen('demo.csv','w');
%  >> for a=1:5, for b=1:3, fprintf(fid,'%d,',demo_dat(a,b)); end; fprintf(fid,'\n'); end;
%  >> fclose(fid);
%  >> csv_matrix=csv2cell('demo.csv',',');
%
% Author:
%  David Groppe
%  9/3/2011
%  Mehtalab
%  Cushing Neuroscience Institutes
%  New Hyde Park, New York

function csv_matrix=csv2cell(csv_fname,delimiter)

if nargin<1,
    help csv2cell
end

if nargin<2,
    delimiter=',';
else
    if ~ischar(delimiter)
        error('Specified delimiter needs to be a character.');
    elseif length(delimiter)>1
        error('Delimiter needs to be a single character.');
    end
end

[fid, msg]=fopen(csv_fname,'r');
if fid==-1,
   error('Cannot open %s because: %s.\n',csv_fname,msg);
end

csv_matrix=cell(1,1);
row_ct=1;
col_ct=1;
while 1
    tline = fgetl(fid);
    if ~ischar(tline), break, end

    while ~isempty(tline),
        [t, tline]=parse_by_char(tline,delimiter);
        csv_matrix{row_ct,col_ct}=t;
        col_ct=col_ct+1;
    end
    col_ct=1;
    row_ct=row_ct+1;
end
fclose(fid);

function [pre, post]=parse_by_char(str,delimiter)

char_ids=find(str==delimiter);
if isempty(char_ids),
    pre=str;
    post=[];
else
    pre=str(1:char_ids(1)-1);
    post=str(char_ids(1)+1:end);
end

repo: mcbench-benchmarks