text:
function test_withinScansMotionComp
%Test test_withinScansMotionComp as well as the creation of a new data type
%
%   test_withinScansMotionComp()
%
%
% Tests: intiHiddenInplane, test_withinScansMotionComp, dtGet, motionCompSelScan
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_withinScansMotionComp()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2014
%
%

% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','prfInplane');

% This is the validation file
stored = mrtGetValididationData('withinScansMotionComp');

%
% dtNum = 2;
% stored.name = dtGet(dataTYPES(dtNum),'Name');
% stored.annotation = dtGet(dataTYPES(dtNum),'Annotation', 1);
% stored.nFrames = dtGet(dataTYPES(dtNum),'nFrames', 1);
% stored.framePeriod = dtGet(dataTYPES(dtNum),'Frame Period', 1);
% stored.numSlices = length(dtGet(dataTYPES(dtNum),'slices', 1));
% stored.numScans = dtGet(dataTYPES(dtNum),'N Scans');
% stored.motionEstimates = motionEstimates;
% save(vFile, '-struct',  'stored')

% Relative tolerance of assertAlmostEqual
relTol = 1e-10;

% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% Get data structure:
vw = initHiddenInplane;
mrGlobals;

typeName        = 'WitinScansMotionComp';
baseScan        = 1;
baseFrame       = 1;
scansToCorrect  = 1;
nSmooth         = 1;

% compute the within scan motion, and return the view structure
vw = motionCompSelScan(vw, typeName, scansToCorrect, baseFrame, nSmooth, baseScan);

% close the graph window created by motionCompSelScan
close(gcf)

%Now that we have created the necessary data, we can do the actual tests.
dtNum  = length(dataTYPES);
vw     = viewSet(vw, 'current dt', dtNum);
nScans = viewGet(vw, 'num scans');

% get the within scan motion
motionEstimates = dtGet(dataTYPES(dtNum), 'within scan motion', scansToCorrect);

% Now check that the motion estimates agree with stored estimates
assertVectorsAlmostEqual(motionEstimates, stored.motionEstimates, 'relative', relTol);

% Check that the new dataTYPE got the proper name
assertEqual(typeName, dtGet(dataTYPES(dtNum),'Name'));

% Check that the new dataTYPE got the proper annoation
assertEqual(stored.annotation, dtGet(dataTYPES(dtNum),'Annotation', nScans));

% Check that the new dataTYPE got the proper numbr of frames
assertEqual(stored.nFrames, dtGet(dataTYPES(dtNum),'nFrames', nScans));

% Check that the new dataTYPE got the proper TR
assertEqual(stored.framePeriod, dtGet(dataTYPES(dtNum),'Frame Period', nScans));

% Check that the new dataTYPE got the proper number of slices
assertEqual(stored.numSlices, length(dtGet(dataTYPES(dtNum),'slices', nScans)));

% Check that the new dataTYPE got the proper number of scans
assertEqual(stored.numScans, dtGet(dataTYPES(dtNum),'N Scans'));

% clean up vistadata repository because this test script wrote new data
% test_CleanUpSVN
mrvCleanWorkspace;

cd(curDir);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pascalScore = computePascalScore(bb1,bb2)
%compute the Pascal score of the bb1, bb2 (intersection/union)

intersectionArea = computeIntersectionArea(bb1,bb2);
pascalScore = intersectionArea/(computeArea(bb1)+computeArea(bb2)-intersectionArea);
return

author: Cloud-CV

repo: object-proposals


text:
function [dat] = bz_scrubTracking(dat)
% USAGE
%  [dat] = bz_scrubTracking(dat)
%
% INPUTS
%   dat - struct created by calling dat = importdata(file.csv) on a csv
%         file exported from a .tak file
%
% OUTPUTS
%   dat - the same struct 1) outliers removed
%                         2) median filtering
%                         3) interpolation over missing frames
%
%
%
% this functions takes the unedited *.csv exported from Motive, smoothes and
% resorts columns to match the position tracking formatting standards
%
% David Tingley, 2017
%
%TODO
% check that units are in meters...

% find extreme outliers outside of tracking volume
dat.data((abs(dat.data(:,7))>mean(dat.data(:,7))+std(dat.data(:,7))*5),7)=nan;
dat.data((abs(dat.data(:,8))>mean(dat.data(:,8))+std(dat.data(:,8))*5),8)=nan;
dat.data((abs(dat.data(:,9))>mean(dat.data(:,9))+std(dat.data(:,9))*5),9)=nan;

% find poorly tracked frames using the error per marker column
dat.data(find(dat.data(:,10)>mean(dat.data(:,10))+std(dat.data(:,10))*10),:)=nan;

% find dropped frames or large translational shifts and call them nan's here
f = find(diff(dat.data(:,8))==0 | abs(diff(dat.data(:,8)))> mean(std(abs(diff(dat.data(:,8)))))+ std(abs(diff(dat.data(:,8)))) * 10);
dat.data(f+1,8)=nan;
f = find(diff(dat.data(:,9))==0 | abs(diff(dat.data(:,9)))> mean(std(abs(diff(dat.data(:,9)))))+std(abs(diff(dat.data(:,9)))) * 10);
dat.data(f+1,9)=nan;
f = find(diff(dat.data(:,7))==0 | abs(diff(dat.data(:,7)))> mean(std(abs(diff(dat.data(:,7)))))+std(abs(diff(dat.data(:,7)))) * 10);
dat.data(f+1,7)=nan;

% median filtering
dat.data(:,7) = medfilt1(dat.data(:,7),12);
dat.data(:,8) = medfilt1(dat.data(:,8),12);
dat.data(:,9) = medfilt1(dat.data(:,9),12);

% pattern matching interpolation
dat.data(:,7) = fillmissing(dat.data(:,7),'pchip');
dat.data(:,8) = fillmissing(dat.data(:,8),'pchip');
dat.data(:,9) = fillmissing(dat.data(:,9),'pchip');

%  now change column order to new standards
% columns = [7 9 8 3 5 4 6 10 1 2]; % optitrack is Y-oriented instead of Z
% dat.data = dat.data(:,columns);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function CHECK_FILE_EXIST(filename)

if exist(filename, 'file') == 0
  error('%s does not exist', filename);
end

end

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
function st_io_test02 ( )

%*****************************************************************************80
%
%% ST_IO_TEST02 tests ST_HEADER_READ, ST_DATA_READ.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2014
%
%  Author:
%
%    John Burkardt
%
  input_filename = 'kershaw.st';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_IO_TEST02\n' );
  fprintf ( 1, '  ST_HEADER_READ reads the header from an ST file.\n' );
  fprintf ( 1, '  ST_DATA_READ reads the data from an ST file.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data from "%s".\n', input_filename );

  [ i_min, i_max, j_min, j_max, m, n, nst ] = st_header_read ( input_filename );

  st_header_print ( i_min, i_max, j_min, j_max, m, n, nst );

  [ ist, jst, ast ] = st_data_read ( input_filename, m, n, nst );

  st_print ( m, n, nst, ist, jst, ast, '  Sparse Triplet (ST) data from file' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = mrdivide(A, B)
%/   Right matrix divide for a BNDFUN.
%   A/B divides the BNDFUN A by a scalar B. More generally, it gives the
%   least-squares solution (with respect to the continuous L^2 norm) to X*B = A
%   when either A or B is a BNDFUN.  Note that in the latter case, formally it
%   is X.' that is returned, as BNDFUN objects are always columns.
%
% See also QR, RDIVIDE, MLDIVIDE.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Case B is a matrix: X*B = A  ==> X = A/B = (Q*R)/B = Q*(R/B)
%
% Case A is a matrix: X*B = X*(Q*R) = A ==> X = (A/R)*Q' ==> X' = Q*(A/R)'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( (size(B, 2) ~= size(A, 2)) && ~(isa(B, 'double') && isscalar(B)) )
    error('CHEBFUN:BNDFUN:mrdivide:size', 'Matrix dimensions must agree.');

elseif ( isa(B, 'double') )  % BNDFUN / double

    if ( isscalar(B) )
        % Scalar case is easy:
        X = A;                              % Copy A to X
        X.onefun = X.onefun/B;              % mrdivide of the onefun
    else
        % Call MRDIVIDE at the ONEFUN level:
        X = A;
        X.onefun = (A.onefun/B);

%         % Alternatively, we could call QR() at the BNDFUN level
%         % For matrix case, we do least squares via QR:
%         [Q, R] = qr(A, 0);
%         X = Q*(R/B);
    end

elseif ( isa(A, 'double') )  % double / BNDFUN

    % Call MRDIVIDE at the ONEFUN level:
    X = B;
    X.onefun = (A/B.onefun);
    X = X/(.5*diff(B.domain));

%     % Alternatively we could call QR() at the BNDFUN level:
%     % Do least squares via QR:
%     [Q, R] = qr(B, 0);
%     % Return the transpose for the output.
%     X = Q*(A/R).';

elseif ( isa(B, 'bndfun') && isa(A, 'bndfun') )
    error('CHEBFUN:BNDFUN:mrdivide:bndfunDivBndfun', ...
        'Use ./ to divide BNDFUN by a BNDFUN.');

else
    error('CHEBFUN:BNDFUN:mrdivide:badArg', '%s/%s is not well-defined.', ...
        class(A), class(B));

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [A,Pt] = mci_lds_dfdx (x,u,P,M)
% Jacobian for linear system, dx/dt=Ax, with constrained connectivity
% FORMAT [A,Pt] = mci_lds_dfdx (x,u,P,M)
%
% x     State vector
% u     input
% P     parameters (vectorised)
% M     model structure
%
% A     f=Ax
% Pt    Parameters (transformed from latent pars)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_lds_dfdx.m 6548 2015-09-11 12:39:47Z will $

[Pt,a,b] = mci_lds_lat2par (P,M);
A=diag(a);

Nb=length(b);
for k=1:Nb,
    i=M.Aconn(k,1);
    j=M.Aconn(k,2);
    A(i,j)=b(k);
end

author: spm


text:
function SIM_GEN_f_FuncGen(plant)
    fileID = fopen('./funcgen/SIM_GEN_f.m','w');

    fprintf(fileID, 'function f = SIM_GEN_f(u,x,p)\n');
    fprintf(fileID, '   f = fSim_Wrapper(u,x,p);\n');
    fprintf(fileID, 'end');

    fclose(fileID);

end

author: deng-haoyang

repo: ParNMPC


text:
%% Grain Tutorial
%
%%
% The following script is a quick guide through the grain reconstruction
% capabilities of MTEX. It uses the same data set as in the corresponding
% publication
% <https://www.researchgate.net/publication/51806709_Grain_detection_from_2d_and_3d_EBSD_data-Specification_of_the_MTEX_algorithm
% Grain detection from 2d and 3d EBSD data>. Data courtesy of Daniel Rutte
% and Bret Hacker, Stanford.

mtexdata mylonite

% plot a phase map
plot(ebsd)

%%
% The phase map shows a multi-phase rock specimen with Andesina, Quartz,
% Biotite and Orthoclase. Lets restrict it to a smaller region of interest.
% The rectangle is defined by [xmin, ymin, xmax-xmin, ymax-ymin].

region = [19000 1500 4000 1500];
% overlay the selected region on the phase map
rectangle('position',region,'edgecolor','k','linewidth',2)

%%
% Now copy the EBSD data within the selected rectangle to a new variable

ebsd_region = ebsd(inpolygon(ebsd,region))

%% Grain Reconstruction
% Next we reconstruct the grains and grain boundaries in the region of
% interest, using a 15 degree orientation change threshold.

grains = calcGrains(ebsd_region,'angle',15*degree)

% plot a phase map of the region of interest
plot(ebsd_region)

% overlay the grain boundaries
hold on
plot(grains.boundary,'color','k','linewidth',1.5)
hold off

%%
% We may also visualize the different quarz orientations together with the
% grain boundaries.

% plot a phase map of three of the phases based on the grains data
plot(grains({'Andesina','Biotite','Orthoclase'}),'FaceAlpha',0.4)

hold on
% add the quarz orientations as ipf map based on EBSD data
plot(ebsd_region('Quartz'),ebsd_region('Quartz').orientations)

% plot grain boundaries so that those in the Quartz are shown
plot(grains.boundary,'color','black');
legend off
hold off

%%
% For the map created, most of the phases are coloured based on where they
% exist, while only the Quartz phase is colored according to the
% orientation. The quartz orientations are colured using the following ipf
% color key

close all
ipfKey = ipfColorKey(ebsd_region('Quartz'));
plot(ipfKey)

%%
% Alternatively, we may colorize each quarz grain according to its mean
% orientation.  Again, the other phases are colured based on where they
% exist.

plot(grains({'Andesina','Biotite','Orthoclase'}),'FaceAlpha',0.4)
hold on
plot(grains('Quartz'),grains('Quartz').meanOrientation)
legend off

%% Highlight specific boundaries
% We can create a phase map with certain grain boundaries highlighted.  In
% this case, we highlight where adjacent grains of Andesina and Orthoclase
% have a misorientation with rotational axis close to the c-axis.

close all
% copy all boundaries between Andesina Orthoclase to a new variable
AOboundary = grains.boundary('Andesina','Orthoclase');
% copy the misorientation angle of this boundary in radians to a new variable.
angle = AOboundary.misorientation.angle;

plot(grains,'FaceAlpha',0.4)
hold on
% highlight boundaries where the angle between the Andesina and Orthoclase phase is over 160 degrees
plot(AOboundary(angle>160*degree),'linewidth',2,'linecolor','red')
hold off

%%
% We can also represent the angular misorientation data between these two
% phases as a histogram.

figure;histogram(angle./degree)
xlabel('angle in degrees of boundary segment')
ylabel('count of boundary segments')
title('angular relationships between Andesina and Orthoclase')

author: mtex-toolbox

repo: mtex

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


text:
function fns_elec_write(pnt, vsize, dimpos, elecfile)
%
% This function write the electrodes locations on disk and other
% information
%   1. pnt is the location of the electrodes [NX3]
%   3. voxel_sizes is the dimension of a voxel in the cartesian coordinates in mm [3X1]
%   4. nodes_sizes is the dimension of the dipoles grid [NX NY NZ]

% $Copyright (C) 2010 by Hung Dang$

hdf5write(elecfile, '/electrodes/locations', pnt);

hdf5write(elecfile, '/electrodes/gridlocs', int32(pnt), ...
          'WriteMode', 'append');       % Assume the electrodes
                                        % locations are mm and the
                                        % voxel sizes is 1mm x 1mm
                                        % x 1mm.

hdf5write(elecfile, '/electrodes/voxel_sizes', vsize, ...
          'WriteMode','append');

hdf5write(elecfile, '/electrodes/node_sizes', int32(dimpos), ...
          'WriteMode', 'append');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pt_start(mCatalog, hFigure, bMap, rContainer, sName)
    % Starts the probabilistic forecast testing.
    %
    % pt_start(mCatalog, hFigure, bMap, rContainer)
    %
    % The function invokes the parameter dialog, calls the
    %   calculation function and invokes the dialog for displaying the results
    %
    % Input parameters:
    %   mCatalog      Earthquake catalog to use for the testing
    %   hFigure       Handle of figure where the user should select the grid (i.e. the seismicity map)
    %   bMap          Map/cross-section switch. If the testing is carried out on a map set bMap = true,
    %                 on a cross-section set bMap = false
    %   rContainer    Just a container (structure proposed) to store any variables that should be available in the
    %                 probabilistic forecast testing code
    %
    % Danijel Schorlemmer
    % June 4, 2002

    report_this_filefun();

    % Launch GUI
    hMenuFig = pt_options(bMap);

    % Set up parameter struct and store input parameters
    params.mCatalog = mCatalog;
    params.sCatalogName = sName;
    params.bMap = bMap;
    params.rContainer = rContainer;

    % Analyze Output
    if ~ishandle(hMenuFig)
        answer = 0;
    else
        handles = guidata(hMenuFig);
        answer = handles.answer;
        % OK pressed
        if answer == 1
            % Get parameter values from the gui
            % ---------------------------------
            % Hypothesis
            if get(handles.radHNumber, 'Value') == 1
                rOptionsH.nGriddingMode = 0;
                rOptionsH.nNumberEvents = str2double(get(handles.txtHNumber, 'String'));
                rOptionsH.fRadius = 0;
                rOptionsH.fSizeRectX = 0;
                rOptionsH.fSizeRectY = 0;
                rOptionsH.nMinimumNumber = 0;
                rOptionsH.fMaximumRadius = str2double(get(handles.txtHMaximumRadius, 'String'));
            elseif get(handles.radHRadius, 'Value') == 1
                rOptionsH.nGriddingMode = 1;
                rOptionsH.nNumberEvents = 0;
                rOptionsH.fRadius = str2double(get(handles.txtHRadius, 'String'));
                rOptionsH.fSizeRectX = 0;
                rOptionsH.fSizeRectY = 0;
                rOptionsH.nMinimumNumber = str2double(get(handles.txtHMinimumNumber, 'String'));
                rOptionsH.fMaximumRadius = 0;
            else % rectangle
                rOptionsH.nGriddingMode = 2;
                rOptionsH.nNumberEvents = 0;
                rOptionsH.fRadius = 0;
                rOptionsH.fSizeRectX = str2double(get(handles.txtHSizeRectX, 'String'));
                rOptionsH.fSizeRectY = str2double(get(handles.txtHSizeRectY, 'String'));
                rOptionsH.nMinimumNumber = str2double(get(handles.txtHMinimumNumber, 'String'));
                rOptionsH.fMaximumRadius = 0;
            end
            rOptionsH.nCalculateMC = get(handles.cboHCalculateMC, 'Value');
            if get(handles.radHSpatialB, 'Value') == 1
                rOptionsH.nCalcMode = 0;  % Spatial b-values
            else
                rOptionsH.nCalcMode = 1;  % Overall b-values
            end
            rOptionsH.bCalcBothPeriods = get(handles.chkHCalcBothPeriods, 'Value');

            % Null hypothesis
            if get(handles.radNNumber, 'Value') == 1
                rOptionsN.nGriddingMode = 0;  % Constant number
                rOptionsN.nNumberEvents = str2double(get(handles.txtNNumber, 'String'));
                rOptionsN.fRadius = 0;
                rOptionsN.fSizeRectX = 0;
                rOptionsN.fSizeRectY = 0;
                rOptionsN.nMinimumNumber = 0;
                rOptionsN.fMaximumRadius = str2double(get(handles.txtNMaximumRadius, 'String'));
            elseif get(handles.radNRadius, 'Value') == 1
                rOptionsN.nGriddingMode = 1;  % Constant Radius
                rOptionsN.nNumberEvents = 0;
                rOptionsN.fRadius = str2double(get(handles.txtNRadius, 'String'));
                rOptionsN.fSizeRectX = 0;
                rOptionsN.fSizeRectY = 0;
                rOptionsN.nMinimumNumber = str2double(get(handles.txtNMinimumNumber, 'String'));
                rOptionsN.fMaximumRadius = 0;
            else % rectangle
                rOptionsN.nGriddingMode = 2;  % Rectangles
                rOptionsN.nNumberEvents = 0;
                rOptionsN.fRadius = 0;
                rOptionsN.fSizeRectX = str2double(get(handles.txtNSizeRectX, 'String'));
                rOptionsN.fSizeRectY = str2double(get(handles.txtNSizeRectY, 'String'));
                rOptionsN.nMinimumNumber = str2double(get(handles.txtNMinimumNumber, 'String'));
                rOptionsN.fMaximumRadius = 0;
            end
            rOptionsN.nCalculateMC = get(handles.cboNCalculateMC, 'Value');
            if get(handles.radNSpatialB, 'Value') == 1
                rOptionsN.nCalcMode = 0;  % Spatial b-values
            else
                rOptionsN.nCalcMode = 1;  % Overall b-values
            end
            rOptionsN.bCalcBothPeriods = get(handles.chkNCalcBothPeriods, 'Value');

            % Testing
            if get(handles.radTNumber, 'Value') == 1
                rOptionsT.nGriddingMode = 0;  % Constant number
                rOptionsT.nNumberEvents = str2double(get(handles.txtTNumber, 'String'));
                rOptionsT.fRadius = 0;
                rOptionsT.fSizeRectX = 0;
                rOptionsT.fSizeRectY = 0;
                rOptionsT.nMinimumNumber = 0;
                rOptionsT.fMaximumRadius = str2double(get(handles.txtTMaximumRadius, 'String'));
            elseif get(handles.radTRadius, 'Value') == 1
                rOptionsT.nGriddingMode = 1;  % Constant Radius
                rOptionsT.nNumberEvents = 0;
                rOptionsT.fRadius = str2double(get(handles.txtTRadius, 'String'));
                rOptionsT.fSizeRectX = 0;
                rOptionsT.fSizeRectY = 0;
                rOptionsT.nMinimumNumber = str2double(get(handles.txtTMinimumNumber, 'String'));
                rOptionsT.fMaximumRadius = 0;
            else % rectangle
                rOptionsT.nGriddingMode = 2;  % Rectangles
                rOptionsT.nNumberEvents = 0;
                rOptionsT.fRadius = 0;
                rOptionsT.fSizeRectX = str2double(get(handles.txtTSizeRectX, 'String'));
                rOptionsT.fSizeRectY = str2double(get(handles.txtTSizeRectY, 'String'));
                rOptionsT.nMinimumNumber = str2double(get(handles.txtTMinimumNumber, 'String'));
                rOptionsT.fMaximumRadius = 0;
            end
            rOptionsT.nCalculateMC = [];      % Fill the unused fields
            rOptionsT.nCalcMode = [];         % Fill the unused fields
            rOptionsT.bCalcBothPeriods = false;   % Fill the unused fields
            bSaveParameter = get(handles.chkSaveParameter, 'Value');
            if bSaveParameter
                sSaveParameter = get(handles.lblSaveParameter, 'String');
            end

            % Set up the options array
            params.rOptions = [rOptionsH; rOptionsN; rOptionsT];

            % Options section
            params.bGridEntireArea = get(handles.chkGridEntireArea, 'Value');
            params.fSpacingX = str2double(get(handles.txtSpacingX, 'String'));
            params.fSpacingY = str2double(get(handles.txtSpacingY, 'String'));
            params.fSplitTime = str2double(get(handles.txtSplitTime, 'String'));
            params.bForecastPeriod = get(handles.chkForecastPeriod, 'Value');
            if params.bForecastPeriod
                params.fForecastPeriod = str2double(get(handles.txtForecastPeriod, 'String'));
            else
                params.fForecastPeriod = 0;
            end
            params.bLearningPeriod = get(handles.chkLearningPeriod, 'Value');
            if params.bLearningPeriod
                params.fLearningPeriod = str2double(get(handles.txtLearningPeriod, 'String'));
            else
                params.fLearningPeriod = 0;
            end
            params.bRandomNode = get(handles.chkRandomNode, 'Value');
            if params.bRandomNode
                params.nNumberCalculationNode = str2double(get(handles.txtNumberCalculationNode, 'String'));
            end
            params.bMinMagMc = get(handles.chkMinMag, 'Value');
            if ~params.bMinMagMc
                params.fMinMag = str2double(get(handles.txtMinMag, 'String'));
            else
                params.fMinMag = [];
            end
            params.fMaxMag = str2double(get(handles.txtMaxMag, 'String'));

            % Additional stuff
            params.nTestMethod = get(handles.cboTestMethod, 'Value');
            params.bForceRandomCalculation = false;
            params.sComment = 'First test';

            % Dave Jackson
            params.bSaveRates = true;
            params.sSaveRatesAsciiFilenameH = '/home/danijel/pro/parkfield/Runs/ratesH.txt';
            params.sSaveRatesAsciiFilenameN = '/home/danijel/pro/parkfield/Runs/ratesN.txt';
            params.sSaveRatesMatFilenameH = '/home/danijel/pro/parkfield/Runs/ratesH.mat';
            params.sSaveRatesMatFilenameN = '/home/danijel/pro/parkfield/Runs/ratesN.mat';
            %     delete(params.sSaveRatesAsciiFilenameH);
            %     delete(params.sSaveRatesAsciiFilenameN);
            %     delete(params.sSaveRatesMatFilenameH);
            %     delete(params.sSaveRatesMatFilenameN);
            params.vRatesH = [];
            params.vRatesN = [];

            % Remove figure from memory
            delete(hMenuFig);

            % Select grid
            [params.mPolygon, params.vX, params.vY, params.vUsedNodes] = ex_selectgrid(hFigure, params.fSpacingX, params.fSpacingY, params.bGridEntireArea);

            % Validate polygonsize
            if length(params.vX) < 2  ||  length(params.vY) < 2
                errordlg('Selection is too small. Please select a larger polygon.');
                return;
            end

            % General calculations
            % --------------------

            % Split the catalog
            [params.mLearningCatalog, params.mObservedCatalog, ~, ~, params.fLearningPeriodUsed, params.fObservedPeriodUsed] ...
                = ex_SplitCatalog(params.mCatalog, params.fSplitTime, ...
                params.bLearningPeriod, params.fLearningPeriod, params.bForecastPeriod, params.fForecastPeriod);
            % Create Indices to catalogs
            for nCnt_ = 1:3
                % Create Indices to learning catalogs
                params.rOptions(nCnt_).caLearningNodeIndices = ex_CreateIndexCatalog(params.mLearningCatalog, params.mPolygon, params.bMap, ...
                    params.rOptions(nCnt_).nGriddingMode, params.rOptions(nCnt_).nNumberEvents, params.rOptions(nCnt_).fRadius, ...
                    params.rOptions(nCnt_).fSizeRectX, params.rOptions(nCnt_).fSizeRectY);
                % Create Indices to learning catalogs
                params.rOptions(nCnt_).caObservedNodeIndices = ex_CreateIndexCatalog(params.mObservedCatalog, params.mPolygon, params.bMap, ...
                    params.rOptions(nCnt_).nGriddingMode, params.rOptions(nCnt_).nNumberEvents, params.rOptions(nCnt_).fRadius, ...
                    params.rOptions(nCnt_).fSizeRectX, params.rOptions(nCnt_).fSizeRectY);
            end

            if bSaveParameter
                % Save parameters
                save(sSaveParameter, 'params');
            else
                %       % Determine the overall magnitude of completeness and overall b-value
                %       params.fMcOverall = calc_Mc(params.mCatalog, params.nCalculateMC);
                %       vSel = params.mCatalog.Magnitude >= params.fMcOverall;
                %       [params.fBValueOverall params.fStdDevOverall] = calc_bmemag(params.mCatalog.Magnitude(vSel), 0.1);

                %       % Determine the overall standard deviation of the b-value for the bayesian approach (if used)
                %       if params.nTestMethod == 3
                %         disp(['Calculating overall mean standard deviation']);
                %         params.fStdDevOverall = kj_CalcOverallStdDev(params);
                %         disp(['Standard deviation calculated']);
                %       end

                % Perform the calculation
                [params] = pt_calc(params);

                if params.bSaveRates
                    vRatesH = params.vRatesH;
                    vRatesN = params.vRatesN;
                    save(params.sSaveRatesAsciiFilenameH, 'vRatesH', '-ascii');
                    save(params.sSaveRatesAsciiFilenameN, 'vRatesN', '-ascii');
                    save(params.sSaveRatesMatFilenameH, 'vRatesH');
                    save(params.sSaveRatesMatFilenameN, 'vRatesN');
                end

                % Display the results
                gui_result(params);
            end
        else
            delete(hMenuFig);
        end
    end
end

author: CelsoReyes


text:
fprintf('Compiling mex functions\n');
mex -largeArrayDims maskmult.c
mex -largeArrayDims setsparseentries.c

author: andrewssobral

repo: lrslibrary


text:
% Test file read
% File:  c:\M-files\shortcut_updates\readtxt.m
% 2/11/07
%
clc;clear;
[freq,mag,ph]=textread('c:\Spiceapps\datfiles\Elliptf7.txt','%f %f %f');
%
h=plot(log10(freq),20*log10(mag),'k');
set(h,'LineWidth',2);
hold on
% Plot -210 dB asymptote
h=plot(log10(freq),-210*log10(freq/1000),'k--');
hold off
axis([2 4 -80 40]);
set(h,'LineWidth',2);
grid on
xlabel('Log Freq(Hz)');
ylabel('dBV');
title('7th Order Elliptical LPF');
legend('Output','-210 dB Asymptote',0);
figure
%
h=plot(log10(freq),ph,'k');
grid on
xlabel('Log Freq(Hz)');
ylabel('Deg');
title('Elliptical LPF Phase');
axis([2 4 -200 200]);
set(h,'LineWidth',2);

author: Sable

repo: mcbench-benchmarks


text:
%% Copyright (C) 2014, 2016, 2019 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
%% @defmethod @@sym trace (@var{A})
%% Trace of symbolic matrix.
%%
%% Example:
%% @example
%% @group
%% syms x
%% A = [1 2 x; 3 sym(pi) 4; 13 5 2*x]
%%   @result{} A = (sym 3×3 matrix)
%%       ⎡1   2   x ⎤
%%       ⎢          ⎥
%%       ⎢3   π   4 ⎥
%%       ⎢          ⎥
%%       ⎣13  5  2⋅x⎦
%% trace(A)
%%   @result{} ans = (sym) 2⋅x + 1 + π
%% @end group
%% @end example
%%
%% As an example, we can check that the trace of the product is @emph{not}
%% the product of the traces:
%% @example
%% @group
%% A = sym([1 2; 3 4]);
%% B = sym([pi 3; 1 8]);
%% trace(A*B)
%%   @result{} ans = (sym) π + 43
%% trace(A) * trace(B)
%%   @result{} ans = (sym) 5⋅π + 40
%% @end group
%% @end example
%% However, such a property does hold if we use the Kronecker tensor product
%% (@pxref{@@sym/trace}):
%% @example
%% @group
%% kron(A, B)
%%   @result{} ans = (sym 4×4 matrix)
%%       ⎡ π   3   2⋅π  6 ⎤
%%       ⎢                ⎥
%%       ⎢ 1   8    2   16⎥
%%       ⎢                ⎥
%%       ⎢3⋅π  9   4⋅π  12⎥
%%       ⎢                ⎥
%%       ⎣ 3   24   4   32⎦
%% trace(kron(A, B))
%%   @result{} ans = (sym) 5⋅π + 40
%% trace(A) * trace(B)
%%   @result{} ans = (sym) 5⋅π + 40
%% @end group
%% @end example
%%
%% @seealso{@@sym/det}
%% @end defmethod

function z = trace(x)

  cmd = { 'x, = _ins'
          'if not x.is_Matrix:'
          '    x = sp.Matrix([[x]])'
          'return sp.trace(x),' };

  z = pycall_sympy__ (cmd, x);

end

%!test
%! % scalar
%! syms x
%! assert (isequal (trace(x), x))

%!test
%! syms x
%! A = [x 3; 2*x 5];
%! assert (isequal (trace(A), x + 5))

repo: octsympy


text:
function f = tand(f)
%TAND  Tangent of a CHEBFUN3 in degrees.
%   TAND(F) returns the tangent of a CHEBFUN3 object F in degrees.
%
% See also CHEBFUN3/TAN and CHEBFUN3/COMPOSE.

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

% Empty check:
if ( isempty(f) )
    return
end

op = @(x,y,z) tand(feval(f, x, y, z));      % Resample.
f = chebfun3(op, f.domain);                 % Call constructor.

end


text:
function u2 = fd1d_wave_start ( x_num, x_vec, t, t_delta, alpha, u_x1, u_x2, ...
  ut_t1, u1 )

%*****************************************************************************80
%
%% FD1D_WAVE_START takes the first step for the wave equation.
%
%  Discussion:
%
%    This program solves the 1D wave equation of the form:
%
%      Utt = c^2 Uxx
%
%    over the spatial interval [X1,X2] and time interval [T1,T2],
%    with initial conditions:
%
%      U(T1,X)  = U_T1(X),
%      Ut(T1,X) = UT_T1(X),
%
%    and boundary conditions of Dirichlet type:
%
%      U(T,X1) = U_X1(T),
%      U(T,X2) = U_X2(T).
%
%    The value C represents the propagation speed of waves.
%
%    The program uses the finite difference method, and marches
%    forward in time, solving for all the values of U at the next
%    time step by using the values known at the previous two time steps.
%
%    Central differences may be used to approximate both the time
%    and space derivatives in the original differential equation.
%
%    Thus, assuming we have available the approximated values of U
%    at the current and previous times, we may write a discretized
%    version of the wave equation as follows:
%
%      Uxx(T,X) = ( U(T,   X+dX) - 2 U(T,X) + U(T,   X-dX) ) / dX^2
%      Utt(T,X) = ( U(T+dt,X   ) - 2 U(T,X) + U(T-dt,X   ) ) / dT^2
%
%    If we multiply the first term by C^2 and solve for the single
%    unknown value U(T+dt,X), we have:
%
%      U(T+dT,X) =        (     C^2 * dT^2 / dX^2 ) * U(T,   X+dX)
%                  +  2 * ( 1 - C^2 * dT^2 / dX^2 ) * U(T,   X   )
%                  +      (     C^2 * dT^2 / dX^2 ) * U(T,   X-dX)
%                  -                                  U(T-dT,X   )
%
%    (Equation to advance from time T to time T+dT, except for FIRST step!)
%
%    However, on the very first step, we only have the values of U
%    for the initial time, but not for the previous time step.
%    In that case, we use the initial condition information for dUdT
%    which can be approximated by a central difference that involves
%    U(T+dT,X) and U(T-dT,X):
%
%      dU/dT(T,X) = ( U(T+dT,X) - U(T-dT,X) ) / ( 2 * dT )
%
%    and so we can estimate U(T-dT,X) as
%
%      U(T-dT,X) = U(T+dT,X) - 2 * dT * dU/dT(T,X)
%
%    If we replace the "missing" value of U(T-dT,X) by the known values
%    on the right hand side, we now have U(T+dT,X) on both sides of the
%    equation, so we have to rearrange to get the formula we use
%    for just the first time step:
%
%      U(T+dT,X) =   1/2 * (     C^2 * dT^2 / dX^2 ) * U(T,   X+dX)
%                  +       ( 1 - C^2 * dT^2 / dX^2 ) * U(T,   X   )
%                  + 1/2 * (     C^2 * dT^2 / dX^2 ) * U(T,   X-dX)
%                  +  dT *                         dU/dT(T,   X   )
%
%    (Equation to advance from time T to time T+dT for FIRST step.)
%
%    It should be clear now that the quantity ALPHA = C * DT / DX will affect
%    the stability of the calculation.  If it is greater than 1, then
%    the middle coefficient 1-C^2 DT^2 / DX^2 is negative, and the
%    sum of the magnitudes of the three coefficients becomes unbounded.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Second Edition,
%    Prentice Hall, 1999,
%    ISBN: 0-13-012641-1,
%    LC: QA297.P45.
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes in the X direction.
%
%    Input, real X_VEC(X_NUM), the spatial coordinates of the nodes.
%
%    Input, real T, the time after the first step has been taken.
%    In other words, T = T1 + T_DELTA.
%
%    Input, real T_DELTA, the time step.
%
%    Input, real ALPHA, the stability coefficient, computed by FD1D_WAVE_ALPHA.
%
%    Input, real U_X1(T), U_X2(T), functions for the left and right boundary
%    conditions.
%
%    Input, real UT_T1(X), the function that evaluates dUdT at the initial time.
%
%    Input, real U1(X_NUM), the initial condition.
%
%    Output, real U2(X_NUM), the solution at the first time step.
%
  ut = ut_t1 ( x_num, x_vec );

  u2 = zeros ( 1, x_num );

  u2(1) = u_x1 ( t );
  u2(2:x_num-1) =         alpha^2   * u1(3:x_num) / 2 ...
                  + ( 1 - alpha^2 ) * u1(2:x_num-1) ...
                  +       alpha^2   * u1(1:x_num-2) / 2 ...
                  +        t_delta  * ut(2:x_num-1);
  u2(x_num) = u_x2 ( t );

  return
end


text:
function a = reduce_rows(a);
% function a = reduce_rows(a);
% reduction of inequalities of the form A*x <= b where a = [A,b] due to redundancies
%   Input:  a ... matrix of size m x n
%   Output: a ... matrix of size ... x n
%   (c) Sebastian Siegel, created: 2005/06/08, last modified: 2005/07/06

[m,n] = size(a);
a = sortrows(a);	% sort rows according to 1., 2., 3., ... column (ascending)
			% purpose: "normalize" according to first column later

% before actually reducing similar rows, let's also check if there are similarities that differ by a factor (therefore "normalize")

nneg = sum(a(:,1)<0);	% count negative entries in first column
nzer = sum(a(:,1)==0);	% count zero entries in first column
npos = sum(a(:,1)>0);	% count positive entries in first column

a = [a a]; % purpose: will work on "normalized" a and original a
% divide according to entries in first column (if not zero):
a(1:nneg,1:n) = a(1:nneg,1:n)./-(a(1:nneg,1)*ones(1,n)); % the negative part
a(nneg+nzer+1:m,1:n) = a(nneg+nzer+1:m,1:n)./(a(nneg+nzer+1:m,1)*ones(1,n)); % the positive part
% now a(:,1:n) is "normalized" and a(:,n+1:2*n) is still the original a

% sort according to "normalized" a
a = sortrows(a,1:n);

% delete similar rows (according to a(:,1:n-1)) => (save the min value of b):
nsimilar = sum(sum((a(1:m-1,1:n-1) == a(2:m,1:n-1)),2)==n-1); % ...
% (a(1:m-1,1:n-1) == a(2:m,1:n-1)) ... similar neighbored rows (disregard
%    column n which represents b) will produce [1 1 1 ...], other have entries with '0'
% sum( expression above ,2) ... add entries in a row
% expression above == n-1 ... test if all entries in a row were 1
% sum( expression above ) ... get the number of similar rows
if nsimilar > 0, % true if redundant rows exist
	a(:,n+1:2*n)=a(:,n+1:2*n).*([1;(1-(sum((a(1:m-1,1:n-1) == ...
	    a(2:m,1:n-1)),2)==n-1))]*ones(1,n)); % ...
	% [1;(1-(sum((a(1:m-1,1:n-1) == a(2:m,1:n-1)),2)==n-1))] ... column vector
	%      where each '0' represents a row which is similar to the one above
	% ( expression above *ones(1,n)) ... expand column vector to matrix where
	%      each row consists of the same entries
	% a(:,n+1:2*n).* expression above ... now each redundant row has entries with %      zeros only in the right half of a
	a = sortrows(a(:,n+1:2*n)); % sort rows according to n+1., n+2., n+3., ...
	%      column (ascending) of a and save result to a
	%      (now a has only n columns again)
	[temp,position]= max(sum(a==zeros(m,n),2)==n); % ...
	% find position of the first redundant row
	a = [a(1:position-1,:); a(position+nsimilar:m,:)]; % assemble a without
	%      redundant rows
else
	a = a(:,n+1:2*n); % reduce a to original part (right half)
end


text:
 function [alpha, beta, ok] = nufft_alpha(N, J, K, alpha, beta)
%function [alpha, beta, ok] = nufft_alpha(N, J, K, alpha, beta)
%|
%| Determine alpha and beta, as associated with the scaling factors
%| for min-max interpolation, from arguments.
%|
%| in
%|	alpha	[L,1]	Fourier series coefficients of scaling factors
%|			or a string (see below)
%|	beta		scale gamma=2pi/K by this in Fourier series
%|			typically is K/N (me) or 0.5 (Liu)
%| out
%|	alpha,beta
%|
%| Copyright 2001-3-30, Jeff Fessler, University of Michigan

if nargin < 3, ir_usage, end

if ~isvar('alpha') || isempty(alpha)
	alpha = [1]; % default Fourier series coefficients of scaling factors
end
if ~isvar('beta') || isempty(beta)
	beta = 0.5; % default is Liu version for now
end

if streq(alpha, 'uniform')
	alpha = [1];
	beta = 0.5;
return
end

% see if 'best' alpha is desired
if ischar(alpha)
	if streq(alpha, 'best')
		L = 0;
	elseif streq(alpha, 'best,L=1')
		L = 1;
	elseif streq(alpha, 'best,L=2')
		L = 2;
	else
		fail 'unknown alpha argument'
	end

	[alpha, beta, ok] = nufft_best_alpha(J, L, K/N);
	if ~ok
		warn('optimal alpha unknown for J=%d, K/N=%g, L=%d', J, K/N, L)
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function value = r8_uniform_pdf ( lower, upper, rval )

%*****************************************************************************80
%
%% R8_UNIFORM_PDF evaluates the PDF of a uniform distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real LOWER, UPPER, the lower and upper range limits.
%    LOWER < UPPER.
%
%    Input, real RVAL, the point where the PDF is evaluated.
%
%    Output, real VALUE, the value of the PDF at RVAL.
%
  if ( upper <= lower )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_UNIFORM_PDF - Fatal error!\n' );
    fprintf ( 1, '  For uniform PDF, the lower limit must be \n' );
    fprintf ( 1, '  less than the upper limit.\n' );
    error ( 'R8_UNIFORM_PDF - Fatal error!' );
  end

  if ( rval < lower )
    value = 0.0;
  elseif ( rval <= upper )
    value = 1.0 / ( upper - lower );
  else
    value = 0.0;
  end

  return
end


text:
function sigmaEq = equivalentStress(sigma,varargin)
% von Mises equivalent stress
%
% Input
%  sigma - @stressTensor
%
% Output
%  s - double
%

s = sigma.deviatoricStress;
sigmaEq = sqrt(3/2 * s:s);

% the following code gives the same result
% is more efficient but less readable

% M = sigma.M;
% sigmaEq = sqrt(M(1,1,:).^2 + M(2,2,:).^2 + M(3,3,:).^2 ...
%   - M(1,1,:).*M(2,2,:) - M(2,2,:).*M(3,3,:)- M(1,1,:).*M(3,3,:) ...
%   + 3*(M(1,2).^2 + M(1,3).^2 + M(2,3).^2 ));


text:
function UNew = fluidDirichlet3D(varargin);
% fluidDirichlet3D: solve fluid registraion in 3D with Dirichlet
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin{:});

% multiply F by adjoint of Navier-Lame equations
FNew = adjointNL(F/RegularizerFactor,mu,lambda,0,M,N,P);

% compute sine transform of new force field
FS = discreteSineTransform(FNew,M,N,P);

% construct images of coordinates scaled by pi/(N or M or P)
[a,b,c] = ndgrid(pi*(0:(M-1))/(M-1),pi*(0:(N-1))/(N-1),pi*(0:(P-1))/(P-1));

% construct LHS factor
LHSfactor = mu.*(lambda+2*mu).*(2*cos(a) + 2*cos(b) + 2*cos(c) - 6).^2;

% if gamma is zero, set origin term to 1, as DC term does not matter
LHSfactor(1,1,1) = 1;

% solve for FFT of U
VS = cat(4,FS(:,:,:,1)./LHSfactor,FS(:,:,:,2)./LHSfactor,FS(:,:,:,3)./LHSfactor);

% perform inverse DST
V = discreteSineTransform(VS,M,N,P);

% now perform Euler integration to construct new displacements
UNew = zeros(M,N,P,3);
UNew(:,:,:,1) = (1 - imfilter(V(:,:,:,1),HX,'replicate','same')).*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HY,'replicate','same').*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HZ,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,2) = -imfilter(V(:,:,:,1),HY,'replicate','same').*V(:,:,:,1) + ...
    (1 - imfilter(V(:,:,:,2),HY,'replicate','same')).*V(:,:,:,2) - ...
    imfilter(V(:,:,:,3),HY,'replicate','same').*V(:,:,:,3);
UNew(:,:,:,3) = -imfilter(V(:,:,:,1),HZ,'replicate','same').*V(:,:,:,1) - ...
    imfilter(V(:,:,:,2),HZ,'replicate','same').*V(:,:,:,2) + ...
    (1 - imfilter(V(:,:,:,3),HZ,'replicate','same')).*V(:,:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FS = discreteSineTransform(F,M,N,P);
% compute discrete sine transform of 3-D vector field

% initialize resulting array
FS = F;

% first perform sine transform down columns
len = 2*M-2; ind = 1:M;
for p=1:P
    for n=1:N
        s = fft(FS(:,n,p,:),len,1);
        FS(:,n,p,:) = imag(s(ind,:,:,:));
    end
end
FS = sqrt(2/(M-1))*FS;

% next perform sine transform across rows
len = 2*N-2; ind = 1:N;
for p=1:P
    for m=1:M
        s = fft(FS(m,:,p,:),len,2);
        FS(m,:,p,:) = imag(s(:,ind,:,:));
    end
end
FS = sqrt(2/(N-1))*FS;

% finally perform sine transform across pages
len = 2*P-2; ind = 1:P;
for n=1:N
    for m=1:M
        s = fft(FS(m,n,:,:),len,3);
        FS(m,n,:,:) = imag(s(:,:,ind,:));
    end
end
FS = sqrt(2/(P-1))*FS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FNew = adjointNL(F,mu,lambda,gamma,M,N,P);
% multiply vector field F by adjoint Navier-Lame equations

% initialize FNew
FNew = zeros(M,N,P,3);

% construct filter that implements 3-D Laplacian
L = (lambda+2*mu)*cat(3,[0 0 0;0 1 0;0 0 0],[0 1 0;1 -6 1;0 1 0],[0 0 0;0 1 0;0 0 0]);

% we will need to use L to form two different filters
% L1 = -(lambda+2*mu)*L; L1(2,2,2) = gamma + L1(2,2,2);
% L2 = -mu*L; L2(2,2,2) = gamma + L2(2,2,2);

% construct grad div filters
GD11 = (lambda+mu)*cat(3,zeros(3,3),[0 1 0;0 -2 0;0 1 0],zeros(3,3));
GD22 = ipermute(GD11,[2 1 3]);
GD33 = ipermute(GD11,[3 2 1]);
GD23 = zeros(3,3,3);
GD23(2,1,1) = 1; GD23(2,3,3) = 1; GD23(2,1,3) = -1; GD23(2,3,1) = -1;
GD23 = GD23*(lambda+mu)/4;
GD12 = ipermute(GD23,[3 1 2]);
GD13 = ipermute(GD23,[2 3 1]);

% perform filtering
FNew(:,:,:,1) = imfilter(F(:,:,:,1),L-GD11,'replicate') + ...
    imfilter(F(:,:,:,2),-GD12,'replicate') + ...
    imfilter(F(:,:,:,3),-GD13,'replicate');
FNew(:,:,:,2) = imfilter(F(:,:,:,1),-GD12,'replicate') + ...
    imfilter(F(:,:,:,2),L-GD22,'replicate') + ...
    imfilter(F(:,:,:,3),-GD23,'replicate');
FNew(:,:,:,3) = imfilter(F(:,:,:,1),-GD13,'replicate') + ...
    imfilter(F(:,:,:,2),-GD23,'replicate') + ...
    imfilter(F(:,:,:,3),L-GD33,'replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DU,F,mu,lambda,PixSize,M,N,P,RegularizerFactor,HX,HY,HZ] = parse_inputs(varargin);

% get displacement field and check size
F = varargin{2};
PixSize = varargin{4}(1:3);
M = varargin{5};
N = varargin{6};
P = varargin{7};
mu = varargin{8};
lambda = varargin{9};
RegularizerFactor = varargin{10};
DU = varargin{11};
HX = varargin{12};
HY = varargin{13};
HZ = varargin{14};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
%% RRT for 3 DOF Mobile Robot
%
%

% Do not change anything in rrt.m,
% rrt_star.m and rrt_star_fn.m
%

%

%% Getting started
% Create *main_3dof_mobile_rrt.m* file. You can create m-file with any other name, that is
% perfectly okay, this will not affect to the solution of a path/motion
% planning problem. All the sources could be found in examples/ directory
% of the distribution.
%

%% Step 1: Choosing the map
% Firstly, we should define what map to use. It is done by defining map
% structure the sample code is goes after.
%
%   map = struct('name', 'bench_june1.mat', 'start_point', [-14.5 -7.5], 'goal_point', [10 -0.65]);
%

%%
%
% * *name field* defines the file name in maps/ directory
% * *start_point* defines the initial point of the problem
% * *goal_point* defines the goal point on the map
%

%% Step 2: Setting maximum number of iterations
%
%   max_iter = 20e3;
%
% * *max_iter* variable defines how much iteration should be done to solve
% the path planning problem.

%% Step 3: Do we have to benchmark
%
%   is_benchmark = false;
%
% * *is_benchmark* enables benchmarking. For more details please read the
% sources of rrt.m, rrt_star.m and rrt_star_fn.m

%% Step 4: Setting random seed
%
%   rand_seed = 40;
%
% * *rand_seed* variable is a random seed for random number generator, it
% is used for sampling nodes. It is usually used for benchmarking. We set
% the same map, however the random seed is different for the same map. You
% can use *now* if you don't care about random seed.
%
%   rand_seed = now;
%

%% Step 5: Choosing the class (model) we want
%
%   variant = 'FNSimple3D';
%
% * *variant* defines from what class we should instantiate the object. In
% other words it defines what model we choose for application of RRT.
%
% *FNSimple3D* is a name of a class which contains all necessary methods
% and fields in order to represent simple 3 DOF Mobile Robot model.

%% Step 6: Rapidly-Exploring Random Tree (RRT)
%
%   rrt(map, max_iter, is_benchmark, rand_seed, variant);
%
% Line above runs RRT with given parameters. In addition, *rrt* function
% returns the class object with a certain solution.

%% Sources of *main_3dof_mobile_rrt.m*
% Press
% <matlab:edit('examples/main_3dof_mobile_rrt.m') here>
% to play with example code.
%
%   % 3 DOF mobile robot example.
%   % by Almaskhan Baimyshev
%   % 08/28/2013
%
%   map = struct('name', 'bench_june1.mat', 'start_point', [-14.5 -7.5], 'goal_point', [10 -0.65]);
%   max_iter = 20e3;
%   is_benchmark = false;
%   rand_seed = 40;
%   variant = 'FNSimple3D';
%   result = rrt(map, max_iter, is_benchmark, rand_seed, variant);
%

author: olzhas


function [b,R,t,s] = FitEdges( im,xp,landmarks,shapePC,shapeMU,shapeEV,Ef,Ev,tri,ndims,w_prior, w_edges, w_landmarks,niter )
%FITEDGES Perform morphable model fitting using edges
%   This function initialises by fitting to a sparse set of landmarks, then
%   iteratively fits to edges by finding nearest neighbour between
%   projected model edges and image edges
%
% Inputs:
%   im        - Image to fit to
%   xp        - 2 by nlandmarks matrix containing locations of 2D landmarks
%   landmarks - Indices of vertices corresponding to points in xp
%   shapePC   - 3DMM principal components
%   shapeMU   - 3DMM average shape
%   shapeEV   - standard deviations of each 3DMM dimension
%   Ef        - nedges by 2 matrix storing faces adjacent to each edge
%   Ev        - nedges by 2 matrix storing vertices adjacent to each edge
%   tri       - Face structure of 3DMM mesh
%   ndims     - Number of model dimensions to use
%   numsd     - Hard limit on number of standard deviations each parameter
%               may deviate
%
% Outputs:
%   b         - Estimated shape parameters
%   R,t,s     - Estimated rotation, translation and scale of face

FV.faces = tri;
% PARAMETERS
% Number of model dimensions used for initial landmark fit
%ndims = 40;
% Number of edge fitting iterations (may want to replace with a convergence
% test)
%niter = 1;
% Proportion of nearest-neighbour edge matches used at each iteration, may
% be more sensible to threshold on distance to nearest neighbour
percentile = 1;

% Perform initial landmark-only fit
%w_prior = 0.7;
[b,R,t,s] = FitSingleSOP( xp,shapePC,shapeMU,shapeEV,ndims,landmarks,w_prior );
%b=zeros(20,1);
FV.vertices = reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';

% Display input image
%figure; imshow(im);

% Extract image edges
edges = edge(rgb2gray(im),'canny',0.15);
%figure; imshow(edges);

% Store pixel locations of edges
[r,c]=find(edges);
r = size(edges,1)+1-r;
%bold=b;
for iter=1:niter
    % Compute vertices lying on occluding boundary
    [ occludingVertices ] = occludingBoundaryVertices( FV,Ef,Ev,R );

    % Project occluding boundary vertices
    x2 = R*FV.vertices(occludingVertices,:)';
    x2 = x2(1:2,:);
    x2(1,:) = x2(1,:)+t(1);
    x2(2,:) = x2(2,:)+t(2);
    x2 = x2.*s;

    % Find edge correspondences
    [idx,d] = knnsearch([c r],x2'); %its correct, searching based on columns dont change!

    % Filter edge matches - probably want to do something better here
    sortedd=sort(d);
    threshold = sortedd(round(percentile*length(sortedd)));
    idx = idx(d<threshold);
    occludingVertices = occludingVertices(d<threshold);

    [b,R,t,s] = FitSingleSOP( [c(idx)' xp(1,:); r(idx)' xp(2,:)],shapePC,shapeMU,shapeEV,ndims,[occludingVertices; landmarks],w_edges,w_landmarks,length(landmarks));
    % Note: this is completely refitting from scratch. We could just
    % re-start the nonlinear optimisation using previous estimates as
    % initialisation but this doesn't seem to work well.
    %disp(num2str(norm(b-bold)));
    %differb(iter)=norm(b-bold);
    FV.vertices = reshape(shapePC(:,1:ndims)*b+shapeMU,3,size(shapePC,1)/3)';
    %bold=b;
end

end


text:
function J = evalObjectiveFCN(u,x,xref,Q,R,Ru)
% [J,Js,Ju]
%% Cost function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%   xref:   state references, constant from time k+1 to k+N
%   u0:     previous controller output at time k-1
%
% Output:
%   J:      objective function cost
%

%% Nonlinear MPC design parameters

%% Cost Calculation
% Set initial cost and input
N = size(x,2);
J = zeros(N,1);
u0 = 0;

% Loop through each time step
for ct=1:N

    % Accumulate state tracking cost from x(k+1) to x(k+N)
    J(ct) = (x(:,ct)-xref(:,ct))'*Q*(x(:,ct)-xref(:,ct));

    % Accumulate MV rate of change cost from u(k) to u(k+N-1)
    if ct==1
        J(ct) = J(ct) + (u(:,ct)-u0)'*R*(u(:,ct)-u0) + u(:,ct)'*Ru*u(:,ct);
    else
        J(ct) = J(ct) + (u(:,ct)-u(:,ct-1))'*R*(u(:,ct)-u(:,ct-1)) + u(:,ct)'*Ru*u(:,ct);
    end
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function varargout = cos(varargin)
%COS   Cosine of a DISKFUN.
%   COS(F) returns the cosine of F.
%
% See also DISKFUN/COSH and DISKFUN/SIN

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

[varargout{1:nargout}] = cos@separableApprox(varargin{:});

end


text:
function [ node_num, element_num, hole_num ] = example1_q4_mesh_size ( )

%*****************************************************************************80
%
%% EXAMPLE1_Q4_MESH_SIZE sets sizes for example #1 Q4 mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer ELEMENT_NUM, the number of elements.
%
%    Output, integer HOLE_NUM, the number of holes.
%
  element_num = 17;
  hole_num = 1;
  node_num = 28;

  return
end


text:
function kern = diagKernParamInit(kern)

% DIAGKERNPARAMINIT DIAG kernel parameter initialisation.
% The diag covariance function takes a one dimensional input and outputs a diagonal noise that is provided by an exponentiated and scaled version of the input.
%
% k(x_i, x_j) = delta_ij sigma2 exp(x_i)
%
% The only parameter is sigma2, the process variance (kern.variance).
%
% SEEALSO : whiteKernParamInit
%
% FORMAT
% DESC initialises the diagonal noise covariance function
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

kern.variance = exp(-2);
kern.nParams = 1;

kern.transforms.index = 1;
kern.transforms.type = optimiDefaultConstraint('positive');
kern.trans = optimiDefaultConstraint('positive');

kern.isStationary = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% sub2indn
% Below is a demonstration of the features of the |sub2indn| function

%%
clear; close all; clc;

%% Syntax
% |ind = sub2indn(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%  Figure 3.14     Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig3_14.m
clf;
num=[2 1];
den=[1 3 2];
axis ('square')
pzmap(num,den)
grid


text:
function display_legend(plot_obj)
% Displays the legend on the plot.  Call this function only once
% all curves have been plotted.

if length(plot_obj.legend_strings) > 0
    figure(plot_obj.fh)
    legend(plot_obj.handles_vec,plot_obj.legend_strings)
end
end

author: nesl


text:
function [overlaps] = cat_overlap(varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%% concatenate overlap %%%%%%%%%%%%%%%%%%%%%%%%%%%
% Allows daisy-chaining of multiple registers2p regi variables. NB all must
% have been registered to the same time-point 1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
% -------
% - none: opens a file selection gui, select multiple regi variables

% - cell-array of strings: each cell contains the full file-path to a regi
% variable

% - structure: multiple regi variables saved into one structure with
% fields .rois, .targets and .dat (and their corresponding sub-fields)

% Output:
% -------
% - overlaps.fnames: cell array where each cell contains the name of an
% F...proc.mat file containing the ROIs from which the correspondingly
% indexed columns in the below fields are drawn (i.e. overlaps.fnames{1}
% gives the filename for overlaps.rois.idcs(:,1))

% - overlaps.rois.idcs: n * m matrix where n = number of overlapping ROIs
% and m = timepoints. Each element of a given column can be directly
% indexed into the variables Fcell & FcellNeu of the F...proc.mat file
% dictated by the corresponding cell in overlaps.fnames (see above).

% - overlaps.rois.iscell_idcs: same as overlaps.rois.idcs except these
% indices correspond to the non-zero elements of the iscell classifier. To
% index into Fcell & FcellNeu do:
%       nn_idcs = find(iscell);
%       Fcell(nn_idcs(overlaps.rois.iscell_idcs(:,1)))

% - overlaps.targets.idcs: n * m matrix where n = number of targets user
% imported to registers2p and m = timepoints. Elements are either NaN if
% target didn't overlap with a suite2P ROI, or else the index of the
% overlapping ROI in Fcell & Fcell Neu of F...proc.mat file dictated by the
% corresponding cell in overlaps.fnames (see above).

% - overlaps.targets.iscell_idcs: Same as overlaps.targets.idcs except
% non-NaN elements are indices of non-zero elements of the iscell
% classifier (see overlaps.rois.iscell_idcs for an example of how to index
% this into Fcell & FcellNeu).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
overlaps = [];
f = {};
d = {};
if ~isempty(varargin)
    if iscell(varargin{1})
        if ischar(varargin{1}{1})
            in = varargin{1};
        end
        for i = 1:numel(in)
            [d{i},ff,ee] = fileparts(in{i});
            f{i} = [ff ee];
        end
    elseif isstruct(varargin{1})
        f = varargin{1};
    end
else
    [f,in] = uigetfile('MultiSelect','on');
    if ~iscell(f)
        f = {f};
    end
    for i = 1:numel(f)
        d{i} = in;
    end
end
if iscell(f)
    rois_idcs = cell(numel(f),1);
    rois_iscell_idcs = rois_idcs;
    targets_idcs = rois_idcs;
    targets_iscell_idcs = rois_idcs;
    overlaps.fname = cell(numel(f)+1,1);
    for i = 1:numel(f)
        load([d{i} filesep f{i}])
        rois_idcs{i} = regi.rois.idcs;
        rois_iscell_idcs{i} = regi.rois.iscell_idcs;
        targets_idcs{i} = regi.targets.idcs;
        targets_iscell_idcs{i} = regi.targets.iscell_idcs;
        if i == 1
            [~,overlaps.fname{1}] = fileparts(regi.dat.files(1).filepath);
            [~,overlaps.fname{2}] = fileparts(regi.dat.files(2).filepath);
        else
            [~,overlaps.fname{i+1}] = fileparts(regi.dat.files(2).filepath);
        end
    end
    overlaps.rois.idcs = find_consistent_rois(rois_idcs);
    overlaps.rois.iscell_idcs = find_consistent_rois(rois_iscell_idcs);
    if max(cellfun('isempty',targets_idcs)) == 0
        overlaps.targets.idcs = [targets_idcs{1}(:,1) cell2mat(cellfun(@(x) x(:,2),targets_idcs,'UniformOutput',0)')];
        overlaps.targets.iscell_idcs = [targets_iscell_idcs{1}(:,1) cell2mat(cellfun(@(x) x(:,2),targets_iscell_idcs,'UniformOutput',0)')];
    else
        overlaps.targets.idcs = [];
        overlaps.targets.iscell_idcs = [];
    end
elseif isstruct(f)
    rois_idcs = cell(numel(f),1);
    rois_iscell_idcs = rois_idcs;
    targets_idcs = rois_idcs;
    targets_iscell_idcs = rois_idcs;
    overlaps.fname = cell(numel(f)+1,1);
    for i = 1:numel(f)
        rois_idcs{i} = f(i).rois.idcs;
        rois_iscell_idcs{i} = f(i).rois.iscell_idcs;
        targets_idcs{i} = f(i).targets.idcs;
        targets_iscell_idcs{i} = f(i).targets.iscell_idcs;
        if i == 1
            [~,overlaps.fname{1}] = fileparts(f(i).dat.files(1).filepath);
            [~,overlaps.fname{2}] = fileparts(f(i).dat.files(2).filepath);
        else
            [~,overlaps.fname{i+1}] = fileparts(f(i).dat.files(2).filepath);
        end
    end
    overlaps.rois.idcs = find_consistent_rois(rois_idcs);
    overlaps.rois.iscell_idcs = find_consistent_rois(rois_iscell_idcs);
    if max(cellfun('isempty',targets_idcs)) == 0
        overlaps.targets.idcs = [targets_idcs{1}(:,1) cell2mat(cellfun(@(x) x(:,2),targets_idcs,'UniformOutput',0)')];
        overlaps.targets.iscell_idcs = [targets_iscell_idcs{1}(:,1) cell2mat(cellfun(@(x) x(:,2),targets_iscell_idcs,'UniformOutput',0)')];
    else
        overlaps.targets.idcs = [];
        overlaps.targets.iscell_idcs = [];
    end
else
    warning('Import error - either provide no input and use file-selection gui or provide a cell-array of filepath strings')
end
end

author: cortex-lab


text:
%TEST_bezierFunction.m
%
% This script tests the use of a bezier curve to represent a vector
% function, rather than an arbitrary space curve
%

order = 2;

p = rand(1,order+1);
w = 0.05 + 0.95*rand(1,order+1);

tSpan = [1,4];
t = linspace(tSpan(1),tSpan(2),100);
pGrid = linspace(tSpan(1),tSpan(2),order+1);

tic
x = rationalBezierCurve(p,w,t,tSpan);
toc

figure(1); clf;
plot(t,x); hold on;
scatter(pGrid,p,50*w,'filled');

author: MatthewPeterKelly

repo: dscTutorials


text:
function region = p21_region ( )

%*****************************************************************************80
%
%% P21_REGION returns the name of the integration region for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'sphere';

  return
end

author: johannesgerer


text:
function value = h_03 ( x )

%*****************************************************************************80
%
%% H_03 evaluates x^3+x^2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which F is to be evaluated.
%
%    Output, real VALUE, the value of the function at X.
%
  value = x .* x .* ( x + 1.0 );

  return
end


text:
function [ r, seed ] = r8row_uniform_abvec ( m, n, a, b, seed )

%*****************************************************************************80
%
%% R8ROW_UNIFORM_ABVEC fills an R8ROW with scaled pseudorandom numbers.
%
%  Discussion:
%
%    An R8ROW is an array of R8 values, regarded as a set of row vectors.
%
%    The user specifies a minimum and maximum value for each column.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in
%    the array.
%
%    Input, real A(N), B(N), the lower and upper limits.
%
%    Input/output, integer SEED, the "seed" value, which
%    should NOT be 0.  On output, SEED has been updated.
%
%    Output, real R(M,N), the array of pseudorandom values.
%
  i4_huge = 2147483647;

  for i = 1 : m

    for j = 1 : n

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + i4_huge;
      end

      r(i,j) = a(j) + ( b(j) - a(j) ) * seed * 4.656612875E-10;

    end
  end

  return
end


text:
function stateLogRow = getStateLogRowAtTime(stateLog, time)
    if(time <= min(stateLog(:,1)))
        stateLogRow = stateLog(1,:);
    elseif(time >= max(stateLog(:,1)))
        stateLogRow = stateLog(end,:);
    else
        frameTimeLwr = max(stateLog(stateLog(:,1)<=time,1));
        stateLogRowLwr = stateLog(stateLog(:,1)==frameTimeLwr,:);
        stateLogRowLwr = stateLogRowLwr(1,:);

        frameTimeUpr = min(stateLog(stateLog(:,1)>time,1));
        stateLogRowUpr = stateLog(stateLog(:,1)==frameTimeUpr,:);
        stateLogRowUpr = stateLogRowUpr(1,:);

        interpX = [stateLogRowLwr(1), stateLogRowUpr(1)];
        interpYx = [stateLogRowLwr(2), stateLogRowUpr(2)];
        interpYy = [stateLogRowLwr(3), stateLogRowUpr(3)];
        interpYz = [stateLogRowLwr(4), stateLogRowUpr(4)];
        interpYVx = [stateLogRowLwr(5), stateLogRowUpr(5)];
        interpYVy = [stateLogRowLwr(6), stateLogRowUpr(6)];
        interpYVz = [stateLogRowLwr(7), stateLogRowUpr(7)];

        posX = interp1(interpX,interpYx,time);
        posY = interp1(interpX,interpYy,time);
        posZ = interp1(interpX,interpYz,time);
        velX = interp1(interpX,interpYVx,time);
        velY = interp1(interpX,interpYVy,time);
        velZ = interp1(interpX,interpYVz,time);

        stateLogRow = stateLogRowLwr;
        stateLogRow(2:7) = [posX, posY, posZ, velX, velY, velZ];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef MaF13 < PROBLEM
% <multi/many> <real> <large/none>
% P7

%------------------------------- Reference --------------------------------
% R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin, and X. Yao, A
% benchmark test suite for evolutionary many-objective optimization,
% Complex & Intelligent Systems, 2017, 3(1): 67-81.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = 5; end
            obj.M        = max(obj.M,3);
            obj.lower    = [zeros(1,2),zeros(1,obj.D-2)-2];
            obj.upper    = [ones(1,2),zeros(1,obj.D-2)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            Y = X - 2*repmat(X(:,2),1,D).*sin(2*pi*repmat(X(:,1),1,D)+repmat(1:D,N,1)*pi/D);
            PopObj(:,1) = sin(X(:,1)*pi/2)                   + 2*mean(Y(:,4:3:D).^2,2);
            PopObj(:,2) = cos(X(:,1)*pi/2).*sin(X(:,2)*pi/2) + 2*mean(Y(:,5:3:D).^2,2);
            PopObj(:,3) = cos(X(:,1)*pi/2).*cos(X(:,2)*pi/2) + 2*mean(Y(:,3:3:D).^2,2);
            PopObj(:,4:obj.M) = repmat(PopObj(:,1).^2+PopObj(:,2).^10+PopObj(:,3).^10+2*mean(Y(:,4:D).^2,2),1,obj.M-3);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = R./repmat(sqrt(sum(R.^2,2)),1,3);
            R = [R,repmat(R(:,1).^2+R(:,2).^10+R(:,3).^10,1,obj.M-3)];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 3
                a = linspace(0,pi/2,10)';
                R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% copyright by Jianxiong Xiao http://mit.edu/jxiao
% demo how to sample uniformly on a sphere

%{
Please cite this paper if you use this code in your publication:
J. Xiao, T. Fang, P. Zhao, M. Lhuillier, and L. Quan
Image-based Street-side City Modeling
ACM Transaction on Graphics (TOG), Volume 28, Number 5
Proceedings of ACM SIGGRAPH Asia 2009
%}

clear
clc

tic;
points = icosahedron2sphere(0);
toc;
subplot(1,4,1);
plot3(points(:,1),points(:,2),points(:,3),'.')
title(sprintf('Level %d with %d points',0,size(points,1)))
axis equal
axis tight

tic;
points = icosahedron2sphere(1);
toc;
subplot(1,4,2);
plot3(points(:,1),points(:,2),points(:,3),'.')
title(sprintf('Level %d with %d points',1,size(points,1)))
axis equal
axis tight

tic;
points = icosahedron2sphere(2);
toc;
subplot(1,4,3);
plot3(points(:,1),points(:,2),points(:,3),'.')
title(sprintf('Level %d with %d points',2,size(points,1)))
axis equal
axis tight

tic;
points = icosahedron2sphere(4);
toc;
subplot(1,4,4);
plot3(points(:,1),points(:,2),points(:,3),'.')
title(sprintf('Level %d with %d points',3,size(points,1)))
axis equal
axis tight

author: zouchuhang

repo: LayoutNet


text:
%*****************************************************************************************************
%NAME: esemilogy.m
%AUTHOR: Andri M. Gretarsson
%DATE: 01/21/98
%
%SYNTAX:	esemilogy(X,Y, 'colour')
%
%Note that 'colour' is NOT an optional argument.
%
%This function acts just like the built-in function 'semilogy' but plots error-bars. The error-bars
%are plotted in the colour given by 'colour'.  'X' and 'Y' are Nx2 matrixes, the first column
%representing the values of the coordinate (x or y), the second column representing the uncertainty
%('error') of those values.  'colour' is a one-letter string which must be one of the letters allowed
%in the built-in matlab function 'plot', to specify the plot colour.  Note that the function exits
%with "hold" set to "off".
%
%This function does not print points in addition to the error bars.  Where the error bars cross, is
%the coordinate point.  This means that if both error bars are exceedingly small complared to the
%coordinate values, the mark will be correspondingly small.  In such situations, it may be better to
%use 'semilogy' directly and specify that the error is smaller than the size of the mark.
%
%EXAMPLE:
%
%X=[1.0	0.2
%   2.0	0.2]
%Y=[1.0 	0.25
%   2.0	0.25]
%esemilogy(X,Y,'g')
%
%plots a green cross of width 0.2 and height 0.25 at coordinate (1.0,1.0), and a cross of width 0.2
%and height 0.25 at coordinate( (2.0,2.0), on a linear x, log y scale
%
%LAST MODIFIED:  01/21/98
%*****************************************************************************************************

function esemilogx=semilogxplot(x,y,colourstring)

xvalue=x(:,1);												%For clarity
xerror=x(:,2);
yvalue=y(:,1);
yerror=y(:,2);

semilogy(xvalue-xerror,yvalue-yerror,'w-',xvalue+xerror,yvalue+yerror,'w-'); hold on;
%Sets appropriate axes but otherwise invisible on a white background.

for i=1:length(xvalue)
   semilogy([xvalue(i)-xerror(i) xvalue(i)+xerror(i)], [yvalue(i) yvalue(i)], colourstring);
   semilogy([xvalue(i) xvalue(i)], [yvalue(i)-yerror(i) yvalue(i)+yerror(i)], colourstring);
end
hold off;


text:
function X = slpcarecon(S, Y)
%SLPCARECON Reconstructs the samples in original space
%
% $ Syntax $
%   - Xr = slpcarecon(S, Y)
%
% $ Arguments $
%   - S:        the PCA model struct
%   - Y:        the principal component features
%   - Xr:       the reconstructed samples
%
% $ Description $
%   - Xr = slpcarecon(S, Y) reconstructs the original samples approximately
%     using the principal components Y. If the dimension of Y is less than
%     the subspace dimension, the leading space dimensions will be used.
%
% $ History $
%   - Created by Dahua Lin, on Aug 17, 2006
%   - Modified by Dahua Lin, on Sep 10, 2006
%       - replace sladd by sladdvec to increase efficiency
%

%% parse and verify input

if ~isstruct(S)
    error('sltoolbox:invalidarg', ...
        'S should be a PCA model struct');
end

if ~isnumeric(Y) || ndims(Y) ~= 2
    error('sltoolbox:invalidarg', ...
        'The features Y should be a 2D numeric matrix');
end

dy = size(Y, 1);
if dy > S.feadim
    error('sltoolbox:sizmismatch', ...
        'The feature dimension of Y exceeds the subspace dimension preserved in model');
end

%% reconstruct

if dy == S.feadim
    X = S.P * Y;
else
    X = S.P(:, 1:dy) * Y;
end

X = sladdvec(X, S.vmean, 1);

author: lmthang


text:
function C = xcorrn( A, T, shape )
% n-dimensional cross-correlation.  Generalized version of xcorr2.
%
% For 2 dimensional inputs this function is exactly the same as xcorr2, but
% also works in higher dimensions. Can also be more efficient because it
% performs convolution using convnFast.  Note that xcorr2 is part of the
% 'Signal Processing Toolbox' and may not be available on all systems.
%
% The order of parameters is reversed from normxcorrn. This is to be
% compatible with the matlab functions normxcorr2 and xcorr2 (which take
% parameters in different orders).
%
% USAGE
%  C = xcorrn( A, T, [shape] )
%
% INPUTS
%  A           - first d-dimensional matrix
%  T           - second d-dimensional matrix
%  shape       - ['full'] 'valid', or 'same' (see convn)
%
% OUTPUTS
%  C           - correlation matrix
%
% EXAMPLE
%  T=gaussSmooth(rand(20),2); A=repmat(T,[3 3]);
%  C1=xcorr2(A,T); C2=xcorrn(A,T); C3=rot90(xcorrn(T,A),2);
%  figure(1); im(C1);  figure(2); im(C2);  figure(3); im(C3);
%
% See also XCORR2, NORMXCORRN, XEUCN, CONVNFAST
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.12
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(shape)); shape='full'; end
nd = ndims(A);
if(nd~=ndims(A)); error('xcorrn: T and A must have same ndims'); end;

% flip for conv purposes
if(nd==2); T=rot90(T,2); else for d=1:nd; T=flipdim(T,d); end; end

% convolve [in frequency or spatial domain]
C = convnFast( A, T, shape );

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function combo_test10 ( )

%*****************************************************************************80
%
%% COMBO_TEST10 tests I4VEC_SEARCH_BINARY_D and I4VEC_SORT_INSERT_D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2011
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMBO_TEST10\n' );
  fprintf ( 1, '  Integer vectors:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I4VEC_SORT_INSERT_D descending sorts;\n' );
  fprintf ( 1, '  I4VEC_SEARCH_BINARY_D searches a descending \n' );
  fprintf ( 1, '  sorted vector.\n' );

  a(1:n) = [ 6, 7, 1, 0, 4, 3, 2, 1, 5, 8 ]';

  i4vec_print ( n, a, '  Before descending sort:' );

  a = i4vec_sort_insert_d ( n, a );

  i4vec_print ( n, a, '  After descending sort:' );

  b = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now search for an instance of the value %d\n', b );

  index = i4vec_search_binary_d ( n, a, b );

  fprintf ( 1, '\n' );
  if ( index == 0 )
    fprintf ( 1, '  The value does not occur.\n' );
  else
    fprintf ( 1, '  The value occurs at index = %d\n', index );
  end

  return
end

repo: jburkardt-m


text:
function [ n_data, x, fx ] = i1ml1_values ( n_data )

%*****************************************************************************80
%
%% I1ML1_VALUES returns some values of the I1ML1 function.
%
%  Discussion:
%
%    The function is defined by:
%
%      I1ML1(x) = I1(x) - L1(x)
%
%    I1(x) is the modified Bessel function of the first kind of order 1,
%    L1(x) is the modified Struve function of order 1.
%
%    The data was reported by McLeod.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
     0.97575346155386267134E-03, ...
     0.77609293280609272733E-02, ...
     0.59302966404545373770E-01, ...
     0.20395212276737365307E+00, ...
     0.33839472293667639038E+00, ...
     0.48787706726961324579E+00, ...
     0.59018734196576517506E+00, ...
     0.62604539530312149476E+00, ...
     0.63209315274909764698E+00, ...
     0.63410179313235359215E+00, ...
     0.63417966797578128188E+00, ...
     0.63439268632392089434E+00, ...
     0.63501579073257770690E+00, ...
     0.63559616677359459337E+00, ...
     0.63591001826697110312E+00, ...
     0.63622113181751073643E+00, ...
     0.63636481702133606597E+00, ...
     0.63650653499619902120E+00, ...
     0.63655609126300261851E+00, ...
     0.63657902087183929223E+00 ];

  x_vec = [ ...
       0.0019531250E+00, ...
       0.0156250000E+00, ...
       0.1250000000E+00, ...
       0.5000000000E+00, ...
       1.0000000000E+00, ...
       2.0000000000E+00, ...
       4.0000000000E+00, ...
       8.0000000000E+00, ...
      12.0000000000E+00, ...
      16.0000000000E+00, ...
      16.2500000000E+00, ...
      17.0000000000E+00, ...
      20.0000000000E+00, ...
      25.0000000000E+00, ...
      30.0000000000E+00, ...
      40.0000000000E+00, ...
      50.0000000000E+00, ...
      75.0000000000E+00, ...
     100.0000000000E+00, ...
     125.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function x = inv_triu(U)
% INV_TRIU     Invert upper triangular matrix.

x = solve_triu(U,eye(size(U)));
%x = inv(U);

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function omcleanup (modelname)
%
% Clean up all the compilation stuff...
%
% Feedback/problems: Christian Schaad, ingenieurbuero@christian-schaad.de

delete simo.mos;
delete([modelname,'.makefile']);
delete([modelname,'.log']);
delete([modelname,'.libs']);
delete([modelname,'.c']);
delete([modelname,'.h']);
delete([modelname,'.o']);
delete([modelname,'_functions.c']);
delete([modelname,'_functions.h']);
delete([modelname,'_records.c']);
delete([modelname,'_records.o']);

delete(['_',modelname,'.h']);

delete([modelname,'.cpp']);

author: Sable


text:
function outputImage = flip(this, iDim)
% Flips all elements of a given dimension and updates dimInfo.samplingPoints
%
%   Y = MrDataNd()
%   Y.flip(inputs)
%
% This is a method of class MrDataNd.
%
% IN
%   iDim    index of dimension to be flipped; default: 1
%
% OUT
%
% EXAMPLE
%   flip
%
%   See also MrDataNd perform_unary_operation MrDataNd.fliplr
%   See also MrDataNd.flipud

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-06-15
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 2
    iDim = 1;
end

if ~isnumeric(iDim)
    iDim = this.dimInfo.get_dim_index(iDim);
end

if exist('flip', 'builtin')
    outputImage = this.perform_unary_operation(@(x) flip(x, iDim));
else
    outputImage = this.perform_unary_operation(@(x) flipdim(x, iDim));
end
outputImage.dimInfo.samplingPoints{iDim} = flip(outputImage.dimInfo.samplingPoints{iDim});
end

author: translationalneuromodeling

repo: tapas


text:
%% To-do
% - Delete rhs_shapefunction -> use funs
% - Clean RHSunfitted
% OK! - switch class(mesh.type) @ RHSfactory
% - filter_pde_density using rhs_unfitted (now rhs_shapefunctionfun*
% - filter_pde_levelset name change for unfitted

%% Questions
% - Filter_P1_LevelSet.getP0fromP1() ??
% - Careful: unfittedmesh needed to properly integrate!!!

%% Results
%

%% Long-term
% - Use FeFunctions in TopOpt_Problem
% - PDE belongs to Optimizer, not ShapeFunctional
% - Micro as three elasticity problems

%% Backlog
% - Move Input folder to a separate repository
% - Geometry only in Mesh

% EXTRAS
%  - Investigate: converting data to binary format to save read'n'write
%                 resources for paraview
%  - Tutorial for printing
%  - Study file ouptut size vs time (GiD/Paraview) to see which is better
%    for printing (test + graph)
%  - Check XY component of fgaussfunctions
% - Recuperar gid unfitted mesh photo GiDimagecapturer
%      density --(project)--> unfittedmesh -> innermesh/photo

% - GiDImageCapturer -> some variables should be "user variables", defined
%   once in a centralized file

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [hr] = ns2hr(ns)
% Convert time from nanoseconds to hours.
% Chad Greene 2012
hr = ns*2.777777777778e-13 ;


text:
function coeffsRet=shrinkMultiDimPoly2Fit(coeffs)
%%SHRINKMULTIDIMPOLY2FIT Given a hypermatrix of coefficients of a
%                  multivariate polynomial, many of which might be zero,
%                  return the smallest hypermatrix possible that can
%                  represent the polynomial. This function effectively gets
%                  rid of high-order terms that are all zero.
%
%INPUTS: coeffs A hypermatrix of the coefficients for the multivariate
%               polynomial. These are arranged such that
%               coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%               x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%               number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs). Note that this ordering is the reverse
%               that used in the 1D polyval function that is built into
%               Matlab. The number of elements for each index in coeffs is
%               the maximum order of that dimension +1.
%
%OUTPUTS: coeffsRet The same multivariate polynomial as represented by
%               coeffs, but shrunk as small as possible.
%
%As an example, consider a 2D polynomial:
%440-288*x2+60*x2^2-4*x2^3-110*x1+72*x1*x2-15*x1*x2^2+x1*x2^3
% coeffs=zeros(5,5);
% coeffs(0+1,0+1)=440;
% coeffs(0+1,1+1)=-288;
% coeffs(0+1,2+1)=60;
% coeffs(0+1,3+1)=-4;
% coeffs(1+1,0+1)=-110;
% coeffs(1+1,1+1)=72;
% coeffs(1+1,2+1)=-15;
% coeffs(1+1,3+1)=1;
% coeffsRet=shrinkMultiDimPoly2Fit(coeffs)
%One will find that coeffsRet is only 2X4 in size, versus the 5X5 for the
%original matrix. Both represent the same 2D polynomial.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

dimSizeList=size(coeffs);
numIdx=length(dimSizeList);

newDimSizeList=zeros(numIdx,1);

%Go through to see which if any of the coefficients including other terms
%is nonzero.
idxVec(1:numIdx)={':'};
for curIdx=1:numIdx
    %idxVec(1:(curIdx-1))={':'};
    %idxVec((curIdx+1):end)={':'};

    maxIdx=1;
    for i=1:dimSizeList(curIdx)
        idxVec{curIdx}=i;
        %This type of element selection is faster than loops in Matlab, but
        %actually copying the elements into els just to use the any
        %function is very inefficient.
        els=coeffs(idxVec{:});
        if(any(els(:))~=0)
            maxIdx=i;
        end

    end
    idxVec{curIdx}=':';
    newDimSizeList(curIdx)=maxIdx;
end

%Allocate space for the reduced-size array to return and make it the
%correct size.
totalNewEls=prod(newDimSizeList);
coeffsRet=zeros(totalNewEls,1);
coeffsRet=reshape(coeffsRet,newDimSizeList(:)');
%Now, we must copy all of the elements in coeffs that are being kept into
%the proper positions in coeffsRet.

%This time, we need to use a loop and perform the copy one element at a
%time. The copying is performed by going through the linear indices of the
%destination matrix, translating them into the indices for the source
%matrix, one at a time.
for curEl=1:totalNewEls
    coeffsRet(curEl)=coeffs(nDim2Index(dimSizeList,index2NDim(newDimSizeList,curEl)));
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 r=risk(t,j,varargin)
%RISK Node risk.
%   R=RISK(T) returns an N-element vector R of the risk of the nodes in the
%   tree T, where N is the number of nodes.
%
%   R=RISK(T,J) takes an array J of node numbers and returns the risk
%   values for the specified nodes.
%
%   R=RISK(T,J,'criterion','error') returns risk vector R, where R(J) for
%   node J is the node error E(J) (classification error or mean squared
%   error for regression) weighted by the node probability P(J).
%
%   R=RISK(T,J,'criterion','impurity') computes risk by using the node
%   impurity measure for each node instead of the node error. This option
%   is only valid for classification trees grown using impurity measures
%   such as Gini index or deviance.
%
%   See also CLASSREGTREE, CLASSREGTREE/NODEERR, CLASSREGTREE/NODEPROB.

%   Copyright 2006-20097 The MathWorks, Inc.
%   $Revision: 1.1.8.4 $  $Date: 2009/05/07 18:32:51 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:risk:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

args = {'criterion'};
defs = {'error'};
[~,emsg,crit] = getargs(args,defs,varargin{:});
if ~isempty(emsg)
    error('stats:classregtree:risk:InvalidInput','Invalid input: %s',emsg);
end

if isempty(strmatch(lower(crit),'error')) && isempty(strmatch(lower(crit),'impurity'))
    error('stats:classregtree:risk:InvalidInput',...
        '''crit'' argument must be either ''error'' or ''impurity''.');
end

if strcmpi(crit,'error')
    r = t.nodeprob .* t.nodeerr;
else
    if isempty(t.impurity)
        error('stats:classregtree:risk:InvalidInput',...
        'Node risk cannot be computed using impurity. This is either a regression tree or splits were not found using impurity.');
    end
    r = t.nodeprob .* t.impurity;
end

if nargin>=2
    r = r(j,:);
end
end

author: Cloud-CV

repo: object-proposals


text:
function [C, S] = QLiftDec2MaxMin(X, N)
%-----------------------------------------------------------------------------
% QLiftDec2MaxMin
% Multilevel 2-D decomposition by the lifting scheme and using quincunx grids
%
% The MaxMin scheme has been proposed by Heijmans and Goutsias, see e.g.
%    H.J.A.M. Heijmans, J. Goutsias,
%    Multiresolution signal decomposition schemes.
%    Part 2: morphological wavelets.
%    CWI Report PNA-R9905, Amsterdam, 1999.
%    http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04625D.pdf
%
% Calls for: QLmaxlev,
%            storeQ1001, storeR,
%            getcolor01, getcolor10, getcolor00, getcolor11,
%            putcolor01, putcolor10, putcolor00, putcolor11.
% See also: QLiftRec2MaxMin
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: May 16, 2002.
% (c) 1999-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
%Firstly, check input data
%
if  isempty(X)
  error(' QLiftDec2MaxMin - empty matrix ');
else
  if mod(N, 2) == 1
    error(' QLiftDec2MaxMin - only an even number of levels is accepted ');
  end
  if QLmaxlev(size(X), 'maxmin') < N
    error(' QLiftDec2MaxMin - too many levels requested ');
  end
  if N < 2
    disp([' QLiftDec2MaxMin - WARNING too few levels requested ' ...
          '-> empty decomposition ']);
  end
end
%
%Secondly, start decomposition
%
O = X; % For the sake of efficient use of memory this could be improved upon.
% We descend to coarser grids, integer lev indicates number of scale.
C = []; S = [];

for lev=1:2:N
%
   [nO, mO] = size(O);
   if ( nO < 3 ) || ( mO < 3)
     error(' QLiftDec2MaxMin - too many levels ');
   end
   minO = min(min(O));
   maxO = max(max(O));
   cmin = minO-(maxO-minO);
   cmax = maxO+(maxO-minO);
%
%  The Lifting Scheme proceeds from a rectangular grid
%  towards a quincunx grid.
%
%  Stage: predict
   A00 =getcolor00(O);
   A11 =getcolor11(O);
%  Quincunx grid Q0011 is the union of the values at .00 and .11: "even slots"
%  Quincunx grid Q1001 is the union of the values at .10 and .01: "odd slots"
   Q1001D01 = getcolor01(O) - synA01max(A11, A00, cmin);                % Y1
   Q1001D10 = getcolor10(O) - synA10max(A11, A00, cmin);                % Y1
%  At this point the union (quincunx) of Q1001D01 & Q1001D10
%  contains the DETAILS of O.
%
%  For the inverse transform Q1001D01 and Q1001D10 have to be stored:
   [C, S] = storeQ1001( Q1001D10, Q1001D01, lev, 'd', C, S);
%
%  Stage: update
   Q0011A00 = A00 + ...
        max(zeros(size(A00)), synA00max(Q1001D10, Q1001D01, cmin));     % X1
   clear A00;
   Q0011A11 = A11 + ...
        max(zeros(size(A11)), synA11max(Q1001D10, Q1001D01, cmin));     % X1
   clear A11 Q1001D10 Q1001D01;
%  At this point the union (quincunx) of Q0011A00 & Q0011A11
%  contains the updated APPROXIMATION of O, the DETAILS of O
%  were in the union (quincunx) of Q1001D01 & Q1001D10 (see above).
%
%  The Lifting Scheme proceeds by a subsequent step from quincunx
%  to rectangular grid.
%
%  Q0011 is split into the 11 colour with the "odd slots" and
%  the 00 colour with the "even slots".
%
%  Stage: predict
   DETAIL11 = Q0011A11 - synA11Qmin(Q0011A00, size(Q0011A11), cmax);        % Y2
   clear Q0011A11;
%  Stage: update
   APPROX00 = Q0011A00 + ...
     min(zeros(size(Q0011A00)), synA00Qmin(DETAIL11, size(Q0011A00), cmax));% X2
%
%  DETAIL11 presents the detail gridfunction w.r.t. Q0011
%  APPROX00 now represents the updated version of the approximation of Q0011
%
%  For the inverse transform DETAIL11 has to be stored:
   [C, S] = storeR( DETAIL11, lev+1, 'd', C, S);
   clear Q0011A00 DETAIL11;
%
%  At this point gridfunction DETAIL11 containing the DETAILS has been stored,
%  gridfunction APPROX00 contains the updated APPROXIMATION, on the (down-
%  sampled) rectangular grid and has to be stored as well if at the highest
%  scale.
%  Note that APPROX00 is downsampled onto a rectangular grid with dimensions of
%  half size of the original O.
   if lev+1 >= N
     [C, S] = storeR(APPROX00, lev+1, 'a', C, S);
%    It is obligatory that at least at one scale the Approximation has to be
%    stored or else the scheme cannot be inverted.
     clear APPROX00;
%    In the Lifting Scheme all scales have now been processed!
   else
%    We proceed to the next scale.
     O = APPROX00; clear APPROX00;
   end
end
%-----------------------------------------------------------------------------

repo: mcbench-benchmarks


text:
clear;close all;
im1=imread('page.png');
im2=imread('tshape.png');
bwim1=adaptivethreshold(im1,11,0.03,0);
bwim2=adaptivethreshold(im2,15,0.02,0);
subplot(2,2,1);
imshow(im1);
subplot(2,2,2);
imshow(bwim1);
subplot(2,2,3);
imshow(im2);
subplot(2,2,4);
imshow(bwim2);