text:
function [hdr] = read_ctf_res4(fname)

% READ_CTF_RES4 reads the header in RES4 format from a CTF dataset
%
% Use as
%   [hdr] = read_ctf_res4(filename)
%
% See also READ_CTF_MEG4

% "VSM MedTech Ltd. authorizes the release into public domain under the
% GPL licence of the Matlab source code files "read_ctf_res4.m" and
% "read_ctf_meg4.m" by the authors of said files from the F.C. Donders
% Centre, Nijmegen, The Netherlands."

% Author(s): Jim McKay November 1999
% Last revision: Jim McKay
% Copyright (c) 1999-2000 CTF Systems Inc. All Rights Reserved.
%
% modifications Copyright (C) 2002, Ole Jensen
% modifications Copyright (C) 2003, 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$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read header information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

try
  fid = fopen_or_error(fname,'r','ieee-be');
catch err
  ft_error('Cound not open header file: %s', err.message);
end

% First 8 bytes contain filetype, check is fileformat is correct.
% This function was written for MEG41RS, but also seems to work for some other formats
CTFformat=char(fread(fid,8,'uint8'))';
if ~strcmp(CTFformat(1,1:7),'MEG41RS') && ~strcmp(CTFformat(1,1:7),'MEG42RS') && ~strcmp(CTFformat(1,1:7),'MEG3RES')
  ft_warning('res4 format (%s) is not supported for file %s, trying anyway...', CTFformat(1,1:7), fname);
end

% Read the initial parameters
appName       = char(fread(fid,256,'uint8'))' ;
dataOrigin    = char(fread(fid,256,'uint8'))' ;
dataDescrip   = char(fread(fid,256,'uint8'))' ;
no_trial_avgd = fread(fid,1,'int16')          ;
data_time     = char(fread(fid,255,'uint8'))';
data_date     = char(fread(fid,255,'uint8'))';

fseek(fid,1288,'bof');
% Read the general recording parameters
no_samples  = fread(fid,1,'int32');
no_channels = fread(fid,1,'int16');
fseek(fid,2,'cof');         % hole of 2 bytes due to improper alignment
sample_rate = fread(fid,1,'double');
epoch       = fread(fid,1,'double');
no_trials   = fread(fid,1,'int16');
fseek(fid,2,'cof');         % hole of 2 bytes due to improper alignment
preTrigpts=fread(fid,1,'int32');

fseek(fid,1360,'bof');
% read in the meg4Filesetup structure
run_name     = char(fread(fid,32,'uint8')');
run_title    = char(fread(fid,256,'uint8')');
instruments  = char(fread(fid,32,'uint8')');
coll_desc    = char(fread(fid,32,'uint8')');
subj_id      = char(fread(fid,32,'uint8')');
operator     = char(fread(fid,32,'uint8')');
sensFilename = char(fread(fid,60,'uint8')');

% not nececssary to seek, the file pointer is already at the desired location
% fseek(fid,1836,'bof');

% Read in the run description length
rd_len=fread(fid,1,'uint32');
% Go to the run description and read it in
fseek(fid,1844,'bof');
run_desc=char(fread(fid,rd_len,'uint8')');

% read in the filter information
num_filt=fread(fid,1,'uint16');
for fi=0:(num_filt-1),
  %filt_info=fread(fid,18,'uint8');
  filt_freq =fread(fid,1, 'double');
  filt_class=fread(fid,1, 'uint32');
  filt_type =fread(fid,1, 'uint32');
  num_fparm=fread(fid, 1, 'uint16');
  %num_fparm=filt_info(18);
  if num_fparm ~= 0,
    filt_parm=fread(fid,8*num_fparm,'uint8');
  end % if
end % for fi

% Read in the channel names
for i=1:no_channels,
  temp=fread(fid,32,'uint8')';
  temp(find(temp<32 )) = ' ';       % remove non-printable characters
  temp(find(temp>126)) = ' ';       % remove non-printable characters
  endstr = strfind(temp, '-'); temp(endstr:end) = ' ';  % cut off at '-'
  endstr = strfind(temp, ' '); temp(endstr:end) = ' ';  % cut off at ' '
  chan_name(i,:) = char(temp);      % as char array
  chan_label{i}  = deblank(char(temp)); % as cell-array
end %for

% pre-allocate some memory space
sensGain = zeros([no_channels,1]);
qGain    = zeros([no_channels,1]);
ioGain   = zeros([no_channels,1]);
sensType = zeros([no_channels,1]);

% Read in the sensor information
fp = ftell(fid);
for chan=1:no_channels,
  fread(fid,1,'uint8');         % Read and ignore 1 byte from enum
  sensType(chan)=fread(fid,1,'uint8');  % Read sensor type
  fread(fid,2,'uint8');         % Read and ignore originalRunNum
  fread(fid,4,'uint8');         % Read and ignore coilShape
  sensGain(chan)=fread(fid,1,'double'); % Read sensor gain in Phi0/Tesla
  qGain(chan)=fread(fid,1,'double');        % Read qxx gain (usually 2^20 for Q20)
  ioGain(chan)=fread(fid,1,'double');       % Read i/o gain of special sensors (usually 1.0)
  ioOffset(chan)=fread(fid,1,'double');
  numCoils(chan)=fread(fid,1,'int16');
  grad_order_no(chan)=fread(fid,1,'int16');
  fread(fid,4,'uint8');

  % read the coil positions and orientations
  for i=1:8
    Chan(chan).coil(i).pos = fread(fid,3,'double')';
    fread(fid,1,'double');
    Chan(chan).coil(i).ori = fread(fid,3,'double')';
    fread(fid,3,'double');
  end

  % read the coil positions and orientations in head coordinates(?)
  for i=1:8
    Chan(chan).coilHC(i).pos = fread(fid,3,'double')';
    fread(fid,1,'double');
    Chan(chan).coilHC(i).ori = fread(fid,3,'double')';
    fread(fid,3,'double');
  end

  % jump to the next sensor info record
  fseek(fid, fp+chan*1328, 'bof');
end % for chan

% close the header file
fclose(fid);

% according to Tom Holroyd, the sensor types are
%
% meg channels are 5, refmag 0, refgrad 1, adcs 18.
% UPPT001 is 11
% UTRG001 is 11
% SCLK01 is 17
% STIM is 11
% SCLK01 is 17
% EEG057 is 9
% ADC06 is 18
% ADC07 is 18
% ADC16 is 18
% V0 is 15

% assign all the variables that should be outputted as header information
hdr.Fs          = sample_rate;
hdr.nChans      = no_channels;
hdr.nSamples    = no_samples;
hdr.nSamplesPre = preTrigpts;
hdr.timeVec     = (1:no_samples)/sample_rate - preTrigpts/sample_rate - 1/sample_rate;
hdr.nTrials     = no_trials;
hdr.gainV       = ioGain./(qGain.*sensGain);
hdr.ioGain      = ioGain;
hdr.qGain       = qGain;
hdr.sensGain    = sensGain;
hdr.sensType    = sensType;

hdr.label       = chan_label(:);
hdr.nameALL     = chan_name;
hdr.Chan        = Chan;
% hdr.rowMEG      = rowMEG;
% hdr.rowEEG      = rowEEG;
% hdr.rowTRIG     = rowTRIG;
% hdr.rowREF      = rowREF;
% hdr.nameEEG     = [];
% hdr.nameMEG     = [];
% hdr.nameEOG     = [];
% hdr.trigV       = [];
% hdr.SwapData    = [];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function rmSetDataTYPES(vw, params)
% Write the retinotopic model stimulus parameters to the dataTYPES
%
%   rmSetDataTYPES(vw, params)

mrGlobals;

dt = viewGet(vw, 'dataTypeNumber');

% %This returns an error if the fields already in dt.retinotopyModelParams
% are not exactly the same as those in params.stim
%
% for scan = 1:numel(params.stim)
%     dataTYPES(dt) = dtSet(dataTYPES(dt),...
%         'retinotopyModelParams', params.stim(scan), scan);
% end

dataTYPES(dt) = dtSet(dataTYPES(dt),'retinotopyModelParams', params.stim);

saveSession

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [mMedModF, mStdL, loopout] = brutebootloglike_a2(time_as, time_asf, bootloops,fT1, nMod)
    % BRUTEBOOTLOGLIKE_A2 Bootstrap analysis of Omori parameters calculated by brute force
    % (p1,p2,c1,c2,k1,k2)-pair is mean of the bootstrap values by determining the mean cumulative
    % number modeled a end of the learning period
    % Standard deviations are calculated as the 2nd moment, not to rely fully on normal distributions
    %
    % [mMedModF, mStdL, loopout] = BRUTEBOOTLOGLIKE_A2(time_as, time_asf, bootloops,fT1, nMod);
    % -------------------------------------------------------------------------------
    %
    % Input parameters:
    %   time_as     Delay times [days] of learning period
    %   time_asf    Delay times [days] until end of forecast period
    %   bootloops   Number of bootstraps
    %   fT1         Time of biggest aftershock in learning period
    %   nMod        Model to fit data, three models including a secondary aftershock sequence.
    %               Different models have varying amount of free parameters
    %               before (p1,c1,k1) and after (p2,c2,k2) the aftershock occurence
    %               1: modified Omori law (MOL): 3 free parameters
    %                  p1=p2,c1=c2,k1=k2
    %               2: MOL with one secondary aftershock sequence:4 free parameters
    %                  p1=p2,c1=c2,k1~=k2
    %               3: MOL with one secondary aftershock sequence:5 free parameters
    %                  p1~=p2,c1=c2,k1~=k2
    %               4: MOL with one secondary aftershock sequence:6 free parameters
    %                  p1~=p2,c1~=c2,k1~=k2
    %
    % Output parameters:
    %  mMedModF :  Result matrix including the values for the mean forecast at end of forecast period
    %  mStdL    :  Uncertainties of fit to the data in learning period
    %  loopout     contains all results
    %
    % Samuel Neukomm / S. Wiemer / J. Woessner
    % updated: 05.08.03

    time_as = sort(time_as);
    %bootloops = 50; % number of bootstrap samples
    n = length(time_as);
    loopout = nan(bootloops,9); % 8 from bruteforceloglike_a2, plus variate column.
    % Initialize random seed
    rng('shuffle');

    % i = (1:n)';
    for j = 1:bootloops
        newtas = sort(datasample(time_as, n, 'Replace',true));
        [pv1, pv2, cv1, cv2, kv1, kv2, fAIC, fL] = bruteforceloglike_a2(newtas, fT1, nMod);
        loopout(j,1:8) = [pv1, pv2, cv1, cv2, kv1, kv2, fAIC, fL];
    end

    % New version: Choose mean (p,c,k)-variables by modelling the cumulative number at end of
    % the learning period

    % 2nd moment i.e. Standard deviations
    [pstd1] = std(loopout(:,1),1,'omitnan');
    [pstd2] = std(loopout(:,2),1,'omitnan');
    [cstd1] = std(loopout(:,3),1,'omitnan');
    [cstd2] = std(loopout(:,4),1,'omitnan');
    [kstd1] = std(loopout(:,5),1,'omitnan');
    [kstd2] = std(loopout(:,6),1,'omitnan');

    % Uncertainties of fit
    mStdL = [pstd1 pstd2 cstd1 cstd2 kstd1 kstd2];

    %% Compute best fitting pair of variates
    % TODO (maybe) vectorize this
    n_time_asf = length(time_asf);

    pv1=loopout(:,1);
    pv2=loopout(:,2);
    cv1=loopout(:,3);
    cv2=loopout(:,4);
    kv1=loopout(:,5);
    kv2=loopout(:,6);

    cumnr_model = OmoriModel.doForecast(nMod,time_asf, pv1, cv1, kv1, fT1, kv2, pv2, cv2);
    loopout(:,9) = max(cumnr_model);

    %%
    [Y, in] = sort(loopout(:,9));
    loops = loopout(in,:);
    % % Median values: Old version
    % vMedian = abs(loops(:,9)-median(loops(:,9)));
    % nMedian = (find(vMedian == min(vMedian)));
    %
    % if length(nMedian(:,1)) > 1
    %     nMedian = nMedian(1,1);
    % end
    % pmedian1 = loops(nMedian,1);
    % pmedian2 = loops(nMedian,2);
    % cmedian1 = loops(nMedian,3);
    % cmedian2 = loops(nMedian,4);
    % kmedian1 = loops(nMedian,5);
    % kmedian2 = loops(nMedian,6);
    %
    % mMedModF = [pmedian1, pstd1, pmedian2, pstd2, cmedian1, cstd1, cmedian2, cstd2, kmedian1, kstd1, kmedian2, kstd2];

    % Mean values
    vMean = abs(loops(:,9)-mean(loops(:,9)));
    nMean = (find(vMean == min(vMean)));

    if length(nMean(:,1)) > 1
        nMean = nMean(1,1);
    end
    pMean1 = loops(nMean,1);
    pMean2 = loops(nMean,2);
    cMean1 = loops(nMean,3);
    cMean2 = loops(nMean,4);
    kMean1 = loops(nMean,5);
    kMean2 = loops(nMean,6);

    mMedModF = [pMean1, pstd1, pMean2, pstd2, cMean1, cstd1, cMean2, cstd2, kMean1, kstd1, kMean2, kstd2];

author: CelsoReyes


text:
function sev = read_tdt_sev(filename, dtype, begsample, endsample)

% READ_TDT_SEV
%
% Use as
%   sev = read_tdt_sev(filename, dtype, begsample, endsample)
%
% Note: sev files contain raw broadband data that is streamed to the RS4

fid = fopen(filename, 'rb', 'ieee-le');

if nargin<3
  begsample = 1;
end

if nargin<4
  endsample = inf;
end

switch dtype
  case 0
    fmt   = 'float32';
    wsize = 4;
  case 1
    fmt   = 'int32';
    wsize = 4;
  case 2
    fmt   = 'int16';
    wsize = 2;
  case 3
    fmt   = 'int8';
    wsize = 1;
  case 4
    fmt   = 'double';
    wsize = 8;
  case 5
    error('don''t know what a DFORM_QWORD is');
  otherwise
    error('unknown dtype');
end

fseek(fid, begsample*wsize, 'cof');
sev = fread(fid, endsample-begsample, fmt);
fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function out = extractFIR(x,n)
out = reshape(VBA_vec(x)',n,[]);
out = mean(out,2)';

author: MBB-team

repo: VBA-toolbox


text:
function handle = imageModifyVargplvm(handle, imageValues, imageSize, transpose, negative, ...
		     scale,thresh)

% IMAGEMODIFY Helper code for visualisation of image data.
% FORMAT
% DESC is a helper function for visualising image data using latent
% variable models.
% ARG handle : the handle of the image data.
% ARG imageValues : the values to set the image data to.
% ARG imageSize : the size of the image.
% ARG transpose : whether the resized image needs to be transposed
% (default 1, which is yes).
% ARG negative : whether to display the negative of the image
% (default 0, which is no).
% ARG scale : dummy input, to maintain compatability with
% IMAGEVISUALISE.
% ARG thresh : An array, thresh = [minVal, threshDown, threshUp, maxVal]
% which gives value minVal to all image elements smaller than threshDown
% and value maxVal to all image elements larger than threshUp.
% RETURN handle : a the handle to the image data.
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2006
%
% SEEALSO : imageVisualise, fgplvmResultsDynamic

% SHEFFIELDML

if nargin < 4 || isempty(transpose)
  transpose = 1;
end
if nargin< 5 || isempty(negative)
  negative = 0;
end
if nargin < 6
    thresh = [];
end
if negative
  imageValues = -imageValues;
end
if ~isempty(thresh)
    imageValues(imageValues < thresh(2)) = thresh(1);
    imageValues(imageValues >= thresh(3)) = thresh(4);
end
if transpose
  set(handle, 'CData', reshape(imageValues(1:imageSize(1)*imageSize(2)), imageSize(1), imageSize(2))');
else
  set(handle, 'CData', reshape(imageValues(1:imageSize(1)*imageSize(2)), imageSize(1), imageSize(2)));
end

author: SheffieldML

repo: vargplvm


text:
function im = recomputeImage2(view,clipMode)
%
% im = recomputeImage2(view,clipMode)
%
% Recomputes the image (underlay/anat + overlay) for the
% given view, returning the image. This is different from
% recomputeImage in 3 ways:
%       1) It produces a truecolor image;
%       2) It returns the image rather than the view
%          (I do this because I'm adding the option to
%           work on mosaics of many of these images -- e.g.
%           many inplane slices or flat levels)
%       3) It doesn't take numGrays or numColors as arguments;
%          (since it's true color, they're set to 256 each).
%
%
% djh, sometime in '98
% djh, 2/2001. version 3.0
% ras, 3/2004. An attempt to make the images true-color (3 image planes
% (R,G,B), each ranging from 0-255, to allow for a full 256-value dynamic
% range for both the anatomies and overlays). Will check that this doesn't
% critically slow down updating the screen (with modern processors, this
% seems unlikely).

% Initialize images
anatIm=[];
overlay=[];

% Get cothresh, phWindow, and mapWindow from sliders
cothresh = getCothresh(view);
phWindow = getPhWindow(view);
mapWindow = getMapWindow(view);

% although we accept numGrays and numColors for back-compatibility, now
% that it's truecolor, we'll set these directly:
% numGrays = 64;
% numColors = 64;

numGrays = 256;
numColors = 256;

% Get anatClip from sliders
anatClip = getAnatClip(view);

% Get anatomy image
anatIm = cropCurAnatSlice(view);

% Get overlay
overlay = [];
if ~strcmp(view.ui.displayMode,'anat')
  overlay = cropCurSlice(view,view.ui.displayMode);
end

% Select pixels that satisfy cothresh, phWindow, and mapWindow
pts = [];
if ~isempty(overlay)
  pts = ones(size(overlay));
  curCo=cropCurSlice(view,'co');
  curPh=cropCurSlice(view,'ph');
  curMap=cropCurSlice(view,'map');
  if ~isempty(curCo) & cothresh>0
    ptsCo = curCo > cothresh;
    pts = pts & ptsCo;
  end
  if ~isempty(curPh)
    if diff(phWindow) > 0
      ptsPh = (curPh>=phWindow(1) & curPh<=phWindow(2));
    else
      ptsPh = (curPh>=phWindow(1) | curPh<=phWindow(2));
    end
    pts = pts & ptsPh;
  end
  if strcmp(view.ui.displayMode, 'amp')
    curAmp = cropCurSlice(view, 'amp');
    mnv = min(curAmp(:));
    mxv = max(curAmp(:));
    curMap = (curAmp - mnv) ./ (mxv - mnv);
  end
  if ~isempty(curMap)
    ptsMap = (curMap>=mapWindow(1) & curMap<=mapWindow(2));
    pts = pts & ptsMap;
  end
end

% Rescale anatIm to [1:numGrays], anatClip determines the range
% of anatomy values that gets mapped to the available grayscales.
% If anatClip=[0,1] then there is no clipping and the entire
% range of anatomy values is scaled to the range of available gray
% scales.
minVal = double(min(anatIm(:)));
maxVal = double(max(anatIm(:)));
anatClipMin = min(anatClip)*(maxVal-minVal) + minVal;
anatClipMax = max(anatClip)*(maxVal-minVal) + minVal;
warning off;
anatIm = (rescale2(double(anatIm),[anatClipMin,anatClipMax],[1,numGrays]));
%keyboard

warning backtrace;

% Rescale overlay to [0 numGrays-1]
if ~isempty(overlay)
   if strcmp(clipMode,'auto')
      if ~isempty(find(pts));
         overClipMin = min(overlay(pts));
         overClipMax = max(overlay(pts));
      else
         overClipMin = min(overlay(:));
         overClipMax = max(overlay(:));
      end
   else
      overClipMin = min(clipMode);
      overClipMax = max(clipMode);
   end
   overlay=rescale2(overlay,[overClipMin overClipMax],[0 numGrays-1]);
end

% get, and threshold, the anatomical cmap:
lightRng = [0.6 0.8];
darkRng = [0.2 0.4];
thresh = 0.6;
anatCmap = gray(256);
anatCmap(anatCmap < thresh) = ...
    normalize(anatCmap(anatCmap < thresh),darkRng(1),darkRng(2));
anatCmap(anatCmap >= thresh) = ...
    normalize(anatCmap(anatCmap >= thresh),lightRng(1),lightRng(2));

% % convert into a truecolor RGB image, combining
% % the overlay with the anatomy, if necessary
% % Combine overlay with anatomy image
% im = normalize(repmat(anatIm,[1 1 3]));
%
 if ~isempty(overlay) & ~all(pts==0)
     % for truecolor, we need to get the color map
     % info directly.
     cmapname = eval(['view.ui.' view.ui.displayMode 'Mode.name']);
     cmapname = cmapname(1:end-4); % last 4 chars are 'Cmap'
     cmap = feval(cmapname,0,256);
     cmap = [anatCmap; cmap];
     im = ind2rgb(overlay,cmap);

%     % get a set of R,G,B columns of values for each data point in the
%     % overlay
%     vals = cmap(overlay(pts)+1,:);
%
%     % plug in the overlay into each R/G/B plane of the truecolor image:
%     [indx,indy] = ind2sub(size(anatIm),find(pts));
%     indz = [ones(length(indx),1); 2*ones(length(indx),1); 3*ones(length(indx),1)];
%     indx = repmat(indx,3,1);
%     indy = repmat(indy,3,1);
%     final_ind = sub2ind(size(im),indx,indy,indz);
%    im(final_ind) = vals
else
    im = ind2rgb(anatIm,anatCmap);
end

% 2003.01.10 RFD: the following is no longer necessary- the uint8 data
% can't have any NaNs! Also, it caused problems in matlab versions <6.5.
% 2003.01.23 ARW: But without it, Matlab 6.5 fills in the Nans in the flat map as white.
% Do a version check for now and replace NaNs if >=R13
if (version('-release')>=13)
    indices = find(isnan(im));
    im(indices) = 1;
end

% Finally, set the view.ui.image field
%view.ui.image = uint8(double(im)-1);
view.ui.image = im;

if isempty(overlay)
   view.ui.cbarRange = [];
else
   view.ui.cbarRange = [overClipMin overClipMax];
end

return;


text:
function [in3] = m32in3(m3)
% Convert volume from cubic meters to cubic inches.
% Chad Greene 2012
in3 = m3*61023.744095;

author: Sable

repo: mcbench-benchmarks


text:
function value = r4_csevl ( x, cs, n )

%*****************************************************************************80
%
%% R4_CSEVL evaluates a Chebyshev series.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Broucke,
%    Algorithm 446:
%    Ten Subroutines for the Manipulation of Chebyshev Series,
%    Communications of the ACM,
%    Volume 16, Number 4, April 1973, pages 254-256.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Input, real CS(N), the Chebyshev coefficients.
%
%    Input, integer N, the number of Chebyshev coefficients.
%
%    Output, real VALUE, the Chebyshev series evaluated at X.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_CSEVL - Fatal error!\n' );
    fprintf ( 1, '  Number of terms <= 0.\n' );
    error ( 'R4_CSEVL - Fatal error!' )
  end

  if ( 1000 < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_CSEVL - Fatal error!\n' );
    fprintf ( 1, '  Number of terms > 1000.\n' );
    error ( 'R4_CSEVL - Fatal error!' )
  end

  if ( x < -1.1 || 1.0 < x )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_CSEVL - Fatal error!\n' );
    fprintf ( 1, '  X outside (-1,+1).\n' );
    fprintf ( 1, '  X = %f\n', x );
    error ( 'R4_CSEVL - Fatal error!' )
  end

  b1 = 0.0;
  b0 = 0.0;
  twox = 2.0 * x;

  for i = n : -1 : 1
    b2 = b1;
    b1 = b0;
    b0 = twox * b1 - b2 + cs(i);
  end

  value = 0.5 * ( b0 - b2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@symfun ldivide {(@var{f}, @var{g})}
%% @defopx Operator @@symfun {@var{f} .\ @var{g}} {}
%% Component-wise backslash division of symbolic functions.
%%
%% Simple example:
%% @example
%% @group
%% syms x
%% f(x) = [1 x sin(x)];
%% g(x) = [x x pi];
%% @end group
%%
%% @group
%% h = f .\ g
%%   @result{} h(x) = (symfun)
%%       ⎡        π   ⎤
%%       ⎢x  1  ──────⎥
%%       ⎣      sin(x)⎦
%% @end group
%% @end example
%%
%% @seealso{@@symfun/rdivide}
%% @end defop

function h = ldivide(f, g)
  [vars, s1, s2] = helper_symfun_binops(f, g);
  h = symfun(s1 .\ s2, vars);
end

%!test
%! syms x
%! f(x) = x^2;
%! assert( isa(f .\ f, 'symfun'))
%! assert( isa(f .\ x, 'symfun'))

repo: octsympy


text:
classdef Gaussian2C
%%GAUSSIAN2C A collection of functions for the bivariate Gaussian copulae.
%            The bivariate log-normal and Gaussian copulae have been proven
%            equivalent due to the monotonic relationship between the two
%            distributions [2].
%            Implemented functions include: PDF, CDF, tau, tau2CorrMat.
%
%REFERENCES:
%[1] H. Joe and D. Kurowicka, Dependence Modeling: Vine Copula Handbook.
%    World Scientific, 2011.
%[2] X. Zeng, J. Ren, Z. Wang, S. Marshall, and T. Durrani, "Copulas for
%    statistical signal processing (part i): Extensions and
%    generalization," Signal Processing, vol. 94, pp. 691-702, 2014,
%    ISSN: 0165-1684. DOI: https://doi.org/10.1016/j.sigpro.2013.07.009.
%    [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165168413002880.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    methods(Static)
        function jprob = PDF(u,R)
            %PDF Generates a joint probability for column vector u or
            %    a row vector of joint probabilities for matrix u
            %    with each column as the input vector.
            %
            %INPUTS:
            % u: A 2-by-n matrix where each column represents a
            %    2-dimensional random vector and n is the number of
            %    random vectors. Entries must be in [0,1].
            % R: A 2-by-2 correlation matrix, or a scalar which is the
            %    correlation.
            %
            %OUTPUTS: jprob: A 1-by-n vector of joint probabilities.
            %
            %EXAMPLE 1: Make a contour plot of the density.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            %     for iidx = 1:length(Y)
            %         Z(idx,iidx) = Gaussian2C.PDF([X(idx,iidx);Y(idx,iidx)],0.5);
            %     end
            % end
            % contourf(X,Y,Z)
            %
            %EXAMPLE 2: Generate a joint density plot from random samples
            %           with histograms.
            % rvecs = randn(2,1e4);
            % rho = 0.65;
            % R = chol([1,rho;rho,1],'lower');
            % for idx = 1:size(rvecs,2)
            %     rvecs(:,idx) = R*rvecs(:,idx);
            % end
            % c = corrCoeffSpearman(rvecs(1,:),rvecs(2,:)); %diagonal entries close to rho
            % u1 = GaussianD.CDF(rvecs(1,:));
            % u2 = GaussianD.CDF(rvecs(2,:));
            % pdf = Gaussian2C.PDF([u1;u2],rho).*GaussianD.PDF(rvecs(1,:)).*GaussianD.PDF(rvecs(2,:));
            % fig = jointPlot2D(rvecs(1,:),rvecs(2,:),pdf);
            %
            %EXAMPLE 3: Generate an exact joint density plot.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            % for iidx = 1:length(Y)
            % Z(idx,iidx) = Gaussian2C.PDF([X(idx,iidx);Y(idx,iidx)],0.95);
            % X(idx,iidx) = GumbelD.invCDF(X(idx,iidx),1,2);
            % Y(idx,iidx) = GaussianD.invCDF(Y(idx,iidx));
            % Z(idx,iidx) = Z(idx,iidx).*GumbelD.PDF(X(idx,iidx),1,2).*GaussianD.PDF(Y(idx,iidx));
            % end
            % end
            % surf(X,Y,Z)
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('u','var') || isempty(u)
                jprob = [];
                return
            end
            if ~exist('R','var') || isempty(R)
                rho = 0;
            elseif ~isscalar(R)
                rho = R(2,1);
            else
                rho = R;
            end

            a = sqrt(2)*erfinv(2*u(1,:)-1);
            b = sqrt(2)*erfinv(2*u(2,:)-1);
            c = (0.5*(2*a.*b*rho-rho^2*(a.^2+b.^2))/(1-rho^2));

            jprob = (1/sqrt(1-rho^2))*exp(c);

            %NaN results from subtraction of infinities in the exponential.
            %These should be 0s.
            nanIdxs = isnan(jprob);
            jprob(nanIdxs) = 0;

        end

        function jdist = CDF(u,R)
            %CDF Generates a joint distribution for column vector u or
            %    a row vector of joint probabilities for matrix u
            %    with each column as the input vector.
            %
            %INPUTS:
            % u: A 2-by-n matrix where each column represents a
            %    2-dimensional random vector and n is the number of
            %    random vectors. Entries must be in [0,1].
            % R: A 2-by-2 correlation matrix, or a scalar which is the
            %    correlation.
            %
            %OUTPUTS: jprob: A 1-by-n vector of joint distribution values.
            %
            %Note: This function uses a Monte Carlo method for
            %      approximating the CDF. It will take a few seconds to run
            %      the examples below.
            %
            %EXAMPLE 1: Make a contour plot of the distribution.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            %     for iidx = 1:length(Y)
            %         Z(idx,iidx) = Gaussian2C.CDF([X(idx,iidx);Y(idx,iidx)],0.65);
            %     end
            % end
            % contourf(X,Y,Z)
            %
            %EXAMPLE 2: Generate an exact joint distribution plot.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            %     for iidx = 1:length(Y)
            %         Z(idx,iidx) = Gaussian2C.CDF([X(idx,iidx);Y(idx,iidx)],0.65);
            %         X(idx,iidx) = GumbelD.invCDF(X(idx,iidx),1,2);
            %         Y(idx,iidx) = GaussianD.invCDF(Y(idx,iidx));
            %         Z(idx,iidx) = Z(idx,iidx);
            %     end
            % end
            % tiledlayout(1,2)
            % nexttile
            % surf(X,Y,Z)
            % Zsmooth = zeros(size(Z));
            % Zsmooth(:,end) = Z(:,end);
            % Zsmooth(:,1) = Z(:,1);
            % Zsmooth(end,:) = Z(end,:);
            % Zsmooth(1,:) = Z(1,:);
            % for idx = 2:size(Z,1)-1
            %     for iidx = 2:size(Z,2)-1
            %         Zsmooth(idx,iidx) = sum(Z(idx-1:idx+1,iidx-1:iidx+1),'all')/9;
            %     end
            % end
            % nexttile
            % surf(X,Y,Zsmooth)
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('u','var') || isempty(u)
                jdist = [];
                return
            end
            if ~exist('R','var') || isempty(R)
                rho = 0;
            elseif ~isscalar(R)
                rho = R(2,1);
            else
                rho = R;
            end

            R = [1,rho;rho,1]^2;

            a = sqrt(2)*erfinv(2*u(1,:)-1);
            b = sqrt(2)*erfinv(2*u(2,:)-1);

            %Though this approximation could be used, the function
            %bivarGaussRectangleCDF produces a more exact result. Moreover,
            %it can handle -Inf lower bound, whereas the approximation
            %would have to set some arbitrary llim as the lower bound.
            %  jdist = arrayfun(@(ele1,ele2)GaussianD.integralOverRegion([0;0],R,[llim;llim],[ele1;ele2]),a,b);
            jdist = arrayfun(@(ele1,ele2)bivarGaussRectangleCDF([-Inf;-Inf],[ele1;ele2],[0;0],R),a,b);
        end

        function t = tau(R)
            %%TAU Compute Kendall's tau for the copula.
            %
            %INPUTS:
            % R: A 2-by-2 correlation matrix, or a scalar which is the
            %    correlation.
            %
            %OUTPUTS:
            % t: The scalar value of Kendall's tau.
            %
            %EXAMPLE: Shows Gaussian distribution of independent
            %         correlation parameter estimates using Kendall's
            %         tau.
            % for idx = 1:1e3
            %    x = GaussianD.rand(1e3,[0;0],[1,0.3;0.3,1]);
            %    u = GaussianD.CDF(x);
            %    s1 = GammaD.invCDF(u(1,:),2,1);
            %    s2 = GammaD.invCDF(u(2,:),4,4);
            %    tauEst(idx) = corr(s1',s2','Type','Kendall');
            % end
            % histogram(tauEst)
            % mu = mean(tauEst);
            % sigma = std(tauEst);
            % xline(mu,'k','LineWidth',2,'Label','\mu');
            % xline(mu-sigma,'r--','LineWidth',2,'Label','\mu-\sigma');
            % xline(mu+sigma,'r--','LineWidth',2,'Label','\mu+\sigma');
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('R','var') || isempty(R)
                rho = 0;
            elseif ~isscalar(R)
                rho = R(2,1);
            else
                rho = R;
            end

            t = 2*arcsin(rho)/pi;
        end

        function R = tau2corrMat(t)
            %%TAU2CORRMAT Compute a 2-by-2 correlation matrix estimate
            %             from Kendall's tau.
            %
            %INPUTS:
            % t: Kendall's tau
            %
            %OUTPUTS:
            % R: A 2-by-2 correlation matrix
            %
            %EXAMPLE 1: Generates Gamma distributed marginals correlated
            %           via a Gaussian copula and estimates the correlation
            %           matrix.
            % x = GaussianD.rand(1e4,[0;0],[1,0.3;0.3,1]);
            % u = GaussianD.CDF(x);
            % s1 = GammaD.invCDF(u(1,:),2,1);
            % s2 = GammaD.invCDF(u(2,:),4,4);
            % tau = corr(s1',s2','Type','Kendall');
            % Gaussian2C.tau2corrMat(tau)
            %
            %EXAMPLE 2: Shows Gaussian distribution of independent
            %           correlation parameter estimates using Kendall's
            %           tau.
            % for idx = 1:1e3
            %    x = GaussianD.rand(1e3,[0;0],[1,0.3;0.3,1]);
            %    u = GaussianD.CDF(x);
            %    s1 = GammaD.invCDF(u(1,:),2,1);
            %    s2 = GammaD.invCDF(u(2,:),4,4);
            %    tauEst(idx) = corr(s1',s2','Type','Kendall');
            %    rho = Gaussian2C.tau2corrMat(tauEst(idx));
            %    theta(idx) = rho(1,2);
            % end
            % histogram(theta)
            % mu = mean(theta);
            % sigma = std(theta);
            % xline(mu,'k','LineWidth',2,'Label','\mu');
            % xline(mu-sigma,'r--','LineWidth',2,'Label','\mu-\sigma');
            % xline(mu+sigma,'r--','LineWidth',2,'Label','\mu+\sigma');
            %
            %April 2022 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            rho = sin(t*pi/2);
            R = [1,rho;rho,1];
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end


text:
function [x,y] = swap(x,y)
tempx = x;
x = y;
y = tempx;
clear tempx;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
% PlotSphereIntensity(azimuth, elevation)
% PlotSphereIntensity(azimuth, elevation, intensity)
% h = PlotSphereIntensity(...)
%
% Plots the intensity (as color) of a number of points on a unit sphere.
% Input:
%   azimuth (phi), in degrees
%   elevation (theta), in degrees
%   intensity (optional, if not provided, a green sphere is produced)
%   All inputs must be vectors or matrices of the same size.
%   Data does not have to be evenly spaced. When there aren't enough points
%   to draw a smooth sphere, additional points (with color) are
%   interpolated.
% Output:
%   h - a handle to the patch object
%
% The axes are also plotted:
%   positive x axis is red
%   positive y axis is green
%   positive z axis is blue
%
% Author: David Johnstone (DSTO, WSD, Australia)
%         david.johnstone@dsto.defence.gov.au
%         Ninh Duong (DSTO, WSD, Australia)
%         ninh.duong@dsto.defence.gov.au
% Date: 18 March, 2008

function varargout = PlotSphereIntensity(az, el, varargin)

error(nargchk(2, 3, nargin))
error(nargoutchk(0, 1, nargout))

if (any(size(az) ~= size(el)) || (nargin == 3 && any(size(az) ~= size(varargin{1}))))
    error('input data has inconsistent sizes')
end

az = az(:) * pi / 180;
el = el(:) * pi / 180;

if (nargin == 3)
    c = varargin{1};
    c = c(:);
else
    c = ones(length(az),1);
end

% Construct cartesian points (x, y, z) from the given azimuth and elevation
% data, so that each azimuth elevation data pair are a point on the surface
% of a unit sphere.
x = cos(el) .* cos(az);
y = cos(el) .* sin(az);
z = sin(el);

[x y z c] = Sphericalise(x, y, z, c);

% Create a list of facets.
K = convhulln([x y z]);

% Plot the data
h = trisurf(K, x, y, z, c);
shading interp
axis equal
view([1 1 1])

plotSetting = get(gca, 'NextPlot');
hold on

xAxis = [0 0 0; 10 0 0];
yAxis = [0 0 0; 0 10 0];
zAxis = [0 0 0; 0 0 10];
plotAxis(xAxis, 'r', 'LineWidth', 2)
plotAxis(yAxis, 'g', 'LineWidth', 2)
plotAxis(zAxis, 'b', 'LineWidth', 2)

set(gca, 'NextPlot', plotSetting);

text(1.4,0,0.04,'x', 'FontSize', 16)
text(0,1.4,0.04,'y', 'FontSize', 16)
text(0.04,0.04,1.4,'z', 'FontSize', 16)

if (nargout == 1)
    varargout{1} = h;
end

end

function plotAxis(axis, varargin)
plot3(axis(:,1), axis(:,2), axis(:,3), varargin{:})
end

% [x2 y2 z2 C2] = Sphericalise(x, y, z, C)
%
% Creates a unit sphere with color data linearly interpolated across it
% based on a given set of points.
% Input:
%    x, y and z are equal sized vectors/matrices and describe the
%        points on a sphere.
%    C is the color data, and is the same size as x, y and z.
%    x, y, z and C must all be the same size. They are used as nx1
%        matrices with x(1), y(1), z(1) describing the location of the
%        first point, and C(1) containing the color for it.
% Output:
%   x, y, z and color so that the data now looks like a unit sphere.

% Method:
%   This algorithm works by triangulating the data and splitting any lines
%   that are longer than a certain threshold (chosen as pi/10 as this gives
%   a fairly smooth sphere.) convhulln is used to triangulate the surface
%   of the sphere, and the indices it returns are used to identify the
%   start and end points of the lines. If the line across the surface of
%   the sphere is too long, a new point is created halfway between the
%   endpoints, and it is scaled to have a magnitude of 1 so that it will
%   lie on the surface of the unit sphere. After all the lines have been
%   checked (and split if necessary), this process is repeated if any lines
%   were split. This repeating must be done as the triangulation will
%   introduce new line segments with the new point, and these may be too
%   long.

function [x2 y2 z2 C2] = Sphericalise(x, y, z, C)

splitLine = 1; % set to 1 initially as we want it to do the loop at least
% once, and MATLAB doesn't have a do while loop

x2 = x(:);
y2 = y(:);
z2 = z(:);
C2 = C(:);

while splitLine
    T = convhulln([x2 y2 z2]);

    % It is best to split lines on a per line basis (and not a per triangle
    % basis) as each line is shared by two triangles, so any line that is to be
    % split will be split twice (and the new point will be in the same place)
    lines = unique([T(:,1) T(:,2); T(:,1) T(:,3); T(:,2) T(:,3)], 'rows');

    splitLine = 0;
    for line = lines'

        A = struct('x', x2(line(1)), 'y', y2(line(1)), 'z', z2(line(1)), 'C', C2(line(1)));
        B = struct('x', x2(line(2)), 'y', y2(line(2)), 'z', z2(line(2)), 'C', C2(line(2)));

        % if the length of the geodesic is longer than pi/10
        if acos(dot2(A, B)) > pi/10
            splitLine = 1;

            % add point between A and B

            % It is possible to calculate elevation and azimuth for A and
            % B, and then interpolate this to get elevation and azimuth for
            % P, and finally convert these back to cartesian coordinates,
            % but this 1) is inefficient and 2) doesn't work around the
            % poles (at least if theta and phi are linearly interpolated).

            P.x = (A.x + B.x) / 2;
            P.y = (A.y + B.y) / 2;
            P.z = (A.z + B.z) / 2;
            magnitude = sqrt(P.x^2 + P.y^2 + P.z^2);
            P.x = P.x / magnitude;
            P.y = P.y / magnitude;
            P.z = P.z / magnitude;

            P.C = (A.C + B.C) / 2;

            x2(end+1) = P.x;
            y2(end+1) = P.y;
            z2(end+1) = P.z;
            C2(end+1) = P.C;

        end
    end
end

end

function d = dot2(A, B)
d = dot([A.x A.y A.z], [B.x B.y B.z]);
end


text:
% this script tests 1st Order Statistics features between CERR and pyradiomics on a Wavelet filetered image.
%
% RKP, 03/22/2018

%% Load image
firstOrderParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','tests_for_cerr','test_first_order_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','data_for_cerr_tests','CERR_plans','head_neck_ex1_20may03.mat.bz2');

planC = loadPlanC(cerrFileName,tempdir);
indexS = planC{end};

paramS = getRadiomicsParamTemplate(firstOrderParamFileName);
strNum = getMatchingIndex(paramS.structuresC{1},{planC{indexS.structures}.structureName});
scanNum = getStructureAssociatedScan(strNum,planC);

%wavType = 'coif1';
%scanType = 'Original';
scanType = 'Wavelet';
%dirString = 'HLH';
dirString = paramS.imageType.Wavelets.Direction.val;

%% Calculate features using CERR

firstOrderS = calcGlobalRadiomicsFeatures...
            (scanNum, strNum, paramS, planC);
firstOrderS = firstOrderS.(['Wavelets_coif_1_',dirString]).firstOrderS;
cerrFirstOrderV = [firstOrderS.energy, firstOrderS.totalEnergy, firstOrderS.interQuartileRange, ...
    firstOrderS.kurtosis+3, firstOrderS.max, firstOrderS.mean, firstOrderS.meanAbsDev, ...
    firstOrderS.median, firstOrderS.medianAbsDev, firstOrderS.min, ...
    firstOrderS.P10, firstOrderS.P90, firstOrderS.interQuartileRange, ...
    firstOrderS.robustMeanAbsDev, firstOrderS.rms, firstOrderS.skewness, ...
    firstOrderS.std, firstOrderS.var, firstOrderS.entropy];

% %% Calculate features using pyradiomics
%
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
%     single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% pyradFirstorderNamC = {'Energy', 'TotalEnergy','InterquartileRange','Kurtosis',...
%     'Maximum', 'Mean','MeanAbsoluteDeviation','Median','medianAbsDev',...
%     'Minimum','10Percentile','90Percentile','InterquartileRange',...
%     'RobustMeanAbsoluteDeviation','RootMeanSquared','Skewness',...
%     'StandardDeviation','Variance','Entropy'};
%
% pyradFirstorderNamC = strcat(['wavelet', '_', dirString, '_firstorder_'],pyradFirstorderNamC);
% %pyradFirstorderNamC = strcat(['original', '_firstorder_'],pyradFirstorderNamC);
%
% pyRadFirstOrderV = [];
% for i = 1:length(pyradFirstorderNamC)
%     if isfield(teststruct,pyradFirstorderNamC{i})
%         pyRadFirstOrderV(i) = teststruct.(pyradFirstorderNamC{i});
%     else
%         pyRadFirstOrderV(i) = NaN;
%     end
% end
%
% %% Compare
% diffFirstOrderV = (cerrFirstOrderV - pyRadFirstOrderV) ./ cerrFirstOrderV * 100

%% Compare with previously calculated values of pyradiomics first order features
saved_pyRadFirstOrderV = [87498183.4756395,1049978201.70767,30.4772882461548,16.3069672735214,836.854187011719,-0.191955631352121,48.4171827693622,0.0908048748970032,NaN,-718.366882324219,-73.1680297851562,64.9023208618165,30.4772882461548,16.4187532424825,96.2495122333083,0.546840923205712,NaN,9263.93175818536,3.17432864851443];
diffFirstOrderV = (cerrFirstOrderV - saved_pyRadFirstOrderV) ./ cerrFirstOrderV * 100

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%     Intra Database Experiments     %%%%%
%%%%%     SFA (Li et al. TMM 2018)       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% written by Dingquan Li
% dingquanli@pku.edu.cn
% IDM, SMS, PKU
% Last update: Sept. 15, 2018

clear;clc

caffe_path = '/home/ldq/caffe/matlab/'; % point to the caffe path
addpath(genpath(caffe_path));

database = 'MLIVE2';
% im_dir, im_ids, im_lists, im_names, index, ref_ids, subjective_scores, subjective_scoresSTD
load(['data/' database 'info']); % point to data info file, im_dir should be re-assigned.
% for i = 1:length(im_names) % Update im_lists when im_dir changes
%     im_lists{i} = [im_dir im_names{i}];
% end
% clear i
if ~exist('results','dir')
    mkdir('results');
end
save_path = ['./results/' database '-reproduce']; % point to save path

%% Feature Extraction
layer_names = {'res3d', 'res4f', 'res5c'};
options = 'none';
for k = 1:length(layer_names)
    layer_name = layer_names{k};
    eval(['[feature1.' layer_name ', feature2.' layer_name ', feature3.' layer_name ...
        '] = SFAfeature(im_lists, layer_name, options);']);
end

%%
switch options
    case 'none'
        Rep = 1;
    case 'flipH'
        Rep = 2;
    case 'clipUL'
        T = 6;
        Rep = T^2;
end

N = length(im_lists);
train_ratio = 0.8;
for t = 1 : size(index, 1)
    fprintf('the %d-th iteration\n',t);

    srocc = zeros(length(layer_names), 3, 5);
    for k = 1:length(layer_names)
        layer_name = layer_names{k};
        im_index = index(t,1:ceil(train_ratio*size(index,2)));
        I = buffer(im_index,5);
        for i = 1:5
            train_im_index = I;
            train_im_index(i,:) = [];
            train_im_index(train_im_index==0) = [];
            test_im_index = I(i,:);
            test_im_index(test_im_index==0) = [];
            train_im_index = cell2mat(arrayfun(@(i)find(ref_ids==train_im_index(i))',...
                1:length(train_im_index),'UniformOutput',false));
            test_im_index = cell2mat(arrayfun(@(i)find(ref_ids==test_im_index(i))',...
                1:length(test_im_index),'UniformOutput',false));
            train_labels = subjective_scores(train_im_index); %#ok<NASGU>
            test_labels = subjective_scores(test_im_index);

            p = 10; % number of components;
            eval(['[~,~,~,~,beta1] = plsregress(feature1.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
            eval(['predict_statistics1 = [ones(length(test_im_index(:)),1) feature1.' layer_name '(test_im_index, :)]*beta1;']);
            eval(['[~,~,~,~,beta2] = plsregress(feature2.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
            eval(['predict_statistics2 = [ones(length(test_im_index(:)),1) feature2.' layer_name '(test_im_index, :)]*beta2;']);
            eval(['[~,~,~,~,beta3] = plsregress(feature3.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
            eval(['predict_statistics3 = [ones(length(test_im_index(:)),1) feature3.' layer_name '(test_im_index, :)]*beta3;']);
        %
            predict_statistics = [predict_statistics1 predict_statistics2  predict_statistics3];

            srocc(k,:,i) = corr(predict_statistics, test_labels, 'type', 'Spearman');
        end
    end

    train_im_index = index(t,1:ceil(train_ratio*size(index,2)));
    train_im_index = cell2mat(arrayfun(@(i)find(ref_ids==train_im_index(i))',...
        1:length(train_im_index),'UniformOutput',false));
    test_im_index = index(t,1+ceil(train_ratio*size(index,2)):size(index,2));
    test_im_index = cell2mat(arrayfun(@(i)find(ref_ids==test_im_index(i))',...
        1:length(test_im_index),'UniformOutput',false));
    train_labels = subjective_scores(train_im_index);
    test_labels = subjective_scores(test_im_index);

    Avesrocc = mean(srocc,3);
    [bestK, bestS] = find(Avesrocc==max(Avesrocc(:)),1,'first'); %

    layer_name = layer_names{bestK};
    p = 10; % number of components;
    clear predict_statistics
    eval(['[~,~,~,~,beta1] = plsregress(feature1.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
    eval(['predict_statistics(:, 1) = [ones(length(test_im_index(:)),1) feature1.' layer_name '(test_im_index, :)]*beta1;']);
    eval(['[~,~,~,~,beta2] = plsregress(feature2.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
    eval(['predict_statistics(:, 2) = [ones(length(test_im_index(:)),1) feature2.' layer_name '(test_im_index, :)]*beta2;']);
    eval(['[~,~,~,~,beta3] = plsregress(feature3.' layer_name '(train_im_index(:)+[0:N:(Rep-1)*N],:),repmat(train_labels,Rep,1),p);']);
    eval(['predict_statistics(:, 3) = [ones(length(test_im_index(:)),1) feature3.' layer_name '(test_im_index, :)]*beta3;']);
    w = Avesrocc(bestK,:)'/sum(Avesrocc(bestK,:));
    objective_scores = predict_statistics*w;
    results.srocc(t) = corr(objective_scores, test_labels, 'type', 'Spearman');
    results.plcc(t) = corr(objective_scores, test_labels);
    results.krocc(t) = corr(objective_scores, test_labels, 'type', 'Kendall');
    results.RMSE(t) = sqrt(mean((objective_scores-test_labels).^2));
    results.OR(t) = mean(abs(objective_scores-test_labels)>2*subjective_scoresSTD(test_im_index));
    median(results.srocc)
end
save(save_path);

repo: SFA


text:
classdef ThrottleModelEnum < matlab.mixin.SetGet
    %ThrottleModelEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        PolyModel('Polynominal Model','ThrottlePolyModel',"The vehicle throttle varies as a function of a second degree polynomial.  Can set constant throttle here too.");
        T2WModel('Thrust To Weight Model','T2WThrottleModel',"The vehicle will set the throttle to maintain a fixed thrust to weight ratio.");
    end

    properties
        nameStr char = '';
        classNameStr char = '';
        desc(1,1) string
    end

    methods
        function obj = ThrottleModelEnum(nameStr,classNameStr,desc)
            obj.nameStr = nameStr;
            obj.classNameStr = classNameStr;
            obj.desc = desc;
        end
    end

    methods(Static)
        function [throttleModelNameStrs, m] = getThrottleModelTypeNameStrs()
            [m,~] = enumeration('ThrottleModelEnum');

            throttleModelNameStrs = {};
            for(i=1:length(m)) %#ok<*NO4LP>
                throttleModelNameStrs{end+1} = m(i).nameStr; %#ok<AGROW>
            end
        end

        function ind = getIndOfListboxStrsForThrottleModel(throttleModel)
            [m,~] = enumeration('ThrottleModelEnum');
            inputClass = class(throttleModel);

            ind = -1;
            for(i=1:length(m))
                if(strcmpi(m(i).classNameStr,inputClass))
                    ind = i;
                    break;
                end
            end
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ThrottleModelEnum');
            ind = find(ismember({m.nameStr},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function data = xdf2fieldtrip(filename, varargin)

% XDF2FIELDTRIP reads data from a XDF file with multiple streams. It upsamples the
% data of all streams to the highest sampling rate and concatenates all channels in
% all streams into a raw data structure that is compatible with the output of
% FT_PREPROCESSING.
%
% Use as
%   data = xdf2fieldtrip(filename, ...)
%
% Optional arguments should come in key-value pairs and can include
%   streamindx  = list, indices of the streams to read (default is all)
%
% You can also use the standard procedure with FT_DEFINETRIAL and FT_PREPROCESSING
% for XDF files. This will return (only) the stream with the highest sampling rate,
% which is typically the EEG.
%
% You can use FT_READ_EVENT to read the events from the non-continuous data streams.
% To get them aligned with the samples in one of the specific data streams, you
% should specify the corresponding header structure.
%
% See also FT_PREPROCESSING, FT_DEFINETRIAL, FT_REDEFINETRIAL

% Copyright (C) 2019, 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$

% process the options
streamindx = ft_getopt(varargin, 'streamindx');

% ensure this is on the path
ft_hastoolbox('xdf', 1);

% read all streams
streams = load_xdf(filename);

iscontinuous = false(size(streams));
% figure out which streams contain continuous/regular and discrete/irregular data
for i=1:numel(streams)
  iscontinuous(i) = isfield(streams{i}.info, 'effective_srate');
end

% give some feedback
for i=1:numel(streams)
  if iscontinuous(i)
    ft_info('stream %d contains continuous %s data\n', i, streams{i}.info.name);
  else
    ft_info('stream %d contains non-continuous %s data\n', i, streams{i}.info.name);
  end
end

% select the streams to continue working with
if isempty(streamindx)
  selected = true(size(streams));
else
  selected = false(size(streams));
  selected(streamindx) = true;
end

% discard the non-continuous streams
streams = streams(iscontinuous & selected);

if isempty(streams)
  ft_error('no continuous streams were selected');
end

% convert each continuous stream into a FieldTrip raw data structure
data = cell(size(streams));
for i=1:numel(streams)

  % make a copy for convenience
  stream = streams{i};

  % this section of code is shared with fileio/private/sccn_xdf
  hdr             = [];
  hdr.Fs          = stream.info.effective_srate;
  hdr.nChans      = numel(stream.info.desc.channels.channel);
  hdr.nSamplesPre = 0;
  hdr.nSamples    = length(stream.time_stamps);
  hdr.nTrials     = 1;
  hdr.label       = cell(hdr.nChans, 1);
  hdr.chantype    = cell(hdr.nChans, 1);
  hdr.chanunit    = cell(hdr.nChans, 1);

  prefix = stream.info.name;

  for j=1:hdr.nChans
    hdr.label{j} = [prefix '_' stream.info.desc.channels.channel{j}.label];
    hdr.chantype{j} = stream.info.desc.channels.channel{j}.type;
    hdr.chanunit{j} = stream.info.desc.channels.channel{j}.unit;
  end

  hdr.FirstTimeStamp     = stream.time_stamps(1);
  hdr.TimeStampPerSample = (stream.time_stamps(end)-stream.time_stamps(1)) / (length(stream.time_stamps) - 1);

  % keep the original header details
  hdr.orig = stream.info;

  data{i}.hdr = hdr;
  data{i}.label = hdr.label;
  data{i}.time = {streams{i}.time_stamps};
  data{i}.trial = {streams{i}.time_series};

end % for all continuous streams

% determine the continuous stream with the highest sampling rate
srate = nan(size(streams));
for i=1:numel(streams)
  srate(i) = streams{i}.info.effective_srate;
end
[dum, indx] = max(srate);

if numel(data)>1
  % resample all data structures, except the one with the max sampling rate
  % this will also align the time axes
  for i=1:numel(data)
    if i==indx
      continue
    end

    ft_notice('resampling %s', streams{i}.info.name);
    cfg = [];
    cfg.time = data{indx}.time;
    data{i} = ft_resampledata(cfg, data{i});
  end

  % append all data structures
  data = ft_appenddata([], data{:});
else
  % simply return the first and only one
  data = data{1};
end

author: spm


text:
function bm=EWT_beta(x)

% function used in the construction of Meyer's wavelet
if x<0
    bm=0;
elseif x>1
    bm=1;
else
    bm=x^4*(35-84*x+70*x^2-20*x^3);
end


text:
result = max(pfz(600851475143));
save myResult.mat result
printf("The greates prime factor is %d\n",result)

author: TheAlgorithms

repo: MATLAB-Octave


text:
%% Test the niSelect function
%
% We identify a small region of a white matter mask.  We can use this
% region to limit the computation that we do for tracking and LiFE
%
% RL/BW Vistasoft Team, 2016

%% Go get a white matter mask from the Remote Data client

rdt = RdtClient('vistasoft');
rdt.crp('/vistadata/diffusion/sampleData/dti40/bin');

% This is all of the artifacts
a = rdt.listArtifacts('print',true);

% Write the white matter mask out in the local directory
dFolder = fullfile(vistaRootPath,'local');
fname = rdt.readArtifact(a(7),'destinationFolder',dFolder);

%% Have a look
nii = niftiRead(fname);
montageAll = niftiView(nii);

%% Select out a cube

p.keepLR = 20:30; p.keepPA = 20:30; p.keepIS = 40:50;
niiCube = niftiSelect(nii,p);
montage = niftiView(niiCube);

%% Make an overlay from the
montageOverlay(:,:,1) = montage(:,:);
montageOverlay(:,:,2) = montageAll;
montageOverlay(:,:,3) = montageAll;
mrvNewGraphWin; imagesc(montageOverlay); axis image

%%


text:
function new_doc=publishdepfun(func,opts,add_fun,rm_fun)
%%
% publishdepfun creates a single published html file with called
% functions attached to the end of the root mfile. Hyperlinks will be
% within a "Called Functions" in the Table of Contents area of the
% published document that link to the called function.  Functions that are
% called that are under MATLABROOT are not included in the published
% output.
%
% publishdepfun only goes one level deep....
%
% Example: If fun1 calls fun2, fun2 will be published. If fun2 calls fun3,
%          fun3 will not be published.
%
% Calling: new_doc = publishdepfun(func,opts,add_fun)
%
% Inputs:
%
% * func => character string of root function or script ('Function1').
%           The ".m" is not needed.
% * opts => is the list of options for the publishing. See publish.m help
%           for more information.
% * add_fun => is a cell array of additional functions to be published.
% * rm_fun => is a cell array of functions to remove from publish list.
%
%
% Outputs:
%
% * new_doc => the final html document with attached subfunctions.
%                    Name is the root file name with "_withfuncs" appended.
%
%
%
% Author: Nick Angelini
%
% Version: 2.1
%
% Date: 11/29/2011
%
%
%

%%
% Changelog
%
% Version 2.3
%
% # added ability to cell array of functions to remove from publish list
%
%
% Version 2.2
%
% # added ability to remove mex files from output.
%
%
% Version 2.1
%
% # added support for displaying published html document if no output
% variable is specified.
%
%
% Version 2
%
% # added ability to pass additional functions as a cell array to be published.
% # added '-calltree' opition to the depfun call to find the functions
% passed into other functions like ode's and quad.
%
%
%%
% Known Issues
%
% # Likely dealing with publish options structure isn't honored correctly
% across publishing root and children.
% # Only some basic testing on HTML output, no testing and not expected to
% work for other outputs
% # Can't call publishdepfun on functions that call publish (itself or
% others).  Hit some errors that haven't been explored.
% # Current limitations are if depfun does not return function name it will
% not get published unless the user inputs a secondary list of
% functions. Examples can be found in the help documentation for
% depfun.
% # Currently requires cell blocks to be used, with atleast one block with
% a title other than the overall title of the root file. Need a "Contents"
% section to be created.

%% Error Checking
if nargin<1,error('at least pass in a function'),end
% Default to HTML publishing.
if (nargin < 2)
    opts = 'html';
end
if nargin > 2 && isempty(opts)
    opts = 'html';
end
if nargin < 3
    add_fun = [];
end
if nargin < 4
    rm_fun = [];
end

%% get functions at the top level only
list=depfun(func,'-toponly','-quiet','-calltree');
out=list;

%% Remove hits that are in matlabroot
out(strncmp(matlabroot,list,length(matlabroot)))=[];

%% Remove hits that are probably MEX functions
idx=strfind(out,mexext);
for i=1:length(idx),
    idxs(i) = ~isempty(idx{i});
end
out(idxs)=[];

%% Add user list to functions to be published
out = [out;add_fun];

%% Remove user list
for i=1:length(rm_fun)
    idx=strfind(out,rm_fun{i});
    for i=1:length(idx),
        idxs(i) = ~isempty(idx{i});
    end
    out(idxs)=[];
end

%% Publish main/root/caller function
out_doc{1}=publish(out{1},opts);

%% Publish Dependent Functions
options.evalCode=false;
options.outputDir='';

for i=2:length(out)
    out_doc{i}=publish(out{i},options);
end

%% Initialize variables for combining html files
% check is to find id the correct line to add html code, and check2 is to
% find the correct spot in the line to add code for the "Called Functions"
% heading.

check = '</style>';
check2 = '
';
header = sprintf('
Called Functions

');

%% Collect into one HTML

% get name of main html doc
[pathstr, name, ~] = fileparts(out_doc{1});

% create new html with same file name + _withfuncs label
new_doc=[pathstr filesep name '_withfuncs.html'];
fid_new=fopen(new_doc,'w');

% append html docs together
for i=1:length(out_doc)
    % if first/main html
    if i == 1
        % open first html file. not the _withfuncs label
        fid1 = fopen(out_doc{1});

        % loop through line by line of first html doc
        while 1
            x = fgetl(fid1);
            % break out of loop if end of doc is reached
            if x == -1
                break
            end
            % see if current line matches the check variable
            if ~isempty(regexp(x,check,'ONCE'))
                % find location in check line to place subfunction heading
                id = regexp(x,check2);
                % output current line up to char before check2 + header
                counts = fprintf(fid_new,'%s',[x(1:id-1) header]);

                % loop through out_doc adding hyperlink code for
                % subfunctions
                for j=2:length(out_doc)
                    % get function name
                    [~, mname, ~] = fileparts(out_doc{j});
                    % store function name for hyperlink flag
                    namestore{j} = mname;
                    % create link code and append to current line
                    newLink = sprintf('
%s
',mname,mname);
                    count = fprintf(fid_new,'%s',newLink);
                    % when done with hyperlinks append rest of original
                    % line
                    if j == length(out_doc)
                        count = fprintf(fid_new,'
%s\n',x(id:end));
                    end
                end
                % if not line of interest print the line
            else
                counts = fprintf(fid_new,'%s\n',x);
            end
        end
        % close the original driver html
        fclose(fid1);
        % if not first html doc print the hyperlink flag for subfunction then
        % copy the subfunction html to end of _withfuncs.html
    else
        marker = sprintf('',namestore{i});
        % read in children html publised html documents
        fileStr = fileread(out_doc{i});

        % looking for flags used by function grabcode to allow grabcode to
        % grab only the root function
        matches = regexp(fileStr,'(##### SOURCE BEGIN #####\n.*\n##### SOURCE END #####)','tokens','once');
        newStr = strrep(fileStr,matches{1},'\n');
        count=fprintf(fid_new,'%s\n%s\n',marker,newStr);
    end
end
fclose(fid_new);
if nargout == 0
    web(new_doc)
end


text:
function [str_is_found] = instr(str, pattern)
    % implementation fo contains function, not presents in all matlab
    %
    % SYNTAX:
    %   str_is_found = instr(str, pattern)
    %

    %--- * --. --- --. .--. ... * ---------------------------------------------
    %               ___ ___ ___
    %     __ _ ___ / __| _ | __|
    %    / _` / _ \ (_ |  _|__ \
    %    \__, \___/\___|_| |___/
    %    |___/                    v 1.0RC1
    %
    %--------------------------------------------------------------------------
    %  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
    %  Written by:       Andrea Gatti
    %  Contributors:     Andrea Gatti ...
    %  A list of all the historical goGPS contributors is in CREDITS.nfo
    %--------------------------------------------------------------------------
    %
    %   This program is free software: you can redistribute it and/or modify
    %   it under the terms of the GNU General Public License as published by
    %   the Free Software Foundation, either version 3 of the License, or
    %   (at your option) any later version.
    %
    %   This program is distributed in the hope that it will be useful,
    %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    %   GNU General Public License for more details.
    %
    %   You should have received a copy of the GNU General Public License
    %   along with this program.  If not, see <http://www.gnu.org/licenses/>.
    %
    %--------------------------------------------------------------------------
    % 01100111 01101111 01000111 01010000 01010011
    %--------------------------------------------------------------------------

    str_is_found = ~isempty(strfind(str,pattern)); %#ok<STREMP>
end

author: goGPS-Project


text:
function kern = rbfinfwhiteKernExpandParam(kern, params)

% RBFINFWHITEKERNEXPANDPARAM Create kernel structure from RBF-WHITE kernel's
% parameters (with integration limits between minus infinity and infinity).
% FORMAT
% DESC returns a RBF-WHITE kernel structure filled with the parameters in
% the given vector. This is used as a helper function to enable parameters
% to be optimised in, for example, the NETLAB optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : rbfinfwhiteKernParamInit, rbfinfwhiteKernExtractParam, kernExpandParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

kern.inverseWidth = params(1);
kern.variance = params(2);

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [model,stats] = utl_searchmodel(varargin)
% Find the best predictive model out of a parameterized set, via cross-validation.
% [Model,Stats] = utl_searchmodel(Data, Arguments...)
%
% When a predictive model is learned from some data, several of its key parameters can usually be
% efficiently optimized given the data (namely, if their solutions are available in closed form, or
% if they can be found by solving a well-behaved optimization problem). Some parameters, however,
% usually remain, which can not be optimized by the given learning function. If their optimal values
% are unknown but assumed to lie in a small set of possibilities, a general-purpose approach to find
% them as well, is exhaustive search over all possibilities, coupled with a cross-validation to
% assess the predictive performance for each possibility on the available data.
%
% This function performs this task, and can jointly optimize all unknown parameters of a learning
% function. It requires a data set and (because a cross-validation is implicitly performed) all
% parameters that would be required to apply utl_crossval to the given data. The function whose
% parameters shall be optimized is the 'trainer' (i.e., model learning) function. For any of its
% parameters which has multiple possibilities (over which to optimize), these possibilities can be
% specified using the same syntax as in the grid searching function utl_gridsearch (where the
% default here is to specify them via search() clauses).
%
% In:
%   Data :   some data that can be partitioned using index sets
%
%   Arguments : mandatory arguments:
%               'trainer': training function; receives a partition of the data (as produced by
%                          the partitioner), some further arguments (as specified in args), and
%                          returns a model (of any kind)
%
%               'tester' : testing function; receives a model (as produced by the trainer) and a
%                          partition of the data (as produced by the partitioner), and returns a
%                          prediction for every index of the input data, in a format allowed by
%                          ml_predict
%
%               'args': arguments to the training function, cell array (default: empty)
%                       specified as in utl_gridsearch (format controlled via argform)
%
%               optional arguments (same as in utl_crossval, listed below for convenience)
%               'scheme': cross-validation scheme, can be one of the following: (default: 10)
%                 * 0: skip CV, return NaN and empty statistics
%                 * k: k-fold randomized CV (or, if 0<k<1, k-holdhout CV)
%                 * [r k]: r-time repartitioned k-fold randomized CV (or, if 0<k<1,
%                   k-holdout CV (with k a fraction))
%                 * 'loo': leave-one-out CV
%                 * {'chron', k} or {'block', k}: k-fold chronological/blockwise CV
%                 * {'chron', k, m} or {'block', k, m}: k-fold chronological/blockwise CV with
%                   m indices margin width (between training and test set)
%
%               'partitioner': partitioning function for the data, receives three parameters:
%                              (data, index vector, packed trainer args OR model)
%                               * if the index vector is empty, should return the highest index
%                                 in the data
%                               * otherwise, it should return data subindexed by the index vector
%                              default: provides support for cell arrays, numeric arrays, struct
%                                       arrays and {Data,Target} cell arrays
%                              note: the third parameter is for convenience and may optionally be
%                                    taken into account; trainer args are passed packed into a cell
%                                    array for both index set generation and computation of the
%                                    training partition(s), and the model is passed for the
%                                    computation of testing partitions
%
%               'target': a function to derive the target variable from a partition of the data (as
%                         produced by the partitioner), for evaluation; the allowed format is
%                         anything that may be output by ml_predict default: provides support for
%                         {Data,Target} cell arrays
%
%               'metric': metric to be employed, applied both in each fold and the aggregated data
%                         over all folds
%                          * function handle: a custom, user-supplied loss function; receives target
%                            data in the first argument and prediction data in the second argument;
%                            each can be in any format that can be produced by ml_predict (but can
%                            be expected to be mutually consistent). shall return a real number
%                            indicating the summary metric over all data, and optionally additional
%                            statistics in a struct
%                          * string: use ml_calcloss, with 'metric' determining the loss type
%                          * default: use 'mcr','mse','nll','kld', depending on supplied target &
%                            prediction data formats
%
%               'argform': format of the argument ranges, either 'direct' or 'clauses'
%                          (default: 'clauses')
%                           'direct': search ranges are directly specified as arrays
%                           'clauses': search ranges are specified using the search() clause
%
%               'forcestats': enforce a cross-validation to obtain stats, even when no search is
%                             necessary (default: 0)
%
%               further arguments (same as in utl_crossval)
%               'cache_fold_results' : whether to cache the per-fold results. (default: false)
%
%               'only_cached_results' : load only results that are in the cache. (default: false)
%
%               'tolerate_exceptions' : tolerate exceptions during training step. (default: false)
%
%               further arguments (same as in par_beginschedule, listed below for convenience)
%               'engine_gs': the parallelization engine to be used for the grid search
%                            (default: 'local')
%
%               'engine_ncv': the parallelization engine to be used for the nested cross-validation
%                             (default: 'global')
%
%               'pool'     : node pool to be used for parallelization, when using the BLS scheduler
%                            (default: 'global')
%
%               'policy'   : scheduling policy to be used, when using the BLS scheduler
%                            (default: 'global')
%
% Out:
%   Model : the best model, as determined through cross-validation and and grid search
%           w.r.t. to a measure of the overall performance of the trainer/tester combination
%
%   Stats : additional statistics over all tested argument combinations, as produced by the metric
%
% See also:
%   utl_crossval, utl_gridsearch, utl_nested_crossval
%
% Example:
%   % assuming a feature matrix called trials and a label vector called targets, sized as:
%   %  trials: [NxF] array of training instances
%   %  targets: [Nx1] array of labels
%
%   % find the best-performing SVM classifier on the data (here: ignoring the kernel scale parameter)
%   [model,stats] = utl_searchmodel({trials,targets}, 'args',{{'svm' search(2.^(-5:2:15))}})
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-22
dp;

% parse arguments
opts = arg_define(0:1,varargin, ...
    ... % training data
    arg_norep({'data','Data'},[],[],'Data that can be partitioned using index sets. Such as, for example, {X,y} with X being the [NxF] training data and y being the [Nx1] labels (N=#trials, F=#features).'), ...
    ... % method to use and its arguments
    arg({'args','TrainingArguments'},{},[],'Arguments to the training function. Packed in a cell array. Search ranges are defined according to ArgumentFormat; see also utl_gridsearch for additional clarification on this. Note: if using the default trainer/tester functions, args must at least specify the learning function to be used, optionally followed by arguments to that learning function, e.g. {''lda''} for linear discriminant analysis or {''logreg'',''lambda'',0.1} for logistic regression with ''lambda'',0.1 passed in as user parameters (see ml_train* functions for options).','type','expression'), ...
    arg({'trainer','TrainingFunction'},'@ml_train',[],'Training function. Receives a partition of the data (as produced by the specified partitioner), possibly some further arguments as specified in args, and returns a model (of any kind).','type','expression'), ...
    arg({'tester','PredictionFunction'},'@utl_default_predict',[],'Prediction function. Receives a partition of the data (as produced by the partitioner) and a model (as produced by the trainer), and returns a prediction for every index of the input data, in one of the output formats permitted by ml_predict.','type','expression'), ...
    ... % nested cross-validation parameters
    arg({'scheme','EvaluationScheme'},10,[],'Cross-validation scheme. Defines what parts of the data shall be used for training or testing in each fold. Supports many different formats, see documentation.','type','expression'), ...
    arg({'metric','EvaluationMetric'},'auto',{'auto','mcr','auc','mse','medse','smse','mae','medae','smae','smedae','max','sign','rms','bias','kld','nll','cond_entropy','cross_entropy','f_measure'},'Evaluation metric. Loss measure employed to measure the quality of predictions on a test partition given its known target values; applied both to results in each fold and results aggregated over all folds. Can be either a predefined metric supported by ml_calcloss, or a custom function handle which receives an array of target values in the first argument and an array of predictions in the second argument; each can be in any format that can be produced by ml_predict (but can be expected to be mutually consistent). Shall return a real number indicating the summary metric over all data, and optionally additional statistics in a second output struct.','typecheck',false), ...
    arg({'repeatable','RepeatableResults'},1,[],'Produce repeatable (vs. randomized) results. If nonzero, the value is taken as the random seed to use for the performed calculation. This way, different numbers (aside from 0) give different repeatable runs.'), ...
    ... % support for custom data representations
    arg({'target','TargetFunction'},'@utl_default_target',[],'Target extraction function. A function to derive the target variable from a partition of the data (as produced by the partitioner), for evaluation; the allowed format is anything that may be output by ml_predict.','type','expression'), ...
    arg({'partitioner','PartitioningFunction'},'@utl_default_partitioner',[],'Partitioning function. See documentation for the function contract.','type','expression'), ...
    ... % parallel computing settings
    arg({'engine_gs','ParallelEngineGS'},'global',{'global','local','BLS','Reference','ParallelComputingToolbox'}, 'Parallel engine for grid search. This can either be one of the supported parallel engines (BLS for BCILAB Scheduler, Reference for a local reference implementation, and ParallelComputingToolbox for a PCT-based implementation), or local to skip parallelization altogether, or global to select the currently globally selected setting (in the global tracking variable).'), ...
    arg({'engine_ncv','ParallelEngineNCV'},'global',{'global','local','BLS','Reference','ParallelComputingToolbox'}, 'Parallel engine for cross-validation. This can either be one of the supported parallel engines (BLS for BCILAB Scheduler, Reference for a local reference implementation, and ParallelComputingToolbox for a PCT-based implementation), or local to skip parallelization altogether, or global to select the currently globally selected setting (in the global tracking variable).'), ...
    arg({'pool','WorkerPool'},'global',[], 'Worker pool to use. This is typically a cell array, but can also be the string ''gobal'', which stands for the currently globally set up worker pool (see global tracking variable).','type','expression'), ...
    arg({'policy','ReschedulingPolicy'},'global',[], 'Rescheduling policy. This is the name of the rescheduling policy function that controls if and when tasks are being rescheduled. If set to global, the current global setting will be used.'), ...
    ... % misc arguments
    arg({'argform','ArgumentFormat'},'clauses',{'clauses','direct'},'Argument search range format. If set to clauses, search ranges for individual arguments can be given using the search() expression; if set to direct, a cell array is expected each of whose elements is a cell array that represents a particular parameter set to try for the training function.'), ...
    arg({'forcestats','ForceStatistics'},0,[],'Enforce a cross-validation to obtain stats. Even when no search is required.'), ...
    arg({'cache_fold_results','CacheFoldResults'},false,[],'Whether to cache the per-fold results. This is meant to be used when running very long-running computations on machines that crash frequently enough that partial results need to be saved. In this case, any previously computed results will be loaded from disk.'), ...
    arg({'only_cached_results','OnlyCachedResults'},false,[],'Load only results that are in the cache. This will not run any computations (aside from pre-checks, that can be disabled by setting NoPrechecks to true).'), ...
    arg({'no_prechecks','NoPrechecks'},false,[],'Skip pre-checks that access the data. This can save some time when it would take very long to load the data, especially when performing parallel computation.'), ...
    arg({'tolerate_exceptions','TolerateExceptions'},false,[],'Tolerate exceptions during training. If this happens, folds where the training function yielded errors will be skipped.'), ...
    arg({'collect_models','CollectModels'},false,[],'Collect models per fold. Note that this increases the amount of data returned.'));

% if there are parameters to be searched...
if opts.forcestats || is_needing_search(opts.argform,opts.args)

    % the nested cross-validation (utl_crossval) receives a specific subset of our arguments
    nestedcv_opts = hlp_struct2varargin(opts,'suppress',{'engine_gs','forcestats','args','argform'},'rewrite',{'engine_ncv','engine_cv'});

    % the objective function of our search is utl_crossval with the above-defined options; plus, it
    % has free arguments that are passed into the cross-validation as the 'args' parameter (i.e.,
    % the arguments that it passes into its 'trainer' function)
    objfun = @(varargin) utl_crossval(nestedcv_opts{:}, 'args',varargin);

    % the grid search receives the above-defined objective function and a subset of our arguments
    search_opts = {'argform',opts.argform, 'engine_gs',opts.engine_gs, 'policy',opts.policy, ...
        'pool',opts.pool, 'func',objfun};

    % now execute the grid search over every argument combination specified in opts.args
    [stats.bestidx,stats.inputs,stats.outputs] = utl_gridsearch(search_opts, opts.args{:});

else

    % otherwise produce dummy stats if neither necessary nor forced
    stats = struct('bestidx',1, 'inputs',{{opts.args}}, 'outputs',{{NaN,struct()}});

end

% finally, using the best args, send the input data through the training function & compute a model
bestargs = stats.inputs{stats.bestidx};
model = opts.trainer(opts.data,bestargs{:});

repo: REST


text:
function out = chebcoeffs(f, varargin)
%CHEBCOEFFS   Chebyshev polynomial coefficients of a DELTAFUN.
%   CHEBCOEFFS(F) returns the Chebyshev coefficients of F.FUNPART.

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

% Call CHEBCOEFFS() of the .FUNPART:
out = chebcoeffs(f.funPart, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%*****************************************************************
%% HKMcorr: corrector step for the HKM direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [dX,dy,dZ,resnrm,EinvRc] = HKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
            dX,dZ,coeff,L,X,Z);

    global matfct_options solve_ok

    printlevel = par.printlevel;
%%
    [rhs,EinvRc] = HKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ);
    m = length(rp); ncolU = size(coeff.mat12,2);
    rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
%%
    solve_ok = 1; resnrm = norm(rhs);
    if strcmp(matfct_options,'chol') | strcmp(matfct_options,'spchol') ...
       | strcmp(matfct_options,'ldl') | strcmp(matfct_options,'spldl')
       [xx,resnrm,solve_ok] = symqmr(coeff,rhs,L,[],[],printlevel);
       if (solve_ok<=0) & (printlevel)
          fprintf('\n  warning: symqmr fails: %3.1f.',solve_ok);
       end
    else
       [xx,resnrm,solve_ok] = mybicgstab(coeff,rhs,L,[],[],printlevel);
       if (solve_ok<=0) & (printlevel)
          fprintf('\n  warning: bicgstab fails: %3.1f.',solve_ok);
       end
    end
    if (printlevel>=3); fprintf('%2.0d ',length(resnrm)-1); end
%%
    [dX,dy,dZ] = HKMdirfun(blk,At,par,Rd,EinvRc,X,xx,m);
%%*****************************************************************


text:
function ui = mrViewAddOverlay(ui, mapNum);
%
% ui = mrViewAddOverlay([ui], [mapNum]);
%
% Initialize a new overlay for a mrViewer UI, by
% adding a new entry in the ui.overlays field and
% initializing a new overlay window (if the GUI is
% open).
%
%
% ras 07/05.
if ~exist('ui','var') | isempty(ui), ui = mrViewGet; end
if ishandle(ui), ui = get(ui, 'UserData'); end

% check that there are maps to use for the overlays:
% if not, load one:
if ~isfield(ui,'maps') | isempty(ui.maps)
    ui = mrViewLoad(ui,[],1);
end

% get index of new overlay
if ~isfield(ui,'overlays') | isempty(ui.overlays)
    o = 1;
else
    o = length(ui.overlays)+1;
end

if ~exist('mapNum','var') | isempty(mapNum)
    mapNum = min(o, length(ui.maps));
end

% preferred order of color maps according to overly #
prefOrder = {'hot' 'winter' 'greenred' 'cool' 'autumn' ...
             'jet' 'gray' 'red' 'green' 'blue'};

ui.overlays(o).hide = false; % logical indicator whether to hide overlay

ui.overlays(o).mapNum = mapNum;
for j = 1:3
    ui.overlays(o).thresholds(j).on = 0;
    ui.overlays(o).thresholds(j).mapNum = mapNum;
    ui.overlays(o).thresholds(j).min = 0;
    ui.overlays(o).thresholds(j).max = 1;
    ui.overlays(o).thresholds(j).autoScale = 0;
end
ui.overlays(o).cbar = cbarCreate(prefOrder{o});
ui.overlays(o).clim = 'auto';
ui.overlays(o).subVol = 1; % time / 4th-dimension index

% decide whether to dock based on 'dockFlag' preference variable
if ~ispref('VISTA', 'dockFlag'), setpref('VISTA', 'dockFlag', 0); end
dockFlag = getpref('VISTA', 'dockFlag');

if ishandle(ui.fig)
    ui = mrViewOverlayPanel(ui, dockFlag, o);
end

% Try to be predictive: initialize reasonable thresholds for some maps
if (strncmp(ui.maps(mapNum).dataUnits, '-log', 4) | ...
	strncmp(lower(ui.maps(mapNum).dataUnits), 't', 1)) & ...
        (ui.maps(mapNum).dataRange(1) < 2) & ...
        (ui.maps(mapNum).dataRange(2) > 2)
    % if the map appears to be a contrast map, initialize
	% the lower threshold to 2.0
    set(ui.overlays(o).threshCheck(1), 'Value', 1);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Value', 2);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Visible', 'on');
    mrvSliderSet(ui.overlays(o).threshMax(1), 'Visible', 'on');
    set(ui.overlays(o).threshMap(1), 'Visible', 'on');
    ui.overlays(o).thresholds(1).min = 2;
    ui.overlays(o).thresholds(1).on = 1;

elseif strncmp(ui.maps(mapNum).name, 'Coherence', 9)
	% Coherence map: initial lower thresh is 0.20
    set(ui.overlays(o).threshCheck(1), 'Value', 1);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Value', .2);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Visible', 'on');
    mrvSliderSet(ui.overlays(o).threshMax(1), 'Visible', 'on');
    set(ui.overlays(o).threshMap(1), 'Visible', 'on');
    ui.overlays(o).thresholds(1).min = .2;
    ui.overlays(o).thresholds(1).on = 1;

elseif strncmp(ui.maps(mapNum).name, 'meanMap', 7)
	% mean map: one 'click' up on the slider
	lo = ui.maps(mapNum).dataRange(1);
	hi = ui.maps(mapNum).dataRange(2);
	thresh = lo + .1 * (hi - lo);

    set(ui.overlays(o).threshCheck(1), 'Value', 1);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Value', thresh);
    mrvSliderSet(ui.overlays(o).threshMin(1), 'Visible', 'on');
    mrvSliderSet(ui.overlays(o).threshMax(1), 'Visible', 'on');
    set(ui.overlays(o).threshMap(1), 'Visible', 'on');
    ui.overlays(o).thresholds(1).min = thresh;
    ui.overlays(o).thresholds(1).on = 1;

end

ui = mrViewSetOverlay(ui, ui.panels.overlays(o));
% mrViewRefresh(ui);

return


text:
function H = cross(F, G)
%CROSS   Vector cross product.
%   CROSS(F, G) returns a SPHEREFUNV representing the 3D cross product of
%   the SPHEREFUNV objects F and G.
%
%   See also SPHEREFUNV/DOT.

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

% [TODO]: Implement the following option
%   CROSS(F, G, 'n') returns a SPHEREFUN representing the normal component
%   of the cross product of the SPHEREFUNV objects F and G.
%   For two vector fields tangent to the sphere, the cross product is a
%   a vector that points in the normal direction.  CROSS returns a
%   SPHEREFUN object representing the component in the normal (radial)
%   direction.  Mathematically, this is N dot (F x G), where N is the
%   normal to the sphere.

% Empty check:
if ( isempty(F) || isempty(G) )
    H = spherefunv;
    return
end

% Get the components:
Fc = F.components;
Gc = G.components;

% Do cross:
H = [ Fc{2} .* Gc{3} - Fc{3} .* Gc{2} ; ...
      Fc{3} .* Gc{1} - Fc{1} .* Gc{3} ; ...
      Fc{1} .* Gc{2} - Fc{2} .* Gc{1} ];

end


text:
function [Data]=ReconstructImage(amp,phi,nb)
if (max(size(size(amp)))~=3)
    error('amp and phi must be three dimensional [nf,lat,lon]')
end

[~,ny,nx]=size(amp);

Data=zeros(nb,ny,nx);

h= waitbar(0,'Total Progress:');
WBarOuterPosition=get(h,'OuterPosition');
WBarOuterPosition(2)=WBarOuterPosition(2)-WBarOuterPosition(4);
h2=waitbar(0,'Calculating, please wait ...');
set(h2,'OuterPosition',WBarOuterPosition);

for Sample=1:nx
    waitbar(Sample/nx,h);
    for Line=1:ny
        waitbar(Line/ny,h2,['Line:' num2str(Line) ', Sample:' num2str(Sample)]);
        amp_Pixel=squeeze(amp(:,Line,Sample));
        phi_Pixel=squeeze(phi(:,Line,Sample));
        Data(:,Line,Sample)=ReconHANTSData(amp_Pixel,phi_Pixel,nb);
    end
end
close(h);
close(h2);
end


text:
function [data,net] = prepare_data(P,Y,net,P0,Y0)
% 	Prepare Input Data for the use for NN Training and check for errors
%
% 	Args:
% 		P:		neural network Inputs
% 		Y: 		neural network Targets
% 		P0:		previous input data
% 		Y0:		previous output data
% 		net: 	neural network
% 	Returns:
% 		data:	struct containing data for training or calculating putput

if ~exist('P0', 'var')
    P0 = [];
end
if ~exist('Y0', 'var')
    Y0 = [];
end

%  Ceck if input and output data match structure of NN
if size(P,1) ~= net.nn(1)
    error('Dimension of Input Data P does not match number of inputs of the NN')
end
if size(Y,1) ~= net.nn(end)
    error('Dimension of Output Data Y does not match number of outputs of the NN')
end
if size(Y,2) ~= size(P,2)
    error('Input P and output Y must have same number of datapoints Q')
end

%check if previous data is given
if (~isempty(P0)) && (~isempty(Y0))
    %  Ceck if input and output of prevoius data match structure of NN
    if size(P0,1) ~= net.nn(1)
        error('Dimension of previous Input Data P0 does not match number of inputs of the NN')
    end
    if size(Y0,1) ~= net.nn(end)
        error('Dimension of previous Output Data Y0 does not match number of outputs of the NN')
    end
    if size(Y,2) ~= size(P,2)
    error('Previous Input and output data P0 and Y0 must have same number of datapoints Q0')
    end

    q0 = size(P0,2); %number of prevoius Datapoints given
    a=cell(q0,net.M);  %initialise layer outputs
    for i=1:q0
        for j=1:net.M-1
            a{i,j}=zeros(net.nn(j+1),1); %layer ouputs of hidden layers are unknown -> set to zero
        end
        a{i,net.M}=Y0(:,i)./net.normY; %set layer ouputs of output layer
    end
    %add previous inputs and outputs to inpu/output matrices
    P_ = [P0,P];
    Y_ = [Y0,Y];
else
	%add previous inputs and outputs to inpu/output matrices
    P_ = P;
    Y_ = Y;
    q0=0;
    a={};
end

% normalize
P_norm = P_;
Y_norm = Y_;
if isfield(net, 'normP')==0
    normP = ones(size(P_,1),1);
    for p = 1:size(P_,1)
        normP(p) = max(max(abs(P_(p,:))),1.0);
        P_norm(p,:) = P_(p,:)./normP(p);
    end
    normY = ones(size(Y_,1),1);
    for y = 1:size(Y_,1)
        normY(y) = max(max(abs(Y_(y,:))),1.0);
        Y_norm(y,:) = Y_(y,:)./normY(y);
    end
    net.normP = normP;
    net.normY = normY;
else
    for p = 1:size(P_,1)
        P_norm(p,:) = P_(p,:)./net.normP(p);
    end
    for y = 1:size(Y_,1)
        Y_norm(y,:) = Y_(y,:)./net.normY(y);
    end
end

%Create data dict
data = {};
data.P = P_norm;
data.Y = Y_norm;
data.a = a;
data.q0 = q0;

end

author: yabata

repo: pyrenn

save_path: github-repos/MATLAB/yabata-pyrenn


text:
function box = mergeBoxes(box1, box2)
%MERGEBOXES Merge two boxes, by computing their greatest extent
%
%   BOX = mergeBoxes(BOX1, BOX2);
%
%   Example
%   box1 = [5 20 5 30];
%   box2 = [0 15 0 15];
%   mergeBoxes(box1, box2)
%   ans =
%       0 20 0 30
%
%
%   See also
%   boxes2d, drawBox, intersectBoxes
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% unify sizes of data
if size(box1,1) == 1
    box1 = repmat(box1, size(box2,1), 1);
elseif size(box2, 1) == 1
    box2 = repmat(box2, size(box1,1), 1);
elseif size(box1,1) ~= size(box2,1)
    error('Bad size for inputs');
end

% compute extreme coords
mini = min(box1(:,[1 3]), box2(:,[1 3]));
maxi = max(box1(:,[2 4]), box2(:,[2 4]));

% concatenate result into a new box structure
box = [mini(:,1) maxi(:,1) mini(:,2) maxi(:,2)];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function expanded_array = expand_complex(complex_array)

% Note that if "complex_array" is a column vector, then its real and imaginary
% parts are the row vectors of "expanded_array", because the fastest-varying
% index (i.e., index in the column direction) is assigned for alternating the
% real and imaginary parts.
chkarg(istypeof(complex_array, 'complex'), '"complex_array" should be array with complex elements.');
nD = ndims(complex_array);
expanded_array = cat(nD+1, real(complex_array), imag(complex_array));
expanded_array = permute(expanded_array, [nD+1, 1:nD]);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function output = brtTest( input, brtModel, varargin )

    if isempty(varargin)
        len = length(brtModel)-1;
    else
        len = varargin{1};
    end

    nu = brtModel{end};
    output = brtModel{1};

    for i=2:len
        output = output + nu * regressionTreeTest( input, brtModel{i} );
    end
end


text:
function r8col_unique_count_test ( )

%*****************************************************************************80
%
%% R8COL_UNIQUE_COUNT_TEST tests R8COL_UNIQUE_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2010
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 15;

  a = [ ...
    2.0,  6.0, 10.0; ...
    4.0,  8.0, 12.0; ...
    1.0,  5.0,  9.0; ...
    3.0,  7.0, 11.0; ...
    2.0,  6.0,  0.0; ...
    3.0,  4.0, 18.0; ...
    0.0,  0.0,  0.0; ...
    0.0,  6.0, 10.0; ...
    2.0,  6.0, 10.0; ...
    3.0,  7.0, 11.0; ...
    2.0,  0.0, 10.0; ...
    2.0,  6.0, 10.0; ...
    1.0,  5.0,  9.0; ...
    1.0,  5.0,  9.1; ...
    1.0,  5.1,  9.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8COL_UNIQUE_COUNT_TEST\n' );
  fprintf ( 1, '  For an R8COL;\n' );
  fprintf ( 1, '  R8COL_UNIQUE_COUNT counts unique columns.\n' );
  fprintf ( 1, '\n' );

  r8mat_transpose_print ( m, n, a, '  The R8COL (transposed):' );

  unique_num = r8col_unique_count ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of unique columns is %d\n', unique_num );

  return
end


text:
function [ V ] = R2V( R )
%R2V converts 3x3 rotation matrix into a 1x3 angle-axis vector
%   Inputs -
%   R - a standard 3x3 transformation matrix
%
%   Outputs -
%   V - 1x3 vector of form [rx,ry,rz] where rx,ry,rz is an angle-axis
%   representation of the angle where the unit vector representing the axis
%   has been multipled by the angle of rotation about it

validateattributes(R, {'numeric'},{'size',[3,3]});

R = vrrotmat2vec(double(R));
V = R(1:3)*R(4);
V = V';
end

author: ZacharyTaylor

repo: Camera-to-Arm-Calibration