text:
%  Decode extra NIFTI header information into hdr.extra
%
%  Usage: hdr = extra_nii_hdr(hdr)
%
%  hdr can be obtained from load_nii_hdr
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function hdr = extra_nii_hdr(hdr)

   switch hdr.dime.datatype
   case 1
      extra.NIFTI_DATATYPES = 'DT_BINARY';
   case 2
      extra.NIFTI_DATATYPES = 'DT_UINT8';
   case 4
      extra.NIFTI_DATATYPES = 'DT_INT16';
   case 8
      extra.NIFTI_DATATYPES = 'DT_INT32';
   case 16
      extra.NIFTI_DATATYPES = 'DT_FLOAT32';
   case 32
      extra.NIFTI_DATATYPES = 'DT_COMPLEX64';
   case 64
      extra.NIFTI_DATATYPES = 'DT_FLOAT64';
   case 128
      extra.NIFTI_DATATYPES = 'DT_RGB24';
   case 256
      extra.NIFTI_DATATYPES = 'DT_INT8';
   case 512
      extra.NIFTI_DATATYPES = 'DT_UINT16';
   case 768
      extra.NIFTI_DATATYPES = 'DT_UINT32';
   case 1024
      extra.NIFTI_DATATYPES = 'DT_INT64';
   case 1280
      extra.NIFTI_DATATYPES = 'DT_UINT64';
   case 1536
      extra.NIFTI_DATATYPES = 'DT_FLOAT128';
   case 1792
      extra.NIFTI_DATATYPES = 'DT_COMPLEX128';
   case 2048
      extra.NIFTI_DATATYPES = 'DT_COMPLEX256';
   otherwise
      extra.NIFTI_DATATYPES = 'DT_UNKNOWN';
   end

   switch hdr.dime.intent_code
   case 2
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CORREL';
   case 3
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST';
   case 4
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST';
   case 5
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ZSCORE';
   case 6
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ';
   case 7
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BETA';
   case 8
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BINOM';
   case 9
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GAMMA';
   case 10
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POISSON';
   case 11
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NORMAL';
   case 12
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST_NONC';
   case 13
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ_NONC';
   case 14
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGISTIC';
   case 15
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LAPLACE';
   case 16
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_UNIFORM';
   case 17
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST_NONC';
   case 18
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_WEIBULL';
   case 19
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHI';
   case 20
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_INVGAUSS';
   case 21
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_EXTVAL';
   case 22
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_PVAL';
   case 23
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGPVAL';
   case 24
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOG10PVAL';
   case 1001
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ESTIMATE';
   case 1002
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LABEL';
   case 1003
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NEURONAME';
   case 1004
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GENMATRIX';
   case 1005
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_SYMMATRIX';
   case 1006
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DISPVECT';
   case 1007
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_VECTOR';
   case 1008
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POINTSET';
   case 1009
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TRIANGLE';
   case 1010
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_QUATERNION';
   case 1011
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DIMLESS';
   otherwise
      extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NONE';
   end

   extra.NIFTI_INTENT_NAMES = hdr.hist.intent_name;

   if hdr.hist.sform_code > 0
      switch hdr.hist.sform_code
      case 1
         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
      case 2
         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
      case 3
         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
      case 4
         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_MNI_152';
      otherwise
         extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
      end

      extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
   elseif hdr.hist.qform_code > 0
      extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';

      switch hdr.hist.qform_code
      case 1
         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
      case 2
         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
      case 3
         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
      case 4
         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_MNI_152';
      otherwise
         extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
      end
   else
      extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
      extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
   end

   switch bitand(hdr.dime.xyzt_units, 7)	% mask with 0x07
   case 1
      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_METER';
   case 2
      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MM';	% millimeter
   case 3
      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MICRO';
   otherwise
      extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_UNKNOWN';
   end

   switch bitand(hdr.dime.xyzt_units, 56)	% mask with 0x38
   case 8
      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_SEC';
   case 16
      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_MSEC';
   case 24
      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_USEC';	% microsecond
   otherwise
      extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_UNKNOWN';
   end

   switch hdr.dime.xyzt_units
   case 32
      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_HZ';
   case 40
      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_PPM';	% part per million
   case 48
      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_RADS';	% radians per second
   otherwise
      extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_UNKNOWN';
   end

   %  MRI-specific spatial and temporal information
   %
   dim_info = hdr.hk.dim_info;
   extra.NIFTI_FREQ_DIM = bitand(dim_info, 3);
   extra.NIFTI_PHASE_DIM = bitand(bitshift(dim_info, -2), 3);
   extra.NIFTI_SLICE_DIM = bitand(bitshift(dim_info, -4), 3);

   %  Check slice code
   %
   switch hdr.dime.slice_code
   case 1
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_INC';	% sequential increasing
   case 2
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_DEC';	% sequential decreasing
   case 3
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_INC';	% alternating increasing
   case 4
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_DEC';	% alternating decreasing
   case 5
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_INC2';	% ALT_INC # 2
   case 6
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_DEC2';	% ALT_DEC # 2
   otherwise
      extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_UNKNOWN';
   end

   %  Check NIFTI version
   %
   if     ~isempty(hdr.hist.magic) & strcmp(hdr.hist.magic(1),'n') & ...
	( strcmp(hdr.hist.magic(2),'i') | strcmp(hdr.hist.magic(2),'+') ) & ...
	  str2num(hdr.hist.magic(3)) >= 1 & str2num(hdr.hist.magic(3)) <= 9

      extra.NIFTI_VERSION = str2num(hdr.hist.magic(3));
   else
      extra.NIFTI_VERSION = 0;
   end

   %  Check if data stored in the same file (*.nii) or separate
   %  files (*.hdr/*.img)
   %
   if isempty(hdr.hist.magic)
      extra.NIFTI_ONEFILE = 0;
   else
      extra.NIFTI_ONEFILE = strcmp(hdr.hist.magic(2), '+');
   end

   %  Swap has been taken care of by checking whether sizeof_hdr is
   %  348 (machine is 'ieee-le' or 'ieee-be' etc)
   %
   % extra.NIFTI_NEEDS_SWAP = (hdr.dime.dim(1) < 0 | hdr.dime.dim(1) > 7);

   %  Check NIFTI header struct contains a 5th (vector) dimension
   %
   if hdr.dime.dim(1) > 4 & hdr.dime.dim(6) > 1
      extra.NIFTI_5TH_DIM = hdr.dime.dim(6);
   else
      extra.NIFTI_5TH_DIM = 0;
   end

   hdr.extra = extra;

   return;					% extra_nii_hdr

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% POP_SAVEH - save the EEGLAB session command history stored in ALLCOM
%               or in the 'history' field of the current dataset
%
% Usage:
%   >> pop_saveh( ALLCOM, filename, filepath);
%   >> pop_saveh( EEG.history, filename, filepath);
%
% Inputs:
%   ALLCOM      - cell array of strings containing the EEGLAB command history
%   EEG.history - history field of the current dataset
%   filename    - name of the file to save to (optional, default "eeglabhist.m"
%   filepath    - path of the file to save to (optional, default pwd)
%
% Author: Arnaud Delorme, CNL / Salk Institute, 22 March 2002
%
% See also: EEGH, EEGLAB

% Copyright (C) 22 March 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 03-29-02 added update menu -ad

function com = pop_saveh( allcoms, curfilename, curfilepath);

com = '';
if nargin < 1
	help pop_saveh;
	return;
end

if nargin < 3
	[curfilename, curfilepath] = uiputfile('eeglabhist.m', 'Save the EEGLAB session command history with .m extension -- pop_saveh()');
    drawnow;
	if curfilename == 0 return; end
end

fid = fopen( [ curfilepath '/' curfilename ], 'w');
if fid == -1
    error('pop_saveh(): Cannot open named file');
end;
fprintf(fid, '%% EEGLAB history file generated on the %s\n', date);
fprintf(fid, '%% ------------------------------------------------\n');
if iscell(allcoms)
    disp('Saving the EEGLAB session command history...');
    for index = length(allcoms):-1:1
        fprintf(fid, '%s\n', allcoms{index});
    end
    fprintf(fid, 'eeglab redraw;\n');
else
    disp('Saving the current EEG dataset command history...');
    fprintf(fid, '%s\n', allcoms);
end
fclose(fid);

if iscell(allcoms) && nargout == 1
    com = sprintf('pop_saveh( %s, ''%s'', ''%s'');', inputname(1), curfilename, curfilepath);
else
    com = sprintf('pop_saveh( EEG.history, ''%s'', ''%s'');', curfilename, curfilepath);
end

return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function stim = rmJitterImages(stim, params)
% Offset a retinotopic mapping stimulus to compensate for eye position.
%
%   stim = rmJitterImages(stim, params)
%
% The eye movements are read from a file in stim.jitterFile. This file
% contains two variables, x and y, that indicate the x and y positions of
% the jitter across temporal samples. The number of time samples in the eye
% position must much the number of volume acquisitions (TRs), which are
% ordinarily matched to the number of stimulus frames.
%
% Eye movement position is encoded in visual angle, with one (x, y)
% coordinate pair per stimulus frame (or TR).
%
% 12/2008: JW
%
% When describing the stimulus in the vista UI (Analsysis || Retinotopic
% Model || Set Parameters), any matlab file located in the directory
% 'Stimuli' and containing the string 'jitter' in the name will be
% available in the drop down menu for image jitter.
%
% An example of building a jitter file:
%   A scan had 96 frames + 5 pre-scan (removed) frames.
%   The subject fixated 3 deg to the right and above fixation for the
%   entire scan.
%
%   nFrames = 101;
%   deg = 3;
%   x = ones(1, nFrames) * deg;
%   y = ones(1, nFrames) * deg;
%   save Stimuli/jitter3deg101frames x y;
%
% Example to simulate random eye movements with a standard deviation of 1
% deg:
%
%

%% If no eye movement data, return without doing anything
if ~checkfields(stim, 'jitterFile'), return; end
[p, n, e] = fileparts(stim.jitterFile);
if strcmpi(n,'none'), return; end  % A filename of None returns

%% Build jittered images

% Parse eye eyePosition data
[x, y] = parseEyePositionData(stim);

% Get stimulus mesh grid in visual angle
[m n step] = prfSamplingGrid(params);
nrows = size(m,1); ncols = size(m, 2);

% Convert eye positions from degrees to pixels
%   Note that x must be negated, since a positive shift in x means a
%   rightward eye position (and hence a leftward shift in the image),
%   and leftward image shifts are represented by negative numbers. But
%   y is not negated, since a positive shift in y means an upward eye
%   movement (and hence a downward shift in the image), and downward
%   image shifts are represented by positive numbers.
x = -round(x / step);
y = round(y / step);

% Initialize an image of the correct 2D dimenstions
im = zeros(size(m));
inStimWindow = params.stim(1).instimwindow;

% Jitter stimulus frame-by-frame to compensate for eye position
for f = 1:stim.nFrames
    % Reshape image from 1D to 2D
    im(inStimWindow) = stim.images(:, f);

    % Jitter in opp direction to eye movement

    imShifted = zeros(size(im));

    % We can shift in 4 possible ways:
    if x(f) >= 0
        xShifted = 1:ncols-x(f); xUnshifted =  x(f)+1:ncols;
    else
        xShifted = -x(f)+1:ncols; xUnshifted =  1:ncols+x(f);
    end

    if y(f) >= 0
        yShifted = 1:nrows-y(f); yUnshifted =  y(f)+1:nrows;
    else
        yShifted = -y(f)+1:nrows; yUnshifted =  1:nrows+y(f);
    end

    imShifted(yShifted, xShifted) = im(yUnshifted, xUnshifted);

    % Reshape from 2D to 1D
    stim.images(:, f) = imShifted(inStimWindow);
end

% --------------------------------------------------------------
%% test (un-comment to view jittered images as they are created)
% im = zeros(size(m));
% range = mrvMinmax(stim.images(:));
% for f = 1:stim.nFrames
%     % reshape image from 1D to 2D
%     im(inStimWindow) = stim.images(:, f);
%     figure(99);
%     imagesc(im, range);
%     axis image off;
%     pause(0.05)
% end
%---------------------------------------------------------------

return

function [x, y] = parseEyePositionData(stim)
% Read the stimulated eye positions

if ischar(stim.jitterFile)
    tmp = load(stim.jitterFile);
    x = tmp.x; y = tmp.y;
end

if length(x) ~= length(y), error('x/y eye positions not matched'); end

nFrames = size(stim.images,2);
if length(x) < nFrames
    error('Eye movement length (%d) less than nFrames (%d)\n',length(x),nFrames);
elseif length(x) > nFrames
    fprintf('Eye movement length (%d) too long. Truncating to (%d).\n',length(x),nFrames);
    x = x(1:nFrames);
    y = y(1:nFrames);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function fdparain(gobut)
    % Creates the input window for the parameters of the factal dimension calculation.
    %
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    figure_w_normalized_uicontrolunits('Units','pixel','pos',[200 400 550 200 ],'Name','Parameters','visible','off',...
        'NumberTitle','off','Color',ZG.color_bg,'NextPlot','new'); % was color_fbg
    axis off;

    input1 = uicontrol('Style','popupmenu','Position',[.75 .75 .23 .09],...
        'Units','normalized','String','Automatic Range|Manual Fixed Range|Manual',...
        'Value',1, 'Callback', @callbackfun_001);

    input2 = uicontrol('Style','edit','Position',[.34 .43 .10 .09],...
        'Units','normalized','String',num2str(radm), 'enable', 'off',...
        'Value',1, 'Callback', @callbackfun_002);

    input3 = uicontrol('Style','edit','Position',[.75 .43 .10 .09],...
        'Units','normalized','String',num2str(rasm), 'enable', 'off',...
        'Value',1, 'Callback', @callbackfun_003);

    tx1 = text('Position',[0 .85 0 ], ...
        'FontSize',ZmapGlobal.Data.fontsz.m , 'FontWeight','bold' , 'String',' Distance Range within which D is computed: ');

    tx2 = text('Position',[0 .45 0], ...
        'FontSize',ZmapGlobal.Data.fontsz.m , 'FontWeight','bold' , 'String','Minimum value: ', 'color', 'w');

    tx3 = text('Position',[.52 .45 0], ...
        'FontSize',ZmapGlobal.Data.fontsz.m , 'FontWeight','bold' , 'String','Maximum value: ', 'color', 'w');

    tx4 = text('Position',[.41 .45 0], ...
        'FontSize',ZmapGlobal.Data.fontsz.m , 'FontWeight','bold' , 'String','km', 'color', 'w');

    tx5 = text('Position',[.94 .45 0], ...
        'FontSize',ZmapGlobal.Data.fontsz.m , 'FontWeight','bold' , 'String','km', 'color', 'w');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .20 .15 ],...
        'Units','normalized', 'Callback', @callbackfun_004,'String','Cancel');

    switch (gobut)

        case 1   %Defined in timeplot.m, dorand.m,

            go_button=uicontrol('Style','Pushbutton',...
                'Position',[.20 .05 .20 .15 ],...
                'Units','normalized',...
                'callback',@callbackfun_005,...
                'String','Go');

        case 2  %defined in fdtimin.m

            go_button=uicontrol('Style','Pushbutton',...
                'Position',[.20 .05 .20 .15 ],...
                'Units','normalized',...
                'callback',@callbackfun_006,...
                'String','Go');

        case 3  %defined in Dcross.m

            go_button=uicontrol('Style','Pushbutton',...
                'Position',[.20 .05 .20 .15 ],...
                'Units','normalized',...
                'callback',@callbackfun_007,...
                'String','Go');

    end  %switch(gobut)

    set(gcf,'visible','on');
    watchoff;

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        range=(get(input1,'Value'));
        input1.Value=range;
        actrange(range);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        radm=str2double(input2.String);
        input2.String= num2str(radm);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        rasm=str2double(input3.String);
        input3.String= num2str(rasm);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

        org = [2];
        startfd(2);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

        fdtime;
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

        Dcross('ca');
    end

end

author: CelsoReyes


text:
function output = fourier2crsspctrm(cfg, freq)

% FOURIER2CRSSPCTRM transforms a fourier-containing freq-structure
% into a crsspctrm-containing freq-structure, in which the
% powerspectra are also contained in the cross-spectra, being a
% channelcombination of a channel with itself.
%
% Use as
%   [freq] = fourier2crsspctrm(cfg, freq)
%
% where you have the following configuration options:
%   cfg.channel    = cell-array with selection of channels,
%                    see CHANNELSELECTION for details
%   cfg.channelcmb = cell-array with selection of combinations between
%                    channels, see CHANNELCOMBINATION for details
%   cfg.keeptrials = 'yes' or 'no' (default)
%   cfg.foilim     = 2-element vector defining your frequency limits of
%                    interest. By default the whole frequency range of the
%                    input is taken.
%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~isfield(cfg, 'channel'),     cfg.channel     = {'all'};                       end
if ~isfield(cfg, 'channelcmb'),  cfg.channelcmb  = {};                            end
if ~isfield(cfg, 'foilim'),      cfg.foilim      = [freq.freq(1) freq.freq(end)]; end
if ~isfield(cfg, 'keepfourier'), cfg.keepfourier = 'no';                          end
if ~isfield(cfg, 'feedback'),    cfg.feedback    = 'text';                        end

%select the channels on which the power-spectra will be computed
chn = ft_channelselection(cfg.channel,freq.label);
for j = 1:length(chn)
  chnindx(j,1) = find(strcmp(chn(j), freq.label));
  %chnindx(j,1) = find(strcmp(chn{j}, freq.label));
end

%convert the channelcombinations to indices
chncmb  = ft_channelcombination(cfg.channelcmb, freq.label);
cmbindx = zeros(size(chncmb,1),2);
for j = 1:size(chncmb,1)
  cmbindx(j,1) = find(strcmp(chncmb(j,1), freq.label));
  cmbindx(j,2) = find(strcmp(chncmb(j,2), freq.label));
  %cmbindx(j,1) = find(strcmp(chncmb{j,1}, freq.label));
  %cmbindx(j,2) = find(strcmp(chncmb{j,2}, freq.label));
end

%dimensionality of the input data
Nrpt   = length(freq.cumtapcnt);
Nfrq   = size(freq.fourierspctrm,3);
Ntim   = size(freq.fourierspctrm,4);
Nchn   = length(chnindx);
Ncmb   = size(cmbindx,1);

%%FIXME
%if Ntim>1, ft_error('correct handling of time-frequency data is not yet implemented, no information about tapers is available'); end

%keeping track of the tapers
%in the case of tfr fourier-data cumtapcnt is highly redundant; for each frequency
%the number of tapers is equal, as well as for each trial, thus it is sufficient to
%reduce the original cumtapcnt to a vector of Ntrlx1 containing the number of tapers
cumtapcnt = freq.cumtapcnt;
if ~isempty(strfind(freq.dimord, 'time')),
  %cumtapcnt is NtrlxNfrqxNtim; create just one column-vector
  cumtapcnt = ones(size(cumtapcnt,1),1).*unique(cumtapcnt(~isnan(cumtapcnt(:))));
end
sumtapcnt = cumsum([0; cumtapcnt(:)]);

powspctrm = zeros(Nrpt, Nchn, Nfrq, Ntim);
progress('init', cfg.feedback, 'computing single-trial power-spectral densities');
for j = 1:Nrpt
  progress(j/Nrpt, 'trial %d/%d\n', j, Nrpt);
  tmp1     = freq.fourierspctrm([1+sumtapcnt(j):sumtapcnt(j+1)], chnindx, :, :);
  tmp2     = freq.fourierspctrm([1+sumtapcnt(j):sumtapcnt(j+1)], chnindx, :, :);
  powspctrm(j, :, :, :) = squeeze(sum( tmp1 .* conj(tmp2), 1) ./ size(tmp1,1));
end
progress('close');

crsspctrm = complex(zeros(Nrpt, Ncmb, Nfrq, Ntim), zeros(Nrpt, Ncmb, Nfrq, Ntim));
progress('init', cfg.feedback, 'computing single-trial cross-spectral densities');
for j = 1:Nrpt
  progress(j/Nrpt, 'trial %d/%d\n', j, Nrpt);
  tmp1     = freq.fourierspctrm([1+sumtapcnt(j):sumtapcnt(j+1)], cmbindx(:,1), :, :);
  tmp2     = freq.fourierspctrm([1+sumtapcnt(j):sumtapcnt(j+1)], cmbindx(:,2), :, :);
  crsspctrm(j, :, :, :) = squeeze(sum( tmp1 .* conj(tmp2), 1) ./ size(tmp1,1));
end
progress('close');

output.dimord = freq.dimord;
output.freq   = freq.freq;
output.label  = chn;
output.labelcmb(:,1) = freq.label(cmbindx(:,1));
output.labelcmb(:,2) = freq.label(cmbindx(:,2));
output.cumtapcnt = freq.cumtapcnt;
try, output.grad = freq.grad; end
try, output.time = freq.time; end
output.powspctrm = powspctrm;
output.crsspctrm = crsspctrm;
if strcmp(cfg.keepfourier, 'yes'), output.fourierspctrm = freq.fourierspctrm; end

if isempty(output.crsspctrm), output = rmfield(output, 'crsspctrm'); end
if isempty(output.labelcmb ), output = rmfield(output, 'labelcmb' ); end

% add information about the version of this function to the configuration
cfg.version.name = mfilename('fullpath');
cfg.version.id = '$Id$';

% remember the configuration details of the input data
try, cfg.previous = freq.cfg; end

% remember the exact configuration details in the output
output.cfg = cfg;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% SVARGPLVMPREPAREDATA Load and prepare different kinds of data for the shared variational GPLVM.
% DESC Load and prepare different kinds of data for the shared variational GPLVM.
%
% COPYRIGHT: Andreas C. Damianou, Carl Henrik Ek,  2011
% SEEALSO : demSharedVargplvm1
%
% VARGPLVM

%-- Toy (random) data (samples from GPs)
if isempty(dataSetNames) || (isstr(dataSetNames) && strcmp(dataSetNames,'toy'))
    fprintf('# Creating toy data of type %s...\n', toyDataCreate);
    toyData = 1;
    %toyDataCreate = 'vargplvm'; % options: 'vargplvm', 'human','sampling', 'fols'
    dataType = toyDataCreate;
    switch toyDataCreate
        case 'humanPose'
            load humanPose;
            %ind2 = floor(1:2:size(Y,1));
            %%%__silhouettes are whole images
            if exist('imageSilhouette') && imageSilhouette
                load sil_images
            end
            % Remove the 'drunk' sequence
            Y1 = Y(1:100,:);
            Y2 = Y(337:end,:);
            Y = [Y1; Y2];
            clear('Y1','Y2');
            Z1 = Z(1:100,:);
            Z2 = Z(337:end,:);
            Z = [Z1; Z2];
            clear('Z1','Z2');
            % Subsample
            ind2 = floor(1:2:size(Y,1));
            Ytoy{1} = Y(ind2,:);
            Ytoy{2} = Z(ind2,:);
            dataSetNames={'silhouette', 'pose'};
            % X_init = Xp;
            %mappingKern = {'linard2', 'white'};
            %mappingKern = {'rbfard2', 'white'};
            latentDim = 5; % Anything > 2 and < 10
            %xyzankurAnim(Z_test, 3);
        case 'fols'
            %             addpath('../../ncca/matlab/');
            %             addpath('../../sgplvm/matlab/');
            %             addpath('../../fgplvm(svn)/matlab/');
            %             dem_sgplvm_fols; % This creates the dataset
            %             clear model
            %             clear options
            %alpha = linspace(0,2*pi,100);
            alpha = linspace(0,4*pi,100);
            Z1 = cos(alpha)';
            Z2 = sin(alpha)';
            Z3= (cos(alpha)').^2;
            % Z3 = 2*cos(2*alpha)' + 2*sin(2*alpha)' ; %

            % Scale and center data
            bias_Z1 = mean(Z1);
            Z1 = Z1 - repmat(bias_Z1,size(Z1,1),1);
            scale_Z1 = max(max(abs(Z1)));
            Z1 = Z1 ./scale_Z1;

            bias_Z2 = mean(Z2);
            Z2 = Z2 - repmat(bias_Z2,size(Z2,1),1);
            scale_Z2 = max(max(abs(Z2)));
            Z2 = Z2 ./ scale_Z2;

            Z3 = Z3 - repmat(mean(Z3),size(Z3,1),1);%
            Z3 = Z3 ./ (max(max(abs(Z3))));%

            noiseLevel = 0.3; % Default: 0.1
            % Map 1-D to 10-D and add some noise
            Z2p = Z2*rand(1,10);
            Z2p = Z2p + noiseLevel.*randn(size(Z2p));
            Z1p = Z1*rand(1,10);
            Z1p = Z1p + noiseLevel.*randn(size(Z1p));

            Z3p = Z3*rand(1,10);%
            Z3p = Z3p + noiseLevel.*randn(size(Z3p));%

            % pca(Z2p) % This shows that it is actually a 1-D dataset

            % Y = [Z1p Z2p];
            % We like the numer of latent dims to be 2+numSharedDims, ideally 3. With
            % vargplvm we set Q=6 and expect the other 3 or 4 to be switched off.
            % [U,V] = pca(Y,6);
            % Xp = Y*V;
            % pca(Xp)

            %---
            numSharedDims = 5;
            Z1p(:,1:numSharedDims) = Z3p(:,1:numSharedDims);
            Z2p(:,1:numSharedDims) = Z3p(:,1:numSharedDims);
            bar(pca([Z1p Z2p]))
            % return
            %---
            Ytoy{1} = Z1p;
            Ytoy{2} = Z2p;
            %Ytoy{3}= Z3p;%
            dataSetNames={'fols_cos', 'fols_sin'};
            % X_init = Xp;
            mappingKern = {'linard2', 'white'};
            % mappingKern = {'rbfard2', 'white'};
            latentDim = 6; % Anything > 2 and < 10
            % indPoints = 80;
        case 'vargplvm'
            % In the 'vargplvm' case, we generate X from a big Y (which
            % contains signal from rbf in the first half, and from a matern32
            % in the second half of the dims). Then, we set manually the
            % lengthscales for model_i to favor the i partition of X, and
            % generate a new Y_i from the posterior, hopefully associated with
            % the partition X_i. In the end (after optim), we would expect the lengthscales
            % found to be like the ones we manually set here.

            % Constants
            n = dataToKeep;
            D = 20;
            numSubModels = 2;
            latentDimPerModel = 3;
            numSharedDims = 0; % 2;
            %
            latentDim = latentDimPerModel * numSubModels + numSharedDims;
            x = linspace(-1, 1, n)';
            indPoints = round(2/3 * n);

            % Create toy Y as a combination of an rbf and a matern signal
            krn = 'rbf';
            invWidth = 10;
            kern = kernCreate(x,krn);
            kern.inverseWidth = invWidth;
            Kshared = kernCompute(kern, x);
            Y = zeros(n, D);
            for d=1:round(D/2)
                Y(:,d) = real(gsamp(zeros(1, size(x, 1)), Kshared, 1))';
            end

            krn = 'matern32';
            invWidth = 10;
            kern = kernCreate(x,krn);
            kern.inverseWidth = invWidth;
            Kshared = kernCompute(kern, x);
            for d=round(D/2+1):D
                Y(:,d) = real(gsamp(zeros(1, size(x, 1)), Kshared, 1))';
            end

            % Learn a model and X that fit Y
            options = vargplvmOptions('dtcvar');
            options.kern = 'rbfardjit';
            options.numActive = indPoints;
            options.optimiser = 'scg';
            d = size(Y, 2);
            model = vargplvmCreate(latentDim, d, Y, options);
            model = vargplvmParamInit(model, model.m, model.X);
            model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));
            % Optional____
            model.initVardist = 1;
            model = vargplvmOptimise(model, 1, 100);
            model.initVardist = 0;
            model = vargplvmOptimise(model, 1, 300);
            %_____
            % Create new Y_i by manipulating the ard parameters
            % latentIndices = randperm(size( model.vardist.means,2));
            latentIndices = 1:size(model.vardist.means,2);
            allIndices = 1:size(latentIndices,2);
            startVal = 1;
            endVal = numSharedDims;
            sharedIndices = latentIndices(startVal:endVal);
            dataSetNames = {};
            for i=1:numSubModels
                modelTemp = model;
                startVal = endVal+1;
                endVal = startVal + latentDimPerModel - 1;
                curIndices = [latentIndices(startVal:endVal) sharedIndices];
                switchOffIndices = setdiff(allIndices, curIndices);
                modelTemp.kern.inputScales(switchOffIndices) = 0;
                modelTemp.kern.comp{1}.inputScales(switchOffIndices) = 0;
                modelTemp.kern.comp{1}.inputScales(sharedIndices) = (1/2)* modelTemp.kern.comp{1}.inputScales(sharedIndices); % Optional
                modelTemp.kern.inputScales(sharedIndices) = (1/2)* modelTemp.kern.inputScales(sharedIndices); % Optional
                figure, bar(modelTemp.kern.inputScales)
                Ytoy{i} = vargplvmPosteriorMeanVar(modelTemp, modelTemp.vardist.means);
                dataSetNames = {dataSetNames{:}, ['datasetNo' num2str(i)]};
            end
            numberOfDatasets = length(dataSetNames);
            clear('model','modelTemp','options','x','D','n','Y','kern','Kshared','d');
            %{
        case 'sampling' % DOESN'T WORK
            %-- 2nd way: Generate (sample) X from N(0,I), then arbitrarily split X into
            % X_1, X_2,..., X_n, X_s  and fit a GP to then find the posteriors
            % p(Y_1 | X_1, Xs), ..., p(Y_n | X_n, X_s), ... , p(Y_s | X_s). Then, use
            % these Y's as datasets and see if you can get the X's back!!
            n = dataToKeep; % N
            indPoints = round((2/3)*n);
            subModelsNum = 2;
            numSharedDims = 2;
            Q = subModelsNum * (latentDimPerModel+numSharedDims);
            Xsamp = zeros(n,Q);
            for i=1:n
                Xsamp(i,:) = randn(1,Q);
            end
            latentIndices = randperm(size(Xsamp,2));
            startVal = 1;
            endVal =  numSharedDims;
            sharedIndices = latentIndices(startVal:endVal);
            for i=1:subModelsNum
                startVal = endVal+1;
                endVal = startVal + latentDimPerModel;
                X_i = [Xsamp(startVal:endVal) Xsamp(sharedIndices)];
                % Fit a gp
                % ...
                % Ytoy{i} =
            end
            %}
        case 'human'
            load('../../sgplvm/matlab/nccaDemoData.mat')
            Ytoy{1} = Y_train;
            Ytoy{2} = Z_train;
            numberOfDatasets=2;
        case 'rbfOu'
            %-- For private signals
            dataSetNames = {'rand1','rand2','rand3','rand4'};
            kernels = {'rbf', 'rbf', 'ou','ou'};
            n = dataToKeep; % N
            indPoints = round((2/3)*n);
            dims = {40,40,40,40};
            invWidths = {10,10,15,15};
            numberOfDatasets = length(dataSetNames);
            x = linspace(-1, 1, n)';

            %-- For shared signal
            krn = 'rbf';
            invWidth = 100;
            kern = kernCreate(x,krn);
            kern.inverseWidth = invWidth;
            Kshared = kernCompute(kern, x);
            maxD = max(cell2mat(dims));
            Yshared = zeros(n, maxD);
            for d=1:maxD
                Yshared(:,d) = real(gsamp(zeros(1, size(x, 1)), Kshared, 1))';
            end

            %-- Build the datasets by adding the shared with the private signals
            for i=1:numberOfDatasets
                Ytoy{i} = zeros(n, dims{i});
                krn = kernels{i};
                for d=1:dims{i}
                    kern = kernCreate(x, krn);
                    kern.inverseWidth = invWidths{i};
                    K = kernCompute(kern,x);
                    Ytoy{i}(:,d) = real(gsamp(zeros(1, size(x, 1)), K, 1))';
                end
                % Add the shared signal
                Ytoy{i} = Ytoy{i} + Yshared(:,1:dims{i});
            end

            %-- We can add one extra dataset: the shared one
            Ytoy{i+1} = Yshared;
            numberOfDatasets = numberOfDatasets + 1;
            dataSetNames = {dataSetNames{:}, 'shared'};

            clear('n','K','x','kernels','krn','kern','d');
        otherwise
            %...
    end
    if ~exist('numberOfDatasets')
        numberOfDatasets = length(Ytoy);
    end
    %     for i=1:numberOfDatasets
    %         figure
    %         plot(Ytoy{i})
    %     end
    Yall = Ytoy;
    clear('Ytoy');
else
    toyData = 0;
    if isstruct(dataSetNames)
        for i=1:length(dataSetNames)
            Y{i} = svargplvmLoadData(dataSetNames{i});
        end
        numberOfDatasets = length(dataSetNames);
    else
        [Yall, lbls] = svargplvmLoadData(dataSetNames);
        if ~isempty(lbls)
            height = lbls(1); width = lbls(2);
        end
        numberOfDatasets = length(Yall);
    end
end
%--

if isempty(dataType)
    if length(dataSetNames) ==1
        dataType = dataSetNames{1};
    end
end

author: SheffieldML

repo: vargplvm


text:
function r = ge(p,q)
% this function compares the greater than or equality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVariable')&&isa(q, 'CellVariable'))
    r=p;
    r.value = p.value>=q.value;
elseif isa(p, 'CellVariable')
    r=p;
    r.value = p.value>=q;
else
    r=q;
    r.value = p>=q.value;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function hb_file_write ( output_unit, title, key, totcrd, ptrcrd, indcrd, ...
  valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, ...
  valfmt, rhsfmt, rhstyp, nrhs, nrhsix, colptr, rowind, values, ...
  rhsval, rhsptr, rhsind, rhsvec, guess, exact )

%*****************************************************************************80
%
%% HB_FILE_WRITE writes an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer OUTPUT_UNIT, the unit to which data is written.
%
%    Input, character ( len = 72 ) TITLE, a title for the matrix.
%
%    Input, character ( len = 8 ) KEY, an identifier for the matrix.
%
%    Input, integer TOTCRD, the total number of lines of data.
%
%    Input, integer PTRCRD, the number of input lines for pointers.
%
%    Input, integer INDCRD, the number of input lines for row indices.
%
%    Input, integer VALCRD, the number of input lines for numerical values.
%
%    Input, integer RHSCRD, the number of input lines for right hand sides.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NCOL, the number of columns or elements.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Input, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Input, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer NRHSIX, the number of row indices (set to 0
%    in the case of unassembled matrices.)  Ignored if NRHS = 0.
%
%    Input, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Input, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%
%    Input, real VALUES(NNZERO) or VALUES(NELTVL), the nonzero values
%    of the matrix.
%
%    If RHSTYP(1:1) == 'F':
%
%      Input, real RHSVAL(NROW,NRHS), contains NRHS dense right hand
%      side vectors.
%
%      Input, integer RHSPTR(*), is not used.
%
%      Input, integer RHSIND(*), is not used.
%
%      Input, real RHSVEC(*), is not used.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'A':
%
%      Input, real RHSVAL(NROW,NRHS), is not used.
%
%      Input, integer RHSPTR(NRHS+1), RHSPTR(I) points to the location of
%      the first entry of right hand side I in the sparse right hand
%      side vector.
%
%      Input, integer RHSIND(NRHSIX), indicates, for each entry of
%      RHSVAL, the corresponding row index.
%
%      Input, real RHSVEC(NRHSIX), contains the value of the right hand
%      side entries.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'E':
%
%      Input, real RHSVAL(NNZERO,NRHS), contains NRHS unassembled
%      finite element vector right hand sides.
%
%      Input, integer RHSPTR(*), is not used.
%
%      Input, integer RHSIND(*), is not used.
%
%      Input, real RHSVEC(*), is not used.
%
%    Input, real GUESS(NROW,NRHS), the starting guess vectors.
%
%    Input, real EXACT(NROW,NRHS), the exact solution vectors.
%

%
%  Write the header block.
%
  hb_header_write ( output_unit, title, key, totcrd, ptrcrd, indcrd, ...
    valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, ...
    valfmt, rhsfmt, rhstyp, nrhs, nrhsix );
%
%  Write the matrix structure.
%
  hb_structure_write ( output_unit, ncol, mxtype, nnzero, neltvl, ...
    ptrfmt, indfmt, colptr, rowind );
%
%  Write the matrix values.
%
  hb_values_write ( output_unit, valcrd, mxtype, nnzero, neltvl, ...
    valfmt, values );
%
%  Write the right hand sides.
%
  hb_rhs_write ( output_unit, nrow, nnzero, nrhs, nrhsix, ...
    rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp, rhsval, ...
    rhsind, rhsptr, rhsvec );
%
%  Write the starting guesses.
%
  hb_guess_write ( output_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp, ...
    guess );
%
%  Write the exact solutions.
%
  hb_exact_write ( output_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp, ...
    exact );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s = convert(tree,uid,varargin)
% XMLTREE/CONVERT Convert an XML tree into a structure
%
% tree     - XMLTree object
% uid      - uid of the root of the subtree, if provided.
%            Default is root
% s        - converted structure
%__________________________________________________________________________
%
% Convert an XMLTree into a structure, when possible.
% When several identical tags are present, a cell array is used.
% The root tag is not saved in the structure.
% If provided, only the structure corresponding to the subtree defined
% by the uid UID is returned.
%
% Example:
% xml = 'field1<c>field2</c>field3';
% tree = convert(xmltree(xml));
% <=> tree = struct('b',{{'field1', 'field3'}},'c','field2')
%__________________________________________________________________________
% Copyright (C) 2002-2019  https://www.artefact.tk/

% Guillaume Flandin
% $Id: convert.m 7623 2019-06-26 09:34:46Z guillaume $

% Get the root uid of the output structure
if nargin == 1 || isempty(uid)
    % Get the root uid of the XML tree
    uid = root(tree);
end

% Get optional parameters
sub_options({'attributes',false,'str2num',false});
if nargin > 2
    sub_options(varargin);
end

% Build the output structure
s = struct();
s = sub_convert(tree,s,uid,{}); % {get(tree,root_uid,'name')}

%==========================================================================
function s = sub_convert(tree,s,uid,arg)
    type = get(tree,uid,'type');
    switch type
        case 'element'
            child = children(tree,uid);
            l = {};
            ll = {};
            for i=1:length(child)
                if isfield(tree,child(i),'name')
                    ll = [ll, get(tree,child(i),'name')];
                end
            end
            for i=1:length(child)
                if isfield(tree,child(i),'name')
                    name = get(tree,child(i),'name');
                    nboccur = sum(ismember(l,name));
                    nboccur2 = sum(ismember(ll,name));
                    l = [l, name];
                    if nboccur || nboccur2 > 1
                        arg2 = [arg, name, {{nboccur+1}}];
                    else
                        arg2 = [arg, name];
                    end
                else
                    arg2 = arg;
                end
                s = sub_convert(tree,s,child(i),arg2);
            end
            if isempty(child)
                s = sub_setfield(s,arg{:},'');
            end
            if sub_options('attributes')
                attrb = attributes(tree,'get',uid);
                if ~isempty(attrb)
                    arg2 = [arg, 'attributes'];
                    if ~isstruct(attrb), attrb = [attrb{:}]; end
                    try
                        % Saving attributes will work with <c>b</c>
                        % but not with b
                        s = sub_setfield(s,arg2{:},...
                            cell2struct({attrb.val},{attrb.key},2));
                    end
                end
            end
        case 'chardata'
            s = sub_setfield(s,arg{:},get(tree,uid,'value'));
            if sub_options('str2num')
                % Convert strings into their numerical equivalent when possible
                % e.g. string '3.14159' becomes double scalar 3.14159
                v = get(tree,uid,'value');
                cv = str2num(v);
                if isempty(cv)
                    s = sub_setfield(s,arg{:},v);
                else
                    s = sub_setfield(s,arg{:},cv);
                end
            end
        case 'cdata'
            s = sub_setfield(s,arg{:},get(tree,uid,'value'));
        case 'pi'
            % Processing instructions are evaluated if possible
            app = get(tree,uid,'target');
            switch app
                case {'matlab',''}
                    s = sub_setfield(s,arg{:},eval(get(tree,uid,'value')));
                case 'unix'
                    s = sub_setfield(s,arg{:},unix(get(tree,uid,'value')));
                case 'dos'
                    s = sub_setfield(s,arg{:},dos(get(tree,uid,'value')));
                case 'system'
                    s = sub_setfield(s,arg{:},system(get(tree,uid,'value')));
                otherwise
                    try
                        s = sub_setfield(s,arg{:},feval(app,get(tree,uid,'value')));
                    catch
                        warning('[XMLTree] Unknown target application');
                    end
            end
        case 'comment'
            % Comments are ignored
        otherwise
            warning(sprintf('Type %s unknown : not saved',get(tree,uid,'type')));
    end

%==========================================================================
function s = sub_setfield(s,varargin)
% Same as setfield but using '{}' rather than '()'

subs = varargin(1:end-1);
types = repmat({'{}'},1,numel(subs));
types(cellfun(@ischar,subs)) = {'.'};
s = builtin('subsasgn', s, struct('type',types,'subs',subs), varargin{end});

%==========================================================================
function o = sub_options(opt)
persistent opts
if isempty(opts), opts = struct; end
if iscell(opt)
    for i=1:2:numel(opt)
        opts.(lower(opt{i})) = opt{i+1};
    end
else
    o = opts.(lower(opt));
end

author: spm


text:
function  wtree = ind2wtree(pyr, ind)

%this function is called by vifvec.m
% converts the output of Eero Simoncelli's pyramid routines into subbands in a cell array
C=pyr;
S=ind;

offset=0;
numsubs=size(ind,1);
for i=1:numsubs
    wtree{numsubs-i+1}=reshape(C(offset+1:offset+prod(S(i,:))), S(i,1),S(i,2));
    offset=offset+prod(S(i,:));
end

repo: SelfExSR


text:
function [res,ires]=dydt(t,y,yp,rpar)

global DEBUG;

res(1)=yp(1)+rpar(1)*y(1);
res(2)=yp(2)+rpar(2)*y(2);
res(3)=y(3)-y(1)-2.0;

ires = 0;

if ( DEBUG > 9 )
  buf1 = sprintf('%.4e ',y);
  buf2 = sprintf('%.4e ',yp);
  buf3 = sprintf('%.4e ',res);
  disp(sprintf('(Debug) dydt.m: t=%.2e y=%s yp=%s res=%s', ...
               t,buf1,buf2,buf3));
  if ( DEBUG > 99 )
    cbuf = input(sprintf('\tRETURN to continue : '));
  end
end

author: Sable

repo: mcbench-benchmarks


text:
function [PY_x coeff] = bss_tvproj(x,Y,tvshape,tvstep)

% compute the orthogonal projection of x onto the space of shifted windowed versions of the row(s) of Y
%
% Usage: [PY_x coeff] = bss_tvproj(x,Y,tvshape,tvstep)
%
% Input:
%   - x: row vector of length T corresponding to the signal to be projected,
%   - Y: vector or matrix of length T with n_rows rows which windowed rows span the projection space.
%   - tvshape : row vector of length V corresponding to the window applied
%   to Y to define the projection space
%   - tvstep  : number of samples between two adjacent windows
%
% Ouput:
%   - PY_x: row vector of length T containing the orthogonal projection of
%   x onto the range of the shifted windowed versions of the row(s) of Y.
%   - coeff : matrix with n_rows rows and n_frames columns containing the coefficients
%   of the projection
%
% Developers:  - Cedric Fevotte (fevotte@tsi.enst.fr) - Emmanuel Vincent
% (emmanuel.vincent@irisa.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

% 1. Padd the signals with zeros on the left and the right
V        = size(tvshape,2); % the size of the window
NOVERLAP = V-tvstep;        % convert the hop size into a number of overlapping samples
NPadLeft = tvstep;
NPadRight= 2*V;
x = [zeros(size(x,1),NPadLeft) x zeros(size(x,1),NPadRight)];
Y = [zeros(size(Y,1),NPadLeft) Y zeros(size(Y,1),NPadRight)];

% 2. Decompose Y into frames using tvshape as a window

[Y_frames frames_index] = bss_make_frames(Y,tvshape,NOVERLAP);

% 3. Y_frames is a 3-D array of size n_frames x V x n_rows
% we reshape it to get a 2-D array of size (n_frames x n_rows) x V

[n_frames V1 n_rows]= size(Y_frames); % The frames are Y_frames(f,:,n)
Y_frames = permute(Y_frames,[2 1 3]);
% Y_frames is now V x n_frames x nrows, the frames are Y_frames(:,f,n)

Y_frames = reshape(Y_frames,V,n_frames*n_rows);
% Y_frames is now V x (n_frames x n_rows), the frames are Y_frames(:,k)
% with k=(n-1)*n_frames+f

Y_frames = permute(Y_frames,[2 1]);
% Y_frames is now (n_frames x n_rows) x V, the frames are Y_frames(k,:)
% with k=(n-1)*n_frames+f

% 4. Compute the inner products between x and the frames of the row(s) of Y

x_frames  = bss_make_frames(x,ones(1,V),NOVERLAP);
% x_frames is an array n_frames x V x 1 which frames are x_frames(n,:)

ip = zeros(n_frames*n_rows,1);
for n=1:n_rows % loop on rows of Y
    idxrange = ((n-1)*n_frames)+(1:n_frames); % the index of the frames in Y_frames corresponding to the selected row
    ip(idxrange) = sum(x_frames .* Y_frames(idxrange,:),2);
end

% 5. Compute the Gram matrix, which is square of size (n_frames x n_rows) x
% (n_frames x n_rows)
Gram = sparse(n_frames*n_rows,n_frames*n_rows);

for f=1:n_frames
    for f1 = 1:n_frames
        % locate the range of the intersection between frames
        first = max(frames_index(f),frames_index(f1));
        last  = min(frames_index(f),frames_index(f1))+V-1;
        % if the intersection is non empty, fill it
        if (first <= last)
            nrange  = (first:last)-frames_index(f)+1;
            nrange1 = (first:last)-frames_index(f1)+1;
            idxrange  = f  + n_frames*((1:n_rows)-1);
            idxrange1 = f1 + n_frames*((1:n_rows)-1);
            Gram(idxrange,idxrange1) = Y_frames(idxrange,nrange)*Y_frames(idxrange1,nrange1)';
        end
    end
end

% 6. Apply the inverse of the Gram matrix to ip to get coeff
%coeff = full(Gram)\ip;
coeff = pinv(full(Gram))*ip;
% 7.  Reconstruct using coefficients alpha
PY_x = zeros(size(x));

for f=1:n_frames
    idxrange = f + n_frames*((1:n_rows)-1);
    trange   = frames_index(f)+(0:(V-1));
    PY_x(trange) = PY_x(trange) + coeff(idxrange)'*Y_frames(idxrange,:);
end

% 8. Truncate the projection to the support of the signal
PY_x = PY_x(:,(NPadLeft+1):(end-NPadRight));

% 9. The coefficients are a n_frames*n_rows column vector with entries
% coeff(f+(n-1)*n_frames). We reshape it
coeff = reshape(coeff,n_frames,n_rows);
coeff = permute(coeff,[1 2]);

author: posenhuang

repo: deeplearningsourceseparation


text:
%                        experimentImage2D.m
%
% Create a Fourier measurement operator and measurements for reconstructing
% an image.
%
% Inputs:
% numMasks: number of random octanary Fourier masks.
% imagePath: The path of the image to be recovered
%
% Outputs:
%  A     : A function handle: n1*n2 x 1 -> n1*n2*numMasks x 1. It returns
%          A*x.
%  At    : The transpose of A
%  b0    : A n1*n2*numMasks x 1 real, non-negative vector consists of the
%          measurements abs(A*x).
%  Xt    : The true signal - a vectorization of the image
%
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

function [A, At, b0, Xt, plotter] = experimentImage2D(numMasks, imagePath)
    %%  Read Image
    % Below Xt is n1 x n2 x 3; i.e. we have three n1 x n2 images, one for each of the 3 color channels
    image = rgb2gray(imread([imagePath]));
    image = double(image);
    dims = size(image);
    L = numMasks;

    %% Make octanary masks and linear sampling operators
    % Each mask has iid entries following the octanary pattern.
    % Entries have the form b = b1*b2 , where b1 is sampled
    % from {-1, 1, -i, i} with equal probability 1/4, and b2 from
    % { sqrt(1/2), sqrt(3)} with probability 4/5 and 1/5 respectively.
    b1 = [-1;1;-1i;1i];
    b2 = [repmat(sqrt(0.5),4,1); sqrt(3)];
    % Masks has size [n1,n2,L]
    masks = b1(randi(4,[dims,L])) .* b2(randi(5,[dims,L])); % Storage for L masks, each of dim n1 x n2

    %% Make linear operators that act on a vectorized image
    A = @(x) fourierMeasurementOperator(x, masks, dims);
    At = @(y) transposeOperator(y, masks, dims);

    Xt = image(:);
    b0 = abs(A(Xt(:)));

    % Set up plotting
    subplot(1,2,1);
    imagesc(image);
    title('original');
    subplot(1,2,2);

    plotter = @(x) imagesc(reshape(abs(x),dims));

end

function y = fourierMeasurementOperator(x, masks, dims)
    x = reshape(x, dims);   % image comes in as a vector.  Reshape to rectangle
    [n1,n2] = size(x);
    L = size(masks,3);              % get number of masks

    % Compute measurements
    copies = repmat(x,[1,1,L]);
    y = fft2(masks.*copies);
    y = y(:);
end

function x = transposeOperator(y, masks, dims)
    n1 = dims(1);
    n2 = dims(2);
    L = size(masks,3);              % get number of masks
    y = reshape(y, [n1,n2,L]);   % image comes in as a vector.  Reshape to rectangle

    x = n1*n2*ifft2(y).*conj(masks);
    x = sum(x,3);
    x = x(:);
end

author: tomgoldstein

repo: phasepack-matlab


text:
% This subroutine assigns creates a grid with
% spacing dx,dy (in degreees). The size will
% be selected interactiVELY. The bvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/95

disp ('This is /src/stressdepth_ratio.m');
%% disp  ('PLEASE SET my_dir TO YOUR WORKING DIRECTORY!!!!!');

global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 0.1;
    dy = 0.1 ;
    ni = 1000;
    Nmin = 50;
    stan2 = nan;
    stan = nan;
    prf = nan;
    av = nan;
    mid_point = 5;
    top_zonet = 0;
    top_zoneb = 5;
    bot_zonet = 7;
    bot_zoneb = 15;
    topstan2 = nan;
    botstan2 = nan;
    topstan = nan;
    botstan = nan;
    topav = nan;
    botav = nan;
    use_old_win = 1;
    lab1 = 'b-value-depth-ratio:';

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Depth Ratio Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ wex+200 wey-200 650 300]);
    axis off
    labelList2=[' Automatic Mcomp (max curvature) | Fixed Mc (Mc = Mmin) | Automatic Mcomp (90% probability) | Automatic Mcomp (95% probability) | Best (?) combination (Mc95 - Mc90 - max curvature)'];

    labelPos = [0.2 0.60 0.6  0.08];
    hndl2=uicontrol(...
        'Style','popup',...
        'Position',labelPos,...
        'Units','normalized',...
        'String',labelList2,...
        'Callback','inb2 =get(hndl2,''Value''); ');

    % set(hndl2,'value',5);

    % creates a dialog box to input grid parameters
    %
    %

    % mid_point_field=uicontrol('Style','edit',...
    %     'Position',[.47 .80 .12 .08],...
    %     'Units','normalized','String',num2str(mid_point),...
    %     'Callback','mid_point=str2double(get(mid_point_field,''String'')); set(mid_point_field,''String'',num2str(mid_point));');

    top_zonet_field=uicontrol('Style','edit',...
        'Position',[.36 .80 .06 .06],...
        'Units','normalized','String',num2str(top_zonet),...
        'Callback','top_zonet=str2double(get(top_zonet_field,''String'')); set(top_zonet_field,''String'',num2str(top_zonet));');
    top_zoneb_field=uicontrol('Style','edit',...
        'Position',[.36 .74 .06 .06],...
        'Units','normalized','String',num2str(top_zoneb),...
        'Callback','top_zoneb=str2double(get(top_zoneb_field,''String'')); set(top_zoneb_field,''String'',num2str(top_zoneb));');

    bot_zonet_field=uicontrol('Style','edit',...
        'Position',[.78 .80 .06 .06],...
        'Units','normalized','String',num2str(bot_zonet),...
        'Callback','bot_zonet=str2double(get(bot_zonet_field,''String'')); set(bot_zonet_field,''String'',num2str(bot_zonet));');
    bot_zoneb_field=uicontrol('Style','edit',...
        'Position',[.78 .74 .06 .06],...
        'Units','normalized','String',num2str(bot_zoneb),...
        'Callback','bot_zoneb=str2double(get(bot_zoneb_field,''String'')); set(bot_zoneb_field,''String'',num2str(bot_zoneb));');

    freq_field=uicontrol('Style','edit',...
        'Position',[.30 .50 .12 .08],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));set(tgl2,''value'',0); set(tgl1,''value'',1)');

    freq_field0=uicontrol('Style','edit',...
        'Position',[.70 .50 .12 .08],...
        'Units','normalized','String',num2str(ra),...
        'Callback','ra=str2double(get(freq_field0,''String'')); set(freq_field0,''String'',num2str(ra)) ; set(tgl2,''value'',1); set(tgl1,''value'',0)');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.30 .40 .12 .08],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.30 .30 .12 .08],...
        'Units','normalized','String',num2str(dy),...
        'Callback','dy=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dy));');

    tgl1 = uicontrol('Style','checkbox',...
        'string','Number of Events:',...
        'Position',[.09 .50 .2 .08], 'Callback','set(tgl2,''value'',0)',...
        'Units','normalized');

    set(tgl1,'value',1);

    tgl2 =  uicontrol('Style','checkbox',...
        'string','OR: Constant Radius',...
        'Position',[.47 .50 .2 .08], 'Callback','set(tgl1,''value'',0)',...
        'Units','normalized');

    freq_field4=uicontrol('Style','edit',...
        'Position',[.30 .20 .12 .08],...
        'Units','normalized','String',num2str(Nmin),...
        'Callback','Nmin=str2double(get(freq_field4,''String'')); % set(freq_field4,''String'',num2str(Nmin));');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback',' inb1 =get(hndl2,''Value'');tgl1 =get(tgl1,''Value'');tgl2 =get(tgl2,''Value'');close,sel =''ca'', stressdepth_ratio',...
        'String','Go');
    Nmin;
    text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.20 .75 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Please choose an Mc estimation option   ');

    mid_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.24 1.0 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Depth limits for depth ratio calculation  ');

    top_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.10 .85 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Top and bottom for TOP zone(km):');

    bot_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.40 .85 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Top and bottom for BOTTOM zone(km):');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.64 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.40 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in x (dx) in deg:');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.29 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in y (dy) in deg:');

    txt7 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.17 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Min. No. of events > Mc:');

    set(gcf,'visible','on');
    watchoff

end   % if nargin ==0

% for plotting number of events used (ni) in view_bdepth
ni_plot = ni;

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'

    [file1,path1] = uiputfile([my_dir fs '*.mat'], 'Grid Datafile Name?') ;

    selgp
    itotal = length(newgri(:,1));
    if length(gx) < 4  ||  length(gy) < 4
        errordlg('Selection too small! (Dx and Dy are in degreees! ');
        return
    end

    zmap_message_center.set_info(' ','Running bdepth_ratio... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    bvg = [];
    allcount = 0.;
    nobv = 0;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','b-value grid - percent done');;
    drawnow

    % sort by depth

    % [s,is] = sort(a.Depth);
    % adepth = a(is(:,1),:);

    % find row index of ratio midpoint
    l = a.Depth >= top_zonet & a.Depth <  top_zoneb;
    top_zone = a.subset(l);

    l = a.Depth >= bot_zonet & a.Depth <  bot_zoneb;
    bot_zone = a.subset(l);

    %
    % overall b-value
    [bv magco stan av me mer me2,  pr] =  bvalca3(top_zone,inb1,inb2);
    tbo1 = bv; tno1 = length(top_zone(:,1));

    [bv magco stan av me mer me2,  pr] =  bvalca3(bot_zone,inb1,inb2);
    bbo1 = bv; bno1 = length(bot_zone(:,1));

    depth_ratio = tbo1/bbo1;

    disp(depth_ratio);
    hodis = [hodi '/stinvers'];
    do = ['cd  ' hodis ]; eval(do)

    % loop over all points
    for i= 1:length(newgri(:,1))
        x = newgri(i,1);y = newgri(i,2);
        allcount = allcount + 1.;
        i2 = i2+1;

        % calculate distance from center point and sort wrt distance
        l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2) ;
        [s,is] = sort(l);
        b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

        if tgl1 == 0   % take point within r
            l3 = l <= ra;
            b = a.subset(l3);      % new data per grid point (b) is sorted in distanc  (from center point)
            rd = ra;
        else
            % take first ni points
            b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
            l2 = sort(l); rd = l2(ni);

        end

        %estimate the completeness and b-value
        newt2 = b;

        % sort by depth

        l = b(:,7) >= top_zonet & b(:,7) <  top_zoneb;
        topb = b(l,:);
        per_in_top = (length(topb)/length(b))*100.0;
        l = b(:,7) >= bot_zonet & b(:,7) <  bot_zoneb;
        botb = b(l,:);
        per_in_bot = (length(botb)/length(b))*100.0;

        if length(topb) >= Nmin & length(botb) >= Nmin

            tmpi = [topb(:,10:12)];
            fid = fopen('data2','w');
            str = ['Inversion data'];str = str';

            fprintf(fid,'%s  \n',str');
            fprintf(fid,'%7.3f  %7.3f  %7.3f\n',tmpi');

            fclose(fid);
            delete data2.slboot Xtemp.slboot

            unix(' slfast data2 ');
            load data2.slboot
            dtop = data2;

            tmpi = [botb(:,10:12)];
            fid = fopen('data2','w');
            str = ['Inversion data'];str = str';

            fprintf(fid,'%s  \n',str');
            fprintf(fid,'%7.3f  %7.3f  %7.3f\n',tmpi');

            fclose(fid);
            delete data2.slboot Xtemp.slboot

            unix(' slfast data2 ');
            load data2.slboot
            dbot = data2;

            ltopb = length(topb);
            lbotb = length(botb);
            bvg = [bvg ; dtop(2,2:7) dtop(1,1) dbot(2,2:7) dbot(1,1) ltopb lbotb ];

            waitbar(allcount/itotal)
        end

    end

    % plot the results
    % old and re3 (initially ) is the b-value matrix
    %
    normlap2=ones(length(tmpgri(:,1)),1)*nan;
    normlap2(ll)= bvg(:,1);
    re3=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,5);
    r=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,6);
    meg=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,2);
    old1=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,7);
    avm=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,8);
    Prmap=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,9);
    top_b=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,10);
    bottom_b=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,11);
    per_top=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,12);
    per_bot=reshape(normlap2,length(yvect),length(xvect));

    %  normlap2(ll)= bvg(:,13);
    %  ltopb=reshape(normlap2,length(yvect),length(xvect));

    % normlap2(ll)= bvg(:,14);
    % lbotb=reshape(normlap2,length(yvect),length(xvect));

    old = re3;

    % View the b-value map
    view_bdepth

end   % if sel = na

% Load exist b-grid
if sel == 'lo'
    [file1,path1] = uigetfile(['*.mat'],'b-value gridfile');
    if length(path1) > 1
        think
        load([path1 file1])
        normlap2=ones(length(tmpgri(:,1)),1)*nan;

        normlap2(ll)= bvg(:,1);
        re3=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,5);
        r=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,6);
        meg=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,2);
        old1=reshape(normlap2,length(yvect),length(xvect));

        %  normlap2(ll)= bvg(:,7);
        %  pro=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,7);
        avm=reshape(normlap2,length(yvect),length(xvect));

        %  normlap2(ll)= bvg(:,9);
        % stanm=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,8);
        Prmap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,9);
        top_b=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,10);
        bottom_b=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,11);
        per_top=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,12);
        per_bot=reshape(normlap2,length(yvect),length(xvect));

        %    normlap2(ll)= bvg(:,13);
        %    ltopb=reshape(normlap2,length(yvect),length(xvect));

        %   normlap2(ll)= bvg(:,14);
        %  lbotb=reshape(normlap2,length(yvect),length(xvect));

        old = re3;

        view_bdepth
    else
        return
    end
end


text:
function [hx, hy, lx, ly] = Q0011gridfdims(F00, F11)
%------------------------------------------------------------------------------
%
% This function determines the dimensions of quincunx gridfunction {F00 U F11}.
%
% [hx, hy] = Q0011gridfdims(F00, F11) yields meshwidths in x- and y-dir.
%
% [hx, hy, lx, ly] = Q0011gridfdims(F00, F11) yields meshwidths
%                    in x- and y-dir
%                    and dimensions of domain of {F00 U F11}.
%
%         x m
%     o---->
%     |
%   y |
%     v
%   n
%
% See also: gridfdims, Q1001gridfdims
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: February 1, 2001.
% (c) 1999-2001 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[n00, m00]=size(F00);
[n11, m11]=size(F11);
if n11 > n00
  error(' Q0011gridfdims - dimensions do not match for n11 > n00 ')
end
if m11 > m00
  error(' Q0011gridfdims - dimensions do not match for m11 > m00 ')
end
if n00 > n11+1
  error(' Q0011gridfdims - dimensions do not match for n00 > n11+1 ')
end
if m00 > m11+1
  error(' Q0011gridfdims - dimensions do not match for m00 > m11+1 ')
end
%
% m11 <= m00 <= m11+1 is satisfied
% n11 <= n00 <= n11+1 is satisfied
%
%------------------------------------------------------------------------------
n = n00 + n11;
m = m00 + m11;
if nargout == 2
  if m<n
    hx = 1.0/m;
    hy = hx;
  else
    hy = 1.0/n;
    hx = hy;
  end
elseif nargout == 4
  if m<n
    lx = 1.0;
    ly = (n*lx)/m;
    hx = lx/m;
    hy = hx;
  else
    ly = 1.0;
    lx = (m*ly)/n;
    hy = ly/n;
    hx = hy;
  end
else
  error(' Q0011gridfdims - wrong number of output arguments ')
end
%------------------------------------------------------------------------------


text:
function [segclassifier] = mcmcTrainSegmentClassifier2(features, labels, weights, maxdata, classparams)

if exist('classparams') && ~isempty(classparams)
    nnodes = classparams(1);
    ntrees = classparams(2);
    stopval = classparams(3);
else
    ntrees = 20;
    nnodes = 8;
    stopval = 0;
end

if ~exist('maxdata') || isempty(maxdata)
    maxdata = 25000;
end

[data, lab, w] = formatData(features, labels, weights, maxdata);

segclassifier = train_boosted_dt_2c(data, [], lab, ntrees, nnodes, stopval, w);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [data, lab, w] = formatData(features, labels, weights, maxdata)
% concatenate data
%disp(num2str(maxdata))
nimages = numel(features);

[tmp, nvars] = size(features{1});

% count segments
tseg = 0;
nseg = 0;
for f = 1:nimages
    nseg = nseg + sum((labels{f}~=0) & (weights{f}>0));
    tseg = tseg + numel(labels{f});
end
%disp(num2str([nseg tseg]))

data = zeros(nseg, nvars);
lab = zeros(nseg, 1);
w = zeros(nseg, 1);

% concatenate data
vc = 0;
for f = 1:nimages
    ind = find((labels{f}~=0) & (weights{f}>0));
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    lab(vc+1:vc+numel(ind)) = (labels{f}(ind)>0)*2-1;
    w(vc+1:vc+numel(ind)) = weights{f}(ind); % weight according to area in image
    vc = vc + numel(ind);
end

if nseg > maxdata
    rind = randperm(nseg);
    rind = rind(1:maxdata);
    data = data(rind, :);
    lab = lab(rind);
    w = w(rind);
end

w = w / sum(w);

author: zouchuhang

repo: LayoutNet


text:
function tagS = multi_frame_module_tags
%"multi_frame_module_tags"
%   Return the tags used to represent a multiframe study as specified by
%   C.7.6.6 in PS3.3 of 2006 DICOM specification.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%JRA 07/10/06
%
%Usage:
%   tagS = multi_frame_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
%template = tagS;

%Add tags based on PS3.3 attribute lists.

%Number of Frames
tagS(end+1) = struct('tag', '00280008', 'tagdec', 2621448, 'type', '1C', 'children', []);

%Frame Increment Pointer
tagS(end+1) = struct('tag', '00280009', 'tagdec', 2621449, 'type', '1C', 'children', []);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function export_VPSC(ori,filename,varargin)
% export individual orientations to the VPSC format
%
% Syntax
%   export_VPSC(ori,'file.txt')
%   export_VPSC(ori,'file.txt','weights',weights)
%
% Input
%  ori      - individual @orientation to be exported
%  filename - name of the ascii file
%  weights  - list weights with the same size as the orientations
%
% See also
% quaternion/export SO3Fun/export_VPSC

% allocate memory
d = zeros(length(ori),4);

% add Euler angles
d(:,1:3) = ori.Euler(varargin{:});
if ~check_option(varargin,{'radians','radiant','radiand'})
  d = d ./ degree;
end

w = get_option(varargin,'weights',ones(size(ori)));

% add weight
d(:,4) = w./ sum(w);

fid = efopen(filename,'w');

% header
% fourth line has to include the angle convention (B for Bunge)
% and the total number of grains in the phase
fprintf(fid,'\n\n\nB %d\n',length(ori));

fprintf(fid,'%7.2f %7.2f %7.2f %11.7f\n',d');

fclose(fid);

author: mtex-toolbox

repo: mtex

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


text:
function dt = getMetadata(fileName)
	c = benchmarkPaths(0);
	dt = load(sprintf('%s/metadata/%s.mat', c.benchmarkDataDir, fileName));
end

author: s-gupta

repo: rgbd

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


text:
function [ spath ] = tokens2path( tokens )
%TOKENS2PATH Summary of this function goes here
%   Detailed explanation goes here

    spath = tokens{1};

    for i = 2 : numel(tokens)
        spath = [spath '/' tokens{i}];
    end;

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

function X = construct_initial_guess(L, F, r, n)
% Basicially only the first micro-step of ALS

X = TTeMPS_rand( r, n );
X = 1/norm(X) * X;

d = X.order;
n = X.size;

X = orthogonalize(X, 1);
Fi = contract( X, F, 1 );
sz = [X.rank(1), X.size(1), X.rank(2)];

[left, right] = Afun_prepare( L, X, 1 );
B1 =  prepare_precond( L.A{1}, X, 1 );

Ui = pcg( @(y) Afun( L, y, 1, sz, left, right), ...
         Fi(:), ...
         1e-10, 1000, ...
         @(y) apply_precond( L.A{1}, B1, y, sz ), [],...
         X.U{1}(:) );

X.U{1} = reshape( Ui, sz );

X = orth_at( X, 1, 'left', true );

end

function [left, right] = Afun_prepare( A, x, idx )
    y = A.apply(x);
    if idx == 1
        right = innerprod( x, y, 'RL', idx+1 );
        left = [];
    elseif idx == x.order
        left = innerprod( x, y, 'LR', idx-1 );
        right = [];
    else
        left = innerprod( x, y, 'LR', idx-1 );
        right = innerprod( x, y, 'RL', idx+1 );
    end
end

function res = Afun( A, U, idx, sz, left, right )

    V = reshape( U, sz );
    V = A.apply( V, idx );

    if idx == 1
        tmp = tensorprod_ttemps( V, right, 3 );
    elseif idx == A.order
        tmp = tensorprod_ttemps( V, left, 1 );
    else
        tmp = tensorprod_ttemps( V, right, 3);
        tmp = tensorprod_ttemps( tmp, left, 1);
    end

    res = tmp(:);
end

function B1 = prepare_precond( L0, X, idx )

    if idx == 1
        B1 = [];
        return
    end

    n = size(L0, 1);
    r = X.rank;

    X1 = matricize( X.U{1}, 2);
    Y = X;
    Y.U{1} = tensorize( L0*X1, 2, [r(1), n(1), r(2)] );
    B1 = innerprod( X, Y, 'LR', idx-1);
end

function res = apply_precond( L0, B1, rhs, sz )

    n = size(L0, 1);
    rhs = reshape( rhs, sz );
    if isempty(B1) %idx == 1
        res = L0 \ unfold( rhs, 'left' );
        res = reshape( res, sz );
    else
        res = B1 \ unfold(rhs, 'right');
        res = reshape( res, sz );
    end
    res = res(:);
end

author: NicolasBoumal

repo: manopt


text:
function p = emptyNumericalModel;
p.F_struc = [];
p.K.f = 0;
p.K.l = 0;
p.K.q = 0;
p.K.s = 0;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [net,stats] = cnn_train_dag(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = @solver.sgd; % Empty array - optimised SGD solver
[opts, varargin] = vl_argparse(opts, varargin);
if isempty(opts.solver)
  opts.solverOpts.momentum = 0.9;
else
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver - a function handle with two outputs expected.');
  % A call without any input arg - def opts
  opts.solverOpts = opts.solver();
end

opts.saveSolverState = true ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  state.momentum = cellfun(@gpuArray, state.momentum, 'uniformoutput', false) ;
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, batchSize, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  state.solverState = cellfun(@gather, state.solverState, 'uniformoutput', false) ;
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if isempty(params.solver)
        if isempty(state.solverState{p})
          state.solverState{p} = zeros(size(parDer), 'like', parDer);
        end

        state.solverState{p} = vl_taccum(...
          params.solverOpts.momentum,  state.solverState{p}, ...
          - (1 / batchSize), parDer) ;
        net.params(p).value = vl_taccum(...
          (1 - thisLR * thisDecay / (1 - params.solverOpts.momentum)),  ...
          net.params(p).value, ...
          thisLR, state.solverState{p}) ;
      else
        grad = (1 / batchSize) * parDer + thisDecay * net.params(p).value;
        % call solver function to update weights
        [net.params(p).value, state.solverState{p}] = ...
          params.solver(net.params(p).value, state.solverState{p}, ...
          grad, params.solverOpts, thisLR) ;
      end

    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tflow vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

repo: CREST-Release


text:
function geometry_test044 ( )

%*****************************************************************************80
%
%% TEST044 tests SEGMENTS_INT_1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 7;

  q1_test = [ -1.0, 3.0, 1.0,  0.5, 0.25, 0.5, 2.0 ];
  q2_test = [  1.0, 2.0, 2.0, -3.0, 0.50, 0.5, 2.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST044\n' );
  fprintf ( 1, '  SEGMENTS_INT_1D determines the intersection [R1,R2]\n' );
  fprintf ( 1, '  of line segments [P1,P2] and [Q1,Q2] in 1D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  DIST is negative for overlap,\n' );
  fprintf ( 1, '  0 for point intersection,\n' );
  fprintf ( 1, '  positive if there is no overlap.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '  Test      P1        P2        Q1        Q2        R1        R2        DIST\n' );
  fprintf ( 1, '\n' );

  p1 = -1.0;
  p2 = 1.0;

  for test = 1 : test_num

    q1 = q1_test(test);
    q2 = q2_test(test);

    [ dist, r1, r2 ] = segments_int_1d ( p1, p2, q1, q2 );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %4d  %8f  %8f  %8f  %8f  %8f  %8f  %8f\n', ...
      test, p1, p2, q1, q2, r1, r2, dist );

  end

  return
end