text:
function F = retrieveR(level, o, H, L)
%------------------------------------------------------------------------------
%
% This function extracts a two-dimensional gridfunction F from H.
% F is supposed to be uniquely determined by the integer level and
% character o describing its type. F is defined on a rectangular
% domain.
% This function is a two-dimensional lifting scheme utility.
%
% F = 2D gridfunction of coefficients extracted from H.
%
% level = integer designated as the level of F.
%
% o = character, should be either 'a' or 'd',
%     describing the type of F:
%     'a' relates to approximation (coefficients) and
%     'd' relates to detail (coefficients)
%
% L = 2D integer array of bookkeeping, see function storeR.
%
% H = 1D array that functions as storage (heap). The coefficients of F
%     are extracted from H as a result of calling retrieveR.
%
% See also: storeR, retrieveQ
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: May 5, 2002.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
if isempty(L)
  if ~isempty(H)
    disp([' retrieveR - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveR - books empty but heap is not ')
  else
    F = [];
  end
else
  if isempty(H)
    disp([' retrieveR - ERROR at type ' o ' with level ' num2str(level)]);
    error(' retrieveR - books not empty but heap is ')
  else
    [nL, mL] = size(L);
    if mL ~= 6
      disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
      error(' retrieveR - books do not fit ')
    else
      jL = -1;
      j = 1;
      foundit = 0;
      while j <= nL && ~foundit
         if L(j, 1) == level
           switch o
               case 'a' , if L(j, 3) == 0 && L(j, 2) == 0
                            jL = j;
                            foundit = 1;
                          end
               case 'd' , if L(j, 3) == 1 && L(j, 2) == 0
                            jL = j;
                            foundit = 1;
                          end
               otherwise
               disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
               error(' retrieveR - unknown type of coefficients ')
           end
         end
         j = j + 1;
      end
      if jL == -1
        disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
        error(' retrieveR - no such coefficients in heap')
      else
        nF = L(jL, 4);
        if nF < 1
          disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
          error(' retrieveR - unvalid 1st dimension of target')
        end
        mF = L(jL, 5);
        if mF < 1
          disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
          error(' retrieveR - unvalid 2nd dimension of target')
        end
        heaptr = L(jL, 6);
        if heaptr < 2
          disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
          error(' retrieveR - bookkeeping error ')
        else
          [nH, mH] = size(H);
          if mH ~=1
            disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
            error(' retrieveR - heap should be column vector')
          else
            beginptr = heaptr-nF*mF;
            if heaptr-1 > nH
              disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
              error(' retrieveR - heap not that large')
            elseif beginptr < 1
              disp([' retrieveR - ERROR at type ' o ' level ' num2str(level)]);
              error(' retrieveR - heap not that large, dimensions?')
            else
              F = reshape(H(beginptr:(heaptr-1)), nF, mF);
            end
          end
        end
      end
    end
  end
end
%------------------------------------------------------------------------------

author: Sable


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

addpath(genpath('...\cst api'));

cst = actxserver('CSTStudio.application');
mws = cst.invoke('NewMWS');

Geometry = 'mm';
Frequency = 'GHz';
Time = 'ns';
TemperatureUnit = 'Kelvin';
Voltage = 'V';
Current = 'A';
Resistance = 'Ohm';
Conductance  = 'S';
Capacitance = 'PikoF';
Inductance = 'NanoH';

CstDefineUnits(mws,Geometry, Frequency, Time, TemperatureUnit, Voltage, Current, Resistance, Conductance, Capacitance, Inductance)
CstMeshInitiator(mws)

CstDefineFrequencyRange(mws,0.5,4)
minfrequency = 0.5;
x='expanded open';
CstDefineOpenBoundary(mws,minfrequency,x,x,x,x,x,x)
XminSpace = 0;
XmaxSpace = 0;
YminSpace = 0;
YmaxSpace = 0;
ZminSpace = 0;
ZmaxSpace = 0;
CstDefineBackroundMaterial(mws,XminSpace,XmaxSpace, YminSpace, YmaxSpace, ZminSpace, ZmaxSpace)

Name = 'Groundplane';
component = 'component1';
material = 'PEC';
Xrange = [-40 40];
Yrange = [-40 40];
Zrange = [0 2];

Cstbrick(mws, Name, component, material, Xrange, Yrange, Zrange)

Name = 'Conductor';
OuterRadius = 2;
InnerRadius = 0;
Xcenter = 0;
Ycenter = 0;
Zrange = [8 38];

Cstcylinder(mws, Name, component, material, 'Z', OuterRadius, InnerRadius, Xcenter, Ycenter, Zrange)

PortNumber = 1;
SetP1 = [0 0 2];
SetP2 = [0 0 8];

CstDiscretePort(mws,PortNumber,SetP1,SetP2)

monitorindex = 2.25;
CstDefineEfieldMonitor(mws,strcat('e-field', num2str(monitorindex)), monitorindex)
CstDefineFarfieldMonitor(mws,strcat('Farfield',num2str(monitorindex)), monitorindex)
CstSaveProject(mws)

CstDefineTimedomainSolver(mws,-30)

exportpath = '...\cst api\NameofyourFile';%NameofyourFile should be the same with the name of the file of line 90 bellow

format = 'db';

CstExportTouchstone(mws,exportpath, format)
data = read(rfdata.data,'NameofyourFile.s1p');
freq = data.Freq;
s_params = extract(data,'S_PARAMETERS',50);
spar = s_params(1,:);
s11 = 20*log10(sqrt(real(spar).^2 + imag(spar).^2));
plot(freq,s11); grid on;

repo: CST-MATLAB-API


text:
function [X]=tt_minres_selfprec2(A,  eps, varargin)
%Computation of the approximate TT-matrix inverse using self-prec method
%   [X]=TT_MINRES_SELFPREC2(A,EPS,OPTIONS) Computation of the approximate
%   TT-matrix inverse using Saad self-prec method. Options are provided
%   in form 'PropertyName1',PropertyValue1,'PropertyName2',PropertyValue2
%   and so on. The parameters are set to default (in brackets in the
%   following) The list of option names and default values are:
%       o matvec - type of the local matvec [ {mm+compr} | mmk2 ]
%       o max_rank - maximal TT-rank bound [1000]
%       o prec_type - left or right inversion [ {left} | right ]
%       o maxit - maximal number of iterations [10]
%       o tol - the requested inversion tolerance [EPS]
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

% matvec='mmk2';
matvec='mm+compr';
max_rank=1000;
prec_type='left';
tol=eps;
maxit=10;

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'matvec'
            matvec=varargin{i+1};
        case 'max_rank'
            max_rank=lower(varargin{i+1});
        case 'prec_type'
            prec_type=varargin{i+1};
        case 'maxit'
            maxit=varargin{i+1};
        otherwise
            error('Unrecognized option: %s\n',varargin{i});
    end
end

%d=max(size(A));
d=ndims(A);
ns=A.n;
%ns=tt_size(A);

I=tt_eye(ns, d); I=tt_matrix(I);
if (strcmp(prec_type, 'left'))
   % X=tt_scal(tt_transp(A), 1e-15);
   X=(1e-15)*conj(A');
else
  %  X=tt_scal(I, 1e-150);
   X=1e-150*I;
end;
normf=norm(I);
%normf=sqrt(tt_dot(tt_mat_to_vec(I), tt_mat_to_vec(I)));
err=1;
err_old=2;
sp=0;

for iout=1:maxit
    if (strcmp(matvec, 'mmk2'))
        resid=tt_mmk2(A, X, eps, max_rank);
    else
        %resid=tt_mm(A,X);
        resid=round(A*X,eps,max_rank*2);
        %resid=tt_mat_compr(resid, eps,max_rank*2);
    end;
    %resid=ttm_add(I, tt_scal(resid, -1));
    resid=I-resid;

    if (iout>1)||(strcmp(prec_type, 'left'))
        if (strcmp(matvec, 'mmk2'))
            Xresid=tt_mmk2(X, resid, eps,max_rank);
        else
            %Xresid=tt_mm(X,resid);
            Xresid=round(X*resid,eps,max_rank*2);
            %Xresid=tt_mat_compr(Xresid, eps,max_rank*2);

        end;
    else
        Xresid=resid;
    end;

    max_xrrank=max(rank(Xresid));

    if (strcmp(prec_type, 'left')) % we use left preconditioner (igmres)
        if (strcmp(matvec, 'mmk2'))
            w=tt_mmk2(A, Xresid, eps,max_rank);
            w=tt_mmk2(X, w, eps,max_rank);
        else
            w=round(A*Xresid,eps,max_rank);
            w=round(X*w,eps);
            %w=tt_mm(A,Xresid);
            %w=tt_mat_compr(w, eps,max_rank);
            %w=tt_mm(X,w);
%             w=tt_mat_compr(w, eps);
        end;
        %w = tt_mat_to_vec(w);
        %beta=sqrt(tt_dot(tt_mat_to_vec(Xresid),tt_mat_to_vec(Xresid)));
        beta=norm(Xresid);
        wr=dot(w,Xresid);
        %wr = tt_dot(w, tt_mat_to_vec(Xresid));
        ww=dot(w,w);

        H=zeros(2,1);
        H(1,1)=wr/(beta^2);
        H(2,1)=sqrt(ww/(beta^2)-wr^2/(beta^4));
        rhs=H(1,1)*beta;

        y = rhs/(H'*H);
        err = err*norm(H*y-[beta; 0])/beta;
        y=y/beta;

%         y = tt_dot(w,w);
%         if (y~=0)
%             y = tt_dot(w, tt_mat_to_vec(Xresid))/y;
%         else
%             y=1;
%         end;
    else % Use right preconditioner and fgmres
        if (strcmp(matvec, 'mmk2'))
            w=tt_mmk2(A, Xresid, eps,max_rank);
        else
            w=A*Xresid;
            w=round(w,eps);
%             w=tt_mat_compr(w, eps);
        end;
        w=tt_tensor(w);
        %w = tt_mat_to_vec(w);
        beta=norm(resid);
        %beta=sqrt(tt_dot(tt_mat_to_vec(resid),tt_mat_to_vec(resid)));
        wr=dot(w,tt_tensor(resid));
        ww=dot(w,w);
        %wr = tt_dot(w, tt_mat_to_vec(resid));
        %ww = tt_dot(w,w);

        H=zeros(2,1);
        H(1,1)=wr/(beta^2);
        H(2,1)=sqrt(ww/(beta^2)-wr^2/(beta^4));
        rhs=H(1,1)*beta;

        y = rhs/(H'*H);
        err = norm(H*y-[beta; 0])/normf;
        y=y/beta;

%         y2=ww;
%         if (y~=0)
%             y2 = wr/y2;
%         else
%             y2=1;
%         end;
    end;

    max_wrank=max(rank(w));

    %Xresid=tt_mat_to_vec(Xresid);
    Xresid=tt_tensor(Xresid);
    X=tt_tensor(X);
    %X=tt_mat_to_vec(X);
    %if (err<err_old)
        %X = tt_axpy(1, X, y, Xresid, eps, max_rank);
        X=round(X+y*Xresid,eps,max_rank);

        %err_old=err;
    %else
    %    sp=sp+1;
    %end;
    %max_xrank=max(tt_ranks(X));
    max_xrank=max(rank(X));
    %X = tt_vec_to_mat(X, ns, ns);
    X=tt_matrix(X,ns,ns);
%     if (strcmp(matvec, 'mmk2'))
%         resid=tt_mmk2(A, X, eps, max_rank);
%     else
%         resid=tt_mm(A, X);
% %         resid=tt_mat_compr(resid, eps);
%     end;
%     resid=ttm_add(I, tt_scal(resid, -1));
%     resid=tt_mat_to_vec(resid);
%
%     err_real=sqrt(tt_dot(resid,resid))/normf;
    err_real=0;
    fprintf('iter [%d], res_real=%3.3e, x_rank=%d, Xr_rank=%d, w_rank=%d, err=%3.3e\n', iout, err_real, max_xrank, max_xrrank, max_wrank, err);
    if (err<tol)||(sp>0)
        break;
    end

end
%keyboard;

author: oseledets

repo: TT-Toolbox


text:
  function b = ir_mid3_pic1(x, file, varargin)
%|function b = ir_mid3_pic1(x, file, varargin)
%|
%| write 'mid3' type picture to .png file
%|
%| in
%|	x	[nx ny nz]	3D image
%|	file	char		file name
%|
%| option
%|	'clim'			display limit, default [800 1200]
%|
%| out
%|	b	[nx+nz ny+nz]	2D 'mid3' image, uint8
%|
%| 2012-10-01, Jeff Fessler

if nargin < 1, ir_usage, end
if streq(x, 'test'), ir_mid3_pic1_test, return, end

arg.clim = [800 1200];
arg = vararg_pair(arg, varargin);

y = jf_mip3(x, 'type', 'mid');

clim = arg.clim;
b = 255 * (y - clim(1)) / (clim(2) - clim(1));
b = min(b,255);
b = max(b,0);
b = uint8(b);

if exist(file, 'file')
	warn('file "%s" exists already.', file)
	yn = input('over-write? y|n ', 's');
	if ~streq(yn, 'y', 1)
		error('file name')
	end
end
imwrite(b', file); % trick: transpose

% ir_mid3_pic1_test
function ir_mid3_pic1_test

ig = image_geom('nx', 64, 'nz', 30, 'dx', 1, 'dz', 2);
x = ellipsoid_im(ig, [], 'hu_scale', 1000);

file = 'test.png';
b = ir_mid3_pic1(x, file, 'clim', [950 1050]);

%h = im(y, arg.clim);
%z = get(h, 'cdata');
%minmax(b)

tmp = imread(file)'; % trick
jf_equal(tmp, b)

% !xv test.png &

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = process_matlab_eval( varargin )
% PROCESS_MATLAB_EVAL: Evaluate a matlab script.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Run Matlab command';
    sProcess.FileTag     = 'matlab';
    sProcess.Category    = 'Filter';
    sProcess.SubGroup    = 'Pre-process';
    sProcess.Index       = 80;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/TutUserProcess#Alternative';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = -1;
    sProcess.processDim  = [];
    % === Matlab command
    sProcess.options.matlab.Comment = 'Type your Matlab script below: ';
    sProcess.options.matlab.Type    = 'textarea';
    sProcess.options.matlab.Value   = ['% Available variables: Data, TimeVector' 10 10 ...
                                       'Data = Data;' 10];
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data', 'raw'};
end

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

%% ===== RUN =====
function sInput = Run(sProcess, sInput) %#ok<DEFNU>
    % Available variables
    Data = sInput.A;
    TimeVector = sInput.TimeVector;
    Comment = sInput.Comment;
    % Evaluate Matlab code
    try
        eval(sProcess.options.matlab.Value);
    catch
        % Get error message
        e = lasterror();
        if ~isempty(e)
            % Remove HTML tags
            errMsg = str_striptag(e.message);
            % Strip "Error:"
            errMsg = strrep(errMsg, 'Error: ', '');
            % Strip first line
            iNewLine = find(errMsg == 10);
            if ~isempty(iNewLine)
                errMsg = errMsg(iNewLine(1)+1:end);
            end
            errMsg = ['Error while executing Matlab command:' 10 errMsg];
        else
            errMsg = 'Error while executing Matlab command.';
        end
        % Report error
        bst_report('Error', sProcess, sInput, errMsg);
        sInput = [];
        return;
    end
    % Check if time was modified
    if (length(TimeVector) ~= length(sInput.TimeVector))
        % Not allowed for only a few sensors
        if strcmpi(sInput.FileType, 'data') && isfield(sProcess.options, 'sensortypes') && isfield(sProcess.options.sensortypes, 'Value') && ~isempty(sProcess.options.sensortypes.Value)
            sInput = [];
            bst_report('Error', sProcess, sInput, ...
                ['This process can modify the time definition only if all the channels are processed at once.' 10 ...
                 'Either leave the option "Sensor types or names" empty to process all the channels at once,' 10 ...
                 'or keep the number of time points of the original file.']);
            return;
        % Not allowed for raw files
        elseif strcmpi(sInput.FileType, 'raw')
            sInput = [];
            bst_report('Error', sProcess, sInput, 'This process cannot modify the time definition of a continuous file.');
            return;
        end
    end
    % Report results
    sInput.A = Data;
    sInput.TimeVector = TimeVector;
    sInput.Comment = Comment;
    % Do not keep the Std field in the output
    if isfield(sInput, 'Std') && ~isempty(sInput.Std)
        sInput.Std = [];
    end
end

author: brainstorm-tools


text:
function K=interpc(coef,lat)

j=fix(lat/15);
if j<1
    K=coef(1); return;
end

if j>4
    K=coef(5); return;
end

K=coef(j)*(1-lat/15+j)+coef(j+1)*(lat/15-j);

return

repo: GINav


text:
function B = colFirstNonZero(A)
    % Get first non zero element of each column of a matrix
    % SYNTAX:
    % B = colFirstNonZero(A)
    [~, c] = max( A ~=0, [], 1 );
    d = (0 : (size(A, 2) - 1)) * size(A, 1) + c;
    B = A(d);
end

author: goGPS-Project


text:
function [min_val, max_val] = arange(X)
%ARANGE Returns the range (min and max) of an entire array.
% Usage:
%   R = arange(X)
%   [min_val, max_val] = arange(X)
%
% See also: alims, amin, amax

min_val = min(X(:));
max_val = max(X(:));

if nargout < 2
    min_val = [min_val, max_val];
end
end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
classdef ProductComputerForTensor < ProductComputer

    properties
    end

    methods (Access = public)

        function obj = ProductComputerForTensor(C,e)
            obj.generate(C,e)
        end
    end

    methods (Access = protected)

        function computeProduct(obj)
            C = obj.fourthOrder.getValue();
            e = obj.secondOrder.getValue();
            d = obj.secondOrder.getDimension();
            s = zeros(size(e));
            for i = 1:d
                for j = 1:d
                    for k = 1:d
                        for l = 1:d
                            cijkl = C(i,j,k,l);
                            ekl = e(k,l);
                            s(i,j) = s(i,j) + cijkl*ekl;
                        end
                    end
                end
            end
            obj.secondOrderOut.setValue(s)
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function db=lpcar2db(ar,np)
%LPCAR2DB LPC: Convert AR coefs to power spectrum in dB DB=(AR)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcar2db.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
if nargin<2 np=p1-1; end
ff=rfft(ar.',2*np+2).';
db=-10*log10(real(ff.*conj(ff)));
if nargout==0
   plot((0:np+1)/(2*np+2),db.');
   ylabel('dB');
end

author: decouples


text:
classdef SuperEllipseParamsRelator < handle

   methods (Access = public, Static)

       function mx = mx(xi,rho,q)
            c = SuperEllipseParamsRelator.c();
            n = (1-rho).*tan(xi);
            d = c(q);
            mx = sqrt(n./d);
       end

       function my = my(xi,rho,q)
            c = SuperEllipseParamsRelator.c();
            n = (1-rho);
            d = c(q).*tan(xi);
            my = sqrt(n./d);
       end

       function rho = rho(mx,my,q)
            c = SuperEllipseParamsRelator.c();
            rho =  1 - c(q).*mx.*my;
       end

       function xi = xi(mx,my)
           xi = atan(mx./my);
       end

       function c = c()
          c = @(q) gamma(1 + 1./q).^2./gamma(1 + 2./q);
       end

       function my = myFromMxAndTxi(mx,xi)
          my = mx./tan(xi);
       end

       function mx = mxFromMyAndTxi(my,xi)
          mx = my.*tan(xi);
       end

       function rho = rhoFromMxAndTxi(mx,xi,q)
           sE  = SuperEllipseParamsRelator();
           my  = sE.myFromMxAndTxi(mx,xi);
           rho = sE.rho(mx,my,q);
       end

       function rho = rhoFromMyAndTxi(my,xi,q)
           sE  = SuperEllipseParamsRelator();
           mx  = sE.mxFromMyAndTxi(my,xi);
           rho = sE.rho(mx,my,q);
       end

       function xi = xiFromMxAndRho(mx,rho,q)
           sE  = SuperEllipseParamsRelator();
           c  = sE.c();
           xi = atan(c(q).*mx.^2./(1-rho));
       end

       function xi = xiFromMyAndRho(my,rho,q)
           sE  = SuperEllipseParamsRelator();
           c  = sE.c();
           xi = atan((1-rho)./(c(q).*my.^2));
       end

       function xiUB = xiUB(rho)
           sE  = SuperEllipseParamsRelator();
           qMax = 32;
           qMin = 2;
           mxUB = 0.99;
           myLB = 0.01;
           xi1 = sE.xiFromMxAndRho(mxUB,rho,qMax);
           xi2 = sE.xiFromMyAndRho(myLB,rho,qMin);
           xiUB = min(xi1,xi2);
       end

       function xiLB = xiLB(rho)
           sE  = SuperEllipseParamsRelator();
           qMax = 32;
           qMin = 2;
           mxLB = 0.01;
           myUB = 0.99;
           xi1 = sE.xiFromMxAndRho(mxLB,rho,qMin);
           xi2 = sE.xiFromMyAndRho(myUB,rho,qMax);
           xiLB = min(xi1,xi2);
       end

       function [mXb,mYb] = updateMxMyWithOtherQ(mXa,mYa,q,qNew)
            sE   = SuperEllipseParamsRelator;
            rho  = sE.rho(mXa,mYa,q);
            xi   = sE.xi(mXa,mYa);
            mXb  = sE.mx(xi,rho,qNew);
            mYb  = sE.my(xi,rho,qNew);
       end

   end

end


text:
function F = diff(F, dim, n)
%DIFF   Componentwise derivative of a DISKFUNV.
%   DIFF(F) is the derivative of each component of F in
%   x-direction.
%
%   DIFF(F, DIM) is the first derivative of F along the
%   dimension DIM.
%     DIM = 1 (default) is the derivative in the x-direction.
%     DIM = 2 is the derivative in the y-direction.
%
%   DIFF(F, DIM, N) is the Nth derivative each component of F
%   along the dimension specified.
%
% See also CURL and DIV

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

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

% Default to x-derivative:
if ( nargin == 1 || isempty(dim) )
    dim = 1;
end

% Default to the first derivative:
if ( nargin < 3 )
    n = 1;
end

% Differentiate each component:
F.components{1} = diff(F.components{1}, dim, n);
F.components{2} = diff(F.components{2}, dim, n);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [y,am]=anabpsk(N,Ncomp,f0);
%ANABPSK Binary Phase Shift Keying (BPSK) signal.
% 	[Y,AM]=ANABPSK(N,NCOMP,F0) returns a succession of complex
%	sinusoids of NCOMP points each, with a normalized frequency F0 and
%	an amplitude equal to -1 or +1, according to a discrete	uniform
%	law. Such signal is only 'quasi'-analytic.
%
%	N     : number of points
%	NCOMP : number of points of each component (default: N/5)
% 	F0    : normalized frequency.              (default: 0.25)
% 	Y     : signal
% 	AM    : resulting amplitude modulation     (optional).
%
%	Example :
%	 [signal,am]=anabpsk(300,30,0.1); clf; figure(gcf);
%  	 subplot(211); plot(real(signal)); subplot(212); plot(am);
%
%	See also ANAFSK, ANAQPSK, ANAASK.

%	O. Lemoine - June 1995, F. Auger - August 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('The number of parameters must be at least 1.');
elseif (nargin == 1),
 Ncomp=round(N/5); f0=0.25;
elseif (nargin == 2),
 f0=0.25;
end;

if (N <= 0),
 error('The signal length N must be strictly positive' );
elseif (f0<0)|(f0>0.5),
 error('f0 must be between 0 and 0.5');
end;

MatlabVersion=version; MatlabVersion=str2num(MatlabVersion(1));
if (MatlabVersion==4), rand('uniform'); end;

m=ceil(N/Ncomp);
jumps=2.0*round(rand(m,1))-1.0;
am=kron(jumps,ones(Ncomp,1)); am=am(1:N,1);
y=am.*fmconst(N,f0,1);

author: HeLiangHIT


text:
function res = UpSum(p)
%UPSUM        Rounded upwards result of sum(p)
%
%   resD = DownSum(p)
%
%On return, resD is sum(p) rounded upwards, also in the presence
%  of underflow. Input vector p may be single or double precision.
%
%Adapted from Algorithm 7.5 from
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%Requires (4m+4k+4)n flops for m,k executions of the repeat-until loops
%  in the calls of TransformK, respectively.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
% modified 05/09/09     S.M. Rump  rounding to nearest, complex input
%

  if ~isreal(p)
    res = complex(UpSum(real(p)),UpSum(imag(p)));
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if isa(p,'double')
    nmax = 2^26-2;                          % nmax = 67,108,864
  else
    nmax = 2^12-2;                          % nmax = 4,094
  end

  if length(p)<=nmax                        % not huge dimension
    [res,R,p,sigma,Ms] = TransformK(p,0);   % s-res = R+sum(p)
    delta = TransformK(p,R,sigma,Ms);       % delta faithful rounding of s-res
  else                                      % huge dimension
    [res,tau1,tau2,tau,p] = AccSumHugeN(p); % s=tau1+tau2+sum(tau)+sum(p)
    delta = AccSumHugeN([tau1;tau2;tau(:);p(:);-res]);
  end
  if delta>0                                % s > res
    res = succ(res);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function reset_SPMcfg()
% resets columns in SPMcfg by removing all non-intercept nuisance covariates.
% runs on the SPMcfg.mat file in the current directory

    if exist([pwd filesep 'SPMcfg.mat']) == 2
        load SPMcfg
    else
        warning(['No SPMcfg.mat file found in current directory:' pwd])
    end

    nsess=length(Sess);
    xX.X(:,xX.iB(1:end - nsess)) = [];

    if length(xX.Xnames) > size(xX.X,2)
        xX.Xnames(xX.iB(1:end - nsess)) = [];
    end

    xX.iB = size(xX.X,2) - nsess + 1: size(xX.X,2);

    F_iX0.iX0 = xX.iB;
    save SPMcfg F_iX0 SPMid Sess VY xGX xM xX xsDes

    fprintf(1,'\n SPMcfg xX modified to remove nuisance covariates and saved.\n')

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function rphiList=Cart2Hypersphere(x)
%%HYPERSPHERE2CART Given a point in n-dimensional Cartesian coordinates
%                  convert the point to n-dimensional spherical coordinates
%                  (hyperspherical coordinates.
%
%INPUTS: x An nXN matrix of the values converted into n-dimensional
%          Cartesian coordinates. n>=2.
%
%OUTPUTS: rPhiList The nXN values convert from Cartesian coordinates to
%                 hyperspherical coordinates. Each value is a range r
%                 (0<=r<Inf) followed by (n-1) angular values. The first
%                 (n-2) values are bound from 0 to pi and the final one can
%                 span the entire circle (-pi to pi).
%
%Hyperspherical coordinates are derived in [1]. The definition here is
%nearly identical, except the final two coordinates are switched (the final
%term is all sines). Hyperspherical coordinates arise in tracking, for
%example, when generating false measurements in gates in simulations, such
%as in [2]. Note that in [2], yet another definition of the system is used.
%
%The definition here is
% x(1)=cos(phi(1));
% x(2)=sin(phi(1))*cos(phi(2));
% x(3)=sin(phi(1))*sin(phi(2))*cos(phi(3));
% ...
% x(n-1)=sin(phi(1))*...*sin(phi(n-2))*cos(phi(n-1));
% x(n)=sin(phi(1))*...*sin(phi(n-2))*sin(phi(n-1));
%
%The definition of the angles differs in 3D from those use in the
%Cart2Sphere function.
%
%REFERENCES:
%[1] L. E. Blumenson, "A derivation of n-dimensional spherical
%    coordinates," The American Mathematical Monthly, no. 1, pp. 63-66,
%    Jan. 1960.
%[2] H. Sun and M. Farooq, "Note on the generation of random points
%    uniformly distributed in hyper-ellipsoids," in Proceedings of the
%    Fifth International Conference on Information Fusion, Annapolis, MD,
%    8-11 Jul. 2002, pp. 489-496.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(x);

phiList=zeros(n-1,1);
r=norm(x);

for curDim=1:(n-1)
    normVal=norm(x(curDim:n));

    if(normVal==0)
        if(x(curDim)>=0)
            phiList(curDim)=0;
        else
            phiList(curDim)=pi;
        end
    else
        phiList(curDim)=acos(x(curDim)/normVal);
    end
end

if(x(n)<0)
    phiList(n-1)=2*pi-phiList(n-1);
end
phiList(n-1)=wrapRange(phiList(n-1),-pi,pi);

rphiList=[r;phiList];
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 yp = p30_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P30_FUN evaluates the function for problem P30.
%
%  Discussion:
%
%    1 equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  csum = 0.0;
  for i = 1 : 19
    csum = csum + ( i )^(4.0/3.0);
  end
  pprime = 0.0;
  for i = 1 : 19
    pprime = pprime + ( 4.0 / 3.0 ) * r8_sign ( t - i ) ...
      * ( abs ( t - i ) )^(1.0/3.0);
  end
  yp(1) = pprime * y(1) / csum;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef form_acp < mp.form_ac
%MP.FORM_ACP  MATPOWER Formulation class for AC polar voltage formulations
%   Each concrete Network Model Element class must inherit, at least
%   indirectly, from both MP.NM_ELEMENT and MP.FORM.
%
%   Subclass of MP.FORM_AC.
%   MP.FORM provides properties and methods related to the specific
%   formulation (e.g. DC version, AC polar power version, etc.)
%
%   Properties
%       (model parameters inherited from MP.FORM_AC)
%
%   Methods
%       form_name() - returns string w/name of formulation ('AC-polar formulation')
%       form_tag() - returns string w/short label for formulation ('acp')

%   MATPOWER
%   Copyright (c) 2019-2020, 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 name = form_name(obj)
            name = 'AC-polar';
        end

        function tag = form_tag(obj)
            tag = 'acp';
        end

        function vtypes = model_vvars(obj)
            vtypes = {'va', 'vm'};
        end

        function [Iva, Ivm] = port_inj_current_jac(obj, ...
                n, v_, Y, M, invdiagvic, diagSlincJ)
            % [Iva, Ivm] = obj.port_inj_current_jac(...)

            %% intermediate terms
            diagv = sparse(1:n, 1:n, v_, n, n);
            C = invdiagvic * (diagSlincJ - conj(M * diagv));
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);

%             %% linear current term
%             Ivm = Y * diagv;
%             Iva = 1j * Ivm;
%             Ivm = Ivm * D;
%
%             %% + current from linear power term
%             Iva = Iva + 1j * C;
%             Ivm = Ivm - C * D;

            A = Y * diagv;
            Iva = 1j * (A + C);
            Ivm = (A - C) * D;
        end

        function [Ivava, Ivavm, Ivmvm] = port_inj_current_hess_v(obj, x_, lam, v_, z_, diaginvic, Y, M, diagSlincJ, dlamJ)
            % [Ivava, Ivavm, Ivmvm] = obj.port_inj_current_hess_v(x_, lam)
            % [Ivava, Ivavm, Ivmvm] = obj.port_inj_current_hess_v(x_, lam, sysx)
            % [Ivava, Ivavm, Ivmvm] = obj.port_inj_current_hess_v(x_, lam, sysx, idx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, v_, z_, diaginvic, Y, M, diagSlincJ, dlamJ)

            if nargin < 10
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                [Y, M, N, s] = obj.get_params(idx, {'Y', 'M', 'N', 's'});
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                %% compute linear power injections
                if isempty(z_)
                    Slin = M*v_ + s;
                else
                    Slin = M*v_ + N*z_ + s;
                end

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diaginvic = sparse(1:ni, 1:ni, 1 ./ conj(vi_), ni, ni);
                if isempty(idx)     %% all ports
                    diagSlincJ = sparse(1:n, 1:n, conj(Slin), n, n);
                    dlamJ      = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    diagSlincJ = sparse(1:ni, idx, conj(Slin), ni, n);
                    dlamJ      = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            diagv  = sparse(1:n, 1:n, v_, n, n);
            A = diaginvic * diagSlincJ;
            B = diaginvic * conj(M);
            % B2 = B * conj(diagv);
            % C = A - B2;
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);
            % E = diaginvic * conj(N);
            F = dlamJ.';
            G = F * B * conj(diagv);    %% F * B2
            dBtlam = sparse(1:n, 1:n, B.' * lam, n, n);
            H = dBtlam * conj(diagv);
            K = (F * A).';
            GG = G + G.';
            LL = GG - H - K;
            MM = D * (2*K - GG) * D;

            %% linear current term
            dYtlam = sparse(1:n, 1:n, Y.' * lam, n, n);
            Ivava = -dYtlam * diagv;
            Ivavm = -j * Ivava * D;

            %% + current from linear power term
            Ivava = Ivava + LL;
            Ivavm = Ivavm + 1j * LL.' * D;
            Ivmvm = MM;
        end

        function [Ivazr, Ivazi, Ivmzr, Ivmzi] = port_inj_current_hess_vz(obj, x_, lam, v_, z_, diaginvic, N, dlamJ)
            % [Ivazr, Ivazi, Ivmzr, Ivmzi] = obj.port_inj_current_hess_vz(x_, lam)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, sysx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, sysx, idx)
            % [...] = obj.port_inj_current_hess_vz(x_, lam, v_, z_, diaginvic, N, dlamJ)

            if nargin < 8
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                N = obj.get_params(idx, 'N');
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diaginvic = sparse(1:ni, 1:ni, 1 ./ conj(vi_), ni, ni);
                if isempty(idx)     %% all ports
                    dlamJ = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    dlamJ = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);
            E = diaginvic * conj(N);
            NN = dlamJ.' * E;

            %% current from linear power term
            Ivazr = 1j * NN;
            Ivazi = NN;
            Ivmzr = -D * NN;
            Ivmzi = -1j * Ivmzr;
        end

        function [Sva, Svm] = port_inj_power_jac(obj, ...
                n, v_, Y, M, diagv, diagvi, diagIlincJ)
            % [Sva, Svm] = obj.port_inj_power_jac(...)

            %% intermediate terms
            A = diagvi * diagIlincJ;
%             B = diagvi * conj(Y);
%             C = B * conj(diagv);
            C = diagvi * conj(Y * diagv);
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);

%             %% linear power term
%             Svm = M * diagv;
%             Sva = 1j * Svm;
%             Svm = Svm * D;
%
%             %% + power from linear current term
%             Sva = Sva + 1j * (A - C);
%             Svm = Svm + (A + C) * D;

            Svm = M * diagv + A;
            Sva = 1j * (Svm - C);
            Svm = (Svm + C) * D;
        end

        function [Svava, Svavm, Svmvm] = port_inj_power_hess_v(obj, x_, lam, v_, z_, diagvi, Y, M, diagIlincJ, dlamJ)
            % [Svava, Svavm, Svmvm] = obj.port_inj_power_hess_v(x_, lam)
            % [Svava, Svavm, Svmvm] = obj.port_inj_power_hess_v(x_, lam, sysx)
            % [Svava, Svavm, Svmvm] = obj.port_inj_power_hess_v(x_, lam, sysx, idx)
            % [...] = obj.port_inj_power_hess_v(x_, lam, v_, z_, diagvi, Y, M, diagIlincJ, dlamJ)

            if nargin < 10
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                [Y, L, M, i] = obj.get_params(idx, {'Y', 'L', 'M', 'i'});
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                %% compute linear current injections
                if isempty(z_)
                    Ilin = Y*v_ + i;
                else
                    Ilin = Y*v_ + L*z_ + i;
                end

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diagvi = sparse(1:ni, 1:ni, vi_, ni, ni);
                if isempty(idx)     %% all ports
                    diagIlincJ = sparse(1:n, 1:n, conj(Ilin), n, n);
                    dlamJ      = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    diagIlincJ = sparse(1:ni, idx, conj(Ilin), ni, n);
                    dlamJ      = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            diagv  = sparse(1:n, 1:n, v_, n, n);

            A = diagvi * diagIlincJ;
            B = diagvi * conj(Y);
            C = B * conj(diagv);
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);
            % E = diagvi * conj(L);
            F = dlamJ.';
            G = F * C;
            dBtlam = sparse(1:n, 1:n, B.' * lam, n, n);
            H = conj(diagv) * ((F*B).' - dBtlam);
            K = F * (C - A);

            %% linear power term
            dMtlam = sparse(1:n, 1:n, M.' * lam, n, n);
            Svava = -dMtlam * diagv;
            Svavm = -j * Svava * D;

            %% + power from linear current term
            Svava = Svava + H + K;
            Svavm = Svavm + 1j * (H - K).' * D;
            Svmvm = D * (G + G.') * D;
        end

        function [Svazr, Svazi, Svmzr, Svmzi] = port_inj_power_hess_vz(obj, x_, lam, v_, z_, diagvi, L, dlamJ)
            % [Svazr, Svazi, Svmzr, Svmzi] = obj.port_inj_power_hess_vz(x_, lam)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, sysx)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, sysx, idx)
            % [...] = obj.port_inj_power_hess_vz(x_, lam, v_, z_, diagvi, L, dlamJ)

            if nargin < 8
                sysx = v_;
                if nargin < 5
                    idx = []
                else
                    idx = z_;
                end
                L = obj.get_params(idx, 'L');
                [v_, z_, vi_] = obj.x2vz(x_, sysx, idx);

                n  = length(v_);    %% number of all port voltages
                ni = length(vi_);   %% number of selected port voltages
                diagvi = sparse(1:ni, 1:ni, vi_, ni, ni);
                if isempty(idx)     %% all ports
                    dlamJ = sparse(1:n, 1:n, lam, n, n);
                else                %% selected ports
                    dlamJ = sparse(1:ni, idx, lam, ni, n);
                end
            else
                n  = length(v_);    %% number of all port voltages
                ni = length(lam);
            end

            %% intermediate terms
            D = sparse(1:n, 1:n, 1 ./ abs(v_), n, n);
            E = diagvi * conj(L);
            LL = dlamJ.' * E;
            M = D * LL;

            %% power from linear current term
            Svazr = 1j * LL;
            Svazi = LL;
            Svmzr = M;
            Svmzi = -1j * M;
        end

        function [va, vm] = aux_data_va_vm(obj, ad)
            va = ad.va;
            vm = ad.vm;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function f2 = p01_f2 ( x )

%*****************************************************************************80
%
%% P01_F2 evaluates the second derivative for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  f2 = 2.0;

  return
end


text:
clear all;
close all;
clc;
spx.cluster.ssc.util.bench_subspace_preservation(@ssc_omp, 'ssc_omp');

author: indigits

repo: sparse-plex


text:
function GeometryEstimation(wkdir, dataset, matcher, estimator)
% Matching descriptors and save results
disp('Running FM estimation...');

dataset_dir = [wkdir 'Dataset/' dataset '/'];
matches_dir = [wkdir 'Matches/' dataset '/'];

results_dir = [wkdir 'Results/' dataset '/'];
if exist(results_dir, 'dir') == 0
    mkdir(results_dir);
end

pairs_gts = dlmread([dataset_dir 'pairs_with_gt.txt']);
pairs_which_dataset = importdata([dataset_dir 'pairs_which_dataset.txt']);

pairs = pairs_gts(:,1:2);
l_pairs = pairs(:,1);
r_pairs = pairs(:,2);
F_gts = pairs_gts(:,3:11);

load([matches_dir matcher '.mat']);
Results = Matches;
num_pairs = size(pairs,1);

for idx = 1 : num_pairs
    l = l_pairs(idx);
    r = r_pairs(idx);

    Results{idx}.dataset = dataset;
    Results{idx}.subset = pairs_which_dataset{idx};
    Results{idx}.l = l;
    Results{idx}.r = r;
    Results{idx}.F_gt = reshape(F_gts(idx,:), 3, 3)';

    X_l = Results{idx}.X_l;
    X_r = Results{idx}.X_r;

    F_hat = [];
    inliers = [];
    status = 3; % 0 stands for good, others are bad estimations.

    try
        [F_hat, inliers, status] = estimateFundamentalMatrix(X_l, X_r, 'Method','RANSAC', 'NumTrials', 2000);
    catch
        disp('Estimation Crash');
    end

    Results{idx}.F_hat = F_hat;
    Results{idx}.inliers = inliers;
    Results{idx}.status = status;
end

results_file = [results_dir matcher '-' estimator '.mat'];
save(results_file, 'Results');

disp('Finished.');
end

author: JiawangBian

repo: FM-Bench


text:
f = greasy;
a = 128;
M = 1024;
M2 = floor(M/2) + 1;
gl = M;
L = dgtlength(numel(f),a,M);
g = firwin('blackman',gl);
gd = long2fir(gabdual(g,a,M),gl);
N = L/a;
maxit = 100;

corig = dgtreal(f,{'blackman',gl},a,M);
s = abs(corig) + 1i*zeros(size(corig));

cinPtr = libpointer('doublePtr',complex2interleaved(s));
cout = zeros(2*M2,N);
coutPtr = libpointer('doublePtr',cout);

calllib('libphaseret','phaseret_gla_d',cinPtr,libpointer(),g,L,gl,1,a,M,maxit,coutPtr);

cout2 = interleaved2complex(coutPtr.Value);

frec = idgtreal(cout2,{'dual',{'blackman',gl}},a,M);

s2 = dgtreal(frec,{'blackman',gl},a,M);
magnitudeerrdb(s,s2)

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function pstruct = tapas_hgf_ar1_binary_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

l = (length(pvec)+1)/6;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mu_0      = pvec(1:l);
pstruct.sa_0      = pvec(l+1:2*l);
pstruct.phi       = pvec(2*l+1:3*l);
pstruct.m         = pvec(3*l+1:4*l);
pstruct.ka        = pvec(4*l+1:5*l-1);
pstruct.om        = pvec(5*l:6*l-1);

return;

author: translationalneuromodeling

repo: tapas


text:
function [c, s] = onnls(y, g, lam, shift, win, tol, maxIter, mask, smin)
%% Infer the most likely discretized spike train underlying an AR(2) fluorescence trace
% Solves the sparse non-negative deconvolution problem
%  min 1/2|Ks-y|^2 + lam * |s|_1 subject to s_t = c_t-g c_{t-1} >= 0

%% inputs:
%   y:  T*1 vector, vth dimensional array containing the fluorescence intensities
        %withone entry per time-bin.
%   g:  vector, shape (p,)
%       if p in (1,2): AR coefficients for AR(p) process
%       else: kernel that models the fluorescence implulse response
%   lam:  scalar, sparsity penalty parameter lambda.
%   shift: integer scalar, number of frames by which to shift window from on run of
%       NNLS to the next, default-100
%   win: integer acalar, window size
%   tol: scalar, tolerance parameters
%   maxIter: scalar, maximum number of iterations before termination
%   mask: T * 1 boolean vector, restrict potential spike times
%   smin: scalar, minimum spike size
%% outputs
%   c: T*1 vector, the inferred denoised fluorescence signal at each time-bin.
%   s: T*1 vector, discetized deconvolved neural activity (spikes)

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% input arguments
T = length(y);
y = reshape(y, [], 1);

if ~exist('lam', 'var') || isempty(lam)
    lam = 0;
end
if ~exist('shift', 'var') || isempty(shift)
    shift = 100;
end
if ~exist('win', 'var') || isempty(win)
    win = 200;
end

if ~exist('tol', 'var') || isempty(tol)
    tol = 1e-9;
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = [];
end
if ~exist('mask', 'var') || isempty(mask)
    mask = true(T,1);
end
if ~exist('smin', 'var') || isempty(smin)
    smin = 0;
end
%% get the response kernel
w = win;
K = zeros(w);
[u, t] = meshgrid(1:w, 1:w);
ind = 1+t-u;
if length(g)==1
    h = exp(log(g)*(0:(w-1)));
elseif length(g)==2
    temp = roots([1, -g(1), -g(2)]);
    d = max(temp);
    r = min(temp);
    h = (exp(log(d)*(1:w)) - exp(log(r)*(1:w))) / (d-r); % convolution kernel
else
    h = g;
end
K(ind>0) = h(ind(ind>0));   % convolution matrix
KK = K'*K;

%% initialization
a = sum(inv(K));
yp = y - lam * a(1);
yp((end-w+1):end) = yp((end-w+1):end) - lam * a';
s = zeros(T, 1);
c = zeros(T, 1);

%% run online deconvolution
t = 1;
yp0 = yp;
while t <= T-w+1
    ind = t:(t+w-1);
    s(ind) = nnls(KK, K'*yp(ind), s(ind), tol, maxIter, mask(ind));
    yp(ind) = yp(ind) - K(:, 1:shift)*s(t:(t+shift-1));
    c(ind) = c(ind) + K(:, 1:shift)*s(t:(t+shift-1));
    t = t + shift;
end
s(t:T) = nnls(KK((t+w-T):w, (t+w-T):w), K(1:(T-t+1), 1:(T-t+1))'*yp(t:T), ...
    s(t:T), tol, maxIter, mask(t:T));
c(t:T) = c(t:T) + K((t+w-T):w, (t+w-T):w) * s(t:T);

%% running thresholded version
if smin>0
    yp = yp0;
    c = zeros(size(c));
    mask = (s>smin/100);
    t = 1;
    while t <= T-w+1
        ind = t:(t+w-1);
        s(ind) = nnls(KK, K'*yp(ind), s(ind), tol, maxIter, mask(ind), smin);
        yp(ind) = yp(ind) - K(:, 1:shift)*s(t:(t+shift-1));
        c(ind) = c(ind) + K(:, 1:shift)*s(t:(t+shift-1));
        t = t + shift;
    end
    s(t:T) = nnls(KK((t+w-T):w, (t+w-T):w), K(1:(T-t+1), 1:(T-t+1))'*yp(t:T), ...
        s(t:T), tol, maxIter, mask(t:T), smin);
    c(t:T) = c(t:T) + K((t+w-T):w, (t+w-T):w) * s(t:T);
end

function s = nnls(KK, Ky, s, tol, maxIter, mask, smin)
%% fast algorithm for solving nonnegativity constrained least squared
% problem minize norm(y-K*s, 2), s.t. s>=0.

%% inputs:
%   KK: p x p matrix, K'*K
%   Ky: n x 1 vector, K'*y
%   s: p x 1 vector, warm started s
%   tol: scalar, smallest nonzero values
%   maxIter: scalar, maximum nonzero values
%   mask: p x 1 vector, mask to restrict potential spike times considered
%   smin: scala, minimize size of the spike

%% outputs:
%   s: p x 1 vector, solution

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging
% Bro R & Jong S, Journal of Chemometrics 1997, A FAST NON-NEGATIVITY-CONSTRAINED LEAST SQUARES ALGORITHM

%% input arguments
if ~exist('mask', 'var')||isempty(mask)
    mask = true(size(KK,1), 1);
elseif any(mask)
    KK = KK(mask, mask);
    Ky = Ky(mask);
else
    s = double(mask);
    return;
end

p = size(KK,2);      % dimension of s
if ~exist('smin', 'var') || isempty(smin)
    smin = 0;
end
vth = ones(p,1)*smin;
if ~exist('s', 'var') || isempty(s)
    s = zeros(p, 1);
    l = Ky;
    Pset = false(p,1);
else
    s = s(mask);
    s = s - smin;
    Pset = (s>0);
    s(~Pset) = 0;
    l = Ky - KK*s - KK * Pset*smin;
end
if ~exist('tol', 'var') || isempty(tol)
    tol = 1e-9;
end
if ~exist('maxIter', 'var') || isempty(maxIter)
    maxIter = p;
end

% outer loop: loop for passive set
[~, ind] = max(l);
for miter =1:maxIter
    % remove element from the active set
    Pset(ind) = true;

    % solve unconstrained least squares over the passive set
    try
        mu = KK(Pset, Pset) \ (Ky(Pset)-KK(Pset, Pset)*vth(Pset));
    catch % sigular issue
        mu = (KK(Pset, Pset) + tol*eye(sum(Pset))) \ (Ky(Pset)...
            -KK(Pset, Pset)*vth(Pset));
    end

    % inner loop: correct nonnegativity violations
    while any(mu<=0)
        temp = s(Pset) ./ (s(Pset)-mu);
        a = min(temp(mu<0));
        s(Pset) = s(Pset) + a*(mu-s(Pset));
        Pset(s<=tol) = false;

        % solve unconstrained least squares over the passive set
        try
            mu = KK(Pset, Pset) \ (Ky(Pset)-KK(Pset, Pset)*vth(Pset));
        catch % sigular issue
            mu = (KK(Pset, Pset) + tol*eye(sum(Pset))) \ (Ky(Pset)...
                -KK(Pset, Pset)*vth(Pset));
        end
    end

    s(Pset) = mu;
    l = Ky - KK*s - KK*Pset*smin;
        % at least one iteration
    [lmax, ind] = max(l);
    if lmax < tol         % no more passive set
        break;
    end
end
s(Pset) = mu + smin;
temp = double(mask);
temp(mask) = s;
s = temp;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [warp_pts center]= getLKcorner(warp_p, sz)

template_nx = sz(2);
template_ny = sz(1);

tmplt_pts= [1, 1;
        1, template_ny;
        template_nx, template_ny;
        template_nx, 1]';

if size(warp_p,1)==2
    M = [warp_p; 0 0 1];
    M(1,1) = M(1,1) + 1;
    M(2,2) = M(2,2) + 1;
else
    M = warp_p;
%     M(1,3) = M(1,3) + 1;
%     M(2,3) = M(2,3) + 1;
end

warp_pts = M * [tmplt_pts; ones(1, size(tmplt_pts,2))];

c = [(1+template_nx)/2; (1+template_ny)/2; 1];

center = M * c;

warp_pts = warp_pts(1:2,:);

center = center(1:2)';

author: flyers

repo: drone-tracking


text:
function [shortestPath, totalCost] = dijkstra(netCostMatrix, s, d)
%==============================================================
% shortestPath: the list of nodes in the shortestPath from source to destination;
% totalCost: the total cost of the  shortestPath;
% farthestNode: the farthest node to reach for each node after performing the routing;
% n: the number of nodes in the network;
% s: source node index;
% d: destination node index;
%==============================================================
%  Code by:
% ++by Xiaodong Wang
% ++23 Jul 2004 (Updated 29 Jul 2004)
% ++http://www.mathworks.com/matlabcentral/fileexchange/5550-dijkstra-shortest-path-routing
% Modifications (simplifications) by Meral Shirazipour 9 Dec 2009
%==============================================================
n = size(netCostMatrix,1);
for i = 1:n
    % initialize the farthest node to be itself;
    farthestPrevHop(i) = i; % used to compute the RTS/CTS range;
    farthestNextHop(i) = i;
end

% all the nodes are un-visited;
visited(1:n) = false;

distance(1:n) = inf;    % it stores the shortest distance between each node and the source node;
parent(1:n) = 0;

distance(s) = 0;
for i = 1:(n-1),
    temp = [];
    for h = 1:n,
         if ~visited(h)  % in the tree;
             temp=[temp distance(h)];
         else
             temp=[temp inf];
         end
     end;
     [t, u] = min(temp);      % it starts from node with the shortest distance to the source;
     visited(u) = true;         % mark it as visited;
     for v = 1:n,                % for each neighbors of node u;
         if ( ( netCostMatrix(u, v) + distance(u)) < distance(v) )
             distance(v) = distance(u) + netCostMatrix(u, v);   % update the shortest distance when a shorter shortestPath is found;
             parent(v) = u;     % update its parent;
         end;
     end;
end;

shortestPath = [];
if parent(d) ~= 0   % if there is a shortestPath!
    t = d;
    shortestPath = [d];
    while t ~= s
        p = parent(t);
        shortestPath = [p shortestPath];

        if netCostMatrix(t, farthestPrevHop(t)) < netCostMatrix(t, p)
            farthestPrevHop(t) = p;
        end;
        if netCostMatrix(p, farthestNextHop(p)) < netCostMatrix(p, t)
            farthestNextHop(p) = t;
        end;

        t = p;
    end;
end;

totalCost = distance(d);

%return;

author: Sable

repo: mcbench-benchmarks


text:
function [ alpha ] = LiMapS(s, D, DINV, gamma, maxIter)
%
% [ alpha ] = LIMAPS(s, D, DINV, gamma)
%
%   Returns the sparsest vector alpha which satisfies underdetermined
%   system of equations  D*alpha=s, using  Lipschitzian Mapping for Sparsity
%   Recovery algorithm.
%   The dictionary matrix  D should have more columns than rows. The
%   number of rows of dictionary D  must be equal to the length of the
%   column observation vector(or matrix) s.
%
%     The observations signal(s) s and the dictionary D necessarily, must
%   be provided by the user. The other parameters have defult values, or
%   may be provided by the user.
%
%   This algorithm work fine with matrix too, returning the sparsest matrix
%   alpha that satisfy the underdetermined system of equations D*alpha=s.
%   Then the dictionary D is of size [n x m], s must be of size [ s x T ] and
%   as consiquence the coefficients matirx alpha is of size [ m x T]
%
%   s:
%       The observation vector (matrix)
%
%   D:
%       The dictionary
%
%   DINV (optional):
%       The pseudo-inverse of the dictionary matrix D
%
%   gamma (optional):
%       The increase factor (default 1.01)
%
%   maxIter (optional):
%       Maximum number of iterations (default 1000)
%
% Authors: Alessandro Adamo and Giuliano Grossi
% Version: 1.0
% Last modified: 1 Aug. 2011.
%
% WEB PAGE:
% ------------------
%    http://dalab.dsi.unimi.it/limaps
%
% HISTORY:
%--------------
% Version 1.0: 1 Aug. 2011
%        first official version.
%

if(nargin<2||nargin>5)
    error('Wrong parameters number');
end
if(nargin==2)
    DINV=pinv(D);
    gamma=1.01;
    maxIter=1000;
end
if(nargin==3)
    gamma=1.01;
    maxIter=1000;
end
if(nargin==4)
    maxIter=1000;
end

%% --  INITIALIZATION  --
alpha = DINV*s;

lambda = 1/max(abs(alpha(:)));
epsilon=1e-5;                    % stopping criteria
alpha_min=1e-2;

%% --  CORE  --
for i=1:maxIter

    alphaold=alpha;
    % apply sparsity constraction mapping: increase sparsity
    beta = alpha.*(1-exp(-lambda.*abs(alpha)));

    beta(abs(beta)<alpha_min)=0;

    % apply the orthogonal projection
    alpha = beta-DINV*(D*beta-s);

    % update the lambda coefficient
    lambda = lambda*gamma;

    % check the stopping criteria
    if (norm(alpha(:)-alphaold(:))<epsilon)
        break;
    end

    if (mod(i,100))
        if(sum(abs(1-exp(-lambda*abs(alpha)))>1/lambda)>size(D,1))
           break;
        end
    end

end

%% -- REFINEMENT --

% threshold the coefficients
alpha(abs(alpha)<alpha_min) = 0;

% final refinements of the solution
for i=1:size(alpha,2);
    alpha_ind = alpha(:,i)~=0;
    D1 = D(:,alpha_ind);
    alpha(alpha_ind,i) = alpha(alpha_ind,i) - pinv(D1)*(D1*alpha(alpha_ind,i)-s(:,i));
end

end


text:
function events = in_events_micromed(sFile, EventFile)
% IN_EVENTS_MICROMED: Open a Micromed .EVT file.
%
% OUTPUT:
%    - events(i): array of structures with following fields (one structure per event type)

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

% Initialize returned value
events = [];
% Read XML file
sXml = in_xml(EventFile);
if ~isfield(sXml, 'EventFile') || ~isfield(sXml.EventFile, 'Events') || ~isfield(sXml.EventFile.Events, 'Event') || isempty(sXml.EventFile.Events.Event)
    return;
end

% Get recordings start date (timestamp)
acq = sFile.header.acquisition;
start_time = posixtime(datetime(acq.year, acq.month, acq.day, acq.hour, acq.min, acq.sec, 0));

% Parse structure
sEvents = sXml.EventFile.Events.Event;
Markers = cell(length(sEvents), 3);
for iMrk = 1:length(sEvents)
    % Get time from the beginning of the file
    Markers{iMrk,1} = [posixtime(datetime(sEvents(iMrk).Begin.text, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSS')) - start_time; ...
                       posixtime(datetime(sEvents(iMrk).End.text, 'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSS')) - start_time];
    % TEXT: Interpret as event ename if the string is short
    if (length(sEvents(iMrk).Text.text) < 50)
        Markers{iMrk,2} = sEvents(iMrk).Text.text;
        if isfield(sEvents(iMrk), 'ExtraText') && isfield(sEvents(iMrk).ExtraText, 'text') && ~isempty(sEvents(iMrk).ExtraText.text)
            Markers{iMrk,3} = sEvents(iMrk).ExtraText.text;
        end
    % Otherwise, consider it as a note for an event EVT
    else
        Markers{iMrk,2} = 'EVT';
        Markers{iMrk,3} = sEvents(iMrk).Text.text;
    end
end

% List of events
uniqueEvt = unique(Markers(:,2)');
% Initialize returned structure
events = repmat(db_template('event'), 0);
% Create events list
for iEvt = 1:length(uniqueEvt)
    % Find all the occurrences of event #iEvt
    iOcc = find(strcmpi(Markers(:,2)', uniqueEvt{iEvt}));
    % Compute samples from the file start
    samples = round([Markers{iOcc,1}] .* sFile.prop.sfreq);
    notes = Markers(iOcc,3)';
    % Detect if simple events (duration <= 2 samples)
    if ~any(samples(2,:) - samples(1,:) > 2)
        samples(2,:) = [];
    end
    % Adjust if there are segments in the file
    if isfield(sFile.header, 'segment') && ~isempty(sFile.header.segment)
        for iOcc = 1:size(samples,2)
            iSeg = find(samples(1,iOcc) >= [sFile.header.segment.time], 1, 'last');
            if ~isempty(iSeg)
                samples(:,iOcc) = ...
                    samples(:,iOcc) - sFile.header.segment(iSeg).time ...                    % Offset from beginning of segment (in samples)
                    + sFile.header.segment(iSeg).sample - sFile.header.segment(1).sample ... % Offset from beginning of file
                    + round(sFile.prop.times(1) .* sFile.prop.sfreq);                        % Start of the file
            end
        end
    end
    % Detect events outside of the file definition
    iOut = find((samples(end,:) < sFile.prop.times(1) .* sFile.prop.sfreq) | (samples(1,:) > sFile.prop.times(2) .* sFile.prop.sfreq));
    if ~isempty(iOut)
        disp([sprintf('BST> Event "%s": %d occurrences outside of the recordings [%1.3fs,%1.3fs]:   ', uniqueEvt{iEvt}, length(iOut), sFile.prop.times(1), sFile.prop.times(2)), ...
              sprintf('%1.3fs ', samples(1,iOut) ./ sFile.prop.sfreq)]);
        samples(:,iOut) = [];
        notes(iOut) = [];
    end
    if isempty(samples)
        continue;
    end
    % Add event structure
    iNew = length(events) + 1;
    events(iNew).label      = uniqueEvt{iEvt};
    events(iNew).epochs     = ones(1, size(samples,2));
    events(iNew).times      = samples ./ sFile.prop.sfreq;
    events(iNew).reactTimes = [];
    events(iNew).select     = 1;
    events(iNew).channels   = [];
    events(iNew).notes      = notes;
end

author: brainstorm-tools


text:
function [x,w] = GaussHermite2(N)
%
% [X W] = HERMQUAD(N)
%
% Find the Gauss-Hermite abscissae and weights.
%
% Arguments:
%  N - The number of abscissae and weights to return.
%
% Return Values:
%  X - A column vector containing the abscissae.
%  W - A column vector containing the corresponding weights.
%
% Gauss-Hermite quadrature approximates definite integrals of the form
%
%     \int^{-\infty}_{\infty} dx W(x) f(x)
%
% where
%
%     W(x) = \exp( - x^2 )
%
% with the sum
%
%     \sum_{n=1}^{N} w_{n} f(x_{n}).
%
% This function returns the set of abscissae and weights
%
%     {x_{n}, w_{n}}^{N}_{n=1}
%
% for performing this calculation given N, the number of abscissae.
% These abscissae correspond to the zeros of the Nth Hermite
% polynomial.  It can be shown that such integration is exact when f(x)
% is a polynomial of maximum order 2N-1.
%
% The procedure in this calculation is taken more or less directly from
%
% @BOOK{ press-etal-1992a,
%	 AUTHOR   = { Press, William  H.   and
%                     Flannery, Brian  P.  and
%                     Teukolsky, Saul  A.  and
%                     Vetterling, William  T. },
%       ISBN      = {0521431085},
%	MONTH     = {October},
%	PUBLISHER = {{Cambridge University Press}},
%	TITLE     = {Numerical Recipes in C : The Art of Scientific Computing},
%	YEAR      = {1992}
%      }
%

% precision
EPS = 3.0e-14;

% 1/\pi^{1/4}
PIM4 = 0.7511255444649425;

% maximum number of loops
MAXIT = 10;

% allocate the return values
x = zeros([N 1]);
w = zeros([N 1]);

for i=1:(N+1)/2

  % good guesses at initial values for specific roots
  if i == 1
    z = sqrt(2.0*N+1.0) - 1.85575*((2.0*N+1)^(-0.16667));
  elseif i == 2
    z = z - (1.14 * N^0.426 / z);
  elseif i == 3
    z = 1.86 * z - 0.86 * x(1);
  elseif i == 4
    z = 1.91 * z - 0.91 * x(2);
  else
    z = 2.0*z - x(i-2);
  end

  for iter=1:MAXIT+1
    p1 = PIM4;
    p2 = 0.0;

    for j=1:N
      p3 = p2;
      p2 = p1;
      p1 = z * sqrt(2.0/j) * p2 - sqrt((j-1.0)/j) * p3;
    end

    % the derivative
    pp = sqrt(2.0*N) * p2;

    % newton step
    z1 = z;
    z  = z1 - p1/pp;

    if abs(z-z1) <= EPS
      break;
    end
  end

  if iter == MAXIT+1
    fprintf('Too many iterations in hermquad.\n');
  end

  x(i)     = z;
  x(N+1-i) = -z;
  w(i)     = 2.0/(pp*pp);
  w(N+1-i) = w(i);

end

repo: Aero-matlab


text:
function h = p32_h ( n, x )

%*****************************************************************************80
%
%% P32_H evaluates the Hessian for problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  m = 7;
  a = [
    4.0, 4.0, 4.0, 4.0; ...
    1.0, 1.0, 1.0, 1.0; ...
    8.0, 8.0, 8.0, 8.0; ...
    6.0, 6.0, 6.0, 6.0; ...
    3.0, 7.0, 3.0, 7.0; ...
    2.0, 9.0, 2.0, 9.0; ...
    5.0, 5.0, 3.0, 3.0 ]';

  c = [ 0.1, 0.2, 0.2, 0.4, 0.6, 0.6, 0.3 ]';

  for ii = 1 : n
    for jj = 1 : n
      for j = 1 : m
        d = c(j) + sum ( ( x(1:n) - a(1:n,j) ).^2 );
        h(ii,jj) = h(ii,jj) - 8.0 * ( x(ii) - a(ii,j) ) * ( x(jj) - a(jj,j) ) / d^3;
        if ( ii == jj )
          h(ii,jj) = h(ii,jj) + 2.0 / d^2;
        end
      end
    end
  end

  return
end


text:
%Copyright 2013 The MathWorks, Inc.

% These are the files we want to build.
files = {
    'applyScaleFactors'

    }';

% Check that we are running Windows.
assert ( ispc() )

% Check if we are running 32 or 64 bit MATLAB.
bittedness = 32;
if ( strcmp( computer, 'PCWIN64' ) )
	bittedness = 64;
end

% Hard link to Visual Studion (Express) compiler binaries.
VCpath = '"C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\bin"';

% Link to the MATLAB include files to get the type declarations right.
incpath = [ '"' matlabroot() '/extern/include"' ];

for ii=1:numel(files)
	cmd = sprintf( 'nvcc --machine %d -ptx -Xptxas=-v -arch sm_20 -ccbin %s -I %s %s.cu --output-file %s.%s', bittedness, VCpath, ...
          incpath , ...
          files{ii}, files{ii}, parallel.gpu.ptxext );
    fprintf( 'Compiling %s... ', files{ii} );
    system( cmd );
    fprintf( '%s done.\n', repmat(' ',[1,max(0,20-length(files{ii}))]) );
end


text:
function opti_package
% OPTI_PACKAGE  Generates GitHub releases .zip files

% Get Solver Names (ones we want to upload)
s = lower(optiSolver('all'));
removeIdx = strcmp(s,'baron') | strcmp(s,'cplex') | strcmp(s,'gmatlab') | strcmp(s, 'matlab') | strcmp(s,'mosek') | strcmp(s,'scip') | strcmp(s,'sedumi');
solvers = s(~removeIdx);

% Utility Names
utilities = {'asl', 'coinR', 'coinW', 'mklJac', 'rmathlib'};

% Build zip file list
files = [solvers utilities];
for i = 1:length(files)
    files{i} = which([files{i} '.' mexext]);
end

% Zip the files
fprintf('Compressing MEX Files...');
name = ['optiMEXFiles_' mexext '_' sprintf('%.0f_%.0f',optiver,mod(optiver*100,100))];
zip(name,files);
rehash;
rehash;
fprintf('Done!\n');

% Update the contents file
UpdateContentsFile('OPTI',optiver,'Utilities/opti');

% Copy latest SCIP to download folder
copyfile(which(['scip.' mexext]), '../../../Inverse Problem/OPTI Academic Solvers');

function UpdateContentsFile(name,tbxver,contentsFile)
%Update contents file description line

if(~isempty(strfind(contentsFile,'Contents.m')))
    p = contentsFile;
elseif(contentsFile(end) == filesep)
    p = [contentsFile 'Contents.m'];
else
    p = [contentsFile filesep 'Contents.m'];
end
s = fileread(p);
%Find initial date
idx = strfind(s,'(C)');
idxe = strfind(s(idx:end),'-');
stDate = strtrim(s(idx+3:idx+idxe-2));
mver = ver('matlab');
%Regenerate Contents File
fid = fopen(p,'w');
if(fid)
    try
        fprintf(fid,'%% %s Toolbox\n',upper(name));
        fprintf(fid,'%% Version %.2f %s %s\n',tbxver,mver.Release,datestr(now,1));
        fprintf(fid,'%% Copyright (C) %s-%s Jonathan Currie (Inverse Problem Ltd)\n',stDate,datestr(now,10));
        fprintf(fid,'%% License: https://inverseproblem.co.nz/OPTI/index.php/DL/License\n');
    catch ME
        fclose(fid);
        rethrow(ME);
    end
end
fclose(fid);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [P_st, s] = flicker_sim(u, fs, f_line)
% flicker_sim - Flickermeter Simulator according IEC 61000-4-15
%
% [P_st, s, cpf] = flicker_sim(u, fs, f_line)
%
% This function implements a flickermeter according [1] and [2].
% Requires MATLAB with Signal Procesing Toolbox installed or Octave.
% For more information refer to [3].
%
% Inputs:
%   u:      vector of voltage samples
%   fs:     sampling frequency of u in Hz (should be >= 2000)
%   f_line: line frequency in Hz (must be 50 or 60 Hz)
%
% Outputs:
%   P_st:   short term flicker
%   s:      instantaneous flicker severity
%===============================================================================
% References:
% [1] IEC 61000-4-15, Electromagnetic compatibility (EMC), Testing and
%     measurement techniques, Flickermeter, Edition 1.1, 2003-02
% [2] Wilhelm Mombauer: "Messung von Spannungsschwankungen und Flickern mit
%     dem IEC-Flickermeter", ISBN 3-8007-2525-8, VDE-Verlag
% [3] http://www.solcept.ch/en/embedded-tools/flickersim.html
%===============================================================================

%% Configuration

SHOW_TIME_SIGNALS           = 0; % enable to plot internal signals of flickermeter
SHOW_CUMULATIVE_PROBABILITY = 0; % enable to plot the statistical evaluation of the instantaneous flicker severity
SHOW_FILTER_RESPONSES       = 0; % enable to plot the filter responses of all internal filter stages (for model verification)

IS_OCTAVE = exist('OCTAVE_VERSION') ~= 0;

%% Check inputs

if (nargin ~= 3)
  error('Invalid number of arguments');
end

if (~isvector(u))
  error('First input argument must be a vector');
end
% convert to row vector if needed
u = reshape(u, 1, length(u));

if ((f_line ~= 50) && (f_line ~= 60))
  error('Line frequency must be 50 or 60 Hz');
end

if (fs < 2000)
  warning('Sampling frequency should be >= 2000 Hz');
end

%% Block 1: Input voltage adaptor

% remove DC component
u = u - mean(u);

% normalize input (scale with peak-amplitude value)
u_rms = sqrt(mean(u.^2));
u = u / (u_rms * sqrt(2));

%% Block 2: Quadratic demodulator

u_0 = u .^ 2;

%% Block 3: Bandpass and weighting filter

% bandpass filter

HIGHPASS_ORDER  = 1;
HIGHPASS_CUTOFF = 0.05;

LOWPASS_ORDER = 6;
if (f_line == 50)
  LOWPASS_CUTOFF = 35;
end
if (f_line == 60)
  LOWPASS_CUTOFF = 42;
end

% subtract DC component to limit filter transients at start of simulation
u_0_ac = u_0 - mean(u_0);

[b_hp, a_hp] = butter(HIGHPASS_ORDER, HIGHPASS_CUTOFF / (fs / 2), 'high');
u_hp = filter(b_hp, a_hp, u_0_ac);

% smooth start of signal to avoid filter transient at start of simulation
smooth_limit = min(round(fs / 10), length(u_hp));
u_hp(1 : smooth_limit) = u_hp(1 : smooth_limit) .* linspace(0, 1, smooth_limit);

[b_bw, a_bw] = butter(LOWPASS_ORDER, LOWPASS_CUTOFF / (fs / 2), 'low');
u_bw = filter(b_bw, a_bw, u_hp);

% weighting filter

if (f_line == 50)
  K = 1.74802;
  LAMBDA = 2 * pi * 4.05981;
  OMEGA1 = 2 * pi * 9.15494;
  OMEGA2 = 2 * pi * 2.27979;
  OMEGA3 = 2 * pi * 1.22535;
  OMEGA4 = 2 * pi * 21.9;
end
if (f_line == 60)
  K = 1.6357;
  LAMBDA = 2 * pi * 4.167375;
  OMEGA1 = 2 * pi * 9.077169;
  OMEGA2 = 2 * pi * 2.939902;
  OMEGA3 = 2 * pi * 1.394468;
  OMEGA4 = 2 * pi * 17.31512;
end

num1 = [K * OMEGA1, 0];
den1 = [1, 2 * LAMBDA, OMEGA1.^2];
num2 = [1 / OMEGA2, 1];
den2 = [1 / (OMEGA3 * OMEGA4), 1 / OMEGA3 + 1 / OMEGA4, 1];
if (IS_OCTAVE)
  [b_w, a_w] = bilinear(conv(num1, num2), conv(den1, den2), 1 / fs);
else
  [b_w, a_w] = bilinear(conv(num1, num2), conv(den1, den2), fs);
end

u_w = filter(b_w, a_w, u_bw);

%% Block 4: Squaring and smoothing

LOWPASS_2_ORDER  = 1;
LOWPASS_2_CUTOFF = 1 / (2 * pi * 300e-3);  % time constant 300 msec
SCALING_FACTOR   = 1238400;  % scaling of output to perceptibility scale  (according [2])

u_q = u_w .^ 2;

[b_lp, a_lp] = butter(LOWPASS_2_ORDER, LOWPASS_2_CUTOFF / (fs / 2), 'low');
s = SCALING_FACTOR * filter(b_lp, a_lp, u_q);

%% Block 5: Statistical evaluation

NUMOF_CLASSES = 10000;

[bin_cnt, cpf.magnitude] = hist(s, NUMOF_CLASSES);
cpf.cum_probability = 100 * (1 - cumsum(bin_cnt) / sum(bin_cnt));

p_50s = mean([get_percentile(cpf, 30), get_percentile(cpf, 50), get_percentile(cpf, 80)]);
p_10s = mean([get_percentile(cpf, 6), get_percentile(cpf, 8), ...
  get_percentile(cpf, 10), get_percentile(cpf, 13), get_percentile(cpf, 17)]);
p_3s = mean([get_percentile(cpf, 2.2), get_percentile(cpf, 3), get_percentile(cpf, 4)]);
p_1s = mean([get_percentile(cpf, 0.7), get_percentile(cpf, 1), get_percentile(cpf, 1.5)]);
p_0_1 = get_percentile(cpf, 0.1);

P_st = sqrt(0.0314 * p_0_1 + 0.0525 * p_1s + 0.0657 * p_3s + ...
  0.28 * p_10s + 0.08 * p_50s);

%% Optional graphical output

% time signals
if (SHOW_TIME_SIGNALS)
  t = 0 : 1 / fs : (length(u) - 1) / fs;
  filter_len = round(10 / 1000 * fs);
  u_0_m = filter(ones(1, filter_len) / filter_len * 2, 1, u_0);

  figure
  clf
  subplot(2, 2, 1)
  hold on
  plot(t, u, 'b');
  plot(t, u_0, 'm');
  plot(t, u_hp, 'r');
  plot(t, u_0_m, 'c');
  hold off
  legend('u', 'u_0', 'u_h_p');
  grid on
  subplot(2, 2, 2)
  hold on
  plot(t, u_bw, 'b');
  plot(t, u_w, 'm');
  legend('u_b_w', 'u_w');
  hold off
  grid on
  subplot(2, 2, 3)
  plot(t, u_q, 'b');
  legend('u_q');
  grid on
  subplot(2, 2, 4)
  plot(t, s, 'b');
  legend('s');
  grid on
end

% cumulative probability function
if (SHOW_CUMULATIVE_PROBABILITY)
  figure
  clf
  plot(cpf.magnitude, cpf.cum_probability);
  grid
end

% frequency responses of filters
if (SHOW_FILTER_RESPONSES)
  [h_hp, f] = freqz(b_hp, a_hp, 4096, fs);
  [h_bw, f] = freqz(b_bw, a_bw, 4096, fs);
  [h_w, f] = freqz(b_w, a_w, 4096, fs);
  [h_lp, f] = freqz(b_lp, a_lp, 4096, fs);

  figure
  clf
  subplot(2, 1, 1)
  hold on
  plot(f, abs(h_hp), 'b')
  plot(f, abs(h_bw), 'r')
  plot(f, abs(h_w), 'g')
  plot(f, abs(h_lp), 'm')
  hold off
  grid
  axis([0, 35, 0, 1]);

  subplot(2, 1, 2)
  hold on
  plot(f, 180 / pi * unwrap(angle(h_hp)), 'b')
  plot(f, 180 / pi * unwrap(angle(h_bw)), 'r')
  plot(f, 180 / pi * unwrap(angle(h_w)), 'g')
  plot(f, 180 / pi * unwrap(angle(h_lp)), 'm')
  hold off
  grid
  axis([0, 35, -200, 300]);
end

end  % end of function flicker_sim

%% Subfunction: get_percentile
function val = get_percentile(cpf, limit)
  [dummy, idx] = min(abs(cpf.cum_probability - limit));
  val = cpf.magnitude(idx);
end


text:
function title = p10_title ( )

%*****************************************************************************80
%
%% P10_TITLE returns the title of problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Y = 2 + 5*X + 10*N(0,1).';

  return
end


text:
function view_rccross_a2(lab1,valueMap)
    % view_rccross_a2 plots ratechanges and p values calculated
    % derived with rc_cross_a2.m or other similar values as a color map.
    % needs valueMap, gx, gy
    %
    % define size of the plot etc.

    report_this_filefun();

    myFigName='RC-Cross-section';
    myFigFinder=@() findobj('Type','Figure','-and','Name',myFigName);

    % This is the info window text
    %
    ttlStr='The b and p -Value Map Window                 ';
    hlpStr1zmap= ...
        ['                                                '
        ' This window displays b-values and p-values     '
        ' using a color code.                            '
        ' Some of the menu-bar options are               '
        ' described below:                               '
        '                                                '
        ' Threshold: You can set the maximum size that   '
        '   a volume is allowed to have in order to be   '
        '   displayed in the map. Therefore, areas with  '
        '   a low seismicity rate are not displayed.     '
        '   edit the size (in km) and click the mouse    '
        '   outside the edit window.                     '
        'FixAx: You can chose the minimum and maximum    '
        '        values of the color-legend used.        '
        'Polygon: You can select earthquakes in a        '
        ' polygon either by entering the coordinates or  '
        ' defining the corners with the mouse            '];
    hlpStr2zmap= ...
        ['                                                '
        'Circle: Select earthquakes in a circular volume:'
        '      Ni, the number of selected earthquakes can'
        '      be edited in the upper right corner of the'
        '      window.                                   '
        ' Refresh Window: Redraws the figure, erases     '
        '       selected events.                         '

        ' zoom: Selecting Axis -> zoom on allows you to  '
        '       zoom into a region. Click and drag with  '
        '       the left mouse button. type <help zoom>  '
        '       for details.                             '
        ' Aspect: select one of the aspect ratio options '
        ' Text: You can select text items by clicking.The'
        '       selected text can be rotated, moved, you '
        '       can change the font size etc.            '
        '       Double click on text allows editing it.  '
        '                                                '
        '                                                '];

    % Set up the Seismicity Map window Enviroment
    %

    % Find out if figure already exists
    %
    hRccross=myFigFinder();

    if isempty(hRccross)
        hRccross = figure_w_normalized_uicontrolunits( ...
            'Name',myFigName,...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        %lab1 = 'p-value:';
        create_my_menu();

        %valueMap = pvalg;
        ZG.tresh_km = nan; re4 = valueMap;

        colormap(jet)
        ZG.tresh_km = nan; minpe = nan; Mmin = nan; minsd = nan;
    end   % This is the end of the figure setup.

    % Now lets plot the color-map!
    %
    figure(hRccross);
    delete(findobj(hOmoricross,'Type','axes'))
    %delete(sizmap)
    ax=gca;
    reset(ax) % automatically sets NextPlot to replace
    watchon;
    set(ax,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','SortMethod','childorder')

    rect = [0.18,  0.10, 0.7, 0.75];
    rect1 = rect;

    % find max and min of data for automatic scaling
    ZG.maxc = max(valueMap(:));
    ZG.maxc = fix(ZG.maxc)+1;
    ZG.minc = min(valueMap(:));
    ZG.minc = fix(ZG.minc)-1;

    re4 = valueMap;

    % plot image
    orient landscape

    %Plots re4, which contains the filtered values.
    axes('position',rect)
    set(gca,'NextPlot','add')
    pco1 = pcolor(gx,gy,re4);

    axis([ min(gx) max(gx) min(gy) max(gy)])
    axis image
    set(gca,'NextPlot','add')
    shading(ZG.shading_style);

    % make the scaling for the recurrence time map reasonable
    if lab1(1) =='T'
        re = valueMap(~isnan(valueMap));
        caxis([min(re) 5*min(re)]);
    end

    fix_caxis.ApplyIfFrozen(gca);

    title([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
        'Color','r','FontWeight','bold')

    xlabel('Distance [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Depth [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

    % plot overlay
    %
    set(gca,'NextPlot','add')
    zmap_update_displays();
    globalcatalog=ZG.primeCatalog;
    ploeq = plot(globalcatalog.X, globalcatalog.Y,'k.');
    set(ploeq,'Tag','eq_plot','MarkerSize',ZG.ms6,'Marker',ty,'Color',ZG.someColor,'Visible','on')

    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    h1 = gca;
    hzma = gca;

    % Create a colorbar
    %
    h5 = colorbar('horiz');
    set(h5,'Pos',[0.35 0.05 0.4 0.02],...
        'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s,'TickDir','out')

    rect = [0.00,  0.0, 1 1];
    axes('position',rect)
    axis('off')
    %  Text Object Creation
    txt1 = text(...
        'Units','normalized',...
        'Position',[ 0.33 0.06 0 ],...
        'HorizontalAlignment','right',...
        'FontSize',ZmapGlobal.Data.fontsz.s,....
        'FontWeight','bold',...
        'String',lab1);

    % Make the figure visible
    %
    set(gca,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','TickDir','out')
    figure(hRccross);
    axes(h1)
    watchoff(hRccross)

    %% ui functions
    function create_my_menu()
        add_menu_divider();
        add_symbol_menu('eq_plot');

        options = uimenu('Label',' Analyze ');
        uimenu(options,'Label','Refresh ','MenuSelectedFcn',@callbackfun_001)
        uimenu(options,'Label','Select EQ in Circle - Constant R',...
            'MenuSelectedFcn',@callbackfun_002)
        uimenu(options,'Label','Select EQ with const. number',...
            'MenuSelectedFcn',@callbackfun_003)

        op1 = uimenu('Label',' Maps ');

        %Menu for adjusting several parameters.
        adjmenu =  uimenu(op1,'Label','Adjust Map Display Parameters'),...
            uimenu(adjmenu,'Label','Adjust Mmin cut',...
            'MenuSelectedFcn', @(~,~)cb_viewagain('mag') )
        uimenu(adjmenu,'Label','Adjust Rmax cut',...
            'MenuSelectedFcn', @(~,~)cb_viewagain('rmax') )
        uimenu(adjmenu,'Label','Adjust goodness of fit cut',...
            'MenuSelectedFcn', @(~,~)cb_viewagain('gofi') )
        uimenu(adjmenu,'Label','Adjust p-value st. dev. cut',...
            'MenuSelectedFcn', @(~,~)cb_viewagain('pstdc') )

        uimenu(op1,'Label','Relative rate change (bootstrap)',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('Sigma') )
        uimenu(op1,'Label','Model',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('Model') )
        uimenu(op1,'Label','KS-Test',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('Rejection') )
        uimenu(op1,'Label','KS-Test Statistic',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('KS distance') )
        uimenu(op1,'Label','KS-Test p-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('KS-Test p-value') )
        uimenu(op1,'Label','RMS of fit',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('RMS') )
        uimenu(op1,'Label','Resolution Map (Number of events)',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('Number of events') )
        uimenu(op1,'Label','Resolution Map (Radii)',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('Radius / [km]') )
        uimenu(op1,'Label','p-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('p-value') )
        uimenu(op1,'Label','p-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('p-valstd') )
        uimenu(op1,'Label','c-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('c-value') )
        uimenu(op1,'Label','c-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('c-valuestd') )
        uimenu(op1,'Label','k-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('k-value') )
        uimenu(op1,'Label','k-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('k-valuestd') )
        uimenu(op1,'Label','p2-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('p2-value') )
        uimenu(op1,'Label','p-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('p2-valuestd') )
        uimenu(op1,'Label','c2-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('c2-value') )
        uimenu(op1,'Label','c2-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('c2-valuestd') )
        uimenu(op1,'Label','k2-value',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('k2-value') )
        uimenu(op1,'Label','k2-value standard deviation',...
            'MenuSelectedFcn', @(~,~)cb_viewagainOtherdata('k2-valuestd') )

        add_display_menu(1);
    end

    %% callback functions

    function callbackfun_001(mysrc,myevt)
        view_rccross_a2(lab1,valueMap);
    end

    function callbackfun_002(mysrc,myevt)
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        met = 'ra';
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        plot_circbootfit_a2;
        watchoff(hRccross);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG=ZmapGlobal.Data;
        ZG.hold_state2=true;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=true;
        plot_constnrbootfit_a2;
        watchoff(hRccross);
    end

    function cb_viewagain(asel)
        % set asel, adju2, then view
        ZG=ZmapGlobal.Data;
        adju2;
        view_rccross_a2(lab1,ZG.valueMap);
    end

    function cb_viewagainOtherdata(label)
        % set asel, adju2, then view
        ZG=ZmapGlobal.Data;
        switch label
            case 'Sigma'
                valueMap=mRelchange; %8
            case 'Model'
                valueMap=mMod; %9
            case 'Rejection'
                valueMap=mKstestH; %10
            case 'KS distance'
                valueMap=mKsstat; %11
            case 'KS-Test p-value'
                valueMap=mKsp; %12
            case 'RMS'
                valueMap=mRMS; %13
            case 'Number of events'
                valueMap=mNumevents; %14
            case 'Radius / [km]'
                valueMap=vRadiusRes; %15
            case 'p-value'
                valueMap=mPval;
            case 'p-valstd'
                valueMap=mPvalstd;
            case 'c-value'
                valueMap=mCval;
            case 'c-valstd'
                valueMap=mCvalstd;
            case 'k-value'
                valueMap=mKval;
            case 'k-valstd'
                valueMap=mKvalstd;
            case 'p2-value'
                valueMap=mPval2;
            case 'p2-valstd'
                valueMap=mPvalstd2;
            case 'c2-value'
                valueMap=mCval2;
            case 'c2-valstd'
                valueMap=mCvalstd2;
            case 'k2-value'
                valueMap=mKval2;
            case 'k2-valstd'
                valueMap=mKvalstd2;
            otherwise
                error('unknown choice');
        end
        ZG.valueMap=valueMap;
        view_rccross_a2(label,ZG.valueMap);
    end

end

author: CelsoReyes


text:
function [obX,obY,obSize,obSizeOrig,sagPts,sagSlice,lp,obPts,obSlice] = ...
              mrReloadParams(lp,curInplane,obXM,obYM,sagSize,numSlices,volume,cTheta,aTheta,curSag,reflections,scaleFac)

% PURPOSE: Implements the reloading of variables from AlignParams
% AUTHOR: SPG
% DATE: 11.21.97
% NOTES:

if ~isempty(obXM)
    obX = obXM(curInplane,:);
    obY = obYM(curInplane,:);
    lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);
else
    obX = [0,0];
    obY = [0,0];
end

[obSize,obSizeOrig] = mrFindObSize(obX,obY,sagSize,numSlices,curInplane);

[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [triggerLinesGood, triggerCounts, triggersList, fs] = sqdgettriggers(sqdfile, triggerLinesPotential, triggerHighThresh, triggerLowThresh)

%SQDGETTRIGGERS Detects trigger lines and triggers on those lines from an MEG sqd file.
%
%    [triggerLines, triggerCounts, triggersList, fs] = SQDGETTRIGGERS(sqdfile, triggerLinesPotential, triggerHighThresh, triggerLowThresh);
%
%    OUTPUT variables:
%    triggerLines: which lines were actually used as trigger lines
%    triggerCounts: how many triggers were detected on each trigger line
%    triggersList: a cell array (each cell corresponds to a
%       different trigger line) in which each cell is a list of the trigger
%       events. The trigger events are labeled by sample number not time. To
%       convert into time, divide the sample numbers by the sample
%       frequency.
%    fs: the sample frequency, in Hz.
%
%    INPUT variables:
%    sqdfile: the MEG file holding the trigger (and MEG) data
%    triggerLinesPotential: a list of which lines should be inspected for
%       triggers. The default is [160:191], i.e. all channels not containing
%       magnetic field data. This function is much more efficient, though, if
%       the triggerLinesPotential is the list of channels actually used as
%       trigger lines. Because of the way the sqd file is read, however, the
%       fastest method is to set triggerLinesPotential to be list of all
%       channels from the lowest used to the highest used.
%    triggerHighThresh (rarely needed): the threshold, in Volts, above
%       which a trigger line goes high. This value is nominally slightly
%       larger than 5 V, so the default value is 5 V.
%    triggerHighThresh (rarely needed): the threshold, in Volts, below
%       which a trigger line is considered to be low. This
%       value is nominally 0 V, so the default value is 0.2 V.

%
%    Version 0.9 beta 2
%    9 January 2010
%
%    by Jonathan Z. Simon

% Create triggerLinesPotential if not supplied or if equal to []. The
%    The default is all high channels, starting at 160, i.e. 160:191.
%    The numbering convention is that the first MEG channel is 0 (not 1).
if ~exist('triggerLinesPotential','var')
    triggerLinesPotential = [];
end
if isempty(triggerLinesPotential)
    triggerLinesPotential = 160:191;
end

% Create triggerHighThresh if not supplied or if equal to [].
if ~exist('triggerHighThresh','var')
    triggerHighThresh = [];
end
if isempty(triggerHighThresh)
    triggerHighThresh = 5; % mV
end

% Create triggerLowThresh if not supplied or if equal to [].
if ~exist('triggerLowThresh','var')
    triggerLowThresh = [];
end
if isempty(triggerLowThresh)
    triggerLowThresh = 0.2; % mV
end

triggerHighThresh_mV = triggerHighThresh*1000; % mV
triggerLowThresh_mV = triggerLowThresh*1000; % mV
clear triggerHighThresh triggerLowThresh

info = sqdread(sqdfile,'info');

fs = info.SampleRate; % Hz
triggerData = sqdread(sqdfile,'Channels',triggerLinesPotential);  % mV

% First find channels whose maximum is at least the "high" threshold and
% whose minimum is at most the "low" threshold.
triggerLinesTest = ( (max(triggerData) > triggerHighThresh_mV) & (min(triggerData) < triggerLowThresh_mV) );
triggerDataProbable = triggerData(:,triggerLinesTest);
triggerLinesProbable = triggerLinesPotential(triggerLinesTest);
clear triggerData

% (Not sure this is necessary.)
% Next, make sure that the time spent at rest is at least 100 times more
% than the time spent at the set value.
triggersNearHigh = mean(triggerDataProbable > triggerHighThresh_mV);
triggersNearLow = mean(triggerDataProbable < triggerLowThresh_mV);

triggerLinesHealthyHighRest =  (triggersNearHigh > 25 * triggersNearLow);
triggerDataHighRestGood = triggerDataProbable(:,triggerLinesHealthyHighRest);
triggerLinesHighRestGood = triggerLinesProbable(triggerLinesHealthyHighRest);
triggerLinesHighRestGoodCount = length(triggerLinesHighRestGood);

triggerLinesHealthyLowRest =  (triggersNearLow > 25 * triggersNearHigh);
triggerDataLowRestGood = triggerDataProbable(:,triggerLinesHealthyLowRest);
triggerLinesLowRestGood = triggerLinesProbable(triggerLinesHealthyLowRest);
triggerLinesLowRestGoodCount = length(triggerLinesLowRestGood);

if triggerLinesHighRestGoodCount > triggerLinesLowRestGoodCount
    triggerDataGood = triggerDataHighRestGood;
    triggerLinesGood = triggerLinesHighRestGood;
    triggerLinesGoodCount = triggerLinesHighRestGoodCount;
else
    triggerDataGood = triggerDataLowRestGood;
    triggerLinesGood = triggerLinesLowRestGood;
    triggerLinesGoodCount = triggerLinesLowRestGoodCount;
end

clear triggerDataProbable

% Inititialize outputs that hold trigger info
triggerCounts = nan(triggerLinesGoodCount,1);
triggersList = cell(triggerLinesGoodCount,1);

for iTrigger = 1:triggerLinesGoodCount
    if triggerLinesHighRestGoodCount > triggerLinesLowRestGoodCount
        triggerSet = find(triggerDataGood(:,iTrigger) < triggerLowThresh_mV);
    else
        triggerSet = find(triggerDataGood(:,iTrigger) > triggerHighThresh_mV);
    end
    % the next line uses find and diff to discard samples of triggerSet
    % that immediately follow any other values that cross
    % the "set" threshold.
    triggerSamples = (triggerSet([1; find(diff(triggerSet) > 1)+1]));
    triggerCounts(iTrigger) = length(triggerSamples);
    triggersList{iTrigger} = triggerSamples;
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip