text:
function imgGlare = ComputeGlareImage( img, PSF, hot_pixels_pos, C)
%
%       imgGlare = ComputeGlareImage( img, PSF )
%
%       This function computes the glare image of an input HDR image given
%       a point spread function (PSF) as a kernel.
%
%        Input:
%           -img: an HDR image
%           -PSF: a point spread function stored as a kernel
%           -hot_pixels_pos: hot pixels' coordinates
%
%        Output:
%           -imgGlare: the estimated glare in img.
%
%     Copyright (C) 2014  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%

m = size(hot_pixels_pos, 2);

imgGlare = zeros(size(img));

[r,c,col] = size(img);

[X, Y] = meshgrid(1:c, 1:r);

for i=1:m
    x_p = hot_pixels_pos(1, i);
    y_p = hot_pixels_pos(2, i);

    r = max(sqrt((X-x_p).^2 + (Y-y_p).^2), 2);
    value = C(1) + C(2)./r + C(3)./(r.^2) + C(4)./(r.^3);

    tmp_glare = zeros(size(img));

    for j=1:col
        tmp_glare(:,:,j) = value * img(y_p, x_p, j);
    end

    imgGlare = imgGlare + tmp_glare;

end

%
% hot_pixels_col = zeros(m, size(img, 3));
% for i=1:m
%     hot_pixels_col(i, :) = img(hot_pixels_pos(2, i), hot_pixels_pos(1, i), :);
% end
%
%
% PSF_col = zeros(size(PSF,1), size(PSF, 2), size(img, 3));
% for i=1:size(img, 3)
%     PSF_col(:,:,i) = PSF;
% end
%[imgGlare, ~] = imSplat(size(img, 1), size(img, 2), PSF_col, hot_pixels_pos, hot_pixels_col);

%compensation
while(1)
    ind = find(imgGlare > img);

    if(isempty(ind))
        break
    end

    scale = img(ind(1)) / imgGlare(ind(1));

    imgGlare = imgGlare * scale;
end

end

author: banterle


text:
function quality_test11 ( dim_num, n, z )

%*****************************************************************************80
%
%% TEST11 tests POINTSET_SPACING.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n');
  fprintf ( 1, 'TEST11\n');
  fprintf ( 1, '  POINTSET_SPACING computes pointset spacing parameters.\n');

  gamma = pointset_spacing ( dim_num, n, z );

  gamma_min = min ( gamma(1:n) );
  gamma_max = max ( gamma(1:n) );

  gamma_ave = sum ( gamma(1:n) ) / n;

  if ( 1 < n )
    gamma_std = sqrt ( sum ( ( gamma(1:n) - gamma_ave ).^2 )  / ( n - 1 ) );
  else
    gamma_std = 0.0;
  end

  fprintf ( 1, '\n');
  fprintf ( 1, '  Minimum spacing          GAMMA_MIN = %f\n', gamma_min );
  fprintf ( 1, '  Average spacing          GAMMA_AVE = %f\n', gamma_ave );
  fprintf ( 1, '  Maximum spacing          GAMMA_MAX = %f\n', gamma_max );
  fprintf ( 1, '  Spacing standard dev     GAMMA_STD = %f\n', gamma_std );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function FbVisualize( FB, show )
% Used to visualize a series of 1D/2D/3D filters.
%
% For 1D and 2D filterabnks also shows the Fourier spectra of the filters.
%
% USAGE
%  FbVisualize( FB, [show] )
%
% INPUTS
%  FB      - filter bank to visualize (either 2D, 3D, or 4D array)
%  show    - [1] figure to use for display
%
% OUTPUTS
%
% EXAMPLE
%  FB=FbMake(1,1,0);  FbVisualize( FB, 1 );  %1D
%  load FbDoG.mat;    FbVisualize( FB, 2 );  %2D
%  FB=FbMake(3,1,0);  FbVisualize( FB, 3 );  %3D
%
% See also FBAPPLY2D, FILTERVISUALIZE
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<2 || isempty(show) ); show=1; end
if( show<=0); return; end;

% get Fourier Spectra for 1D and 2D filterbanks
nd = ndims(FB)-1;
if( nd==1 || nd==2 )
  FBF=zeros(size(FB));
  if( nd==1 )
    for n=1:size(FB,1);  FBF(n,:)=abs(fftshift(fft(FB(n,:)))); end
  else
    for n=1:size(FB,3);  FBF(:,:,n)=abs(fftshift(fft2(FB(:,:,n)))); end
  end
end

% display
figure(show); clf;
if( nd==1 )
  r = (size(FB,2)-1)/2;
  subplot(1,3,1); plot( -r:r, FB );
  subplot(1,3,2); plot( (-r:r)/(2*r+1), FBF );
  subplot(1,3,3); stem( (-r:r)/(2*r+1), max(FBF,[],1) );

elseif( nd==2 )
  subplot(1,3,1); montage2(FB);  title('filter bank');
  subplot(1,3,2); montage2(FBF);
  title('filter bank fft');
  subplot(1,3,3); im(sum(FBF,3));  title('filter bank fft coverage');

elseif( nd==3 )
  n = size(FB,4); nn = ceil( sqrt(n) ); mm = ceil( n/nn );
  for i=1:n
    subplot(nn,mm,i);
    filterVisualize( FB(:,:,:,i), 0 );
  end
end

author: zouchuhang

repo: LayoutNet


text:
function mrSetRadioButton(buttonNum,buttonIdList)
%mrSetRadioButton(buttonNum,buttonIdList)
%
%Sets (highlights) one radio button and turns other off
%
%bttonNum:  button number in list to turn on
%buttonIdList:  array of handles to radio buttons.

%7/30/97 gmb  Wrote it.

for i=1:length(buttonIdList)
  if buttonNum == i
    set(buttonIdList(i),'Value',1);
  else
    set(buttonIdList(i),'Value',0);
  end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function opt = ft_setopt(opt, key, val)

% FT_SETOPT assigns a value to an configuration structure or to a cell-array
% with key-value pairs. It will overwrite the option if already present, or
% append the option if not present.
%
% Use as
%   s = ft_setopt(s, key, val)
% where s is a structure or a cell-array.
%
% See also FT_GETOPT, FT_CHECKOPT

% Copyright (C) 2011-2012, Robert Oostenveld
%
% 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 isa(opt, 'struct') || isa(opt, 'config')

  % just replace or add the option
  opt.(key) = val;

elseif isa(opt, 'cell')
  % determine whether the key already exists
  fn = opt(1:2:end);
  sel = find(strcmp(key, fn));
  if isempty(sel)
    % append it
    opt{end+1} = key;
    opt{end+1} = val;
  elseif length(sel)==1
    % replace the current value
    keyindex = 2*sel-1;
    valindex = keyindex+1;
    opt{valindex} = val;
  elseif length(sel)>1
    % first remove all occurences
    keyindex = 2*sel-1;
    valindex = keyindex+1;
    opt([keyindex valindex]) = [];
    % then append it
    opt{end+1} = key;
    opt{end+1} = val;
  end
end % isstruct or iscell

author: spm


text:
function Fitness = CalFitness(PopObj,PopCon)
% Calculate the fitness of each solution

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Kangjia Qiao

    N = size(PopObj,1);
    if nargin == 1
        CV = zeros(N,1);
    else
        CV = sum(max(0,PopCon),2);
    end

    %% Detect the dominance relation between each two solutions
    Dominate = false(N);
    for i = 1 : N-1
        for j = i+1 : N
            if CV(i) < CV(j)
                Dominate(i,j) = true;
            elseif CV(i) > CV(j)
                Dominate(j,i) = true;
            else
                k = any(PopObj(i,:)<PopObj(j,:)) - any(PopObj(i,:)>PopObj(j,:));
                if k == 1
                    Dominate(i,j) = true;
                elseif k == -1
                    Dominate(j,i) = true;
                end
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)
    R = zeros(1,N);
    for i = 1 : N
        R(i) = sum(S(Dominate(:,i)));
    end

    %% Calculate D(i)
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    D = 1./(Distance(:,floor(sqrt(N)))+2);

    %% Calculate the fitnesses
    Fitness = R + D';
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = grad( f )
%GRAD   Numerical gradient of a SEPARABLEAPPROX.
%   [FX FY] = GRAD(F) returns the numerical gradient of the SEPARABLEAPPROX F, where FX
%   is the derivative of F in the x direction and FY is the derivative of F in
%   the y direction. Both derivatives are returned as SEPARABLEAPPROX objects.
%
%   G = GRAD(F) returns a CHEBFUN2V which represents
%
%            G = (F_x ; F_y )
%
%  This command is shorthand for GRADIENT(F).
%
% See also GRADIENT.

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

% Call GRADIENT:
if ( nargout <= 1 )
    out = gradient( f );
    varargout = { out };
else
    [fx, fy] = gradient( f );
    varargout = {fx, fy};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function A = makehatch(hatch)
%MAKEHATCH Predefined hatch patterns
%  MAKEHATCH(HATCH) returns a matrix with the hatch pattern for HATCH
%   according to the following table:
%      HATCH        pattern
%     -------      ---------
%        /          right-slanted lines
%        \          left-slanted lines
%        |          vertical lines
%        -          horizontal lines
%        +          crossing vertical and horizontal lines
%        x          criss-crossing lines
%        .          single dots
%
%  See also: APPLYHATCH

%  By Ben Hinkle, bhinkle@mathworks.com
%  This code is in the public domain.

n = 6;
A=zeros(n);
switch (hatch)
case '/'
  A = fliplr(eye(n));
case '\'
  A = eye(n);
case '|'
  A(:,1) = 1;
case '-'
  A(1,:) = 1;
case '+'
  A(:,1) = 1;
  A(1,:) = 1;
case 'x'
  A = eye(n) | fliplr(diag(ones(n-1,1),-1));
case '.'
  A(1:2,1:2)=1;
otherwise
  error(['Undefined hatch pattern "' hatch '".']);
end

author: jacoxu


text:
function makecomment(filename)
% makecomment(filename)
% Forces user to make commetn and stores it into the FILENAME
% If the file already exist, it will append the comment to the end of the
% file...
comment = input('Make some comments (line string):\n','s');
fid = fopen(filename,'a');
fprintf(fid, '\nComments (%s):\n%s',datestr(now),comment);
fclose(fid);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = process_evt_detect_ecg( varargin )
% PROCESS_EVT_DETECT_ECG: Detect heartbeats in a continuous file, and create set of events called "cardiac"

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, Elizabeth Bock, 2011-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Detect heartbeats';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 43;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsDetect#Detection:_Heartbeats';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Channel name
    sProcess.options.channelname.Comment = 'Channel name: ';
    sProcess.options.channelname.Type    = 'channelname';
    sProcess.options.channelname.Value   = '';
    % Channel name comment
    sProcess.options.channelhelp.Comment = 'You can use the montage syntax here: "ch1, -ch2"';
    sProcess.options.channelhelp.Type    = 'label';
    % Channel name comment
    sProcess.options.channelhelp.Comment = 'You can use the montage syntax here: "ch1, -ch2"';
    sProcess.options.channelhelp.Type    = 'label';
    % Time window
    sProcess.options.timewindow.Comment = 'Time window:';
    sProcess.options.timewindow.Type    = 'timewindow';
    sProcess.options.timewindow.Value   = [];
    % Event name
    sProcess.options.eventname.Comment = 'Event name: ';
    sProcess.options.eventname.Type    = 'text';
    sProcess.options.eventname.Value   = 'cardiac';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Pre-defined set of properties
    %sProcess.options.eventname.Value = 'cardiac';
    sProcess.options.threshold.Value = {4, 'x std', 2};
    sProcess.options.bandpass.Value  = {[10 40], 'Hz', 2};
    sProcess.options.blanking.Value  = {0.5, 'ms', []};
    sProcess.options.isnoisecheck.Value = 1;
    sProcess.options.isclassify.Value   = 1;
    % Call the generic version of the event detection
    OutputFiles = process_evt_detect('Run', sProcess, sInputs);
end

author: brainstorm-tools


text:
function fx = quartic ( x )

%*****************************************************************************80
%
%% QUARTIC evaluates a function defined by a sum of fourth powers.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    R ONeill,
%    Algorithm AS 47:
%    Function Minimization Using a Simplex Procedure,
%    Applied Statistics,
%    Volume 20, Number 3, 1971, pages 338-345.
%
%  Parameters:
%
%    Input, real X(10), the argument.
%
%    Output, real FX, the value of the function.
%
  fx = sum ( x(1:10).^4 );

  return
end


text:
function scale = p02_scale ( neqn )

%*****************************************************************************80
%
%% P02_SCALE returns scale factors for problem p02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end


text:
function [Opt, R, E] = RetroTS(SN)
%    [Opt, OptR, OptE] = RetroTS(Opt)
%This function creates slice-based regressors for regressing out
% components of heart rate, respiration and respiration volume per time.
%
%  Opt is the options structure with the following fields
%     Mandatory:
%     ----------
%     Respfile: Respiration data file
%     Cardfile: Cardiac data file
%     PhysFS: Physioliogical signal sampling frequency in Hz.
%     Nslices: Number of slices
%     VolTR: Volume TR in seconds
%     Optional:
%     ---------
%     Prefix: Prefix of output file
%     SliceOffset: Vector of slice acquisition time offsets in seconds.
%                  (default is equivalent of alt+z)
%     RVTshifts: Vector of shifts in seconds of RVT signal.
%                (default is [0:5:20])
%     RespCutoffFreq: Cut off frequency in Hz for respiratory lowpass filter
%                     (default 3 Hz)
%     CardCutoffFreq: Cut off frequency in Hz for cardiac lowpass filter
%                     (default 3 Hz)
%     ResamKernel: Resampling kernel.
%                 (default is 'linear', see help interp1 for more options)
%     FIROrder: Order of FIR filter. (default is 40)
%     Quiet: [1]/0  flag. (defaut is 1) Show talkative progress as the program runs
%     Demo: [1]/0 flag. (default is 0)
%     RVT_out: [1]/0 flag for writing RVT regressors
%     Card_out: [1]/0 flag for writing Card regressors
%     Resp_out: [1]/0 flag for writing Resp regressors
%     SliceOrder: ['alt+z']/'alt-z'/'seq+z'/'seq-z'/'Custom'/filename.1D
%                 Slice timing information in seconds. The default is
%                 alt+z. See 3dTshift help for more info. 'Custom' allows
%                 the program to use the values stored in the
%                 Opt.SliceOffset array. If a value is placed into the
%                 SliceOrder field other than these, it is assumed to be
%                 the name of a 1D / text file containing the times for
%                 each slice (also in seconds).
%
%Example:
%
%  Opt.Respfile = 'Resp_epiRT_scan_14.dat'
%  Opt.Cardfile = 'ECG_epiRT_scan_14.dat'
%  Opt.VolTR = 2
%  Opt.Nslices = 20;
%  Opt.PhysFS = 1./0.02;   %20 msec sampling period, 50 samples per second
%  RetroTS(Opt);
%

% Input Mode 2 (for testing purposes only):
%  Opt: Scan number for file triplet to be processed.
%      Files called Resp*SN*, ECG*SN*, and scan*SN* are presumed to
%      exist in the directory from which RetroTS is called.
%      Many parameters' value are hard coded to defaults
%
% Output:
%  Opt: Structure of options including default settings.
%

% This option is not to be used because window width calculations do not use it
%     ResampFS: Frequency of resampled signal (default is same as PhysFS)

%Implementation Notes:
%%%%%%%%%%%%%%%%%%%%%%
% The script is intended as a prototype for development in C or Python
% The important routines are:
%    hilbert: Easily implemented with fft and ifft
%    interp1: A table lookup interpolation
%    fir: Tool for designing filters (we can just take it's coefficients)
%    filter: function to apply fir filter parameters (easy)
%
% All of the above can be easily implemented in C. However, I find it
% very useful to be able to plot the various steps in the process as we
% will undoubtedly face problems in the future. So I would vote for
% Python, assuming library vintage is not an issue.
%

if (nargin < 1),
   fprintf(2,'Need some input.\n');
   return;
end

R = struct([]);
E = struct([]);

if (~isstruct(SN)), %mode 1, toy mode
   iscan = 12;

   lll = zglobb({ sprintf('Resp*%d*',iscan),...
                  sprintf('ECG*%d*',iscan),...
                  sprintf('scan*%d*', iscan)});

   %Get some info from header file and set params
   f = fopen(lll(3).name, 'r');
   s = fscanf(f,'%c');
   fclose(f);
   ns = length(s);
   pat = 'RT Physio:\W*sampling\W*';
   Opt.PhysFS = 1000/str2num(strtok(s(regexp(s,pat,'end'):ns)));
   Opt.Nslices = 20;
   Opt.VolTR = 2;
   Opt.SliceMajor = 1;
   Opt.ResampFS = Opt.PhysFS; %upsampling frequency in Hz
   Opt.Quiet = 1;
   Opt.ResamKernel = 'linear';   %resampling filter for envelopes and phase
   Opt.FIROrder = 40;  %order of fir filter
   Opt.RVTshifts = [0:5:20];  %shifts, in seconds, applied to RVT curve
   Opt.Demo = 0;
   Opt.zerophaseoffset = 0;
   Opt.fcutoff = 3; %cut off frequency for filter
   Opt.RespCutoffFreq = 3;
   Opt.CardCutoffFreq = 3;
   Opt.Respfile = lll(1).name;
   Opt.Cardfile = lll(1).name;
   Opt.SliceOffset = ...
      [0:Opt.VolTR./Opt.Nslices:Opt.VolTR-Opt.VolTR./Opt.Nslices];
   Opt.Prefix = sprintf('%d',iscan);
   Opt.SepDups = 0;
   clear ('s');
   clear ('SN');
else,
   Opt = SN; clear ('SN');
   Opt.err = 1; Opt.zerophaseoffset = 0;
   if ( (~isfield(Opt,'Respfile') | isempty(Opt.Respfile))),
      Opt.Respfile = '';
      Opt.Resp_out = 0;
      Opt.RVT_out = 0;
   end
   if ( (~isfield(Opt,'Cardfile') | isempty(Opt.Cardfile))),
      Opt.Cardfile = '';
      Opt.Card_out = 0;
   end
   if ( (~isfield(Opt,'Respfile') | isempty(Opt.Respfile)) & (~isfield(Opt,'Cardfile') | isempty(Opt.Cardfile))),
      fprintf(2,'No Respfile or Cardfile\n');
      return;
   end
   if ( ~isfield(Opt,'PhysFS') | isempty(Opt.PhysFS)),
      fprintf(2,'Missing field PhysFS\n');
      return;
   end
   if ( ~isfield(Opt,'SliceMajor') | isempty(Opt.SliceMajor)),
      Opt.SliceMajor = 1;
   end
   if ( ~isfield(Opt,'Nslices') | isempty(Opt.Nslices)),
      fprintf(2,'Missing field Nslices\n');
      return;
   end
   if ( ~isfield(Opt,'VolTR') | isempty(Opt.VolTR)),
      fprintf(2,'Missing field VolTR\n');
      return;
   end
   if ( ~isfield(Opt,'RVTshifts') | isempty(Opt.RVTshifts)),
      Opt.RVTshifts=[0:5:20];
   end
   if ( ~isfield(Opt,'ResampFS') | isempty(Opt.ResampFS)),
      Opt.ResampFS=Opt.PhysFS;
   end
   if ( ~isfield(Opt,'RespCutoffFreq') | isempty(Opt.RespCutoffFreq)),
      Opt.RespCutoffFreq=3;
   end
   if ( ~isfield(Opt,'CardCutoffFreq') | isempty(Opt.CardCutoffFreq)),
      Opt.CardCutoffFreq=3;
   end
   if ( ~isfield(Opt,'ResamKernel') | isempty(Opt.ResamKernel)),
      Opt.ResamKernel='linear';
   end

   if ( ~isfield(Opt,'FIROrder') | isempty(Opt.FIROrder)),
      Opt.FIROrder=40;
   end
   if ( ~isfield(Opt,'Quiet') | isempty(Opt.Quiet)),
      Opt.Quiet=1;
   end
   if ( ~isfield(Opt,'Demo') | isempty(Opt.Demo)),
      Opt.Demo=0;
   end
   if ( ~isfield(Opt,'Prefix') | isempty(Opt.Prefix)),
      Opt.Prefix = 'oba';
   end
   if ( ~isfield(Opt,'Resp_out') | isempty(Opt.Resp_out)),
      Opt.Resp_out = 1;
   end
   if ( ~isfield(Opt,'Card_out') | isempty(Opt.Card_out)),
      Opt.Card_out = 1;
   end
   if ( ~isfield(Opt,'RVT_out') | isempty(Opt.RVT_out)),
      Opt.RVT_out = 1;
   end
   if ( ~isfield(Opt,'SepDups') | isempty(Opt.SepDups)),
      Opt.SepDups = 0;
   end

   dtt = Opt.VolTR/Opt.Nslices; tt = 0.0;

   % & ~isfield(Opt, 'SliceOffset')
   % & (Opt.SliceOrder ~= 'alt+z')

      % default slice offset times are for alt+z (alternating slice timing)
   if ( ~isfield(Opt,'SliceOffset') | isempty(Opt.SliceOffset))
      Opt.SliceOffset=zeros(Opt.Nslices,1);
   end
   if(~isfield(Opt,'SliceOrder'))
      Opt.SliceOrder = 'alt+z'
   end

   if (isfield(Opt,'SliceOrder'))
      Opt.SliceOffset=zeros(Opt.Nslices,1);
      if(strcmpi(Opt.SliceOrder,'alt+z'))
         for (i=1:2:Opt.Nslices),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
         for (i=2:2:Opt.Nslices),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
      elseif(strcmpi(Opt.SliceOrder, 'alt+z2'))
         for (i=2:2:Opt.Nslices),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
         for (i=1:2:Opt.Nslices),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
      elseif(strcmpi(Opt.SliceOrder, 'seq+z'))
         for (i=1:1:Opt.Nslices),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
      elseif(strcmpi(Opt.SliceOrder,'seq-z'))
         for (i=Opt.Nslices:-1:1),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
      elseif(strcmpi(Opt.SliceOrder,'alt-z'))
         for (i=Opt.Nslices:-2:1),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
         for (i=Opt.Nslices-1:-2:1),
            Opt.SliceOffset(i) = tt; tt = tt+dtt;
         end
      elseif(strcmpi(Opt.SliceOrder,'Custom'))
          % timing already set in Opt.SliceOffset, do nothing
      else
         % read in time offsets from a file (SliceOrder is actually a
         % filename)
         readopt.verb = 0;
         [err, Opt.SliceOffset] = Read_1D(Opt.SliceOrder,readopt);
         if(length(Opt.SliceOffset)~=Opt.Nslices)
            fprintf('Could not read enough slice offsets from file');
            exit(1);
         end
      end
   end
   if(~Opt.Quiet)
      fprintf('Slice timing:'); Opt.SliceOffset
   end
   if ( ~isfield(Opt,'ShowGraphs') | isempty(Opt.ShowGraphs)),
      Opt.ShowGraphs = 1; % show graphs by default
   end
end

if (Opt.SepDups),
   fprintf(1,'WARNING: SepDups should not be used\n');
   fprintf(1,'         It is kept in the code for debugging\n');
   fprintf(1,'         purposes.\n');
end

%create option copy for each type of signal
   OptR = Opt;
      OptR.fcutoff = Opt.RespCutoffFreq;
      OptR.AmpPhase = 1;   %amplitude based phase for respiration
      %OptR.as_percover = 50; %percent overlap of windows for fft
      %OptR.as_windwidth = 0; %window width in seconds for fft, 0 for full window
      %OptR.as_fftwin = 0 ; %1 == hamming window. 0 == no windowing
   OptE = Opt;
      OptE.fcutoff = Opt.CardCutoffFreq;
      OptE.AmpPhase = 0;   %time based phase for cardiac signal

%Get the peaks for R and E
if (~isempty(Opt.Respfile)),
   [R,e]= PeakFinder(Opt.Respfile,OptR);
   if (e), fprintf(2,'Died in PeakFinder\n'); return; end
else
   R = struct([]);
end
if (~isempty(Opt.Cardfile)),
   [E,e] = PeakFinder(Opt.Cardfile,OptE);
   if (e), fprintf(2,'Died in PeakFinder\n'); return; end
else
   E = struct([]);
end
%get the phase
if (~isempty(R)),
   fprintf(2, 'Estimating phase for R\n');
   R = PhaseEstimator(R,OptR);
end
if (~isempty(E)),
   fprintf(2, 'Estimating phase for E\n');
   E = PhaseEstimator(E,OptE);
end

%Now do the RVT for Respiration
if (~isempty(R)),
   fprintf(2,'Computing RVT from peaks\n');
   R = RVT_from_PeakFinder(R, OptR);
end

%Show some results
if(Opt.ShowGraphs)
   if (~isempty(R)),
      fprintf(2, 'Showing RVT Peaks for R\n');
      Show_RVT_Peak(R,1);
   end
end

if (0),
   %write retroicor regressors
   for (i=1:1:Opt.Nslices),
      fname = sprintf('%s.RetroCard.slc%02d.1D', Opt.Prefix, i);
      wryte3(E.phz_slc_reg(:,:,i), fname, 1);
      fname = sprintf('%s.RetroResp.slc%02d.1D', Opt.Prefix, i);
      wryte3(R.phz_slc_reg(:,:,i), fname, 1);
   end

   %and write the RVT puppy, plus or minus a few seconds delay
   fname = sprintf('%s.RetroRVT.1D', Opt.Prefix);
   wryte3(R.RVTRS_slc, fname, 1);
end

%also generate files as 3dREMLfit likes them
nn = 0;
nRv = 0; nRp = 0; nE = 0;
if (~isempty(R)),
   nn = length(R.tst);
   nRp = size(R.phz_slc_reg,2);
   nRv = size(R.RVTRS_slc,2);
end
if (~isempty(E)), %must have E
   nn = length(E.tst); %ok to overwrite length(R.tst), should be same.
   nE = size(E.phz_slc_reg,2);
end

if ( ~Opt.Card_out & ~Opt.Resp_out & ~Opt.RVT_out ),
   fprintf(2, 'Options Card_out, Resp_out, and RVT_out all 0.\nNo output required.\n');
   return;
end
Opt.RemlOut = zeros(  nn,...
                  Opt.Nslices .* ...
                     (  (Opt.RVT_out~=0) .*nRv + ...
                        (Opt.Resp_out~=0).*nRp + ...
                        (Opt.Card_out~=0).*nE ) );
cnt = 0;
head = sprintf([ '# <RetroTSout\n',...
                  '# ni_type = "%d*double"\n'...
                  '# ni_dimen = "%d"\n'...
                  '# ColumnLabels = "'],...
                  size(Opt.RemlOut,2), size(Opt.RemlOut,1) );
tail = sprintf('"\n# >\n');
tailclose = sprintf('# </RetroTSout>\n');

label = head;

if (Opt.SliceMajor == 0), %old approach, not handy for 3dREMLfit
   %RVT
   if (Opt.RVT_out),
      for (j=1:1:size(R.RVTRS_slc,2)),
         for (i=1:1:Opt.Nslices),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = R.RVTRS_slc(:,j); %same for each slice
            label = sprintf('%s s%d.RVT%d ;', label, i-1, j-1);
          end
      end
   end
   %Resp
   if (Opt.Resp_out),
      for (j=1:1:size(R.phz_slc_reg,2)),
         for (i=1:1:Opt.Nslices),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = R.phz_slc_reg(:,j,i);
            label = sprintf('%s s%d.Resp%d ;', label, i-1, j-1);
         end
      end
   end
   %Card
   if (Opt.Card_out),
      for (j=1:1:size(E.phz_slc_reg,2)),
         for (i=1:1:Opt.Nslices),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = E.phz_slc_reg(:,j,i);
            label = sprintf('%s s%d.Card%d ;', label, i-1, j-1);
         end
      end
   end
   fid = fopen(sprintf('%s.retrots.1D', Opt.Prefix),'w');
else
   for (i=1:1:Opt.Nslices),
      if (Opt.RVT_out),
         %RVT
         for (j=1:1:size(R.RVTRS_slc,2)),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = R.RVTRS_slc(:,j); %same regressor for each slice
            label = sprintf('%s s%d.RVT%d ;', label, i-1, j-1);
         end
      end
      if (Opt.Resp_out),
         %Resp
         for (j=1:1:size(R.phz_slc_reg,2)),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = R.phz_slc_reg(:,j,i);
            label = sprintf('%s s%d.Resp%d ;', label, i-1, j-1);
         end
      end
      if (Opt.Card_out),
         %Card
         for (j=1:1:size(E.phz_slc_reg,2)),
            cnt = cnt + 1;
            Opt.RemlOut(:,cnt) = E.phz_slc_reg(:,j,i);
            label = sprintf('%s s%d.Card%d ;', label, i-1, j-1);
         end
      end
   end
   fid = fopen(sprintf('%s.slibase.1D', Opt.Prefix),'w');
end

%remove very last ';'
label = label(1:end-1);

fprintf(fid,'%s',label);
fprintf(fid,'%s ',tail);
for(i=1:1:size(Opt.RemlOut,1)),
   fprintf(fid,'%g ', Opt.RemlOut(i,:));
   fprintf(fid,'\n ');
end
fprintf(fid,'%s',tailclose);
fclose(fid);

Opt.err = 0;

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ierror = part_rsf_check ( n, npart, a )

%*****************************************************************************80
%
%% PART_RSF_CHECK checks a reverse standard form partition of an integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPART, the number of parts of the partition.
%    1 <= NPART <= N.
%
%    Input, integer A(NPART), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.  The entries must be in ASCENDING order.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    -1, N is illegal.
%    -2, NPART is illegal.
%    -3, the entries do not add up to N.
%    I, the I-th entry of A is illegal.
%
  ierror = 0;

  if ( n < 1 )
    ierror = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PART_RSF_CHECK - Fatal error!\n' );
    fprintf ( 1, '  N < 1.\n' );
    error ( 'PART_RSF_CHECK - Fatal error!\n' );
  end

  if ( npart < 1 || n < npart )
    ierror = -2;
    return
  end
%
%  Every entry must lie between 1 and N.
%
  for i = 1 : npart
    if ( a(i) < 1 || n < a(i) )
      ierror = i;
      return
    end
  end
%
%  The entries must be in ascending order.
%
  for i = 2 : npart
    if ( a(i) < a(i-1) )
      ierror = i;
      return
    end
  end
%
%  The entries must add up to N.
%
  if ( sum ( a(1:npart) ) ~= n )
    ierror = -3;
  end

  return
end


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

function particle_system = smc_resample(particle_system)
% SMC_RESAMPLE: Function resamples binary particles (or models) based on
% the effective sample size at the current iteration.

% Extract inputs from particle_system
models     = particle_system.models;
model_val  = particle_system.model_val;
weights    = particle_system.weights;

% Generate new particles and data
[new_models, new_idx] = binary_resampling(models, weights);
n_models = size(new_models,1);

% Assign uniform weights
new_weights  = 1/n_models*ones(n_models,1);

% Re-order model data
model_val  = model_val(new_idx);

% Save models and weights
particle_system.models     = new_models;
particle_system.model_val  = model_val;
particle_system.weights    = new_weights;

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


function y = cf_RN_NIG( u,r,T,alpha,beta,delta)
%
%
asq = alpha^2;
bsq = beta^2;
temp = sqrt(asq-bsq);
y = -delta*(sqrt(asq - (beta +1i*u).^2) - temp);  %Psi_s
RNmu = r + delta*(sqrt(asq - (beta+1)^2)-temp);
y = exp(T*(1i*u*RNmu + y));

end


text:
function [ t_start, y_start ] = p28_start ( neqn )

%*****************************************************************************80
%
%% P28_START returns the starting point for problem p28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = [ 0.0; 0.0 ];

  return
end


text:
function spectrumscope_ex3
% SPECTRUMSCOPE_EX3   Spectrum scope example, using the Data Acquisition Toolbox
%
% See comments on bottom of file for extra notes on this application

%    Scott Hirsch 2-25-04
%    shirsch@mathworks.com
%    Copyright 1998-2004 The MathWorks, Inc.

%%
% Create analoginput object (Data Acquisition Toolbox)
daqreset
ai = analoginput('winsound');
chan = addchannel(ai,1:2);

%%
% Configure analoginput object
Fs = 44100;     % Sample Rate
set(ai,'SampleRate',Fs);
set(ai,'SamplesPerTrigger',4096);

%%
% Initialize scope
Nfft = ai.SamplesPerTrigger/2;                  % FFT Length.
NTraces = length(chan);                         % Number of traces on scope
hSpectrumScope = spectrumscope(Fs,Nfft,NTraces);% Initialize spectrum scope

%%
% Configure acquisition to update scope
set(ai,'TriggerRepeat',inf)
set(ai,'TriggerFcn',{@updateplot,hSpectrumScope,Nfft});

% Put a stop button on the scope
% This is available on MATLAB Central
if exist('daqstopbutton','file')
    daqstopbutton(gcf,ai);
else
    start(ai);
    disp('Type "stop(daqfind)" to stop the application.  Then, go to MATLAB Central and download daqstopbutton.m!');
end;

function updateplot(ai,event,hSpectrumScope,Nfft)
d = peekdata(ai,Nfft);
spectrumscope(hSpectrumScope,d);

%%
% Comments:
%
% For an application that simple gives me some real-time visualization of
% the data that I am acquiring, I prefer to use PEEKDATA instead of
% GETDATA.  Since GETDATA always grabs the oldest data in the buffer, there
% can easily be a significant lag from when an event occurs to when the
% response to the event is shown on the screen. PEEKDATA grabs the most
% recently acquired data, so it makes for much more responsive
% applications.  PEEKDATA is also non-blocking, so it won't slow down your
% application while waiting for data.  This means that you need to make
% sure there is enough data in the buffer before you call PEEKDATA.  This
% is why I set the amount of data I grab per call to peekdata (Nfft) to
% something less than the amount of data acquired per trigger.
%
% There are a bunch of different callback functions for running continuous
% applications.  I happen to like using Triggers, but you could just as
% easily use TimerFcn or SamplesAcquiredFcn

author: Sable

repo: mcbench-benchmarks


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

disp('Ex0501: Perform FOE Alignment')

SpharmMatDir = '../';
CodeDir = [SpharmMatDir 'code'];
addpath(CodeDir);

DataFolder = 'data';

%% (1) Input data directory
InDataDir = [SpharmMatDir DataFolder '/Ex0501/hip04_des'];

%% (2) List input file names
inFiles = dir([InDataDir '/*_des.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [InDataDir '/' inFiles(i).name];
end

%% (3) Display input objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'none';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

%% (4) Perform FOE Alignment
    % Available values for CPoint- 'x';'y';'z'
confs.CPoint = 'y';
    % Available values for NPole- 'x';'y';'z'
confs.NPole = 'z';
confs.MaxSPHARMDegree = 12;
confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0501/hip05_reg'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

outNames = SpharmMatAlignment(confs, inNames, 'AligFOE');

%% (5) Display output objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'quad64';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'adc_paramap';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = 1;
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, outNames, CodeDir);

dispConfs.Degree = 12;
SpharmMatUtilDisplayObjs(dispConfs, outNames, CodeDir);

rmpath(CodeDir);

clear;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function epsilon = CauchyStrain(F)
% Cauchy strain or small strain
%
% Syntax
%   epsilon = CauchyStrain(F)
%
% Input
%  F - @deformationTensor
%
% Output
%  epsilon - @strainTensor
%

epsilon = strainTensor(F.sym - tensor(eye(3),'rank',2));

author: mtex-toolbox

repo: mtex

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


text:
function [paramsOut] = postProcessGUI(hObject, callbackdata, inputData) %#ok<INUSL>
title ='Post process parameters';
defaultStruct = inputData.userData.postProcess;
mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
userdata = get(mainFigure, 'UserData');
if isempty(userdata) || ~isfield(userdata, 'postProcess')
    paramsOut = struct();
else
    paramsOut = userdata.postProcess;
end
[defaultStruct, errors] = checkStructureDefaults(paramsOut, defaultStruct);

if ~isempty(errors)
    warning('postProcessGUI:bad parameters', getMessageString(errors)); %#ok<CTPCT>
end
if defaultStruct.cleanupReference.value
    cleanUpCheckValue = 1;
else
    cleanUpCheckValue = 0;
end
if defaultStruct.keepFiltered.value
    filteredCheckValue = 1;
else
    filteredCheckValue = 0;
end
if defaultStruct.removeInterpolatedChannels.value
    removeCheckValue = 1;
else
    removeCheckValue = 0;
end
closeOpenWindows(title);
geometry = {[1, 1], [1, 1], [1, 1]};
geomvert = [];
uilist = {{'style', 'text', 'string', 'Clean up reporting fields', 'TooltipString', ''}...
    {'style', 'checkbox', 'tag', 'cleanupReference', 'Value', cleanUpCheckValue}...
    {'style', 'text', 'string', 'Keep filtered', 'TooltipString', ''}...
    {'style', 'checkbox', 'tag', 'keepFiltered', 'Value', filteredCheckValue}...
    {'style', 'text',  'string', 'Remove bad interpolated channels'}...
    {'style', 'checkbox', 'tag', 'removeInterpolatedChannels', 'Value', removeCheckValue}};
[~, ~, ~, paramsOut] = inputgui('geometry', geometry, ...
    'geomvert', geomvert, 'uilist', uilist, 'title', title, ...
    'helpcom', 'pophelp(''pop_prepPipeline'')');
if (isempty(paramsOut))  % Cancelled out
    return;
end
mainFigure = findobj('Type', 'Figure', '-and', 'Name', inputData.name);
userdata = get(mainFigure, 'UserData');
if paramsOut.cleanupReference
    userdata.postProcess.cleanupReference = true;
else
    userdata.postProcess.cleanupReference = false;
end
if paramsOut.keepFiltered
    userdata.postProcess.keepFiltered = true;
else
    userdata.postProcess.keepFiltered = false;
end
if paramsOut.removeInterpolatedChannels
    userdata.postProcess.removeInterpolatedChannels = true;
else
    userdata.postProcess.removeInterpolatedChannels = false;
end
set(mainFigure, 'UserData', userdata);
end

author: VisLab

repo: EEG-Clean-Tools


text:
% This demo Will show how to reconstruct real data.
%
% We will reconstruct the so-called SophiaBeads dataset, where a bunch of
% "balls" are put in a Nikon ?-CT machine.
%
% The data can be downloaded at:
% https://zenodo.org/record/16474
%
% And the codes can be downloaded at:
% http://zenodo.cern.ch/record/16539
%
% The SophiaBeads codes and dataset ARE NOT part of TIGRE toolbox, and are
% not licensed under the BSD license, be aware.
%
% Additionally, if you use the dataset or the codes, they have each of them
% their own bilbiography, that you need to reference, as per their own
% requirements.

% cite as:

% citeme('SophiaBeads data')
% citeme('SophiaBeads code')

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------

%% How to run this code:
%
% 1.- Download both the code and the data (any size you prefer).
%
% 2.- Open SophiaBeads.m in the codes and replace pathname and filename (in
%     lines 13-14 by the data path and name you chosen
%
% 3.- In line 27 add the followinh code:     break;
%
% 4.- Run SophiaBeads.m
%
% 5.- Without clearing the variables, come to this code and run it.
%
%

% from sophiabeads to TIGRE
geo.DSD = geom.d_sd;
geo.DSO = abs(geom.source.x);

% Detector parameters
geo.nDetector=[geom.dets.ny;geom.dets.nz];
geo.dDetector=[mean(diff(geom.dets.y)); mean(diff(geom.dets.z))];
geo.sDetector=geo.nDetector.*geo.dDetector;

% Image parameters
geo.nVoxel=geom.voxels.';  % this is the datasheets standard size, but its
% very big size and most ocmputers will not be able to run it. We redefine
% it in the next line, but feel free to use the original
geo.nVoxel=[500 500 200].';
geo.sVoxel=[geom.voxel_size.*geom.voxels].';
geo.dVoxel=geo.sVoxel./geo.nVoxel;

% Offsets
geo.offOrigin=[0,0,0].';
geo.offDetector=[0; 0];

% Auxiliary
geo.accuracy=0.5;
geo.mode='cone';
geo.COR=geom.source.y;

% Angles
angles=-geom.angles.';

data=reshape(data,1564,1564,size(angles,2));
data=permute(data,[2 1 3]);

%% Now use any algorithms you want

% example:
% beadsSART=SART(data,geo,angles,15);
beadsCGLS=CGLS(data,geo,angles,15);
ImgASD_POCS_20=OS_ASD_POCS(data,geo,angles,50,...
'TViter',2,'maxL2err',0.5,'alpha',0.002,'alpha_red',0.2,... % these are very important
'lambda',1,'lambda_red',0.99,'Ratio',1,'Verbose',true);
% plotImg([beadsSART beadsCGLS],'Dim','Z','Step',1)

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
% WAITBAR Display wait bar.
%    H = WAITBAR(X,'message', property, value, property, value, ...)
%    creates and displays a waitbar of fractional length X.  The
%    handle to the waitbar figure is returned in H.
%    X should be between 0 and 1.  Optional arguments property and
%    value allow to set corresponding waitbar figure properties.
%    Property can also be an action keyword 'CreateCancelBtn', in
%    which case a cancel button will be added to the figure, and
%    the passed value string will be executed upon clicking on the
%    cancel button or the close figure button.
%
%    WAITBAR(X) will set the length of the current bar to the fractional
%    length X.
%
%    WAITBAR(X,H) will set the length of the bar in waitbar H
%    to the fractional length X.
%
%    WAITBAR(X,H,'message') will update the message text in
%    the waitbar figure, in addition to setting the fractional
%    length to X.
%
%    WAITBAR is typically used inside a FOR loop that performs a
%    lengthy computation.
%
%    Example:
%        h = waitbar(0,'Please wait...');
%        for i=1:1000,
%            % computation here %
%            waitbar(i/1000,h)
%        end
%
%    See also DIALOG, MSGBOX.
%
%    Reference page in Doc Center
%       doc waitbar
%
%

author: jmaih


text:
function F = ml_Fmeasure(cluster_labels, class_labels)
%ML_FMEASURE Compute the Fmeasure for clustering with external validity
%
%
%   http://www.uni-weimar.de/medien/webis/teaching/lecturenotes/machine-learning/unit-en-cluster-analysis-evaluation.pdf
%   slide number 23-26
%
%   Macro average F-measure
%
%   Fm is an N x M matrix (N number of classes, M number of clusters) where
%   Fm(i,j) is the F-measure of the jth cluster computed with respect to
%   the ith class
%
%
%   input -----------------------------------------------------------------
%
%       - cluster_labels : N x 1 vector of cluster labels given by the clustering algorithm
%       (N is the size of data)
%
%       - class_labels : N x 1 vector of class labels (ground truth)
%
%   output ----------------------------------------------------------------
%
%       - F = macro averaged F-measure

nbOfClusters = length(unique(cluster_labels));
nbOfClasses  = length(unique(class_labels));
sizeClasses  = zeros(nbOfClasses,1);

Fm = zeros(nbOfClasses,nbOfClusters);

unique_cluster_labels = unique(cluster_labels);
unique_class_labels   = unique(class_labels);

% for i = 1:nbOfClasses
%     for j = 1:nbOfClusters
%         N=0;
%         for k = 1:length(cluster_labels)
%             if(cluster_labels(k) == j && class_labels(k) == i)
%                 N = N+1;
%             end
%         end
%         Precision = N/length(find(cluster_labels==j));
%         Recall = N/length(find(class_labels==i));
%         Fm(i,j) = 2*Precision*Recall/(Precision+Recall);
%     end
% end

for i = 1:nbOfClasses
    for j = 1:nbOfClusters
        N=0;
        for k = 1:length(cluster_labels)
            if(cluster_labels(k) == unique_cluster_labels(j) && class_labels(k) == unique_class_labels(i))
                N = N+1;
            end
        end
        Precision = N/length(find(cluster_labels==unique_cluster_labels(j)));
        Recall = N/length(find(class_labels==unique_class_labels(i)));
        Fm(i,j) = 2*Precision*Recall/(Precision+Recall);
    end
end

for m = 1:nbOfClasses
    sizeClasses(m,:) = length(find(class_labels==unique_class_labels(m)));
end

F = sum(max(Fm,[],2).*sizeClasses)/sum(sizeClasses);

end

author: epfl-lasa


text:
%
function a = fitellipse(X,Y)

% FITELLIPSE  Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE(X,Y) returns the parameters of the best-fit
%        ellipse to 2D points (X,Y).
%        The returned vector A contains the center, radii, and orientation
%        of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)
%
% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
%
%  @Article{Fitzgibbon99,
%   author = "Fitzgibbon, A.~W.and Pilu, M. and Fisher, R.~B.",
%   title = "Direct least-squares fitting of ellipses",
%   journal = pami,
%   year = 1999,
%   volume = 21,
%   number = 5,
%   month = may,
%   pages = "476--480"
%  }
%
% This is a more bulletproof version than that in the paper, incorporating
% scaling to reduce roundoff error, correction of behaviour when the input
% data are on a perfect hyperbola, and returns the geometric parameters
% of the ellipse, rather than the coefficients of the quadratic form.
%
%  Example:  Run fitellipse without any arguments to get a demo
if nargin == 0
  % Create an ellipse
  t = linspace(0,2);

  Rx = 300;
  Ry = 200;
  Cx = 250;
  Cy = 150;
  Rotation = .4; % Radians

  NoiseLevel = .5; % Will add Gaussian noise of this std.dev. to points

  x = Rx * cos(t);
  y = Ry * sin(t);
  nx = x*cos(Rotation)-y*sin(Rotation) + Cx + randn(size(t))*NoiseLevel;
  ny = x*sin(Rotation)+y*cos(Rotation) + Cy + randn(size(t))*NoiseLevel;

  % Clear figure
  clf
  % Draw it
  plot(nx,ny,'o');
  % Show the window
  figure(gcf)
  % Fit it
  params = fitellipse(nx,ny);
  % Note it may return (Rotation - pi/2) and swapped radii, this is fine.
  Given = round([Cx Cy Rx Ry Rotation*180]);
  Returned = round(params.*[1 1 1 1 180]);

  % Draw the returned ellipse
  t = linspace(0,pi*2);
  x = params(3) * cos(t);
  y = params(4) * sin(t);
  nx = x*cos(params(5))-y*sin(params(5)) + params(1);
  ny = x*sin(params(5))+y*cos(params(5)) + params(2);
  hold on
  plot(nx,ny,'r-')

  return
end

% normalize data
mx = mean(X);
my = mean(Y);
sx = (max(X)-min(X))/2;
sy = (max(Y)-min(Y))/2;

x = (X-mx)/sx;
y = (Y-my)/sy;

% Force to column vectors
x = x(:);
y = y(:);

% Build design matrix
D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

% Build scatter matrix
S = D'*D;

% Build 6x6 constraint matrix
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

% Solve eigensystem
if 0
  % Old way, numerically unstable if not implemented in matlab
  [gevec, geval] = eig(S,C);

  % Find the negative eigenvalue
  I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));

  % Extract eigenvector corresponding to negative eigenvalue
  A = real(gevec(:,I));
else
  % New way, numerically stabler in C [gevec, geval] = eig(S,C);

  % Break into blocks
  tmpA = S(1:3,1:3);
  tmpB = S(1:3,4:6);
  tmpC = S(4:6,4:6);
  tmpD = C(1:3,1:3);
  tmpE = inv(tmpC)*tmpB';

  % Return if tmpE contains nan
  if sum(sum(isnan(tmpE))) ~= 0
      a = [];
      return;
  end
  [evec_x, eval_x] = eig(inv(tmpD) * (tmpA - tmpB*tmpE));

  % Find the positive (as det(tmpD) < 0) eigenvalue
  I = find(real(diag(eval_x)) < 1e-8 & ~isinf(diag(eval_x)));

  % Extract eigenvector corresponding to negative eigenvalue
  A = real(evec_x(:,I));

  % Recover the bottom half...
  evec_y = -tmpE * A;
  A = [A; evec_y];
end

% unnormalize
par = [
  A(1)*sy*sy,   ...
      A(2)*sx*sy,   ...
      A(3)*sx*sx,   ...
      -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
      -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
      A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
      - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
      + A(6)*sx*sx*sy*sy   ...
      ]';

% Convert to geometric radii, and centers

thetarad = 0.5*atan2(par(2),par(1) - par(3));
cost = cos(thetarad);
sint = sin(thetarad);
sin_squared = sint.*sint;
cos_squared = cost.*cost;
cos_sin = sint .* cost;

Ao = par(6);
Au =   par(4) .* cost + par(5) .* sint;
Av = - par(4) .* sint + par(5) .* cost;
Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;
Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;

% ROTATED = [Ao Au Av Auu Avv]

tuCentre = - Au./(2.*Auu);
tvCentre = - Av./(2.*Avv);
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

uCentre = tuCentre .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;

Ru = -wCentre./Auu;
Rv = -wCentre./Avv;

Ru = sqrt(abs(Ru)).*sign(Ru);
Rv = sqrt(abs(Rv)).*sign(Rv);

a = [uCentre, vCentre, Ru, Rv, thetarad];

repo: Fall-Detection


text:
function h = p08_h ( n, x )

%*****************************************************************************80
%
%% P08_H evaluates the Hessian for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  ap = 0.00001;

  t1 = - 0.25 + sum ( x(1:n).^2 );

  d1 = 2.0 * ap;
  th = 4.0 * t1;

  for i = 1 : n
    h(i,i) = d1 + th + 8.0 * x(i)^2;
    for j = 1 : i - 1
      h(i,j) = 8.0 * x(i) * x(j);
    end
  end

  for i = 1 : n
    for j = i + 1 : n
      h(i,j) = h(j,i);
    end
  end

  return
end

author: johannesgerer


text:
function v = setLayerInputs(obj, layer, inputs)
%SETLAYERINPUTS  Set or change the inputs to a layer
%   Example: NET.SETLAYERINPUTS('layerName', {'input1', 'input2', ...})

v = [] ;
l = obj.getLayerIndex(layer) ;
for input = inputs
  v(end+1) = obj.addVar(char(input)) ;
end
obj.layers(l).inputs = inputs ;
obj.rebuild() ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function[I]=besselitilde(varargin)
%BESSELITILDE  I-type Bessel function after factoring off exponential growth.
%
%   ITILDE=BESSELITILDE(NU,Z) returns the modified Bessel function of the
%   first kind of order NU at argument Z, after factoring off the
%   asymptotic behavior of EXP(Z).
%
%   BESSELITILDE is useful for products of modified Bessel functions in
%   which the exponential behaviors cancel, but that cannot be evaluated
%   directly because of numerical overflow.
%
%   BESSELITILDE(NU,Z,N) use a summation truncated at N terms. The default
%   behavior uses N=30 and is highly accurate.
%
%   See 10.40.1 of https://dlmf.nist.gov/10.40.
%
%   This is low-level code used by WINDTRANS using an algorithm described
%   in Lilly and Elipot (2021).
%
%   See also BESSELKTILDE.
%
%   'besselitilde --t' runs a test.
%
%   Usage: I=besselitilde(nu,z);
%          I=besselitilde(nu,z,nterms);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2019--2021 J.M. Lilly --- type 'help jlab_license' for details

if strcmp(varargin{1}, '--t')
    besselitilde_test,return
end

nu=varargin{1};
z=varargin{2};
nterms=30;
if nargin==3
    nterms=varargin{3};
end

sizez=size(z);
z=z(:);
%let's sum over the third dimensions
zk=vrep(-z,nterms,2);  %this makes the alternating sign
zk(:,1)=1;
zk=cumprod(zk,2);

k=(0:nterms-1)';
ak=(4*nu.^2-(2*k-1).^2);
ak(1)=1;
ak=cumprod(ak,1);
ak=frac(ak,factorial(k).*8.^k);
if any(~isfinite(ak))
    ak=ak(1:find(~isfinite(ak),1,'first')-1);
end
zk=zk(:,1:length(ak));

I=frac(1,sqrt(2*z*pi)).*((1./zk)*ak);
I=reshape(I,sizez);
%nterms

% %let's sum over the third dimensions
% zk=vrep(-z,nterms,3);  %this makes the alternating sign
% zk(:,:,1)=1;
% zk=cumprod(zk,3);
%
% k=(0:nterms-1)';
% ak=(4*nu.^2-(2*k-1).^2);
% ak(1)=1;
% ak=cumprod(ak,1);
% ak=frac(ak,factorial(k).*8.^k);
% if any(~isfinite(ak))
%     ak=ak(1:find(~isfinite(ak),1,'first')-1);
% end
% zk=zk(:,:,1:length(ak));
%
% ak=vrep(permute(ak,[3 2 1]),size(z),[1 2]);
%
% I=frac(1,sqrt(2*z*pi)).*sum(frac(ak,zk),3);

function[]=besselitilde_test

for s=[1 -1]
    z=sqrt(s*1i)*[23:0.01:100]';
    %z=[20:0.01:100]';
    [bi0,bi]=vzeros(length(z),2);
    for i=1:2
        bi0(:,i)=exp(-z).*besseli(i-1,z);
        bi(:,i)=besselitilde(i-1,z);
    end

    if s==1
        reporttest('BESSELITILDE for z with phase of pi/4',allall(abs((bi0-bi)./bi)<1e-14))
    else
        reporttest('BESSELITILDE for z with phase of -pi/4',allall(abs((bi0-bi)./bi)<1e-14))
    end

	bi=vzeros(length(z),2);
    for i=1:2
        bi(:,i)=besselitilde(i-1,z,2);
    end

    if s==1
        reporttest('BESSELITILDE 2-term for z with phase of pi/4, order 0 and 1',allall(abs((bi0(:,1:2)-bi)./bi)<1e-3))
    else
        reporttest('BESSELITILDE 2-term for z with phase of -pi/4, order 0 and 1',allall(abs((bi0(:,1:2)-bi)./bi)<1e-3))
    end

    bi0=sqrt(frac(1,2*pi*z)).*(1+frac(1,8*z));
    bi1=sqrt(frac(1,2*pi*z)).*(1-frac(3,8*z));

    if s==1
        reporttest('BESSELITILDE 2-term vs. analytic for z with phase of pi/4',allall(abs(([bi0 bi1]-bi)./bi)<1e-15))
    else
        reporttest('BESSELITILDE 2-term  vs. analytic for z with phase of -pi/4',allall(abs(([bi0 bi1]-bi)./bi)<1e-15))
    end
end

% figure,plot(abs(z),abs(bi-bi0)),ylog

% nu=5;
% z=[0:0.001:100]';
% tic;bi(:,1)=besselitilde(nu,z);toc
% tic;bi(:,2)=sqrt(2*z*pi).*exp(-z).*besseli(nu,z);toc
% figure,plot(z,log10(abs((bi(:,1)-bi(:,2))./bi(:,2))))
% %wow, really excellent
%
%
% tic;bi0=sqrt(2*z*pi).*exp(-z).*besseli(nu,z);toc
% for i=1:50
%     bi(:,i)=besselitilde(nu,z,i);
% end
% bi0=vrep(bi0,size(bi,2),2);
% plot(z,log10(abs((bi-bi0)./bi0)))
%
%

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout = dcm2metaimage( str, res, scale, crop, ext, file )
% DCM2METAIMAGE  Read a batch of DICOM files, collate them and save as a
% single MetaImage volume (one .mha and one .raw file)
%
% DCM2METAIMAGE(STR)
%
%   STR is a string with the path and output mha header file name. For
%   example,
%
%     STR='/home/john/data/study01/study01.mha'; % linux
%     STR='C:\data\study01\study01.mha';         % windows
%
%   It is assumed that the MetaImage file (.mha) will be created from the
%   DICOM files found in the target directory.
%
% DCM2METAIMAGE(STR, RES, SCALE, CROP, EXT, FILE)
%
%   RES is a 3-vector with the pixel size in the x-, y- and z-coordinates.
%   If RES is not provided, the function tries to read the values from the
%   DICOM header. If the resolution values cannot be read, by default
%   RES=[1 1 1].
%
%   SCALE is a factor to reduce the size of each frame (the number of
%   frames doesn't change). By default, SCALE=1.
%
%   CROP is a vector with values [xmin, xmax, ymin, ymax, zmin, zmax].
%   These values allow to crop the data volume. It follows the same
%   convention as Seg3D [1]. For example, if xmin=0, xmax=3, then voxels 1,
%   2 and 3 are selected.
%   Note that the output MetaImage file will have an offset header so that
%   the data retains its true coordinates.
%
%   EXT is a string with the file extension of the images. By default,
%   EXT='tif'.
%
%   FILE is the list of image files. By default, all the files in the
%   target directory are read, but FILE allows to change the system listing
%   order, select a subset of files, etc. This is useful, for example, when
%   the system listing order does not correspond to the order in which you
%   want to collate the slices. FILE is expected to have the same struct
%   format as the output of function DIR().
%
%   If you have filenames like
%
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.10.2009.10.23.16.02.38.562500.9699412.DCM
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.11.2009.10.23.16.02.38.562500.9699451.DCM
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.1.2009.10.23.16.02.38.562500.9699056.DCM
%    MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.12.2009.10.23.16.02.38.562500.9699490.DCM
%
%   where the slice order is given by the field between the "6" and the
%   "2009", then you can use Gerardus function sortdirbynum() to order the
%   slices correctly. Example:
%
%   >> str1 = 'MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.*.DCM';
%   >> str2 = 'MISAS_DT_22_10_2009.MR.HEAD_GENERAL.6.mha';
%   >> file = sortdirbynum( str1, 5, '.' );
%   >> dcm2metaimage( str2, [], [], [], 'DCM', file );
%
% IM = DCM2METAIMAGE(...);
%
%   IM is the image volume.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2009-2010 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error( nargchk( 1, 6, nargin, 'struct' ) );
error( nargoutchk( 0, 1, nargout, 'struct' ) );

% defaults
% if ( nargin < 2 || isempty( res ) )
%     res = [ 1, 1, 1 ];
% end
% "res" default set below
if ( nargin < 3 || isempty( scale ) )
    scale = 1;
end
if ( nargin < 4 || isempty( crop ) )
    crop = [];
end
if ( nargin < 5 || isempty( ext ) )
    ext = 'tif';
end
if ( nargin < 6 || isempty( file ) )
    file = [];
end

% check arguments
if ( scale <= 0 || scale > 1 )
    error( 'SCALE must be in (0, 1]' );
end

% get path to the data files
[ dirdata, name ] = fileparts( str );
if (isempty(dirdata))
    dirdata = '.';
end

% get list of image files
if isempty( file )
    file = dir( [ dirdata filesep '*.' ext ] );
end

% if resolution is not provided, try to read it from DICOM file
if ( nargin < 2 || isempty( res ) )
    % get DICOM header
    info = dicominfo( [ dirdata filesep file( 1 ).name ] );

    % if resolution is in the DICOM header, read it...
    if ( isfield( info, 'SliceThickness' ) && isfield( info, 'PixelSpacing' ) )
        res = [ info.PixelSpacing' info.SliceThickness ];
    else % ... if not, then fall back to default resolution
        res = [ 1, 1, 1 ];
    end
end

% adjust the image resolution to the scaling factor
res( 1:2 ) = res( 1:2 ) / scale;

% load first slice
frame = dicomread( [ dirdata filesep file( 1 ).name ] );

% remove alpha channel or other layers if present
frame = frame( :, :, 1 );

% init volume to load data
im = zeros( [ size( frame )*scale length( file ) ], class( frame ) );

% load slices of whole volume
im( :, :, 1 ) = imresize( frame, scale, 'bilinear' );
for I = 2:length( file )
    frame = dicomread( [ dirdata filesep file( I ).name ] );
    frame = frame( :, :, 1 );
    % resize frames
    im( :, :, I ) = imresize( frame, scale, 'bilinear' );
end

% crop volume (this is less efficient than initialising the volume as
% cropped, and then cropping each frame, but it makes the code simpler)
%
% The "+1" is necessary to follow the Seg3D convention
if ~isempty( crop )
    im = im( ...
        crop( 1 )+1:crop( 2 ), ...
        crop( 3 )+1:crop( 4 ), ...
        crop( 5 )+1:crop( 6 ) ...
        );
end

% compute offset in the coordinates of the cropped volume
if isempty( crop )
    offset = [ 0.0 0.0 0.0 ];
else
    offset = crop( [ 1 3 5 ] ) .* res;
end

% write mha file
writemetaimagefile([dirdata filesep name '.mha'], im, res, offset);

% avoid outputing the image volume unless the user has requested it
% explicitly
if ( nargout == 0 )
    varargout = [];
else
    varargout = { im };
end


text:
function [dat, dimord] = read_shm_data(hdr, chanindx, begtrial, endtrial)

% READ_SHM_DATA reads the data in real-time from shared memory
% this is a helper function for READ_DATA

% Copyright (C) 2007-2009, Robert Oostenveld
%
% 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$

% this persistent variable is used for caching the 600 packets
% which inproves the throughput when reading overlapping data segments
persistent ctf_shm

% read the data from shared memory, first the meta information only
[msgType, msgId, sampleNumber, numSamples, numChannels] = read_ctf_shm;

if isempty(ctf_shm)
  ctf_shm.msgType      = nan(size(msgType));
  ctf_shm.msgId        = nan(size(msgId));
  ctf_shm.sampleNumber = nan(size(sampleNumber));
  ctf_shm.numSamples   = nan(size(numSamples));
  ctf_shm.numChannels  = nan(size(numChannels));
  ctf_shm.data = {};
  ctf_shm.hit   = 0;
  ctf_shm.mis   = 0;
end

% there seems to be a bug in Acq, causing the messageId to wrap around
% hence it cannot be used as index into the packets, so construct a new trial numbering vector
trlNum = nan(size(msgId));
trlNum(msgType==1) = sampleNumber(msgType==1)./numSamples(msgType==1) + 1;

% allocate memory for the data, fill with NaNs
dat = nan(length(chanindx), hdr.nSamples, endtrial-begtrial+1);

% determine which trials/packets to read
sel = find((trlNum>=begtrial) & (trlNum<=endtrial));

% this is for calibrating the integer values to physical values
if isfield(hdr, 'orig') && isfield(hdr.orig, 'gainV')
  gain = diag(hdr.orig.gainV(chanindx));
elseif isfield(hdr, 'orig') && isfield(hdr.orig, 'res4')
  gain = diag([hdr.orig.res4.senres(chanindx).gain]);
end

if length(chanindx)>1
  % this speeds up the multiplication, but would result in a sparse matrix when nchans=1
  gain = sparse(gain);
end

for i=1:length(sel)
  % nchan = numChannels(i);
  % nsmp  = numSamples(i);
  nchan = hdr.nChans;
  nsmp  = hdr.nSamples;
  % buf = read_ctf_shm(sel(i));
  % buf = buf(1:nchan*nsmp);
  if all([msgType(sel(i)) msgId(sel(i)) sampleNumber(sel(i)) numSamples(sel(i)) numChannels(sel(i))] ==  [ctf_shm.msgType(sel(i)) ctf_shm.msgId(sel(i)) ctf_shm.sampleNumber(sel(i)) ctf_shm.numSamples(sel(i)) ctf_shm.numChannels(sel(i))])
    % get the data from the cache
    buf = ctf_shm.data{sel(i)};
    ctf_shm.hit = ctf_shm.hit+1;
  else
    % read the data from shared memory, update the cache
    buf = read_ctf_shm(sel(i),nchan*nsmp);
    ctf_shm.data{sel(i)}         = buf;
    ctf_shm.msgType(sel(i))      = msgType(sel(i));
    ctf_shm.msgId(sel(i))        = msgId(sel(i));
    ctf_shm.sampleNumber(sel(i)) = sampleNumber(sel(i));
    ctf_shm.numSamples(sel(i))   = numSamples(sel(i));
    ctf_shm.numChannels(sel(i))  = numChannels(sel(i));
    ctf_shm.mis                  = ctf_shm.mis+1;
  end
  buf = reshape(buf, nchan, nsmp);
  thistrial = trlNum(sel(i));
  % apply calibration to the selected channels
  dat(:,:,thistrial-begtrial+1) = gain * double(buf(chanindx,:));
end

% if any(isnan(dat(:)))
%   ft_warning('data has been padded with NaNs');
%   fprintf('trials present   = %d - %d\n', min(trlNum), max(trlNum));
%   fprintf('trials requested = %d - %d\n', begtrial, endtrial);
% end

dimord = 'chans_samples_trials';


text:
function radian = toRadian(degree)
% degree to radian
radian = degree/180*pi;

author: xuedidi


text:
function [ori,varargout] = loadOrientation(fname,varargin)

warning('loadOrientation is depreciated. Please use instead orientation.load');
[ori,varargout{1:nargout-1}] = orientation.load(fname,varargin{:});


text:
function distance = compute_interfiber_distances_multires(fgFile, maxDistanceOfInterest,  range1start, range1end, range2start,  range2end)
%Multiresolution approach: first approximation with mass_center distances,
%for the closest fibers pointwise average distance is computed.

%First we compute center-of-mass interfiber distances for fg.fibers structures in infile.mat
%Then we threshold these distances to retain only those that are smaller
%than 20mm. For those interfiber parwise_dist distances are computed.
%save distance matrix into outfile. Range1 ([Nfrst Nlasst]) and range2 specify groups of fibers among which the distances should be computed.

%ER 12/2007 Range1/2start/end: supply as a string
%(e.g., '5') -- needed for later to pass bash string arguments, when
%working on the cluster.

%02/2008 Revision:
%1. Omit last 4 parameters if want to have the wholefile processed.
%2. Resampling is no longer performed -- use
%ResampledFiberGroup=resampleFiberGroup(fg, numNodes) and save the
%resampled dataset as save('ResampledDataFile', 'fg', 'versionNum',
%'coordinateSpace')
%3. Added parameter maxDistanceOfInterest: fibers whos mass-center distance
%in mm exceeds this parameters will be considered "infinitely remote" (for the
%purpose off sparse matrix production). A typical value is 40 (in mm) for the first step, 500 for the last step.

epsilon=.000000001;

%Strip off the .mat extension if provided
if fgFile(length(fgFile)-3:length(fgFile))=='.mat'
    fgFile=fgFile(1:length(fgFile)-4);
end
load(fgFile);

method='multiresolution';

%Checking number of arguments: if 1, full range in file is used; if 5, arbitrary ranges for first and second fiber groups are allowed.
if nargin==6
range1start=min(max(str2num(range1start), 1), size(fg.fibers, 1));
range2start=min(max(str2num(range2start), 1), size(fg.fibers, 1)) ;

range1end=max(min(str2num(range1end), size(fg.fibers, 1)), range1start);
range2end=max(min(str2num(range2end), size(fg.fibers, 1)), range2start) ;
samegroupflag=0;

else
    if nargin==2
    samegroupflag=1;
    range1start=1; range2start =1; %By default process full range.
    range1end=size(fg.fibers, 1);
    range2end=size(fg.fibers, 1);
    else
        display('Wrong number of argument supplied');
    end
end

display(['Data: ' fgFile ]);

display(['Distance metric: ' method ]);
display(['Fibergroups analyzed: ' num2str(range1start) ' to ' num2str(range1end) ' and ' num2str(range2start) ' to ' num2str(range2end)]);

fibergroup1=fg.fibers(range1start:range1end);
fibergroup2=fg.fibers(range2start:range2end);
clear fg;

if (range1start==range2start&range1end==range2end)
    display('One fiber group found');

    samegroupflag=1;
    tic; distanceSquared = InterfiberMassCenterDistances(fibergroup1);  toc;

else
    display('Two distinct fiber groups found');
    tic; distanceSquared =  InterfiberMassCenterDistances(fibergroup1, fibergroup2);  toc;

end
display('Mass center distances computed');

%Thresholding distmsr and transform it into the sparse matrix:
distanceSquared(distanceSquared(:)==0)=epsilon; %This is to keep actual zeros in the matrix, as the matrix entries corresponding to "larger-than-threshold" distances will be turned to zeros in the sparse matrix;
distanceSquared(distanceSquared>maxDistanceOfInterest^2)=0;  %Note that maxDistanceOfInterest has to be squared because distanceSquared returned by InterfiberMassCenterDistances is squared!
distanceSquared=sparse(double(distanceSquared));

%Form curve arrays from fg structure
nfibers1=size(fibergroup1, 1);
nfibers2=size(fibergroup2, 1);
npoints=  size(fibergroup1{1}, 2);

curves1=zeros(3, npoints, nfibers1);
%<s>Resample fibers in Fiberr Group 1 using splines</s>
%02/2008: Resampling no longer performed. Will complain if you have not
%resampled data apriori.

for i=1:nfibers1
%   curves1(:, :, i)=dtiFiberResample(fibergroup1{i}, npoints);
   curves1(:, :, i)=fibergroup1{i};
end
if (samegroupflag==1)
    curves2=curves1;
else
     curves2=zeros(3, npoints, nfibers1);
    %<s>Resample fibers in Fiberr Group 2 using splines</s>
    %02/2008: Resampling no longer performed.

    for i=1:nfibers2
       %curves2(:, :, i)=dtiFiberResample(fibergroup2{i}, npoints);
       curves2(:, :, i)=fibergroup2{i};
    end

end
%display('Fibers resampled');

if(samegroupflag==1)
    distanceSquared=tril(distanceSquared); %If the matrix is symmetric, avoid extra computations
end

%Distances for the following pairs of fibers will be fine-tuned:
[i1 j1]=find(distanceSquared~=0);
finetunedElementIndices=find(distanceSquared~=0);

%Computing interfiber distances (pointwise)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
distance=sparse(zeros(size(distanceSquared)));

tic;
for pair=1:size(i1)

distance(i1(pair), j1(pair))=intercurve_dist_pointwise(curves1(:, :, i1(pair)), curves2(:, :, j1(pair)), samegroupflag);
%The new, more "precise" values replace those not precise.
end
toc; %Wow, for 20% of fibers pairs still takes 5 times longer.
display(['Precise distances for pairs of fibers whose mass-center distance did not exceed ' num2str(maxDistanceOfInterest) 'mm computed']);

if(samegroupflag==1)
    distance=distance+distance'-triu(distance); %Since the matrix for "samegroup" case was symmetric, we need to repopulate the upper triangle.
end

% Remember: we have a sparse matrix representation where zeros mean
% "infinite distance" actually (and are not stored). So actual true zero distances should be  made a
% very small number.
distance(nonzeros(finetunedElementIndices.*(distance(finetunedElementIndices)==0)))=epsilon;

outfile=[fgFile 'distances_multires' num2str(range1start) 'to' num2str(range1end) 'vs' num2str(range2start) 'to' num2str(range2end) '.mat'];
[pathstr, name, ext, versn] = fileparts(outfile) ;
if isempty(pathstr)
    pathstr=['.' filesep ];
end

mkdir(pathstr, 'multiresDistances');
outfilefull=fullfile(pathstr, 'multiresDistances', name);
save(outfilefull, 'distance', 'fibergroup1', 'fibergroup2', 'fgFile', 'range1start', 'range2start', 'range1end', 'range2end', 'method');

author: vistalab


text:
function ROI = itkClass2roi(ni, varargin)
% Convert a layer from a nifti segmentation file (created by itkGray) into
% a mrVista Gray ROI
%
% ROI = itkClass2roi(ni, varargin)
%
%   INPUTS:
%       ni: a matrix, a nifti struct, or a path to a nifti file
%       Optional Inputs:
%           spath: a directory in which to save the ROI (if no path, then
%                 don't save; just return the ROI struct)
%           color: a 3-vector or color char (rbgcymkw)
%           comments: string for ROI.comments
%           name: string for ROI.name
%           layer: integer layer number in input matrix to be converted to
%                   ROI (Note that an ITK gray segmentation file can have
%                   any number of layers. We may want to convert only one
%                   of them to an ROI.)
%   OUTPUT:
%       ROI: a mrVista gray-view ROI struct
%
% see nifti2mrVistaAnat.m, itkClasses2roi.m
%
% Example:
%
%   ni =   '3DAnatomy/t1_class.nii';
%   fname = 'LeftWhiteMatter'
%   col = 'm';
%   ROI = itkClass2roi(ni, 'name', fname, 'color', col, 'layer', 2);
%
% April, 2009: JW

error('Obsolete. Use nifti2ROI instead.')

%--------------------------------------------------------------
% variable check
if nargin > 1
    for ii = 1:2:length(varargin)
        switch lower(varargin{ii})
            case {'spath'}
                spath = varargin{ii+1};
            case {'color', 'roicolor'}
                ROI.color = varargin{ii+1} / 255;
            case 'comments'
                ROI.comments = varargin{ii+1};
            case {'name', 'roiname'}
                ROI.name = varargin{ii+1};
            case {'layer'}
                layer = varargin{ii+1};
            case{'roi'}
                ROI = varargin{ii+1};
            otherwise
                warning('Unknown input arg: %s', varargin{ii}); %#ok<WNTAG>
        end
    end
end
%--------------------------------------------------------------

% We assume that the input data structure comes from an ITK-generated
% nifti, which must be permuted for our standard mrVista format
data = nifti2mrVistaAnat(ni);

% Find the indices in the 3d image that have nonzero values
if exist('layer', 'var')
    % if we know the layer, just get voxels with that value
    ROIinds = find(ismember(data,layer));
else
    % otherwise get all the non-zero voxels in the image
    ROIinds = find(data);
end

% Get the XYZ coordinates of these indices
[x y z] = ind2sub(size(data), ROIinds );

% Define a mrVista Gray ROI using these coordinates
ROI.coords   = single([x y z]');
ROI.created  = datestr(now);
ROI.modified = datestr(now);
ROI.viewType = 'Gray';

% Add optional or default ROI fields
if ~isfield(ROI, 'name'),       ROI.name = 'MyROI'; end
if ~isfield(ROI, 'color'),      ROI.color   = 'b'; end
if ~isfield(ROI, 'comments'),   ROI.comments = sprintf('%s: created by anat2roi', ROI.name); end

% If requested, save
if exist('spath', 'var'),
    save(fullfile(spath, ROI.name), 'ROI');
    fprintf('[%s] Saving ROI %s in %s\n', mfilename, ROI.name, spath);
end

% Done
return


text:
function ids = database_ids( database, gt_set )
    switch database,
      case 'nyud4oObj',
        ids = getImageSet(gt_set);
      otherwise,
        index_file = fullfile(mcg_root_dir,'datasets',database, 'gt_sets',[gt_set '.txt']);
        fileID = fopen(index_file);
        ids = textscan(fileID, '%s');
        ids = ids{1};
        fclose(fileID);
    end
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
%% PRIME_BATCH_ITHACA uses the BATCH command to run the PRIME code on the VT ITHACA cluster.
%
%  Discussion:
%
%    The PRIME code is a function, so first we must write a script
%    called PRIME_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_BATCH_ITHACA:\n' );
  fprintf ( 1, '  Run PRIME_SCRIPT on Ithaca.\n' );
%
%  The BATCH command sends the script for execution.
%
  my_job = batch ( 'prime_script', ...
    'Configuration', 'ithaca_2009b', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'prime_fun' }, ...
    'CurrentDirectory', '/home/burkardt/matlab', ...
    'matlabpool', 4 );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  LOAD makes the script's workspace available.
%
%  total = total number of primes.
%
  load ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Total number of primes = %d\n', total );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_BATCH_ITHACA:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

repo: jburkardt-m


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
classdef Matrix < ocl.types.Structure
  %OCLMATRIX Matrix valued structure for variables
  %
  properties
    msize
  end

  methods

    function self = Matrix(size)
      % ocl.types.Matrix(size)
      self.msize = size;
      self.len = prod(size);
    end
    function [N,M] = size(self)
      s = self.msize;
      if nargout>1
        N = s(1);
        M = s(2);
      else
        N = s;
      end
    end
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function bboxeso = bbox_remove_duplicates(bboxes, minSize, maxNum)
% Copyright (C) 2016 Hakan Bilen.
% 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).
bboxeso = cell(size(bboxes));
for i=1:numel(bboxes)
  bbox = bboxes{i};
  % remove small bbox
  isGood = (bbox(:,3)>=bbox(:,1)-1+minSize) & (bbox(:,4)>=bbox(:,2)-1+minSize);
  bbox = bbox(isGood,:);
  % remove duplicate ones
  [dummy, uniqueIdx] = unique(bbox, 'rows', 'first');
  uniqueIdx = sort(uniqueIdx);
  bbox = bbox(uniqueIdx,:);
  % limit number for training
  nB = min(size(bbox,1),maxNum);

  bboxeso{i} = bbox(1:nB,:);
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function net = gpunpak(net, hp)
%GPUNPAK Separates hyperparameter vector into components.
%
%	Description
%	NET = GPUNPAK(NET, HP) takes an Gaussian Process data structure NET
%	and  a hyperparameter vector HP, and returns a Gaussian Process data
%	structure  identical to the input model, except that the covariance
%	bias BIAS, output noise NOISE, the input weight vector INWEIGHTS and
%	the vector of covariance function specific parameters  FPAR have all
%	been set to the corresponding elements of HP.
%
%	See also
%	GP, GPPAK, GPFWD, GPERR, GPGRAD
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check arguments for consistency
errstring = consist(net, 'gp');
if ~isempty(errstring);
  error(errstring);
end
if net.nwts ~= length(hp)
  error('Invalid weight vector length');
end

net.bias = hp(1);
net.noise = hp(2);

% Unpack input weights
mark1 = 2 + net.nin;
net.inweights = hp(3:mark1);

% Unpack function specific parameters
net.fpar = hp(mark1 + 1:size(hp, 2));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%DEMO_MEMORYSAVE  Demonstration of memory save option in GPstuff
%
%  Description
%    This demo consists of various combinations of covariance functions,
%    likelihoods and full/sparse approximations. This demo is intended for
%    showing how to use memory save option in GPstuff. Unfortunately MATLAB
%    doesn't have memory usage monitoring, so the purpose of the demo is only
%    to test that results are the same whether you use memory save or not
%    and to show that the running times with memory saving are little bit longer
%    (more overhead). Memory saving can be enabled with the following command
%
%       gp = gp_set(..., 'savememory', 'on');
%
%  See also
%    GP_SET
%

% Copyright (c) 2012 Ville Tolvanen

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

% Classification data with full GP's

S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];

lik = lik_logit();

pl = prior_t();
pm = prior_sqrtunif();

gpcf = gpcf_neuralnetwork('weightSigma2', [0.9 0.9], 'biasSigma2', 10);
gpcf = gpcf_neuralnetwork(gpcf, 'weightSigma2_prior', pl,'biasSigma2_prior', pm); %

gpcf2 = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 10);
gpcf2 = gpcf_sexp(gpcf2, 'lengthScale_prior', pl,'magnSigma2_prior', pm); %

% GP without memory saving
gp = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-4, 'latent_method', 'EP');
% GP with memory saving option enabled
gp2 = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-4, 'latent_method', 'EP', 'savememory', 'on');

fprintf('Classification (EP), Neuralnetwork covariance function without and with memory saving (optimization and prediction)\n');
opt=optimset('TolX',1e-3,'TolFun',1e-3, 'Display', 'off');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'opt',opt);
[Eft,Varft, lpyt]=gp_pred(gp,x,y,x, 'yt', y);toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'opt',opt);
[Eft2,Varft2, lpyt2]=gp_pred(gp2,x,y,x, 'yt', y);toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft==Eft2)); assert(all(Varft==Varft2)); assert(all(lpyt==lpyt2));

gp = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-6, 'latent_method', 'Laplace');
gp2 = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-6, 'latent_method', 'Laplace', 'savememory', 'on');

fprintf('Classification (Laplace), Squared-Exponential covariance function without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'opt',opt);
[Eft,Varft, lpyt]=gp_pred(gp,x,y,x, 'yt', y);toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'opt',opt);
[Eft2,Varft2, lpyt2]=gp_pred(gp2,x,y,x, 'yt', y);toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft==Eft2)); assert(all(Varft==Varft2)); assert(all(lpyt==lpyt2));

% Regression data with sparse approximations

prevstream=setrandstream();
xx=linspace(1,10,901);
x1=logspace(0,1,100);
x1=round(x1*100)-99;
x=xx(x1)';
y=2*sin(4*x)+0.2*randn(size(x));
xt=[1:0.01:14]';
[n,nin] = size(x);

pl = prior_t('s2', 1);
pm = prior_logunif();
pn = prior_logunif();
gpcfse = gpcf_matern32('lengthScale',0.5,'magnSigma2',1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
lik = lik_gaussian('sigma2', 0.1, 'sigma2_prior', pn);
gp = gp_set('lik', lik, 'cf', gpcfse, 'jitterSigma2', 1e-6);
opt=optimset('TolFun',1e-4,'TolX',1e-4);

fprintf('Regression, FIC GP, Matern-3/2 covariance function, without and with memory saving (optimization and prediction)\n')
Xu=round(10+90*rand(18,1))/10; % Random placement

gp_fic = gp_set(gp, 'type','FIC','X_u',Xu,'infer_params','covariance+likelihood+inducing');
gp_fic2 = gp_set(gp, 'type','FIC','X_u',Xu,'infer_params','covariance+likelihood+inducing', 'savememory', 'on');

opt=optimset('TolFun',1e-4,'TolX',1e-4, 'Display', 'off');

% Optimization and prediction without memory saving
tic,gp_fic=gp_optim(gp_fic,x,y,'opt',opt);
[Eft_fic, Varft_fic] = gp_pred(gp_fic, x, y, xt);
Varft_fic = Varft_fic + gp_fic.lik.sigma2;toc

% Optimization and prediction with memory saving
tic,gp_fic2=gp_optim(gp_fic2,x,y,'opt',opt);
[Eft_fic2, Varft_fic2] = gp_pred(gp_fic2, x, y, xt);
Varft_fic2 = Varft_fic2 + gp_fic.lik.sigma2;toc
setrandstream([],prevstream);

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft_fic==Eft_fic2)); assert(all(Varft_fic==Varft_fic2));

gpcfse = gpcf_sexp('lengthScale',0.5,'magnSigma2',1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gp = gp_set('lik', lik, 'cf', gpcfse, 'jitterSigma2', 1e-6);

gp_fic = gp_set(gp, 'type','VAR','X_u',Xu,'infer_params','covariance+likelihood+inducing');
gp_fic2 = gp_set(gp, 'type','VAR','X_u',Xu,'infer_params','covariance+likelihood+inducing', 'savememory', 'on');

fprintf('Regression, VAR GP, Squared-Exponential covariance function, without and with memory saving (optimization and prediction)\n')

% Optimization and prediction without memory saving
tic,gp_fic=gp_optim(gp_fic,x,y,'opt',opt);
[Eft_fic, Varft_fic] = gp_pred(gp_fic, x, y, xt);
Varft_fic = Varft_fic + gp_fic.lik.sigma2;toc

% Optimization and prediction with memory saving
tic,gp_fic2=gp_optim(gp_fic2,x,y,'opt',opt);
[Eft_fic2, Varft_fic2] = gp_pred(gp_fic2, x, y, xt);
Varft_fic2 = Varft_fic2 + gp_fic.lik.sigma2;toc
setrandstream([],prevstream);

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Eft_fic==Eft_fic2)); assert(all(Varft_fic==Varft_fic2));

% Spatial data with sparse approximations

S = which('demo_spatial1');
data = load(strrep(S,'demo_spatial1.m','demodata/spatial1.txt'));
x = data(1:200,1:2);
ye = data(1:200,3);
y = data(1:200,4);
dims = [1    60     1    35];
[trindex, Xu] = set_PIC(x, dims, 5, 'corners', 0);
[n,nin] = size(x);

pl = prior_t('s2',10);
pm = prior_sqrtunif();

gpcf1 = gpcf_ppcs3('nin',nin,'lengthScale', 1, 'magnSigma2', 0.1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcf2 = gpcf_sexp('lengthScale', 5, 'magnSigma2', 0.05, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
gpcf3 = gpcf_neuralnetwork('weightSigma2', [1 1], 'biasSigma2', 0.05, 'weightSigma2_prior', pl, 'biasSigma2_prior', pm);

lik = lik_negbin();
% GP without memory save
gp = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+inducing');
% GP with memory saving option enabled
gp2 = gp_set('type', 'FIC', 'lik', lik, 'cf', gpcf1, 'X_u', Xu, ...
            'jitterSigma2', 1e-4, 'infer_params', 'covariance+inducing', 'savememory', 'on');
opt=optimset('TolFun',1e-2,'TolX',1e-2, 'Display', 'off');

fprintf('Spatial process (Laplace), FIC GP, PPCS3 covariance function, without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

gp = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-4, 'latent_method', 'EP');
gp2 = gp_set('lik', lik, 'cf', gpcf2, 'jitterSigma2', 1e-4, 'latent_method', 'EP', 'savememory', 'on');

fprintf('Spatial process (EP), FIC GP, Squared-Exponential covariance function, without and with memory saving (optimization and prediction)\n');

% Optimization and prediction without memory saving
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Optimization and prediction with memory saving
tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

gp = gp_set('lik', lik, 'cf', gpcf3, 'jitterSigma2', 1e-4);
gp2 = gp_set('lik', lik, 'cf', gpcf3, 'jitterSigma2', 1e-4, 'savememory', 'on');

fprintf('Spatial process (Laplace), FIC GP, Neuralnetwork covariance function, without and with memory saving (optimization and prediction)\n');
tic,gp=gp_optim(gp,x,y,'z',ye,'opt',opt);
[Ef, Varf] = gp_pred(gp, x, y, x, 'z', ye, 'tstind', [1:n]); toc

tic,gp2=gp_optim(gp2,x,y,'z',ye,'opt',opt);
[Ef2, Varf2] = gp_pred(gp2, x, y, x, 'z', ye, 'tstind', [1:n]); toc

% Check that predictions (and optimization) with and without memory saving are same
assert(all(Ef==Ef2)); assert(all(Varf==Varf2));

author: fieldtrip


text:
function gf=gfobs_L1L2(obsr,obsb,lam)

pi=sdobs(obsr,obsb,1)*lam(1);
pj=sdobs(obsr,obsb,2)*lam(2);

if pi==0||pj==0
    gf=0;
else
    gf=pi-pj;
end

return

repo: GINav


text:
function [p,priors]=create_parameters(start_from_mode)
p=struct();
%--- Targeted steady state values
p.alpha = 0.30; %(1) the average labor income share is 70% (page 15 of paper)
p.r = 1.01; %(2) the average real prime loan rate is 4% per annum (page 15 of paper)
p.ky = 4.6194; %(3) the capital-output ratio is on average 1.15 at the annual frequency (page 15 of paper)
p.ik = 0.2093/4; % (4) the investment-capital ratio is on average 0.209 at the annual frequency
p.qlLeOY = 2.60; %(5) the average land-output ratio is 0.65 at the annual frequency
p.theta_bar = 0.75;  %(6) the average nonfarm and nonfiancial businesses' loan-asset ratio is 0.75 at the annual frequency (page 15 of paper)
p.qlLhOY = 5.8011; %(7) the average housing output ratio is 1.45 at the annual frequency
p.n = 1/4; %(8) average market hours is 25% of time endowment
p.lbar = 1; %arbitrary

priors=struct();
priors.gamma_h          ={ 0.7  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.gamma_e          ={ 0.7  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.omega            ={ 3.0  , 0.1020, 5.994  , 'gamma_pdf(0.9)'    , .00001        , 10          };
priors.g_trans          ={ 0.375, 0.1   , 1.5    , 'gamma_pdf(0.9)'    , .00001        , 10          };
priors.lambda_qbar_trans={ 1.25 , 0.1   , 1.5    , 'gamma_pdf(0.9)'    , .00001        , 10          };
priors.rho_a            ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_z            ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_v            ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_q            ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_mu           ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_phi          ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_psi          ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.rho_xitheta      ={ 0.9  , 0.0256, 0.7761 , 'beta_pdf(0.9)'     , .00000001     , .99999999999};
priors.sig_Eps_a        ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_z        ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_v        ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_q        ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_mu       ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_phi      ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_psi      ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };
priors.sig_Eps_xitheta  ={ 0.01 , 0.0001, 2      , 'inv_gamma_pdf(0.9)', 0.000000000001, 100         };

if start_from_mode
    priors.gamma_e{1}= 0.658435028;
    priors.gamma_h{1}= 0.497555761;
    priors.omega{1}= 0.175347074;
    priors.g_trans{1}= 0.422109209;
    priors.lambda_qbar_trans{1}= 1.21261837;
    priors.rho_a{1}=  0.905471687;
    priors.rho_v{1}= 0.00947706131;	%rho_{nu_z}
    priors.rho_mu{1}= 0.294889827;	   %rho_{nu_q}
    priors.rho_psi{1}=  0.982918278;
    priors.rho_q{1}= 0.561965864;
    priors.rho_xitheta{1}= 0.980427788;
    priors.rho_phi{1}=  0.999758433;
    priors.rho_z{1}= 0.426298781;	 %rho_z
    priors.sig_Eps_z{1}=     0.00418908791;
    priors.sig_Eps_v{1}=   0.00365910563;
    priors.sig_Eps_q{1}=     0.00415238812;
    priors.sig_Eps_mu{1}=   0.00288990925;
    priors.sig_Eps_a{1}=     0.101281626;
    priors.sig_Eps_phi{1}=   0.0462334747;
    priors.sig_Eps_psi{1}=    0.00731583619;
    priors.sig_Eps_xitheta{1}= 0.0111777605;
end

% add the initial conditions of the priors to the parameters
%------------------------------------------------------------
fields=fieldnames(priors);
for ip=1:numel(fields)
    name=fields{ip};
    p.(name)=priors.(name){1};
end

end