text:
function momentMat=quasiMoments2Moments(quasiMomentMat,muN,SigmaN)
%%QUASIMOMENTS2MOMENTS Convert a set of (multivariate) quasi-moments taken
%               with respect to a particular mean vector and covariance
%               matrix into noncentral moments. Quasi-moments are used in
%               generalized Edgeworth series approximations of
%               distributions.
%
%INPUTS: quasiMomentMat A matrix taking n indices, where
%               quasiMomentMat(a1,a2,a3...an) corresponds to the quasi-
%               moment whose multivariate order is given by a1-1,a2-1, etc.
%               It is assumed that quasi moments only up to an order equal
%               to the size of the first dimension of the matrix are
%               available (other entries in quasiMomentMat are ignored).
%               All other dimensions must be at least the same size as the
%               first. Also, the zero-th order moment must be 1 (as the PDF
%               integrates to 1). If a univariate distribution is used,
%               then this is a column vector.
%           muN The nX1 mean vector with respect to with the quasi-moments
%               are computed. If this parameter is omitted or an empty
%               matrix is passed, a zero mean vector is used.
%        SigmaN The nXn covariance matrix with respect to which the
%               quasi-moments are computed. If this parameter is omitted or
%               an empty matrix is passed, the identity matrix is used.
%
%OUTPUTS: momentMat A matrix taking n indices, where
%               momentMat(a1,a2,a3...an) corresponds to the coefficient of
%               a E(x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1)) moment.
%               The order of the moment is a1+a2+..an.
%
%Quasi-moments are introduced in [2] and are defined in terms of a certain
%mean vector and covariance matrix. If mean vector and covariance matrix
%coincide with those of the moments associated with the quasi-moments, then
%the quasi-moments are called proper quasi-moments.
%
%The conversion is based on Equation 2.11 in [1]. The equation is
%m_\alpha=\sum_{\beta}^{\alpha}C_{\alpha}^\beta m^N_{\alpha-\beta} q_\beta
%where q's represent quasi moments, m^N's represent moments of the Normal
%distribution having covariance matrix Sigma, m's are moments of the
%distribution that are to be converted. The indexation is vector based.
%Thus, beta is actually a vector of beta(1)...beta(n) and alpha goes from
%alpha(1)...alpha(n). The sum is actually multiple sums. The
%C_{\alpha}^\beta} is equal to the product of binomial(alpha(i),beta(i))
%over i.
%
%REFERENCES:
%[1] I. E. Poloskov, "CAS Mathematica in random studies," in Proceedings
%    of Computational Science - ICCS 2003, Melbourne, Australia and St.
%    Petersburg, Russia, 2-4 Jun. 2003, pp. 781-790.
%[2] P. I. Kuznetsov, R. L. Stratonovich, and V. I. Tikhonov, "Quasi-moment
%    functions in the theory of random processes," Theory of Probability
%    and its Applications, vol. V, no. 1, pp. 80-97, 1960.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDimList=size(quasiMomentMat);
maxDeg=numDimList(1)-1;
numIdx=length(numDimList);

%If the last dimension is just a singleton dimension, then shrink
%everything by one.
if(numIdx==2&&numDimList(2)==1&&(nargin<3||size(SigmaN,2)==1))
    numIdx=1;
end

%Do the conversion with respect to the standard normal distribution if no
%further details are given.
if(nargin<2||isempty(muN))
    muN=zeros(numIdx,1);
end
if(nargin<3||isempty(SigmaN))
    SigmaN=eye(numIdx,numIdx);
end

%Allocate space to hold a table of binomial values that can be quickly
%looked up in the loops rather than having to compute them again and again.
binomTable=makeBinomTable(maxDeg);

%This is a table that will hold values of moments for the normal
%distribution with zero mean and covariance matrix Sigma. These moments
%play a role in the conversion to quasi-moments as they are related to
%hermite polynomials.
momentMatN=zeros(numDimList);
momentMatN(1)=1;%Zeroth-order moment is always 1.
for curDeg=1:maxDeg
    curPart=getNextMPartition(numIdx+curDeg,numIdx);

    while(~isempty(curPart))
        %Go through all permutation of the orders of each index.
        powTable=genAllMultisetPermutations(curPart-1);
        numVals=size(powTable,2);
        for curVal=1:numVals
            numDerivs=powTable(:,curVal);
            momentVal=GaussianD.momentGenFun(muN,SigmaN,numDerivs);
            curEl=nDim2Index(numDimList,numDerivs+1);
            momentMatN(curEl)=momentVal;
        end

        curPart=getNextMPartition(curPart);
    end
end

%Allocate space for the return values.
momentMat=zeros(numDimList);
%The zeroth-order moment is always 1.
momentMat(1)=1;

%First, fill in the values for the case where only the first dimension has
%nonzero exponents. This is just Equation 11 for the scalar case.
for n=1:maxDeg
    for i=0:n
        momentMat(n+1)=momentMat(n+1)+binomTable(n+1,i+1)*momentMatN(n-i+1)*quasiMomentMat(i+1);
    end
end

for numNonzeroIdx=2:numIdx
    %For a given new index that will be nonzero.
    %We start with the last index (not visited) being nonzero, because we
    %have already considered all of the instances where it is zero.
    for nonzeroDegVal=1:maxDeg
        %curDeg is the degree of all elements except the last one, which is
        %nonzero and which we go through in the outer loop.
        for curDeg=0:(maxDeg-nonzeroDegVal)
            %Fixing the order of highest index, we must now go through all
            %partitions of possible values that the other elements can
            %take.
            curPart=getNextMPartition(numNonzeroIdx-1+curDeg,numNonzeroIdx-1);
            while(~isempty(curPart))
                %curPart-1 is the (unordered) degrees of the indices before
                %the current index (the last nonzero index).
                %For the given partition, we must go through all possible
                %multiset permutations of the values.
                lowerIdxPerms=genAllMultisetPermutations(curPart);
                numPerms=size(lowerIdxPerms,2);

                for curPerm=1:numPerms
                    r=[lowerIdxPerms(:,curPerm)-1;nonzeroDegVal];
                    idxVal=nDim2Index(numDimList(1:numNonzeroIdx),r+1);

                    %Use Equation 2.11 to compute the moment term.
                    momentMat(idxVal)=equation211Update(r,momentMatN,quasiMomentMat,binomTable);
                end
                curPart=getNextMPartition(curPart);
            end
        end
    end
end
end

function sumVal=equation211Update(alpha,momentMatN,quasiMomentMat,binomTable)
%This implements Equation 2.11 in [1] for a given alpha vector.

numIdx=length(alpha);
maxDim=size(quasiMomentMat);
maxDim=maxDim(1:numIdx);%Only the first numIdx indices are considered.

%This stores cumulative values of the product of the binomials going down
%each sum. This way, the values do not have to be constantly recomputed.
binomProdVal=ones(numIdx-1,1);

%This vector holds all of the i values. Initially, all sums start at zero.
iVec=zeros(numIdx,1);

%The loop goes through all of the values in the sums.
curLevel=numIdx;
isAscending=false;
sumVal=0;
while(curLevel>=1)
    if(curLevel==numIdx)
        %Go through the innermost sum
        for i=0:alpha(numIdx)
            binomVal=binomProdVal(numIdx-1)*binomTable(alpha(numIdx)+1,i+1);

            iVec(numIdx)=i;
            kIdxVec=alpha-iVec;

            idx=nDim2Index(maxDim,iVec+1);
            kidx=nDim2Index(maxDim,kIdxVec+1);

            sumVal=sumVal+binomVal*momentMatN(kidx)*quasiMomentMat(idx);
        end
        curLevel=curLevel-1;
        isAscending=true;
        continue;
    elseif(isAscending)
        %Try incrementing the order at this level...
        iVec(curLevel)=iVec(curLevel)+1;
        if(iVec(curLevel)<=alpha(curLevel))
            %If the value is valid, then update the binomial product value
            %and descend to the next level.
            if(curLevel~=1)
                binomProdVal(curLevel)=binomProdVal(curLevel-1)*binomTable(alpha(curLevel)+1,iVec(curLevel)+1);
            else
                binomProdVal(curLevel)=binomTable(alpha(curLevel)+1,iVec(curLevel)+1);
            end
            isAscending=false;
            curLevel=curLevel+1;
            continue;
        else
            %If the value is invalid, then just keep ascending.
            curLevel=curLevel-1;
            continue;
        end
    else%We are descending in the sums here.
        iVec(curLevel)=0;
        binomProdVal(curLevel)=binomProdVal(curLevel-1);
        curLevel=curLevel+1;
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function pasteLossConvertFromClipboard(hFigWithLossConvertAppData)
%pasteLossConvertFromClipboard Summary of this function goes here
%   Detailed explanation goes here
    global GLOBAL_LossConvertClipboard;

    if(~isempty(GLOBAL_LossConvertClipboard))
        setappdata(hFigWithLossConvertAppData,'lossConverts',GLOBAL_LossConvertClipboard);
    else
        GLOBAL_LossConvertClipboard = [];
        warndlg('Could not paste Loss/Convert settings.  Clipboard contents are either empty or invalid.  Try copying mass loss/conversion settings to the clipboard first.','Could not paste loss/convert settings.');
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
clc, close all, clear all

load('teapot'); %loading matrix S
% % Matrix S stores all the control points of all the patches of
% % teapot surface such that
% % S(:,:,:,k) control points of kth patch, where k=1..32
% % Size of S(:,:,:,k) is 4 x 4 x 3, i.e., 16 control points and each
% % control point has three values (x,y,z)

% % S(:,:,1,k): x-coordates of control points of kth patch as 4 x 4 matrix
% % S(:,:,2,k): y-coordates of control points of kth patch as 4 x 4 matrix
% % S(:,:,3,k): z-coordates of control points of kth patch as 4 x 4 matrix
% % ------------------------------------
[r c d np]=size(S);
% % np: number of patches
ni=20; %number of interpolated values between end control points
u=linspace(0,1,ni); v=u;  %uniform parameterization
% % Higher the value of ni smoother the surface but computationally
% % expensive
% % ------------------------------------
% % Cubic Bezier interpolation of control points of each patch
for k=1:np
    Q(:,:,:,k)=bezierpatchinterp(S(:,:,:,k),u,v); %interpolation of kth patch
end
% % ------------------------------------
% % Plotting a signle Bezier Patch in many ways
k=13; %ploting kth patch
plotbezierpatch3D(S(:,:,:,k),Q(:,:,:,k))

% % Plotting Bezier surface (all the patches) in many ways
plotbeziersurface3D(S,Q)
% % --------------------------------
% % This program or any other program(s) supplied with it does not provide any
% % warranty direct or implied.
% % This program is free to use/share for non-commerical purpose only.
% % Kindly reference the author.
% % Author: Dr. Murtaza Khan
% % URL : http://www.linkedin.com/pub/dr-murtaza-khan/19/680/3b3
% % --------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_div( )
% Test divergence

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Check divergence of an empty diskfunv is an empty diskfun.
u = diskfunv;
f = div(u);
pass(1) = isempty(f) & isa(f,'diskfun');

% Check definition:
F = diskfunv(@(x,y) cos(x), @(x,y) sin(y));
divF = diff(F(1), 1,1) + diff(F(2), 2, 1);
pass(2) = ( norm(divF - divergence(F) ) < tol );

% Check divergence of the zero field is zero
f = diskfun(@(x,y) 0*x);
u = diskfunv(f,f);
pass(3) = norm(div(u)) < tol;

% Check divergence of gradient is laplacian:
f = diskfun(@(x,y) cos(x.*y));
pass(4) = ( norm(laplacian( f ) - divergence( gradient( f ) ) ) < tol );

% Check that the div and divergence give the same result.
divergenceu = divergence(u);
divu = div(u);
pass(5) = isequal(divu,divergenceu);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cheby_u_poly_coef_test ( )

%*****************************************************************************80
%
%% CHEBY_U_POLY_COEF_TEST tests CHEBY_U_POLY_COEF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_U_POLY_COEF_TEST\n' );
  fprintf ( 1, '  CHEBY_U_POLY_COEF determines the Chebyshev U \n' );
  fprintf ( 1, '  polynomial coefficients.\n' );

  c = cheby_u_poly_coef ( n );

  for i = 0 : n
    fprintf ( 1, '\n' );
    fprintf ( 1, '  U(%d)\n', i );
    fprintf ( 1, '\n' );
    for j = i : -1 : 0
      if ( j == 0 )
        fprintf ( 1, '    %f\n', c(i+1,j+1) );
      elseif ( j == 1 )
        fprintf ( 1, '    %f * x\n', c(i+1,j+1) );
      else
        fprintf ( 1, '    %f * x^%d\n', c(i+1,j+1), j );
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = cvxdual( prob, name )

y = class( struct( 'problem_', prob, 'name_', name, 'attached_', false ), 'cvxdual', cvxobj );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% pop_biosig() - import data files into EEGLAB using BIOSIG toolbox
%
% Usage:
%   >> OUTEEG = pop_biosig; % pop up window
%   >> OUTEEG = pop_biosig( filename, channels, type);
%
% Inputs:
%   filename - [string] file name
%
% Optional inputs:
%   'channels'   - [integer array] list of channel indices
%   'blockrange' - [min max] integer range of data blocks to import, in seconds.
%                  Entering [0 3] will import the first three blocks of data.
%                  Default is empty -> import all data blocks.
%  'importevent' - ['on'|'off'] import events. Default if 'on'.
%  'importannot' - ['on'|'off'] import annotations (EDF+ only). Default if 'on'
%  'blockepoch'  - ['on'|'off'] force importing continuous data. Default is 'on'
%  'ref'         - [integer] channel index or index(s) for the reference.
%                  Reference channels are not removed from the data,
%                  allowing easy re-referencing. If more than one
%                  channel, data are referenced to the average of the
%                  indexed channels. WARNING! Biosemi Active II data
%                  are recorded reference-free, but LOSE 40 dB of SNR
%                  if no reference is used!. If you do not know which
%                  channel to use, pick one and then re-reference after
%                  the channel locations are read in. {default: none}.
%                  For more information see http://www.biosemi.com/faq/cms&drl.htm
%  'refoptions'  - [Cell] Option for the pop_reref function. Default is to
%                  remove the reference channel if there is one of them and to
%                  keep it if there are several of them from the graphic
%                  interface. From the command line default option is to
%                  keep the reference channel.
%  'rmeventchan' - ['on'|'off'] remove event channel after event
%                  extraction. Default is 'on'.
%  'memorymapped' - ['on'|'off'] import memory mapped file (useful if
%                  encountering memory errors). Default is 'off'.
%
% Outputs:
%   OUTEEG   - EEGLAB data structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Oct. 29, 2003-
%
% Note: BIOSIG toolbox must be installed. Download BIOSIG at
%       http://biosig.sourceforge.net
%       Contact a.schloegl@ieee.org for troubleshooting using BIOSIG.

% Copyright (C) 2003 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, command, dat] = pop_biosig(filename, varargin);
EEG = [];
command = '';

if ~plugin_askinstall('Biosig', 'sopen'), return; end;
biosigpathfirst;

if nargin < 1
	% ask user
	[filename, filepath] = uigetfile('*.*', 'Choose a data file -- pop_biosig()'); %%% this is incorrect in original version!!!!!!!!!!!!!!
    drawnow;

	if filename == 0 return; end;
	filename = [filepath filename];

    % look if MEG
    % -----------
    if length(filepath)>4
        if strcmpi(filepath(end-3:end-1), '.ds'), filename = filepath(1:end-1); end;
    end;

    % open file to get infos
    % ----------------------
    disp('Reading data file header...');
    dat = sopen(filename, 'r', [], 'OVERFLOWDETECTION:OFF');
    if ~isfield(dat, 'NRec')
        error('Unsuported data format');
    end;

    % special BIOSEMI
    % ---------------
    eeglab_options;
    if strcmpi(dat.TYPE, 'BDF')
        disp(upper('We highly recommend that you choose a reference channel IF these are Biosemi data'));
        disp(upper('(e.g., a mastoid or other channel). Otherwise the data will lose 40 dB of SNR!'));
        disp('For more information, see http://www.biosemi.com/faq/cms&drl.htm');
    end;
    uilist = { { 'style' 'text' 'String' 'Channel list (defaut all):' } ...
                 { 'style' 'edit' 'string' '' } ...
                 { 'style' 'text' 'String' [ 'Data range (in seconds) to read (default all [0 ' int2str(dat.NRec) '])' ] } ...
                 { 'style' 'edit' 'string' '' } ...
                 { 'style' 'text' 'String' 'Extract event' } ...
                 { 'style' 'checkbox' 'string' '' 'value' 1 'enable' 'on' } {} ...
                 { 'style' 'text' 'String' 'Import anotations (EDF+ only)' } ...
                 { 'style' 'checkbox' 'string' '' 'value' 1 'enable' 'on' } {} ...
                 { 'style' 'text' 'String' 'Force importing continuous data' 'value' 1} ...
                 { 'style' 'checkbox' 'string' '' 'value' 0 } {} ...
                 { 'style' 'text' 'String' 'Reference chan(s) indices - required for BIOSEMI' } ...
                 { 'style' 'edit' 'string' '' } ...
                  { 'style' 'checkbox' 'String' 'Import as memory mapped file (use if out of memory error)' 'value' option_memmapdata } };
    geom = { [3 1] [3 1] [3 0.35 0.5] [3 0.35 0.5] [3 0.35 0.5] [3 1] [1] };

    result = inputgui( geom, uilist, 'pophelp(''pop_biosig'')', ...
                                 'Load data using BIOSIG -- pop_biosig()');
    if length(result) == 0 return; end;

    % decode GUI params
    % -----------------
    options = {};
    if ~isempty(result{1}), options = { options{:} 'channels'   eval( [ '[' result{1} ']' ] ) }; end;
    if ~isempty(result{2}), options = { options{:} 'blockrange' eval( [ '[' result{2} ']' ] ) }; end;
    if length(result) > 2
        if ~result{3},          options = { options{:} 'importevent'  'off'  }; end;
        if ~result{4},          options = { options{:} 'importannot'  'off'  }; end;
        if  result{5},          options = { options{:} 'blockepoch'   'off' }; end;
        if ~isempty(result{6}), options = { options{:} 'ref'        eval( [ '[' result{6} ']' ] ) }; end;
        if  result{7},          options = { options{:} 'memorymapped' 'on' }; end;
    end;
    if length(eval( [ '[' result{6} ']' ] )) > 1
        options = { options{:} 'refoptions' { 'keepref' 'off' } };
    end;
else
    options = varargin;
end;

% decode imput parameters
% -----------------------
g = finputcheck( options, { 'blockrange'   'integer' [0 Inf]    [];
                            'channels'     'integer' [0 Inf]    [];
                            'ref'          'integer' [0 Inf]    [];
                            'refoptions'   'cell'    {}             { 'keepref' 'on' };
                            'rmeventchan'  'string'  { 'on';'off' } 'on';
                            'importevent'  'string'  { 'on';'off' } 'on';
                            'importannot'  'string'  { 'on';'off' } 'on';
                            'memorymapped' 'string'  { 'on';'off' } 'off';
                            'blockepoch'   'string'  { 'on';'off' } 'off' }, 'pop_biosig');
if isstr(g), error(g); end;

% import data
% -----------
EEG = eeg_emptyset;
[dat DAT interval] = readfile(filename, g.channels, g.blockrange, g.memorymapped);

if strcmpi(g.blockepoch, 'off')
    dat.NRec = 1;
end;

EEG = biosig2eeglab(dat, DAT, interval, g.channels, strcmpi(g.importevent, 'on'));

if strcmpi(g.rmeventchan, 'on') & strcmpi(dat.TYPE, 'BDF') & isfield(dat, 'BDF')
    if size(EEG.data,1) >= dat.BDF.Status.Channel,
        disp('Removing event channel...');
        EEG.data(dat.BDF.Status.Channel,:) = [];
        if ~isempty(EEG.chanlocs) && length(EEG.chanlocs) >= dat.BDF.Status.Channel
            EEG.chanlocs(dat.BDF.Status.Channel) = [];
        end;
    end;
    EEG.nbchan = size(EEG.data,1);
end;

% rerefencing
% -----------
if ~isempty(g.ref)
    disp('Re-referencing...');
    EEG = pop_reref(EEG, g.ref, g.refoptions{:});
%     EEG.data = EEG.data - repmat(mean(EEG.data(g.ref,:),1), [size(EEG.data,1) 1]);
%     if length(g.ref) == size(EEG.data,1)
%         EEG.ref  = 'averef';
%     end;
%     if length(g.ref) == 1
%         disp([ 'Warning: channel ' int2str(g.ref) ' is now zeroed (but still present in the data)' ]);
%     else
%         disp([ 'Warning: data matrix rank has decreased through re-referencing' ]);
%     end;
end;

% test if annotation channel is present
% -------------------------------------
if isfield(dat, 'EDFplus') && strcmpi(g.importannot, 'on')
    tmpfields = fieldnames(dat.EDFplus);
    for ind = 1:length(tmpfields)
        tmpdat = getfield(dat.EDFplus, tmpfields{ind});
        if length(tmpdat) == EEG.pnts
            EEG.data(end+1,:) = tmpdat;
            EEG.nbchan        = EEG.nbchan+1;
            if ~isempty(EEG.chanlocs)
                EEG.chanlocs(end+1).labels = tmpfields{ind};
            end;
        end;
    end;
end;

% convert data to single if necessary
% -----------------------------------
EEG = eeg_checkset(EEG,'makeur');   % Make EEG.urevent field
EEG = eeg_checkset(EEG);

% history
% -------
if isempty(options)
    command = sprintf('EEG = pop_biosig(''%s'');', filename);
else
    command = sprintf('EEG = pop_biosig(''%s'', %s);', filename, vararg2str(options));
end;

% Checking if str2double is on top of the path
biosigpathlast;

% ---------
% read data
% ---------
function [dat DAT interval] = readfile(filename, channels, blockrange, memmapdata);

if isempty(channels), channels = 0; end;
dat = sopen(filename, 'r', channels,'OVERFLOWDETECTION:OFF');

if strcmpi(memmapdata, 'off')
    fprintf('Reading data in %s format...\n', dat.TYPE);

    if ~isempty(blockrange)
        newblockrange    = blockrange;
%         newblockrange    = newblockrange*dat.Dur;
        DAT=sread(dat, newblockrange(2)-newblockrange(1), newblockrange(1));
    else
        DAT=sread(dat, Inf);% this isn't transposed in original!!!!!!!!
        newblockrange    = [];
    end
    sclose(dat);
else
    fprintf('Reading data in %s format (file will be mapped to memory so this may take a while)...\n', dat.TYPE);
    inc = ceil(250000/(dat.NS*dat.SPR)); % 1Mb block

    if isempty(blockrange), blockrange = [0 dat.NRec]; end;
    blockrange(2) = min(blockrange(2), dat.NRec);
    allblocks = [blockrange(1):inc:blockrange(end)];
    count = 1;
    for bind = 1:length(allblocks)-1
        TMPDAT=sread(dat, (allblocks(bind+1)-allblocks(bind))*dat.Dur, allblocks(bind)*dat.Dur);
        if bind == 1
            DAT = mmo([], [size(TMPDAT,2) (allblocks(end)-allblocks(1))*dat.SPR]);
        end;
        DAT(:,count:count+length(TMPDAT)-1) = TMPDAT';
        count = count+length(TMPDAT);
    end;
    sclose(dat);
end;

if ~isempty(blockrange)
     interval(1) = blockrange(1) * dat.SampleRate(1) + 1;
     interval(2) = blockrange(2) * dat.SampleRate(1);
else interval = [];
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [D,L] = spm_opm_create(S)
% Read or simulate magnetometer data and optionally set up forward model
% FORMAT D = spm_opm_create(S)
%   S               - input structure
% Optional fields of S:
% SENSOR LEVEL INFO
%   S.data          - filepath to .bin file        - Default: Simulates data
%   S.channels      - channels.tsv file            - Default: REQUIRED
%   S.fs            - Sampling frequency (Hz)      - Default: REQUIRED if S.meg is empty
%   S.meg           - meg.json file                - Default: REQUIRED if S.fs is empty
%   S.precision     - 'single' or 'double'         - Default: 'single'
% SIMULATION
%   S.wholehead     - whole head coverage flag     - Deafult: 0
%   S.space         - space between sensors(mm)    - Default: 25
%   S.offset        - scalp to sensor distance(mm) - Default: 6.5
%   S.nSamples      - number of samples            - Default: 1000
%   S.Dens          - number of density checks     - Default: 40

% SOURCE LEVEL INFO
%   S.coordsystem   - coordsystem.json file        - Default:
%   S.positions     - positions.tsv file           - Default:
%   S.sMRI          - Filepath to  MRI file        - Default: uses template
%   S.cortex        - Custom cortical mesh         - Default: Use inverse normalised cortical mesh
%   S.scalp         - Custom scalp mesh            - Default: Use inverse normalised scalp mesh
%   S.oskull        - Custom outer skull mesh      - Default: Use inverse normalised outer skull mesh
%   S.iskull        - Custom inner skull mesh      - Default: Use inverse normalised inner skull mesh
%   S.voltype       - Volume conducter Model type  - Default: 'Single Shell'
%   S.meshres       - mesh resolution(1,2,3)       - Default: 1
%   S.lead          - flag to compute lead field   - Default: 0
% Output:
%  D           - MEEG object (also written to disk)
%  L           - Lead field (also written on disk)
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Tim Tierney
% $Id: spm_opm_create.m 7645 2019-07-25 13:58:25Z tim $
spm('FnBanner', mfilename);

%-Set default values
%--------------------------------------------------------------------------
if ~isfield(S, 'voltype'),     S.voltype = 'Single Shell'; end
if ~isfield(S, 'meshres'),     S.meshres = 1; end
if ~isfield(S, 'scalp'),       S.scalp = []; end
if ~isfield(S, 'cortex'),      S.cortex = []; end
if ~isfield(S, 'iskull'),      S.iskull = []; end
if ~isfield(S, 'oskull'),      S.oskull = []; end
if ~isfield(S, 'lead'),        S.lead = 0; end
if ~isfield(S, 'fs'),          S.fs   = 1000; end
if ~isfield(S, 'nSamples'),    S.nSamples   = 1000; end
if ~isfield(S, 'nDens'),       S.nDens   = 40; end
if ~isfield(S, 'space'),       S.space  = 30; end
if ~isfield(S, 'offset'),      S.offset  = 6.5; end
if ~isfield(S, 'data'),        S.data = zeros(1,S.nSamples); end
if ~isfield(S, 'wholehead'),   S.wholehead = 1; end
if ~isfield(S, 'fname'),       S.fname = 'sim_opm'; end
if ~isfield(S, 'precision'),   S.precision = 'single'; end

%- Read Binary File
%----------------------------------------------------------------------
try % to read data
    [direc, dataFile] = fileparts(S.data);
    dat = fopen(S.data);
    S.data = fread(dat,Inf,S.precision,0,'b');
    fclose(dat);
    binData=1;
catch % if not readable check if it is numeric
    if ~isa(S.data,'numeric') % if not numeric throw error
        error('A valid dataest or file was not supplied')
    end
    binData=0;
    direc = pwd();
    dataFile=S.fname;
end
%- identify potential BIDS Files
%----------------------------------------------------------------------
base = strsplit(dataFile,'meg');
chanFile= fullfile(direc,[base{1},'channels.tsv']);
megFile= fullfile(direc,[base{1},'meg.json']);
posFile= spm_select('FPList',direc,[base{1},'positions.tsv']);
coordFile= fullfile(direc,[base{1},'coordsystem.json']);

%- Check for channel Info
%--------------------------------------------------------------------------
try % to load a channels file
    channels = spm_load(S.channels);
catch
    try % to load a BIDS channel file
       channels = spm_load(chanFile);
    catch
        try  % use channel struct if supplied
           channels = S.channels;
        catch % create channel struct
            args=[];
            args.base='Chan';
            args.n= size(S.data,1);
            labs= spm_create_labels(args);
            channels = [];
            channels.name=labs;
            channels.type=repmat({'MEG'},size(S.data,1),1);
            channels.units= repmat({'fT'},size(S.data,1),1);
        end
    end
end

%- reformat data according to channel info
%--------------------------------------------------------------------------
nc = size(channels.name,1);

if binData
    S.data = reshape(S.data,nc,numel(S.data)/nc);
elseif nc~=size(S.data,1)
    error('numer of channels in S.data different to S.channels')
else
    S.data =S.data;
end

%- Check for MEG Info
%--------------------------------------------------------------------------
try % to load a meg file
    meg = spm_load(S.meg);
catch
    try % to load a BIDS meg file
       meg = spm_load(megFile);
    catch
        try % to use meg struct
           meg = S.meg;
        catch
            try % to use S.fs argument to get sampling frequency
                meg =[];
                meg.SamplingFrequency=S.fs;
            catch
                ('A valid meg.json file is required if S.fs is empty');
            end
        end
    end
end

%- Position File check
%----------------------------------------------------------------------
try % to load a channels file
    posOri = spm_load(S.positions);
    positions =1;
catch
    try % to load a BIDS channel file
       posOri = spm_load(posFile);
       positions =1;
    catch
       positions =0;
    end
end

%- Forward model Check
%----------------------------------------------------------------------
subjectSource  = (positions|isfield(S,'space')) & isfield(S,'sMRI');
subjectSensor = ~subjectSource;

if subjectSource
    forward =1;
    template =0;
elseif subjectSensor
    forward =0;
    template =0;
else
    forward =1;
    template =1;
end

%- Create SPM object of simulated or real data
%--------------------------------------------------------------------------
Dtemp = meeg(size(S.data,1),size(S.data,2),size(S.data,3));
Dtemp = fsample(Dtemp,meg.SamplingFrequency);
Dtemp = fname(Dtemp,[dataFile,'.mat']);
Dtemp = path(Dtemp,direc);
Dtemp = chanlabels(Dtemp,1:size(Dtemp,1),channels.name);
Dtemp = units(Dtemp,1:size(Dtemp,1),channels.units);
Dtemp = chantype(Dtemp,1:size(Dtemp,1),channels.type);

%- Overwrite and Save
%--------------------------------------------------------------------------
ma = fullfile(direc,[dataFile,'.mat']);
da = fullfile(direc,[dataFile,'.dat']);

ae = exist(fname(Dtemp),'file')==2;
if(ae)
    delete(ma);
    delete(da);
end
Dtemp.save();
% create data file and insert data
D= blank(Dtemp,[dataFile,'.dat']);
dim=size(D);
D(1:dim(1),1:dim(2),1:dim(3)) = S.data;
D.save();

%- Create Meshes
%--------------------------------------------------------------------------
if forward
    %initially used inverse normalised meshes
    D = spm_eeg_inv_mesh_ui(D,1,S.sMRI,S.meshres);
    save(D);
    % then fill in custom meshes(if they exist)
    args = [];
    args.D = D;
    args.scalp = S.scalp;
    args.cortex = S.cortex;
    args.iskull = S.iskull;
    args.oskull = S.oskull;
    args.template = template;
    D = opm_customMeshes(args);
    save(D);
end

%- Create the Sensor Array
%--------------------------------------------------------------------------
if forward
    try % create positions and orientations
        pos = [posOri.Px,posOri.Py,posOri.Pz];
        ori = [posOri.Ox,posOri.Oy,posOri.Oz];
        cl = posOri.name;
    catch % if no postions and orientations provided then create them
        args = [];
        args.D =D;
        args.offset = S.offset;
        args.space = S.space;
        args.wholehead = S.wholehead;
        args.nDens = S.nDens;
        [pos,ori] = opm_createSensorArray(args);
    end

    nSensors=size(pos,1);
    if nSensors>size(S.data,1) %
        args=[];
        args.base='Chan';
        args.n= nSensors;
        labs= spm_create_labels(args);
        channels = [];
        channels.name=labs;
        channels.type=repmat({'MEG'},nSensors,1);
        channels.units= repmat({'fT'},nSensors,1);
        D = clone(D,fnamedat(D),[nSensors,S.nSamples,1],1);
        D = chanlabels(D,1:size(D,1),channels.name);
        D = units(D,1:size(D,1),channels.units);
        D = chantype(D,1:size(D,1),channels.type);
        cl=chanlabels(D)';
    end
end

%-Place Sensors  in object
%--------------------------------------------------------------------------
if forward
    grad= [];
    grad.label = cl;
    grad.coilpos = pos;
    grad.coilori = ori;
    grad.tra = eye(numel(grad.label));
    grad.chanunit = repmat({'T'}, numel(grad.label), 1);
    grad.chantype= 'MEG';
    grad = ft_datatype_sens(grad, 'amplitude', 'T', 'distance', 'mm');
    D = sensors(D, 'MEG', grad);
    save(D);
end

%- fiducial settings
%--------------------------------------------------------------------------
if subjectSource
    miMat = zeros(3,3);
    fiMat = zeros(3,3);
    fid=[];
    try % to read the coordsystem.json
        coord = spm_load(S.coordsystem);
        fiMat(1,:) = coord.HeadCoilCoordinates.coil1;
        fiMat(2,:) = coord.HeadCoilCoordinates.coil2;
        fiMat(3,:) = coord.HeadCoilCoordinates.coil3;
        miMat(1,:) = coord.AnatomicalLandmarkCoordinates.coil1;
        miMat(2,:) = coord.AnatomicalLandmarkCoordinates.coil2;
        miMat(3,:) = coord.AnatomicalLandmarkCoordinates.coil3;
        fid.fid.label = fieldnames(coord.HeadCoilCoordinates);
        fid.fid.pnt =fiMat;
        fid.pos= []; % headshape field that is left blank (GRB)
        M = fid;
        M.fid.pnt=miMat;
        M.pnt = D.inv{1}.mesh.fid.pnt;
    catch
        try % to find the BIDS coordsystem.json
            coord = spm_load(coordFile);
            fiMat(1,:) = coord.HeadCoilCoordinates.coil1;
            fiMat(2,:) = coord.HeadCoilCoordinates.coil2;
            fiMat(3,:) = coord.HeadCoilCoordinates.coil3;
            miMat(1,:) = coord.AnatomicalLandmarkCoordinates.coil1;
            miMat(2,:) = coord.AnatomicalLandmarkCoordinates.coil2;
            miMat(3,:) = coord.AnatomicalLandmarkCoordinates.coil3;
            fid.fid.label = fieldnames(coord.HeadCoilCoordinates);
            fid.fid.pnt =fiMat;
            fid.pos= []; % headshape field that is left blank (GRB)
            M = fid;
            M.fid.pnt=miMat;
            M.pnt = D.inv{1}.mesh.fid.pnt;
        catch % DEFAULT: transform between fiducials and anatomy is identity
            fid.fid.label = {'nas', 'lpa', 'rpa'}';
            fid.fid.pnt = [0 0 0; -1 0 0; 1 0 0];
            fid.pos= [];
            M = fid;
            M.pnt = D.inv{1}.mesh.fid.pnt;
        end
    end
end

if(template) %make
    fid.fid.label = {'nas', 'lpa', 'rpa'}';
    fid.fid.pnt = D.inv{1}.mesh.fid.fid.pnt(1:3,:);
    fid.pos= []; % headshape field that is left blank (GRB)
    M = fid;
    M.pnt = D.inv{1}.mesh.fid.pnt;
end

%- Coregistration
%--------------------------------------------------------------------------
if(subjectSource)
    D = fiducials(D, fid);
    save(D);
    f=fiducials(D);
    f.pnt =zeros(0,3);
    D = spm_eeg_inv_datareg_ui(D,1,f,M,0);
end
%- 2D view based on mean orientation of sensors
%--------------------------------------------------------------------------
if(forward)
    n1=mean(grad.coilori); n1= n1./sqrt(dot(n1,n1));
    t1=cross(n1,[0 0 1]);
    t2=cross(t1,n1);
    pos2d =zeros(size(grad.coilpos,1),2);
    for i=1:size(grad.coilpos,1)
        pos2d(i,1)=dot(grad.coilpos(i,:),t1);
        pos2d(i,2)=dot(grad.coilpos(i,:),t2);
    end

 nMEG = length(indchantype(D,'MEG'));
    if nMEG~=size(pos2d,1)
        m1 = '2D positions could not be set as there are ';
        m2 =num2str(nMEG);
        m3 = ' channels but only ';
        m4 = num2str(size(pos2d,1));
        m5 =  ' channels with position information.';
        message = [m1,m2,m3,m4,m5];
        warning(message);
    else
        args=[];
        args.D=D;
        args.xy= pos2d';
        args.label=grad.label;
        args.task='setcoor2d';
        D=spm_eeg_prep(args);
        D.save;
    end
end
%- Foward  model specification
%--------------------------------------------------------------------------
if forward
    D.inv{1}.forward.voltype = S.voltype;
    D = spm_eeg_inv_forward(D);
    nverts = length(D.inv{1}.forward.mesh.vert);
    if(S.lead)
    [L,D] = spm_eeg_lgainmat(D,1:nverts);
    end
    spm_eeg_inv_checkforward(D,1,1);
end
save(D);
fprintf('%-40s: %30s\n','Completed',spm('time'));

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Create Sensor Array                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pos,ori] = opm_createSensorArray(S)
% Given a scalp surface even samples the surface with sensors
%
% FORMAT [pos,ori] = spm_opm_createSensorArray(S)
%   S               - input structure
% Fields of S:
%   S.D            - SPM M/EEG object with surface meshes
%   S.offset       - distance to place sensors(mm) from scalp surface
%   S.space        - distance between sensors(mm)
%   S.wholehead    - boolean: Should whole scalp surface should be covered?
%   S.nDens        - number of density optimisations
% _________________________________________________________________________

% Args
%--------------------------------------------------------------------------
D = S.D;
wholehead = S.wholehead;

% Meshes
%--------------------------------------------------------------------------
scalp = gifti(D.inv{1}.mesh.tess_scalp);
cortex = gifti(D.inv{1}.mesh.tess_ctx);
lp = min(cortex.vertices(:,3));

% create convex hull of mesh
%--------------------------------------------------------------------------
[~,V1] = convhull(double(scalp.vertices));

% outward facing vertx normals of convex hull
%--------------------------------------------------------------------------
[Nv,~] = spm_mesh_normals(scalp,true);
ns = Nv;
vs = scalp.vertices;
cog1=mean(vs);

for i = 1:length(ns)
    ad = vs(i,:) + 5*ns(i,:);
    subtrac = vs(i,:) - 5*ns(i,:);

    d1 = sum((cog1 - ad).^2);
    d2 = sum((cog1 - subtrac).^2);

    if(d2>d1)
        ns(i,:) = -ns(i,:);
    else
        ns(i,:) =ns(i,:);
    end
end

% add an offset to the convex hull
%--------------------------------------------------------------------------
offVertices=vs+ns*S.offset;
[~,V2] = convhull(double(offVertices));

% Expand solution space by ratio of Volumes
%--------------------------------------------------------------------------
T= eye(4)*(V2/V1)^(1/3);
T(4,4)=1;
scalp = spm_mesh_transform(scalp,T);

% Translate solution space so centre of gravities overlap
%--------------------------------------------------------------------------

cog2 = mean(scalp.vertices);
T= eye(4);
T(1:3,4)=cog1-cog2+0;
scalp = spm_mesh_transform(scalp,T);

% Create the sensor array
%--------------------------------------------------------------------------
args= [];
args.space=S.space;
args.g=scalp;
args.nDens=S.nDens;
[pos, ~, ~] = spm_mesh_pack_points(args);

% get orientation of scalp
%--------------------------------------------------------------------------
[~,Nf] = spm_mesh_normals(scalp,true);
cog=mean(scalp.vertices);

v=scalp.vertices;
f=scalp.faces;
faceP=zeros(size(f,1),3);
for i =1:length(faceP)
    whichVerts= f(i,:);
    verts= v(whichVerts,:);
    faceP(i,:)=mean(verts);
end

% make all normals point outwards
%--------------------------------------------------------------------------
for i = 1:length(Nf)
    ad = faceP(i,:) + 5*Nf(i,:);
    subtrac = faceP(i,:) - 5*Nf(i,:);

    d1 = sum((cog - ad).^2);
    d2 = sum((cog - subtrac).^2);

    if(d2>d1)
        Nf(i,:) = -Nf(i,:);
    else
        Nf(i,:) =Nf(i,:);
    end
end

% assign orientatation to sensors of closest face normal
%--------------------------------------------------------------------------
  ori=zeros(size(pos));
for i =1:length(ori)
    tmp=pos(i,:);
    di=sqrt(sum(bsxfun(@minus,faceP,tmp).^2,2));
    [~,indmin]=min(di);
    ori(i,:)=-Nf(indmin,:);
end

% Check if wholehead is requested
%--------------------------------------------------------------------------
if(~wholehead)
    C= pos(:,3)>(lp);
    pos= pos(C,:);
    ori = ori(C,:);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Custom Meshes                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function D = opm_customMeshes(S)
% wrapper for adding custom meshes to MEG object
% FORMAT D = spm_opm_customMeshes(S)
%   S               - input structure
% Fields of S:
%   S.D             - Valid MEG object         - Default:
%   S.cortex        - Cortical mesh file       - Default: Use inverse normalised cortical mesh
%   S.scalp         - Scalp mesh file          - Default: Use inverse normalised scalp mesh
%   S.oskull        - Outer skull mesh file    - Default: Use inverse normalised outer skull mesh
%   S.iskull        - Inner skull mesh file    - Default: Use inverse normalised inner skull mesh
%   S.template      - is mesh in MNI space?    - Default: 0
% Output:
%  D           - MEEG object
%--------------------------------------------------------------------------

%- Default values & argument check
%--------------------------------------------------------------------------
if ~isfield(S, 'D'),           error('MEG object needs to be supplied'); end
if ~isfield(S, 'cortex'),      S.cortex = []; end
if ~isfield(S, 'scalp'),       S.scalp = []; end
if ~isfield(S, 'oskull'),      S.oskull = []; end
if ~isfield(S, 'iskull'),      S.iskull = []; end
if ~isfield(S, 'template'),    S.template = 0; end

D = S.D;
if ~isfield(D.inv{1}.mesh,'sMRI')
  error('MEG object needs to be contain inverse normalised meshes already')
end

%- add custom scalp and skull meshes if supplied
%--------------------------------------------------------------------------
if ~isempty(S.scalp)
    D.inv{1}.mesh.tess_scalp = S.scalp;
end

if ~isempty(S.oskull)
    D.inv{1}.mesh.tess_oskull = S.oskull;
end

if ~isempty(S.iskull)
    D.inv{1}.mesh.tess_iskull = S.iskull;
end

%- add custom cortex and replace MNI cortex with warped cortex
%--------------------------------------------------------------------------
if ~isempty(S.cortex)
    D.inv{1}.mesh.tess_ctx = S.cortex;
    if(S.template)
        D.inv{1}.mesh.tess_mni = S.cortex;
    else
        defs.comp{1}.inv.comp{1}.def = {D.inv{1}.mesh.def};
        defs.comp{1}.inv.space = {D.inv{1}.mesh.sMRI};
        defs.out{1}.surf.surface = {D.inv{1}.mesh.tess_ctx};
        defs.out{1}.surf.savedir.savesrc = 1;
        out = spm_deformations(defs);
        D.inv{1}.mesh.tess_mni  = export(gifti(out.surf{1}), 'spm');
    end
end
save(D);

 end

author: spm


text:
function writeTensors(filename, tensors)

% use this function to write .tensor files (features or weights) into Matlab

fin = fopen(filename, 'wb');

for i=1:length(tensors)

    switch tensors(i).type
        case 'half'
            type_id = 0;
            precision = 'uint16';
        case 'float'
            type_id = 1;
            precision = 'single';
        case 'double'
            type_id = 2;
            precision = 'double';
        case 'uint8'
            type_id = 3;
            precision = 'uint8';
        case 'uint16'
            type_id = 4;
            precision = 'uint16';
        case 'uint32'
            type_id = 5;
            precision = 'uint32';
        case 'uint64'
            type_id = 6;
            precision = 'uint64';
        case 'int8'
            type_id = 7;
            precision = 'int8';
        case 'int16'
            type_id = 8;
            precision = 'int16';
        case 'int32'
            type_id = 9;
            precision = 'int32';
        case 'int64'
            type_id = 10;
            precision = 'int64';
        case 'char'
            type_id = 11;
            precision = 'char';
        case 'bool'
            type_id = 12;
            precision = 'uint8';
        otherwise
            throw(MException('writeTensor:UnsupportedFormat','Unsupported Format'))
    end

    if strcmp(tensors(i).type,'half')
        values = float2half(tensors(i).value);
        tensors(i).value = values;
    end

    fwrite(fin, uint8(type_id), 'uint8');
    fwrite(fin, uint32(tensors(i).sizeof), 'uint32');

    str = tensors(i).name;
    fwrite(fin, int32(length(str)), 'int32');
    fwrite(fin, str, 'char*1');

    dim = size(tensors(i).value);

    % matlab is column first, marvin/c++ is row first
    if length(dim)>1
        pvec = 1:length(dim);
        pvec([1 2]) = pvec([2 1]);
        tensors(i).value = permute(tensors(i).value, pvec);
        dim = dim(pvec);
    end

    dim = dim(end:-1:1);
    while tensors(i).dim > numel(dim)
        dim(end+1) = 1;
    end

    fwrite(fin, int32(length(dim)), 'int32');
    for d=1:length(dim)
        fwrite(fin, int32(dim(d)), 'int32');
    end
    fwrite(fin, tensors(i).value, precision);
end

fclose(fin);

author: jianxiongxiao

repo: ProfXkit


text:
function [verbose,verboseI,debug,doPlot,maxFunEvals,maxIter,tolFun,tolX,method,...
    corrections,c1,c2,LS_init,LS,cgSolve,qnUpdate,cgUpdate,initialHessType,...
    HessianModify,Fref,useComplex,numDiff,LS_saveHessianComp,...
    DerivativeCheck,Damped,HvFunc,bbType,cycle,...
    HessianIter,outputFcn,useMex,useNegCurv,precFunc] = ...
    minFunc_processInputOptions(o)

% Constants
SD = 0;
CSD = 1;
BB = 2;
CG = 3;
PCG = 4;
LBFGS = 5;
QNEWTON = 6;
NEWTON0 = 7;
NEWTON = 8;
TENSOR = 9;

verbose = 1;
verboseI= 1;
debug = 0;
doPlot = 0;
method = LBFGS;
cgSolve = 0;

o = toUpper(o);

if isfield(o,'DISPLAY')
    switch(upper(o.DISPLAY))
        case 0
            verbose = 0;
            verboseI = 0;
        case 'FINAL'
            verboseI = 0;
        case 'OFF'
            verbose = 0;
            verboseI = 0;
        case 'NONE'
            verbose = 0;
            verboseI = 0;
        case 'FULL'
            debug = 1;
        case 'EXCESSIVE'
            debug = 1;
            doPlot = 1;
    end
end

LS_init = 0;
c2 = 0.9;
LS = 4;
Fref = 1;
Damped = 0;
HessianIter = 1;
if isfield(o,'METHOD')
    m = upper(o.METHOD);
    switch(m)
        case 'TENSOR'
            method = TENSOR;
        case 'NEWTON'
            method = NEWTON;
        case 'MNEWTON'
            method = NEWTON;
            HessianIter = 5;
        case 'PNEWTON0'
            method = NEWTON0;
            cgSolve = 1;
        case 'NEWTON0'
            method = NEWTON0;
        case 'QNEWTON'
            method = QNEWTON;
            Damped = 1;
        case 'LBFGS'
            method = LBFGS;
        case 'BB'
            method = BB;
            LS = 2;
            Fref = 20;
        case 'PCG'
            method = PCG;
            c2 = 0.2;
            LS_init = 2;
        case 'SCG'
            method = CG;
            c2 = 0.2;
            LS_init = 4;
        case 'CG'
            method = CG;
            c2 = 0.2;
            LS_init = 2;
        case 'CSD'
            method = CSD;
            c2 = 0.2;
            Fref = 10;
            LS_init = 2;
        case 'SD'
            method = SD;
            LS_init = 2;
    end
end

maxFunEvals = getOpt(o,'MAXFUNEVALS',1000);
maxIter = getOpt(o,'MAXITER',500);
tolFun = getOpt(o,'TOLFUN',1e-5);
tolX = getOpt(o,'TOLX',1e-9);
corrections = getOpt(o,'CORR',100);
c1 = getOpt(o,'C1',1e-4);
c2 = getOpt(o,'C2',c2);
LS_init = getOpt(o,'LS_INIT',LS_init);
LS = getOpt(o,'LS',LS);
cgSolve = getOpt(o,'CGSOLVE',cgSolve);
qnUpdate = getOpt(o,'QNUPDATE',3);
cgUpdate = getOpt(o,'CGUPDATE',2);
initialHessType = getOpt(o,'INITIALHESSTYPE',1);
HessianModify = getOpt(o,'HESSIANMODIFY',0);
Fref = getOpt(o,'FREF',Fref);
useComplex = getOpt(o,'USECOMPLEX',0);
numDiff = getOpt(o,'NUMDIFF',0);
LS_saveHessianComp = getOpt(o,'LS_SAVEHESSIANCOMP',1);
DerivativeCheck = getOpt(o,'DERIVATIVECHECK',0);
Damped = getOpt(o,'DAMPED',Damped);
HvFunc = getOpt(o,'HVFUNC',[]);
bbType = getOpt(o,'BBTYPE',0);
cycle = getOpt(o,'CYCLE',3);
HessianIter = getOpt(o,'HESSIANITER',HessianIter);
outputFcn = getOpt(o,'OUTPUTFCN',[]);
useMex = getOpt(o,'USEMEX',1);
useNegCurv = getOpt(o,'USENEGCURV',1);
precFunc = getOpt(o,'PRECFUNC',[]);
end

function [v] = getOpt(options,opt,default)
if isfield(options,opt)
    if ~isempty(getfield(options,opt))
        v = getfield(options,opt);
    else
        v = default;
    end
else
    v = default;
end
end

function [o] = toUpper(o)
if ~isempty(o)
    fn = fieldnames(o);
    for i = 1:length(fn)
        o = setfield(o,upper(fn{i}),getfield(o,fn{i}));
    end
end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function poly = boxToPolygon(box)
%BOXTOPOLYGON Convert a bounding box to a square polygon.
%
%   poly = boxToPolygon(box)
%   Utility function that convert box data in [XMIN XMAX YMIN YMAX] format
%   to polygon data corresponding to the box boundary. The resulting POLY
%   is a 4-by-2 array.
%
%
%   Example
%     box = [ 10 50 20 40];
%     poly = boxToPolygon(box)
%     poly =
%         10    20
%         50    20
%         50    40
%         10    40
%
%   See also
%     boxes2d, polygons2d, boxToRect

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2017-09-10, using Matlab 8.6.0.267246 (R2015b)
% Copyright 2017-2022 INRA - Cepia Software Platform

% extreme coordinates
xmin = box(1);
xmax = box(2);
ymin = box(3);
ymax = box(4);

% convert to polygon
poly = [...
    xmin ymin; ...
    xmax ymin; ...
    xmax ymax; ...
    xmin ymax];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function q=transpose(q)
% transpose array of quaternions

q.a = q.a.';
q.b = q.b.';
q.c = q.c.';
q.d = q.d.';

author: mtex-toolbox

repo: mtex

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


text:
function h = histogram(grains,varargin)
% grain area weighted histogram
%
% Syntax
%   histogram(grains)      % area weighted histogram of the grain area
%
%   prop = grains.equivalentRadius;
%   histogram(grains,prop) % area weighted histogram of property prop
%
%   histogram(grains,prop,n)    % specify number of bins
%   histogram(grains,prop,bins) % specify bin edges
%
%   h = histogram(grains)
%
% Input
%  grains - @grain2d
%  n      - number of bin edges, default ist 15, (number of bins is n-1)
%  bins   - vector of bin edges
%
% Output
%  h - handle to the histogram graphics object
%

[mtexFig,isNew] = newMtexFigure(varargin{:});
mtexFig.keepAspectRatio = false;

% exract area and maybe an aditional property
area = grains.area;
if nargin>1 && isnumeric(varargin{1}) && length(varargin{1}) == length(grains)
  prop = varargin{1};
  varargin(1) = [];
else
    prop = area;
end

% generate bins
if ~isempty(varargin) && isnumeric(varargin{1})
  if numel(varargin{1})==1
    nbins = varargin{1}; %define nbins
  elseif numel(varargin{1})>1
    bins = varargin{1}; %define bin edges automatically
    nbins = numel(varargin{1});
  end
else
  nbins = 15;
end
if ~exist('bins','var'), bins = linspace(0,max(prop)+eps,nbins); end

% loop through all phases
h = [];
for id = grains.indexedPhasesId

  % find for each area the binId
  [~,~,binId] = histcounts(prop(grains.phaseId==id),bins);
  if any(~binId,'all')
    warning([num2str(nnz(~binId)) ' grains outside bin limits not plotted!']);
  end

  % compute the sum of areas belonging to the same bin

  areaPhase = area(grains.phaseId==id);
  cumArea = accumarray(binId(binId>0),areaPhase(binId>0),[length(bins)-1 1],@nansum) ./ sum(area);

  h = [h,optiondraw( histogram('BinEdges',bins,'BinCounts',cumArea,...
    'FaceColor',grains.CSList{id}.color),varargin{:})]; %#ok<AGROW>
  if strcmp(get(h,'DisplayStyle'),'stairs')
    set(h(grains.indexedPhasesId==id),'EdgeColor',grains.CSList{id}.color);
  end
  hold on

end
hold off

% labels and title
if all(prop == area)
  title(mtexFig.gca,'grain size distribution')
  xlabel(mtexFig.gca,'grain area');
end
yticklabels(mtexFig.gca,yticks*100)
ylabel(mtexFig.gca,'relative area (%)')

% legend
min = grains.mineralList(grains.indexedPhasesId);
legend(min{:})

if isNew, mtexFig.drawNow(varargin{:});end

if nargout == 0, clear h; end


text:
% RES = histoMatch(MTX, N, X)
%
% Modify elements of MTX so that normalized histogram matches that
% specified by vectors X and N, where N contains the histogram counts
% and X the histogram bin positions (see histo).

% Eero Simoncelli, 7/96.

function res = histoMatch(mtx, N, X)

if ( exist('histo') == 3 )
  [oN, oX] = histo(mtx(:), size(X(:),1));
else
  [oN, oX] = hist(mtx(:), size(X(:),1));
end

oStep = oX(2) - oX(1);
oC = [0, cumsum(oN)]/sum(oN);
oX = [oX(1)-oStep/2, oX+oStep/2];

N = N(:)';
X = X(:)';
N = N + mean(N)/(1e8);   %% HACK: no empty bins ensures nC strictly monotonic

nStep = X(2) - X(1);
nC = [0, cumsum(N)]/sum(N);
nX = [X(1)-nStep/2, X+nStep/2];

nnX = interp1(nC, nX, oC, 'linear');

if ( exist('pointOp') == 3 )
  res = pointOp(mtx, nnX, oX(1), oStep);
else
  res = reshape(interp1(oX, nnX, mtx(:)),size(mtx,1),size(mtx,2));
end

repo: SelfExSR


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtHmxLowrank.m                               |
%|    #    |   VERSION    : 0.50                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 25.11.2018                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   : Compare compressor with total pivoting        |
%|  `---'  |                                                              |
%+========================================================================+

clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Dimensions
Nx  = 100;
Ny  = 101;
rk  = 10;
tol = 1e-3;

% Compress zeros
M     = zeros(Nx,Ny);
[A,B] = hmxACA(M,tol);
norm(A*B,'inf')

% Compress ones
M     = ones(Nx,Ny);
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress eye
M          = eye(Nx,Ny);
[~,~,flag] = hmxACA(M,tol);
flag

% Compress circular permutation for eye
M = circshift(eye(Nx,Ny),3,2);
[~,~,flag] = hmxACA(M,tol);
flag

% Compress 1-eye
M          = 1-eye(Nx,Ny);
[~,~,flag] = hmxACA(M,tol);
flag

% Compress random
M     = rand(Nx,rk) * rand(rk,Ny);
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress random complex
M     = (-1+2*rand(Nx,rk)+1i*(-1+2*rand(Nx,rk))) * (rand(rk,Ny) + 1i*rand(rk,Ny));
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 1x2 block with zeros
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(Nx,10);
M     = [M Z];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 2x1 block with zeros
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(10,Ny);
M     = [M ; Z];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 2x2 block random
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(Nx,Ny);
M     = [M Z ; Z M];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 4x4 block random
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(Nx,Ny);
M     = [M Z ; Z M];
Z     = zeros(size(M));
M     = [M Z ; Z M];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 2x2 inegal block random (to be hardly improved...)
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(Nx,1);
M     = [M zeros(Nx,1) ; zeros(1,Ny) 1];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Compress 4x4 inegal block random (to be hardly improved...)
M     = rand(Nx,rk) * rand(rk,Ny);
Z     = zeros(Nx,1);
M     = [M zeros(Nx,3) ; zeros(3,Ny) eye(3)];
[A,B] = hmxACA(M,tol);
norm(A*B-M,'inf')/norm(M,'inf')

% Build orthogonal mesh
vtx   = [0 0 0 ; eye(3)];
elt   = [1 3 2 ; 1 2 4];
col   = [1 ; 2];
mesh  = msh(vtx,elt,col);
mesh  = refine(mesh,0.2);
[~,I] = sort(mesh.col);
mesh  = mesh.sub(I);

% Translate
mesh2          = mesh;
mesh2.vtx(:,1) = mesh2.vtx(:,1) + 2;

% Graphical representation
figure
plot(mesh)
hold on
plotNrm(mesh)
plot(mesh2)
axis equal
view(-180,30)

% Build double layer matrix
Gxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',1);
Gxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',1);
Gxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',1);
D      = integral(dom(mesh,3),dom(mesh2,3),fem(mesh,'P0'),Gxy,ntimes(fem(mesh2,'P0')));

% Graphical represenation
figure
imagesc(abs(D))
colorbar

% Compression
[A,B] = hmxACA(D,tol);
norm(A*B-D,'inf')/norm(D,'inf')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
close all

figure, plot(tout, control_results)
title('Joint positions')

figure,
plot(tout, joint1_error), hold
plot(tout, joint2_error)
title('Joint errors')

figure,
plot(tout, motor1_current), hold
plot(tout, motor2_current)
title('Motor current')

figure,
plot(tout, motor1_speed_rpm), hold
plot(tout, motor2_speed_rpm)
title('Motor speed rpm')

repo: ARTE


text:
function drawIdpPnt(MapFig, Lmk, color, MapOpt)

% DRAWIDPPNT  Draw inverse-depth point landmark in MapFig.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

posOffset = [0;0;.2];

% transform to Euclidean
switch Map.type
    case 'ekf'
        [x,P] = propagateUncertainty(       ...
            Map.x(Lmk.state.r),             ...
            Map.P(Lmk.state.r,Lmk.state.r), ...
            @idp2euc);
    case 'graph'
        x = idp2euc(Map.x(Lmk.state.r));
%         P = eye(3); % irrelevant because we will not print ellipses
end

% draw
drawPnt(MapFig.Lmk(Lmk.lmk).mean, x, color.mean)
if MapOpt.showEllip
    drawEllipse(MapFig.Lmk(Lmk.lmk).ellipse, x, P, color.ellip)
end
if MapOpt.showLmkId
    drawLabel  (MapFig.Lmk(Lmk.lmk).label,   x+posOffset, num2str(Lmk.id))
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function F = calc_normalCDF(vX, vXdata)
    % Compute CDF of normal probability distribution for diverse parameters
    % REQUIRES STATISTICS TOOLBOX
    %
    % F = calc_normalCDF(vX, vXdata);
    % ----------------------------------------
    % Compute CDF of normal probability distribution for diverse parameters; use as function in lsqcurvefit (or equivalent) to find
    % best fitting parameters
    %
    % Distribution parameters:
    % vX(1) : mu
    % vX(2) : sigma
    %
    % Author: J. Woessner
    % woessner@seismo.ifg.ethz.ch
    % updated: 11.06.04

    F = normcdf(vXdata,vX(1),vX(2));
    % Check F on NAN
    if isnan(F(1))
        F = zeros(length(F),1);
        warning('Probabilities set to zero')
    end
end

author: CelsoReyes


text:
% This script train an LSTM or DNN to predict complex Fourier transform
% domain temporal filters, which will be used to filter the noisy and
% reverberant Fourier coefficients to produce enhanced speech. The script
% is based on the simulated data of Reverb challenge.
%
% Xiong Xiao, Nanyang Technological University, Singapore
% Last modified: 10 Feb 2016.
%
%clear
function TrainDereverbNet_RegressionComplex(modelType, hiddenLayerSize, hiddenLayerSizeFF, DeltaGenerationType, learning_rate, ...
    step, nSequencePerMinibatch, useFileName, useGPU)

addpath('local', '../beamforming/lib');     % remember to run addMyPath.m in SignalGraph root directory first.
% define SignalGraph settings. Type ParseOptions2() in command line to see the options.
para.IO.nStream = 2;        % two data streams, one is input signal, one is clean target log spectrogram
para.NET.sequential = 1;    % do not randomize at frame level, randomize at sentence level
para.NET.variableLengthMinibatch = 1;   % a minibatch may contain multiple sentences of different lengths
para.NET.nSequencePerMinibatch = nSequencePerMinibatch;    % number of sentences per minibatch
para.NET.maxNumSentInBlock = 100;       % maximum number of sentences in a block
para.NET.L2weight = 3e-4;
para.NET.learning_rate = learning_rate;
para.NET.learning_rate_decay_rate = 0.999;
para.NET.momentum = [0.9];
para.NET.gradientClipThreshold = 1;
para.NET.weight_clip = 10;
para.useGPU     = useGPU;                  % don't use GPU for single sentence minibatch for LSTM. It's even slower than CPU.
para.minItr     = 10;
para.displayGPUstatus = 1;
para.displayInterval = 1;
para.skipInitialEval = 1;
para.displayTag = 'LSTM-Regression';

reverb_root = ChoosePath4OS({'D:/Data/REVERB_Challenge', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});   % you can set two paths, first for windows OS and second for Linux OS.
para.local.wavroot = [reverb_root];
para.local.wsjcam0root = ChoosePath4OS({'D:/Data/wsjcam0', '/media/xiaoxiong/OS/data1/G/wsjcam0_wav'});
para.local.wsjcam0ext = ChoosePath4OS({'wav', 'wv1'});
para.local.useFileName = useFileName;      % if set to 0, load all training data to memory. otherwise, only load file names.
para.local.loadLabel = 0;
para.local.seglen = 100;
para.local.segshift = 100;

para.topology.useChannel = 1;
para.topology.useWav = 1;
para.topology.FilterNetType = modelType;
% There are 257 filters, one for each frequency bin. Each filter has a 2D
% complex-valued weight matrix, with size TempContext x FreqContext.
para.topology.FilterTempContext = 21;   % this is the number of frames used for Fourier coefficient domain filtering.
para.topology.FilterFreqContext = 1;    % this is the number of frequency bins used for filter
para.topology.useCMN = 0;
para.topology.mu_law_factor = 2047;
para.topology.DeltaGenerationType = DeltaGenerationType;
para.topology.MSECostWeightSDA = [1 4.5 10];
para.topology.hiddenLayerSize = hiddenLayerSize;
para.topology.hiddenLayerSizeFF = hiddenLayerSizeFF;
para = ConfigDereverbNet_Regression(para);

[Data_small, para] = LoadParallelWavLabel_Reverb(para, 100, 'train', {'simu'}, {'far'});
[layer, para] = Build_DereverbNet_RegressionComplex(Data_small, para);

% load the training and cv data
[Data_tr, para] = LoadParallelWavLabel_Reverb(para, step, 'train', {'simu'}, {'far'});
[Data_cv, para] = LoadParallelWavLabel_Reverb(para, step, 'dev', {'simu'}, {'far', 'near'});

% generate directory and file names to store the networks.
para.output = 'nnet/DereverbRegressionComplex';
para.output = sprintf('%s.MbSize%d.U%d.%d-%s', para.output, ...
    para.NET.nSequencePerMinibatch, length(Data_tr(1).data),  3*(para.topology.fft_len/2+1), para.topology.RegressionNetType);
for i=1:length(para.topology.hiddenLayerSize)
    para.output = sprintf('%s-%d', para.output, para.topology.hiddenLayerSize(i));
end
if length(para.topology.hiddenLayerSizeFF)>0
    para.output = [para.output '-DNN'];
    for i=1:length(para.topology.hiddenLayerSizeFF)
        para.output = sprintf('%s-%d', para.output, para.topology.hiddenLayerSizeFF(i));
    end
end
para.output = sprintf('%s-%d.L2_%s.LR_%s/nnet', para.output, 3*(para.topology.fft_len/2+1), FormatFloat4Name(para.NET.L2weight),FormatFloat4Name(para.NET.learning_rate));
LOG = [];

% show the configurations
para
para.NET
para.IO
para.topology

% train the network using SGD
trainGraph_SGD(layer, Data_tr, Data_cv, para, LOG);

author: singaxiong

repo: SignalGraph


text:
% Test file for @deltafun/isempty.m

function pass = test_isempty(pref)

if (nargin < 1)
    pref = chebfunpref();
end
%%
d = deltafun();
pass(1) = isempty(d);

f = bndfun([]);
d = deltafun(f);
pass(2) = isempty(d);

d = deltafun(f, []);
pass(3) = isempty(d);

d = deltafun(f, struct('deltaMag', [], 'deltaLoc', []));
pass(4) = isempty(d);

end


text:
function [phd_trimmed, cphd_nb, meta] = ExtractChipPhaseHistory(cphd_filename,SRP_new,img_support_size,channel)
% EXTRACTCHIPPHASEHISTORY extracts user-specified area from phase history
%   data = ExtractChipPhaseHistory(cphd_filename,img_offset,img_width,channel)
%   filters the input phase history to produce a reduced size phase history
%   that supports full resolution imaging of a user specified image chip
%   (slant plane rectangular area of dimensions IMG_SUPPORT_SIZE centered
%   on SRP_NEW.)
%
%   Separable 1-D range and 1-D azimuth processing, range first, azimuth
%   second, because of memory constraints.
%
% ASSUMPTIONS:
%   1) Assumes a spotlight collect, since it is assumed all points on the
%   ground are equally illuminated by all pulses.
%   2) Assumes a (nearly) monostatic collect.
%   3) This code assumes that the full phase history might not fit into
%   memory, but that the range trimmed phase history must be able to fit
%   into memory.
%   4) Also currently assumes that all pulses have common bandwidths and
%   frequency sampling.  If this is not the case, things become messy.
%   Azimuth filtering becomes difficult/impossible for pulses with
%   arbitrarily different frequency content.  Furthermore, its not clear
%   what the Fx0/F_SS of the resulting decimated pulses would be.
%
% INPUTS:
%   cphd_filename    : string : input phase history file
%   SRP_new          : [1x3]  : ECEF coordinates of new scene reference
%                               point.
%   img_support_size : If scalar, this is the radius of a sphere around the
%                      SRP to extract.
%                      If [1x2], this is the image chip extent in meters
%                      [azimuth_extent range_extent]  (in the slant plane).
%   channel          : scalar : Channel in the phase history file to chip.
%                               (Default = 1).
%
% OUTPUTS:
%   phd_trimmed : Phase history trimmed as requested with pulses in columns
%   cphd_nb     : CPHD formatted per pulse narrowband data for phd_trimmed
%   meta        : SICD-like meta structure
%
% Dan Hack, AFIT/ENG, 09AUG11
% Tim Cox, NRL, 19JUL12 - Modifications to bring to CPHD "X" format.
% Wade Schwartzkopf, NGA/IDT, 01DEC12 - Documentation and code clean up.
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%%%%%%%%%%%%%%%%%%
% Constants      %
%%%%%%%%%%%%%%%%%%
c = 299792458; % Speed of light (m/s)
% Algorithm parameters
show_waitbar = true;
pBlockSize = 2000;
verification_pulse = 0; % Pulse number (or set to zero for no verification plots)

%%%%%%%%%%%%%%%%%%
% Pre-processing %
%%%%%%%%%%%%%%%%%%

% Default input parameters
if ~exist('channel','var')
    channel = 1; % Default
end
% Get the narrowband data
reader_obj = open_ph_reader(cphd_filename);
meta = reader_obj.get_meta();
N = double(meta.Data.Channel(channel).NumVectors);
M = double(meta.Data.Channel(channel).NumSamples);  % num_samples per pulse
[ignore, all_nb] = reader_obj.read_cphd(1:N,[],1);
all_nb.ARP = (all_nb.TxPos + all_nb.RcvPos)/2; % Aperture reference position (nearly monostatic assumption)
% This function assumes spotlight data.
if isfield(meta, 'CollectionInfo') && isfield(meta.CollectionInfo, 'RadarMode') && ...
        isfield(meta.CollectionInfo.RadarMode, 'ModeType') && ...
        ~strcmpi(meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
    error('EXTRACTCHIPPHASEHISTORY:NON_SPOTLIGHT_DATA',...
        'ExtractChipPhaseHistory function requires spotlight data.');
end
% Determine some collection parameters
[resolution, extent, delta_azimuth, total_azimuth] = pulse_info_to_resolution_extent(...
    all_nb.ARP([1 end],:) - all_nb.SRPPos([1 end],:), ... % Line-of-sight vector between ARP and SRP
    max(all_nb.SC0  + (all_nb.SCSS * (M-1))),... % Highest frequency sample (which has the least azimuth extent, so most conservative)
    max(all_nb.SCSS),... % Largest sample spacing (which has the least range extent)-- although we don't use it here
    [],... % This parameter should be bandwidth, but we don't need it here, since we aren't using resolution.
    N);

%%%%%%%%%%%%%%%%%%%%%%%
% 1:  Range Windowing %
%%%%%%%%%%%%%%%%%%%%%%%

% Slice out the required range bins from each pulse.

% First we determine how much range distance we need to cut out of each
% pulse.
if isscalar(img_support_size) % A sphere around the SRP, the easiest.
    image_span_distance = img_support_size;
elseif numel(img_support_size)==2
    % Extracting range/azimuth oriented rectangle in the slant plane of
    % dimensions img_support_size ([azimuth, range]).
    %
    % In theory, since the geometry (and potentially frequency content) of
    % each pulse varies, each pulse could require a (potentially
    % drammatically) different set of range bins to exactly span our slant
    % plane rectangle. For the selected range bins to cover all corners of
    % the slant plane rectangle, including geometry and range curvature
    % contraints, the required distance span around the SRP for a single
    % pulse is:
    %    image_span_distance = 2*(sqrt(...
    %         (SRP_range*cos(theta) + (img_support_size(2)/2)).^2 + ...
    %         (SRP_range*abs(sin(theta)) + (img_support_size(1)/2)).^2) - ...
    %         SRP_range)
    % where
    %    theta is the angle of a given pulse off of the reference range
    %       direction for which our range/azimuth slant plane rectangle is
    %       defined.
    %    SRP_range is the distance from the ARP to SRP.
    %
    % Although the above equation may be exact per pulse, this distance
    % likely equates to a different number of range bins for each pulse,
    % and in MATLAB we prefer to work with rectangular arrays.  So for the
    % purposes of simplifying the computation, we pick the minimum distance
    % that spans our slant plane rectangle for all pulses.  This is not
    % optimum from the standpoint of the minimum amount of information
    % required to describe the requested image area, but it is required to
    % keep the computations tractable in MATLAB.
    % For simplicity, we only need to evaluate the above
    % image_span_distance equation with the pulse geometries at the ends of
    % the collect (when the theta angle, and thus the require range, is the
    % greatest) and with the minimum range to SRP for all pulses (since
    % that is when range curvature is at its greatest).
    min_SRP_range = min(sqrt(sum((all_nb.ARP - all_nb.SRPPos).^2,2)));
    image_span_distance = 2*(sqrt(...
        (min_SRP_range*cos(total_azimuth/2) + (img_support_size(2)/2)).^2 + ...
        (min_SRP_range*sin(total_azimuth/2) + (img_support_size(1)/2)).^2) - ...
        min_SRP_range);
    % Do we need some small amount of buffer on top of this?
else % For other areas (or volumes), this computation would have to change.
    error('EXTRACTCHIPPHASEHISTORY:UNRECOGNIZED_IMG_SUPPORT_SIZE',...
        'Unrecognized IMG_SUPPORT_SIZE format.');
end

% Determine range parameters
range_extent = c/(2*max(all_nb.SCSS));  % range profile extent (for pulse with the least extent)
rc_vec = fft_bin_pos_vec(M) * range_extent;  % Distance of each range bin from SRP
rc_ndx = find(abs(rc_vec) <= image_span_distance/2);  % range extraction indices
% These indices extract at least the range width needed to cover the
% requested area for all pulses.

% Allocate storage.  We assume this entire array of range-trimmed phase
% history can fit into memory.
M_new = length(rc_ndx); % Number of frequency samples after decimating
phd_trimmed_rangeonly = zeros(M_new,N);

% This is a vectorized implementation in which a block of pulses is
% processed on each iteration.  We process in blocks because it is likely
% that the entire untrimmed phase history cannot fit into memory.
%
% This loop processes all pulses in the phase history file because we
% assume a spotlight collect where all pulses equally affect our new
% reference point.  For non-spotlight collects, we should limit this loop
% to only pulses for which our new reference point fall in their beam.
% This would further reduce computation.

if show_waitbar
  wb = waitbar(0); tic
end

NumBlocks = ceil(N/pBlockSize); % Determine number of blocks.
for pbDx = 1:NumBlocks
    % Calculate pulse indices for this block
    pDx = ( ((pbDx-1)*pBlockSize+1):min(pbDx*pBlockSize,N) );

    % Read all samples from a block of pulses
    ph = reader_obj.read_cphd(pDx, 1:M, channel);

    % Calculate differential ranges (difference of range to current SRP and
    % range to desired SRP)
    delR = sqrt(sum((all_nb.ARP(pDx,:).' - repmat(SRP_new(:),[1 numel(pDx)])).^2)) - ...
           sqrt(sum((all_nb.ARP(pDx,:).' - all_nb.SRPPos(pDx,:).').^2));

    % Calculate frequency vector for each pulse: Fx0+Fx_SS*(0:M-1)
    f_orig = bsxfun(@plus,all_nb.SC0(pDx).',(0:M-1).'*all_nb.SCSS(pDx).');
    k_orig = 2*pi*f_orig/c;

    % Adjust phase reference (re-mocomp to new scene reference point)
    % This is currently just a linear shift based on difference in range
    % between the old and new reference points.  This is the 99% solution.
    % Might be some other higher-order components like tropospheric
    % adjustment.
    ph = ph.*exp(1i*2*bsxfun(@times,k_orig,delR));

    % Form range profile (IFFT)
    rc = ifft(ph,[],1);

    % Apply range gate and convert back to spatial freq domain
    phd_trimmed_rangeonly(:,pDx) = fft(rc(rc_ndx,:),[],1);

    % Generate plot comparisons for data from a single pulse
    if (verification_pulse ~= 0) && any(pDx == verification_pulse)
        % Recompute so we are using Fx_SS of this specific pulse, not max(Fx_SS)
        range_extent_orig = c/(2*all_nb.SCSS(verification_pulse));
        rc_vec_orig = fft_bin_pos_vec(M) * range_extent_orig;
        % Compute extent and range profile of trimmed pulse
        range_extent_new = c/(2*all_nb.SCSS(verification_pulse)*M/M_new);
        rc_vec_new = fft_bin_pos_vec(M_new) * range_extent_new;
        rc_new = ifft(phd_trimmed_rangeonly(:,verification_pulse));
        % Compare pre- and post- trimming
        figure;
        plot(fftshift(rc_vec_orig), fftshift(abs(rc(:,pDx==verification_pulse))), 'b', ...
            fftshift(rc_vec_new), fftshift(abs(rc_new)), 'r');
        set(gca,'xlim',[min(rc_vec) max(rc_vec)]);
        legend({'Original' 'Windowed'});
        xlabel('Range (m)'); ylabel('Magnitude');
        title(['Range Profile Comparison, Pulse ' num2str(verification_pulse)]);
    end

    if show_waitbar
        % Determine remaining execution time and display
        t_sofar = toc;
        t_est = (t_sofar*NumBlocks/pbDx)-t_sofar;
        wb_message=sprintf('Range Filtering, Time remaining: %s',datestr(datenum(0,0,0,0,0,t_est),13));
        waitbar(pbDx/NumBlocks,wb,wb_message);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2:  Azimuth (Doppler) Windowing %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Due to re-mocomping to new reference point in the range processing step,
% the desired scene center is already at zero Doppler.
%
% Range filtering is exact.  Azimuth filtering (as implemented below) is
% approximate:
%
% 1) Pulses might not have uniform angular sampling. Since the FFT (and
% most filters) assume uniform sampling, this could be a point of concern.
% For now we assume that spacing is "close enough" to uniform, and that any
% distortions from the FFTing of non-uniform data will be somewhat undone
% with the inverse transform.  Also when we interpolate to new ARP
% positions later, we will retain this non-uniformity by interpolating with
% reference to pulse number rather than time.
% 2) Since we are filtering (at least roughly) across frequency sample
% here, there are some "polar formatting" approximations going on here.
%
% Note that if pulses have signficantly different Fx0 and/or Fx_SS, this
% step probably doesn't make any sense.

if show_waitbar
    waitbar(0,wb,'Azimuth Filtering');
end

% Azimuth profile vector
ac_vec = fft_bin_pos_vec(N)*extent(2); % extent(2) is azimuth extent

% Azimuth extraction indices
ac_ndx = find(abs(ac_vec) <= 1.1 * img_support_size(1)/2); % 10% extra buffer since our filtering is inexact

% Number of azimuth frequency samples after decimating
N_new = length(ac_ndx);

if show_waitbar
    waitbar(0.33,wb,'Azimuth Filtering');
end

% Azimuth compression (IFFT)
ac = ifft(phd_trimmed_rangeonly,[],2);

if show_waitbar
    waitbar(0.66,wb,'Azimuth Filtering');
end

% Window and azimuth uncompress (FFT)
phd_trimmed = fft(ac(:,ac_ndx),[],2);

if show_waitbar
    waitbar(1,wb,'Azimuth Filtering');
    pause(0.25);
    close(wb);
end

%%%%%%%%%%%%%%%%%%%%%%%
% Form pulse metadata %
%%%%%%%%%%%%%%%%%%%%%%%

% Generate output data structure.  Generic structure that reflects CPHD
% format metadata.

% Interpolate aperture reference position (ARP) -- assuming that
% interpolating the original trajectory in time will maintain phase
% coherency with respect to decimated pulses.
ARP_interp = interp1(1:N,all_nb.ARP,linspace(1,N,N_new),'spline');
% Time is really irrelevant in CPHD, but we compute it anyway.
time_ref = (all_nb.TxTime + all_nb.RcvTime)/2; % Monostatic approximation
time_interp = interp1(1:N,time_ref,linspace(1,N,N_new).','spline');

cphd_nb.TxTime = time_interp;
cphd_nb.TxPos = ARP_interp;
cphd_nb.RcvTime = time_interp;
cphd_nb.RcvPos = ARP_interp;
cphd_nb.SRPPos = repmat(SRP_new(:).',[N_new 1]);
% Interpolation along Fx probably isn't valid if Fx values jitter.
cphd_nb.SC0 = interp1(1:N,all_nb.SC0,linspace(1,N,N_new).','spline');
cphd_nb.SCSS = interp1(1:N,all_nb.SCSS*(M/M_new),linspace(1,N,N_new).','spline');
cphd_nb.FX1 = interp1(1:N,all_nb.FX1,linspace(1,N,N_new).','spline');
cphd_nb.FX2 = interp1(1:N,all_nb.FX2,linspace(1,N,N_new).','spline');

end

% Compute a vector of the positions of the bins in an FFT'd vector.
% Output is in normalized units (1 is Nyquist).  Need to multiply by extent
% to convert to real units like meters (or Hz).
function bin_pos_norm = fft_bin_pos_vec(Nfft)
    bin_pos_norm = (0:(Nfft-1))/Nfft; % Normalized units (no sampling frequency used)
    ind_to_flip = (floor((Nfft-1)/2)+2):Nfft; % These positions are 0.5 or above.  Should be negative units.
    bin_pos_norm(ind_to_flip) = bin_pos_norm(ind_to_flip) - 1;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function DEM_demo_DFP
% DEM demo for linear deconvolution:  This demo considers the deconvolution
% of the responses of a single-input-multiple output input-state-output
% model (DCM) to disclose the input or causes.  It starts by demonstrating
% Variational filtering with spm_DFP; this is a stochastic filtering scheme
% that propagates particles over a changing variational energy landscape
% such that their sample density can be used to approximate the underlying
% ensemble or conditional density.  We then repeat the inversion using
% spm_DEM (i.e., under a Laplace assumption) which involves integrating the
% path of just one particle (i.e., the mode).

% get a simple convolution model
%==========================================================================
spm_figure('GetWin','DEM');

M        = spm_DEM_M('convolution model');
M(1).V   = exp(8);
M(1).W   = exp(16);
M(1).E.N = 32;

% and generate data
%==========================================================================
N     = 32;                                        % length of data sequence
U     = exp(-((1:N) - N/4).^2/(2*(N/32)^2));       % Gaussian cause
DEM   = spm_DEM_generate(M,U,{},{32 16});

% display
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.pU)

% invert model - VF
%==========================================================================
DEM  = spm_DFP(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');
spm_DEM_qU(DEM.qU,DEM.pU)

% invert model - DEM
%==========================================================================
DEM  = spm_DEM(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.qU,DEM.pU)


text:
classdef math_model_cpf_acpi < mp.math_model_cpf_acp & mp.mm_shared_pfcpf_acpi
%MP.MATH_MODEL_CPF_ACPI  MATPOWER mathematical model for continuation power flow (CPF) problem.
%   ?
%
%   MP.MATH_MODEL_CPF_ACPI ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

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

%     properties
%     end

    methods
        function tag = form_tag(obj)
            tag = 'acpi';
        end

        function name = form_name(obj)
            name = 'AC-polar-current';
        end

        function obj = add_node_balance_constraints(obj, nm, dm, mpopt)
            %% power balance constraints
            ad = obj.aux_data;
            npvq = ad.npv+ad.npq;
            fcn = @(x)node_balance_equations_cpf(obj, x, nm);
            obj.add_nln_constraint({'Irmis', 'Iimis'}, [npvq;npvq], 1, fcn, []);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function str = deltadatenumstr(dt)
% str = deltadatenumstr(dt)
% Nice string representing time interval. Input is in days.

ONEHOUR = 1/24; % days
ONEMINUTE = ONEHOUR/60;
ONESECOND = ONEMINUTE/60;

if dt > 1
  str = sprintf('%d days',round(dt));
elseif dt>ONEHOUR
  str = sprintf('%d hours',round(dt/ONEHOUR));
elseif dt > ONEMINUTE
  str = sprintf('%d minutes',round(dt/ONEMINUTE));
elseif dt > ONESECOND
  str = sprintf('%d seconds',round(dt/ONESECOND));
else
  str = 'less than 1 second';
end

author: kristinbranson

repo: JAABA


text:
function varargout = cluster_tool_thresh(varargin)
% CLUSTER_TOOL_THRESH M-file for cluster_tool_thresh.fig
%      CLUSTER_TOOL_THRESH, by itself, creates a new CLUSTER_TOOL_THRESH or raises the existing
%      singleton*.
%
%      H = CLUSTER_TOOL_THRESH returns the handle to a new CLUSTER_TOOL_THRESH or the handle to
%      the existing singleton*.
%
%      CLUSTER_TOOL_THRESH('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CLUSTER_TOOL_THRESH.M with the given input arguments.
%
%      CLUSTER_TOOL_THRESH('Property','Value',...) creates a new CLUSTER_TOOL_THRESH or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before cluster_tool_thresh_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to cluster_tool_thresh_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help cluster_tool_thresh

% Last Modified by GUIDE v2.5 20-Jul-2006 15:43:35

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @cluster_tool_thresh_OpeningFcn, ...
                   'gui_OutputFcn',  @cluster_tool_thresh_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before cluster_tool_thresh is made visible.
function cluster_tool_thresh_OpeningFcn(hObject, eventdata, handles, varargin)
movegui(hObject,'center');
set(handles.title,'string',cluster_tool_getver);
set(handles.t,'value',1);set(handles.p,'value',0);set(handles.dat,'value',0);
set(handles.t_thresh,'value',1);set(handles.p_thresh,'value',0);
handles.thrshtyp='t';handles.imtype='t';
handles.pwd=pwd;handles.fdir=pwd;
for k=1:size(varargin(:),1)
    if strcmp(varargin{k},'cl'),handles.cl=varargin{k+1};
    elseif strcmp(varargin{k},'fdir'),handles.fdir=varargin{k+1};
    elseif strcmp(varargin{k},'img'),handles.img=varargin{k+1};
        set(handles.imagename,'string',['Image: ' handles.img]);
    end
end
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to cluster_tool_thresh (see VARARGIN)

% Choose default command line output for cluster_tool_thresh
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes cluster_tool_thresh wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = cluster_tool_thresh_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --------------------------------------------------------------------
function imgtype_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to imgtype (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.imtype=get(hObject,'string');
if strcmp(handles.imtype,'data')||strcmp(handles.imtype,'p')
    set(handles.threshtype,'visible','off');set(handles.dftxt,'visible','off');set(handles.df,'visible','off');
else set(handles.threshtype,'visible','on');set(handles.dftxt,'visible','on');set(handles.df,'visible','on');
end
if strcmp(handles.imtype,'p'),set(handles.posneg,'visible','off','value',1);
elseif ~strcmp(handles.thrshtyp,'p'),set(handles.posneg,'visible','on');
end
if ~strcmp(handles.imtype,'p'),set(handles.posneg,'visible','on');end
if strcmp(handles.thrshtyp,'t'),handles.thrshtyp='p';set(handles.p_thresh,'value',1);set(handles.t_thresh,'value',0);end
guidata(handles.figure1,handles);

function thresh_Callback(hObject, eventdata, handles)
try handles.threshold=str2double(get(hObject,'string'));catch set(hObject,'string','');handles.threshold=[];end
guidata(hObject,handles);
% hObject    handle to thresh (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of thresh as text
%        str2double(get(hObject,'String')) returns contents of thresh as a double

% --- Executes during object creation, after setting all properties.
function thresh_CreateFcn(hObject, eventdata, handles)
% hObject    handle to thresh (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function threshtype_SelectionChangeFcn(hObject, eventdata, handles)
if strcmp(handles.imtype,'p')
    set(handles.p_thresh,'value',1);set(handles.t_thresh,'value',0);set(handles.posneg,'visible','off','value',1);
else handles.thrshtyp=get(hObject,'string');if ~strcmp(handles.imtype,'t'),set(handles.posneg,'visible','on');end
end
guidata(handles.figure1,handles);
% hObject    handle to threshtype (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function df_Callback(hObject, eventdata, handles)
handles.dfd=str2double(get(hObject,'string'));
if isnan(handles.dfd),set(hObject,'string','');end
guidata(hObject,handles);
% hObject    handle to df (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of df as text
%        str2double(get(hObject,'String')) returns contents of df as a double

% --- Executes during object creation, after setting all properties.
function df_CreateFcn(hObject, eventdata, handles)
% hObject    handle to df (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function k_Callback(hObject, eventdata, handles)
handles.kd=str2double(get(hObject,'string'));
if isnan(handles.kd),set(hObject,'string','');handles.kd=0;end
guidata(hObject,handles);
% hObject    handle to k (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of k as text
%        str2double(get(hObject,'String')) returns contents of k as a double

% --- Executes during object creation, after setting all properties.
function k_CreateFcn(hObject, eventdata, handles)
% hObject    handle to k (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in imgsel.
function imgsel_Callback(hObject, eventdata, handles)
if isfield(handles,'cl'),handles.append=cluster_tool_callgui('cluster_tool_append');else handles.append=0;end
if isfield(handles,'img'),msgbox('Selected image will be used instead of previously selected image','','warn');end
handles.pwd=pwd;cd(handles.fdir);
[FileName,PathName]=uigetfile('*.img','Select image file.');
cd(handles.pwd);
if isempty(PathName)||isnumeric(PathName),return,end
handles.fdir=PathName;
set(handles.imagename,'string',['Image: ' PathName FileName]);
handles.img=[PathName FileName];
guidata(hObject,handles);
% hObject    handle to imgsel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in done.
function done_Callback(hObject, eventdata, handles)
s='[dat,volInfo,cl]=iimg_threshold(handles.img,''imgtype'',handles.imtype,''thr'',handles.threshold';
if get(handles.posneg,'value')==2&&~strcmp(handles.thrshtyp,'p'),handles.threshold=[-Inf -handles.threshold];
elseif get(handles.posneg,'value')==3,s=[s ',''abs'''];
end
% try
    if ~strcmp(handles.imtype,'data')
        s=[s ',''threshtype'',handles.thrshtyp'];
        if strcmp(handles.thrshtyp,'p')&&strcmp(handles.imtype,'t'),s=[s ',''df'',handles.dfd'];end
    end
    s=[s ',''thr'',handles.threshold'];
    if handles.kd,s=[s ',''k'',handles.kd'];end
    if get(handles.write,'value'),[FileName,PathName]=uiputfile('*.img','Select .img file to write');
        s=[s ',''outnames'',[PathName FileName]'];
    end
    s=[s ');'];
    eval(s)
% catch beep,disp('Error: Not all inputs were specified, please try again.'),return,end
if isfield(handles,'cl'),handles.cl=[handles.cl(:)' cl(:)'];else handles.cl=cl;end
guidata(hObject,handles);
answer=questdlg('Clusters obtained! Return to Acquire/Create Cluster screen?');
if strcmp(answer,'Yes'),ret_Callback(handles.ret,eventdata,handles);end

% hObject    handle to done (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in ret.
function ret_Callback(hObject, eventdata, handles)
if isfield(handles,'cl'),cluster_tool('cl',handles.cl,'fdir',handles.fdir);
else cluster_tool;end
close(handles.figure1);

% hObject    handle to ret (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in write.
function write_Callback(hObject, eventdata, handles)
% hObject    handle to write (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of write

% --- Executes on selection change in posneg.
function posneg_Callback(hObject, eventdata, handles)
% hObject    handle to posneg (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns posneg contents as cell array
%        contents{get(hObject,'Value')} returns selected item from posneg

% --- Executes during object creation, after setting all properties.
function posneg_CreateFcn(hObject, eventdata, handles)
set(hObject,'string','Positive|Negative|Both');
% hObject    handle to posneg (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function imagename_CreateFcn(hObject,eventdata,handles)
set(hObject,'string',['Image: ' handles.img]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore