text:
function mtrResampleSISPathways(fgFilenames,samplerOptsFilename,numSamples,fgResampledFilename,bCalcStats)

% Run metrotrac
% mtrResampleSISPathways({'paths_newpdf.dat'},'met_params.txt',1000000,'res
% amp_newpdf.dat',1);
%% MUST MAKE SURE THIS PROGRAM IS COMPILED AND IN RIGHT LOCATION
if(ispc)
    executable = which('updatestats.exe');
else
    error('Not compiled for Linux.');
end

% Get importance weight file that we need to resample with
for ii = 1:length(fgFilenames)
    statvec_names{ii} = sprintf('temp%d_statvec_iw.dat',ii);
    if( bCalcStats )
        args = sprintf(' -v %s -i %s -o temp.dat -s 4 -sf %s', samplerOptsFilename, fgFilenames{ii}, statvec_names{ii});
        cmd = [executable args];
        disp(cmd); disp('...')
        [s, ret_info] = system(cmd,'-echo');
        disp('Done')
    end
end

% Resample
%resampleSISPathways(numSamples, {'temp_statvec_iw.dat'}, {fgFilename},
%fgResampledFilename);
resampleSISPathways(numSamples, statvec_names, fgFilenames, fgResampledFilename);

if( bCalcStats )
    % Cleanup temp files
    delete('temp.dat');
end

% for ii = 1:length(statvec_names)
%     delete(statvec_names{ii});
% end

return;

function resampleSISPathways(numSamples, weightFilenames, pdbInFilenames, pdbOutFilename)

numFiles = length(weightFilenames);
sizeFiles = [];

% Load multiple importance weight files
iw = [];
for f = 1:numFiles
    iwtemp = mtrLoadStatvec(weightFilenames{f});
    iw = [iw; iwtemp];
    sizeFiles(f) = length(iwtemp);
end
iw(isnan(iw)) = -Inf;
iw = exp(iw);

% Resample pathways according to normalized weight
sample_ind = randsample([1:length(iw)],numSamples,true,iw);
unique_sample_ind = unique(sample_ind);
unique_sample_w = 0;
for i = 1:length(unique_sample_ind)
    unique_sample_w(i) =  sum(sample_ind == unique_sample_ind(i));
end

pdbtemp = [];
% Load the pathways and create resampled database
for f = 1:numFiles
    disp('Loading Pathway Database...');
    pdbtemp{f} = mtrPathwayDatabase();
    pdbtemp{f} = loadDatabase(pdbtemp{f},pdbInFilenames{f});
end

% create new database
pdb = mtrPathwayDatabase();
pdb = pdbtemp{1};
sub_unique_sample_ind = unique_sample_ind( unique_sample_ind <= length(pdbtemp{1}.pathways) & unique_sample_ind > 0);
pdb.pathways = pdbtemp{1}.pathways(sub_unique_sample_ind);

for f = 2:numFiles
    unique_sample_ind = unique_sample_ind - sizeFiles(f-1);
    sub_unique_sample_ind = unique_sample_ind( unique_sample_ind <= length(pdbtemp{f}.pathways) & unique_sample_ind > 0);
    pdb.pathways = [pdb.pathways pdbtemp{f}.pathways(sub_unique_sample_ind)];
end
% Clear statistic headers
pdb.pathway_statistic_headers = [];

% Create statistic header of just new weights
is_luminance_encoding = 1;
is_computed_per_point = 0;
is_viewable_stat = 1;
agg_name = 'resampleWeight';
local_name = 'NA';
uid = '1';
pdb = addStatisticHeader(pdb,agg_name,local_name,is_luminance_encoding, is_computed_per_point, is_viewable_stat, 1);

% Add the weight data to the path
for p = 1:length(pdb.pathways)
    pdb.pathways(p).path_stat_vector = unique_sample_w(p);
end

% Save out resampled database
disp('Saving new resampled database...');
saveDatabase(pdb,pdbOutFilename);

return;

author: vistalab

repo: vistasoft


text:
function tw = trigauss_conversion ( xw, omega )

%*****************************************************************************80
%
%% TRIGAUSS_CONVERSION converts Gauss to trigonometric Gauss quadrature.
%
%  Modified:
%
%    19 May 2013
%
%  Author:
%
%    Gaspare Da Fies, Alvise Sommariva, Marco Vianello
%
%  Parameters:
%
%    Input, real XW(N,2), the points and weights of a Gauss quadrature rule
%    on [-1,+1].
%
%    Input, real OMEGA, the arc angle, 0 < OMEGA <= PI.
%
%    Output, real TW(N,2), the angles and weights for the trigonometic
%    quadrature rule over [-OMEGA,+OMEGA].
%
  tw(:,1) = 2.0 * asin ( sin ( omega / 2.0 ) * xw(:,1) );
  tw(:,2) = xw(:,2);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% use EKF to estimate gyro and accel biases, for
% MEMS IMU grade close to microstrain 3dm gx3-35
% s0 frame of the sensor coincides a earth-fixed n-frame(NED) at epoch t0, n0.

% The state
% \f$ p_S^N \f$ position of the IMU in n0-frame
% \f$ v_N \f$ velocity of the IMU in n0-frame
% \f$ q_N^S \f$ quaternion representing the rotation from n0-frame to IMU frame
% \f$ b_a \f$ accelerometer bias
% \f$ b_g \f$ gyroscope bias

% The error state
% \f$ \delta p_S^N = \hat{p}_S^N - p_S^N \f$
% \f$ \delta v^N = \hat{v}^N - v^N \f$
% \f$ \psi^E \f$ Error in the rotation from n0-frame to IMU frame $C_N^S$
% defined as \f$ C_N^S = (I + \psi \times) \hat{C}_N^S \f$
% \f$ \delta b_a = b_a - \hat{b}_a \f$ accelerometer bias
% \f$ \delta b_g = b_g - \hat{b}_g \f$ gyro bias

classdef EKF_filter_s0frame_bias < handle
    properties (Hidden)
        type = 'ekf';
        tag  = 'EKF_GPS_IMU_GRAVITY';  % ID tag
    end
    % The following properties can be set only by class methods
    properties (SetAccess = private)
        stateTime; % state timestamp
        invalidateIMUerrors; % reset IMU errors when they goes beyond expectation
        imuType; % the type of IMU, e.g., MEMS 3DX GM 3-35, Steval inems, etc.
        imuErrorModel=3; % the IMU error model corresponding to biases
        dt; % IMU sampling interval, unit sec
        rqs02e; % position of s0 in the e frame, qs02e, constant
        rvqs0; % position of current s(k) frame in s0 frame, velocity in s0 frame, q s0 to s(k)
        imuErrors; %ba, bg
        imuOrientSIP=7; % imu orientation covariance start index in s-frame
        imuBiasDriftSIP=10; %imu bias error covariance start index
        covDim=15; % the dimension of P matrix
        velNoiseStd; % noise density driving horizontal velocity
        p_k_k; % the covariance of the entire state vector
    end
    methods
        function filter = EKF_filter_s0frame_bias(options)
            filter.stateTime = options.startTime;
            filter.invalidateIMUerrors=options.InvalidateIMUerrors;
            filter.imuType=options.imutype;
            filter.imuErrorModel=options.imuErrorModel;
            filter.dt=options.dt;

            inillh=options.inillh_ant;
            Ce2n0=llh2dcm_v000(inillh(1:2),[0;1]);
            filter.rqs02e=[ecef2geo_v000(options.inillh_ant,1); rotro2qr(Ce2n0')];
            % assume the IMU is very close to antenna
            filter.rvqs0=zeros(10,1);

            qs2n=quatmult_v001(options.qb2n,rotro2qr(options.Cb2imu),2);
            filter.rvqs0(7:10)=[qs2n(1);-qs2n(2:4)];
            filter.rvqs0(4:6)=options.Vn;
            filter.rvqs0(1:3) = quatrot_v000(qs2n, -options.imu_p_ant, 0);

            filter.imuErrors=options.imuErrors;
            filter.velNoiseStd=options.velNoiseStd;
            filter.p_k_k=zeros(filter.covDim);    % in s0 frame, position, velocity, attitude are of very small variance

            filter.p_k_k(1:3,1:3)=diag([0.2, 0.2, 0.5].^2); %position errors in meter
            filter.p_k_k(4:6,4:6)=diag([0.1, 0.1, 0.2].^2); %vel error in m/s
            filter.p_k_k(filter.imuOrientSIP+(0:2),filter.imuOrientSIP+(0:2))=...
                diag(([options.initAttVar,options.initAttVar,options.initAttVar*3]).^2);
            IMU_ERRDEF=imu_err_defs_v000(options.imutype);
            filter.p_k_k(filter.imuBiasDriftSIP+(0:2),filter.imuBiasDriftSIP+(0:2))=4*eye(3)*IMU_ERRDEF.acc_bias_var;% enlarge initial std by 2
            filter.p_k_k(filter.imuBiasDriftSIP+(3:5),filter.imuBiasDriftSIP+(3:5))=4*eye(3)*IMU_ERRDEF.gyro_bias_var;

            filter.p_k_k=sparse (filter.p_k_k);
        end
        %===============================================================================================
        %-- State transition function
        % propogate state with accelerometer and gyro input at time k-1 to predict state at k,
        % i.e., X(k|k-1), from state at k-1. U1 contains IMU measurement
        % acc, and gyro angular rate, or delta v, and delta theta
        % and 7th row is the previous epoch(k-1) and 8th is the current
        % epoch(k). Optionally, 9-11th rows are preset gravity in s0 frame,
        % 12-14 th row is $\omega_{ie}_{s_0}$, angular rate of local s0-frame w.r.t i-frame
        % represented in s0-frame. We also call s0 frame as w-frame
        % imuaccum records the accumulated delta v and delta angle for two speed
        % covariance update
        function imuaccum= ffun_state(filter, imuaccum, U1, isDelta, isConstantVel)
            if(isDelta)
                dt1=U1(8,end)-U1(7,1);
                gyroinc=sum(U1(4:6,:),2);
                accinc=sum(U1(1:3,:),2);

                %reset the imu errors if estimates diverges
                IMU_ERRDEF=imu_err_defs_v000(filter.imuType);

                if(filter.invalidateIMUerrors)
                    spuracc=find(abs(filter.imuErrors(1:3))>2*IMU_ERRDEF.initacc_bias_err,1);
                    spurgyro=find(abs(filter.imuErrors(4:6))>2*IMU_ERRDEF.initgyro_bias_err,1);
                    if(~isempty(spuracc)||~isempty(spurgyro))
                        disp(['IMU errors diverge at ' num2str(U1(8,end)) '!']);
                        filter.imuErrors=zeros(6,1);
                    end
                end
                angleinc=gyroinc-filter.imuErrors(4:6)*dt1;
                velinc=accinc-filter.imuErrors(1:3)*dt1;
                %imu data accumulator for the covariance update
                imuaccum=imuaccum+[velinc;angleinc];
                gyroinc=angleinc/dt1;
                accinc=velinc/dt1;
            else
                dt1=U1(8,end)-U1(7,1);
                gyroinc=mean(U1(4:6,:),2);
                accinc=mean(U1(1:3,:),2);
                %reset the imu errors if estimates diverges
                IMU_ERRDEF=imu_err_defs_v000(filter.imuType);
                if(filter.invalidateIMUerrors)
                    spuracc=find(abs(filter.imuErrors(1:3))>2*IMU_ERRDEF.initacc_bias_err,1);
                    spurgyro=find(abs(filter.imuErrors(4:6))>2*IMU_ERRDEF.initgyro_bias_err,1);
                    if(~isempty(spuracc)||~isempty(spurgyro))
                        disp(['IMU errors diverge at ' num2str(U1(8,end)) '!']);
                        filter.imuErrors=zeros(6,1);
                    end
                end
                gyroinc=gyroinc-filter.imuErrors(4:6);
                accinc=accinc-filter.imuErrors(1:3);
                %imu data accumulator for the covariance update
                imuaccum=imuaccum+[accinc;gyroinc]*dt1;
            end
            if(size(U1,1)>8)
                gwomegaw=U1(9:end,:);
            else
                %Gravity (most time consuming part of ecef implementations)
                xyz_imu=filter.rqs02e(1:3)+quatrot_v000(filter.rqs02e(4:7),filter.rvqs0(1:3),0);
                Llh=ecef2geo_v000(xyz_imu,0);
                Cn2e=pos2Cne_v000(Llh(1), Llh(2));
                [Rn, Re, gn, sL, cL, WIE_E]=geoparam_v000(Llh);
                gwomegaw=[quatrot_v000(filter.rqs02e(4:7), Cn2e*[0;0;gn], 1);quatrot_v000(filter.rqs02e(4:7),[0;0;WIE_E],1)];
            end
            if(~isConstantVel)
                filter.rvqs0=strapdown_local_quat_bias(filter.rvqs0, filter.rqs02e, accinc, gyroinc, dt1, gwomegaw);
            else
                filter.rvqs0=strapdown_local_quat_constvel(filter.rvqs0, filter.rqs02e, accinc, gyroinc, dt1, gwomegaw);
            end
            filter.stateTime = U1(8);
        end
        function ffun_covariance(filter, imuaccum, covupt_time, curimutime, isConstantVel )
            % propagate the covariance in the local s0 frame
            % the covariance corresponds to states,
            % rs in s0, v s in s0, q s0 2s, ba, bg
            covdt=curimutime-covupt_time;
            if(~isConstantVel)
            [STM, Qd]=sys_local_dcm_bias(filter.rqs02e, filter.rvqs0,...
                imuaccum(1:3)/covdt,imuaccum(4:6)/covdt, covdt,filter.imuType, filter.imuErrorModel);
            else
                [STM, Qd]=sys_local_dcm_constvel(filter.rqs02e, filter.rvqs0,...
                imuaccum(1:3)/covdt,imuaccum(4:6)/covdt, covdt,filter.imuType, filter.velNoiseStd);
            end

            Pvf = STM*filter.p_k_k*STM'+Qd;  % the covariance of the navigation states and imu error terms
            filter.p_k_k=sparse(Pvf);
        end
        %==============================================================================================
        function applied = correctstates(filter, predict, measure, H, R, measurementTime, gatingtest)
            % Suppose the prediction function is h, and the measurement is
            % z, then z = h(p, v, q, ba, bg). Formally,
            % \f$ H_p = - \frac{\partial h(p \oplus \delta)}{\partial \delta} \f$
            % \f$ H_\psi = - \frac{\partial h(C(q) \oplus \psi)}{\partial
            % \psi} \f$. Note in this update function,
            % \f$ p \oplus \delta = p - \delta \f$, and
            % \f$ C(q) \oplus \psi = \exp(\psi) C(q) \f$.

            if nargin < 7
                gatingtest = true;
            end
            if nargin < 6
                measurementTime = filter.stateTime;
            else
                if (measurementTime > filter.stateTime)
                    fprintf('Warn: Measurement time %.6f  >= state time %.6f\n', measurementTime, filter.stateTime);
                end
            end
            p_km1_k=filter.p_k_k;
            inno=predict-measure;
            S=H*p_km1_k*H'+R;
            if gatingtest
                gamma=inno'/S*inno;
                tol = chi2inv(0.99, 3);
                if gamma > tol
                    fprintf(['Discard measurement at %.6f with large ', ...
                             'gamma %.4f > tol %.4f.\n'], measurementTime, gamma, tol);
                    applied = false;
                    return;
                end
            end
            K=p_km1_k*H'/S;
            deltaX=K*inno;

            filter.p_k_k=(eye(size(p_km1_k,1))-K*H)*p_km1_k*(eye(size(p_km1_k,1))-K*H)'+K*R*K';
            % compute updated states
            filter.rvqs0(1:6)=filter.rvqs0(1:6)-deltaX(1:6); % position and velocity
            qst=rvec2quat_v000(deltaX(7:9));
            filter.rvqs0(7:10)=quatmult_v000(qst, filter.rvqs0(7:10));
            filter.imuErrors = filter.imuErrors + deltaX(filter.imuBiasDriftSIP:end);
            applied = true;
        end

        function SaveToFile(filter, ~, preimutime, ffilres)
            formatString = ['%.8f,%.5f,%.5f,%.5f,%.5f,%.5f,%.5f,%.8f,%.8f,%.8f,%.8f,', ...
                '%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.8f,%.8f,%.8f\n'];
            fprintf(ffilres, formatString, [preimutime;filter.rvqs0(1:6); ...
                [-filter.rvqs0(8:10); filter.rvqs0(7)];...
                full(sqrt(diag(filter.p_k_k(1:9,1:9))))]);
        end
        function mag=GetVelocityMag(filter)
            mag=norm(filter.rvqs0(4:6),2);
        end
    end
end

repo: ekfmonoslam


text:
D = dir;

numtests = 0;
results = [];

for di = 1:length(D)
    f = D(di);
    if f.isdir
        continue;
    end

    [pathstr,name,ext] = fileparts(f.name);

	if ~strcmp(ext,'.m')
		continue
    end;

    if strcmp('rtest_all.m', f.name)
        continue;
    end

    if strmatch('rtest',f.name)
        numtests = numtests+1;

        [pathstr,name,ext] = fileparts(f.name);
        rval = eval(name);
        results(numtests) = rval;

        switch(rval)
            case 1
                fprintf('%s passed\n', name);
            case 0
                fprintf('** %s FAILED!!\n', name);
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
if strcmp(computer, 'PCWIN')
  mex midiIn.c  -I../portmidi/pm_common libportmidi_win32.a C:\MinGW\lib\libwinmm.a
  mex midiOut.c -I../portmidi/pm_common libportmidi_win32.a C:\MinGW\lib\libwinmm.a
elseif strcmp(computer, 'PCWIN64')
  mex midiIn.c  -I../portmidi/pm_common libportmidi_win64.a C:\MinGW64\x86_64-w64-mingw32\lib\libwinmm.a
  mex midiOut.c -I../portmidi/pm_common libportmidi_win64.a C:\MinGW64\x86_64-w64-mingw32\lib\libwinmm.a
elseif strcmp(computer, 'GLNX86')
  mex midiIn.c  -I../portmidi/pm_common ../portmidi/libportmidi_linux32.a
  mex midiOut.c -I../portmidi/pm_common ../portmidi/libportmidi_linux32.a
elseif strcmp(computer, 'GLNXA64')
  mex midiIn.c  -I../portmidi/pm_common ../portmidi/libportmidi_linux64.a
  mex midiOut.c -I../portmidi/pm_common ../portmidi/libportmidi_linux64.a
elseif strcmp(computer, 'MACI64')
  % portmidi can be installed from MacPorts with 'sudo port install portmidi' or
  % downloaded and compiled by hand.
  %
  % compiling the mex file requires linking to a number of frameworks. In principle
  % LDFLAGS="$LDFLAGS -framework CoreFoundation -framework CoreAudio -framework CoreMidi -framework Carbon'
  % should copy the existing flags over, but it does not.
  %
  % using the -v option I figured out the difference on the command line with and without the LDFLAG option
  % '-arch x86_64 -Wl,-syslibroot,/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.9.sdk -mmacosx-version-min=10.7 -bundle -Wl,-exported_symbols_list,/Applications/MATLAB_R2014b.app/extern/lib/maci64/mexFunction.map'
  mex midiIn.c  -I/opt/local/include -L/opt/local/lib -lportmidi_s LDFLAGS='-framework CoreFoundation -framework CoreAudio -framework CoreMidi -framework Carbon -arch x86_64 -Wl,-syslibroot,/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.9.sdk -mmacosx-version-min=10.7 -bundle -Wl,-exported_symbols_list,/Applications/MATLAB_R2014b.app/extern/lib/maci64/mexFunction.map'
  mex midiOut.c -I/opt/local/include -L/opt/local/lib -lportmidi_s LDFLAGS='-framework CoreFoundation -framework CoreAudio -framework CoreMidi -framework Carbon -arch x86_64 -Wl,-syslibroot,/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.9.sdk -mmacosx-version-min=10.7 -bundle -Wl,-exported_symbols_list,/Applications/MATLAB_R2014b.app/extern/lib/maci64/mexFunction.map'
else
  ft_error('Unsupported platform');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [hPa] = Torr2hPa(Torr)
% Convert pressure from torr to hectopascals
% Chad Greene 2012
hPa = Torr*1.33322;

author: Sable

repo: mcbench-benchmarks


text:
function hb_structure_write ( output_unit, ncol, mxtype, nnzero, neltvl, ...
  ptrfmt, indfmt, colptr, rowind )

%*****************************************************************************80
%
%% HB_STRUCTURE_WRITE writes the structure of an HB file.
%
%  Discussion:
%
%    The user should already have opened the file, and written
%    the header records.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2004
%
%  Author:
%
%    John Burkardt
%    Xiaoye Li, UC Berkeley.
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, file pointer OUTPUT_UNIT, the unit to which data is written.
%
%    Input, integer NCOL, the number of columns or elements.
%
%    Input, string MXTYPE, the matrix type, up to 3 characters
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, string PTRFMT, the format for writing pointers,
%    up to 16 characters.
%
%    Input, string INDFMT, the format for writing row indices,
%    up to 16 characters.
%
%    Input, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Input, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%

%
%  From the FORTRAN format string PTRFMT, determine the number of
%  integers per line, and their width.
%
  line = ptrfmt(2:length(ptrfmt));
  [ ptr_nperline, count, errmsg, nextindex ] = sscanf ( line, '%d' );
  line = line(nextindex+1:length(line));
  [ ptr_len, count, errmsg, nextindex ] = sscanf ( line, '%d' );
%
%  Column pointers.
%
%  ptrfmt has the form '(nIw)'
%  where n is ptr_nperline
%  and w is ptr_len.
%
  s = int2str ( ptr_len );
  ptr_format = [ '%' s 'd' ];

  for j = 1 : ncol + 1
    fprintf ( output_unit, ptr_format, colptr(j) );
    if ( mod ( j, ptr_nperline ) == 0 )
      fprintf ( output_unit, '\n' );
    end
  end

  if ( mod ( ncol + 1, ptr_nperline ) ~= 0 )
    fprintf ( output_unit, '\n' );
  end
%
%  From the FORTRAN format string INDFMT, determine the number of
%  integers per line, and their width.
%
  line = indfmt(2:length(indfmt));
  [ ind_nperline, count, errmsg, nextindex ] = sscanf ( line, '%d' );
  line = line(nextindex+1:length(line));
  [ ind_len, count, errmsg, nextindex ] = sscanf ( line, '%d' );
%
%  Row indices.
%
  s = int2str ( ind_len );
  ind_format = [ '%' s 'd' ];

  if ( mxtype(3) == 'A' )

    for j = 1 : nnzero
      fprintf ( output_unit, ind_format, rowind(j) );
      if ( mod ( j, ind_nperline ) == 0 )
        fprintf ( output_unit, '\n' );
      end
    end

  elseif ( mxtype(3) == 'E' )

    number = colptr(ncol) - colptr(1);

    for j = 1 : number
      fprintf ( output_unit, ind_format, rowind(j) );
      if ( mod ( j, ind_nperline ) == 0 )
        fprintf ( output_unit, '\n' );
      end
    end

  end

  return
end


text:
clear all; close all; clc

n=200; L=8;
x=linspace(0,L,n);
x1=x(1:100);   % train
x2=x(101:200); % test
n1=length(x1);
n2=length(x2);
ftrain=(x1.^2).';  % train parabola x=[0,4]
ftest=(x2.^2).';   % test parbola x=[4,5]
figure(1), subplot(3,1,1),
plot(x1,ftrain,'r',x2,ftest,'b','Linewidth',[2])
legend('','','Location','Northwest')
legend boxoff

M=30; % number of model terms
Eni=zeros(100,M); Ene=zeros(100,M);
for jj=1:M
    for j=1:jj
        phi_i(:,j)=(x1.').^(j-1); % interpolation key
        phi_e(:,j)=(x2.').^(j-1); % extrapolation key
    end

    f=(x.^2).';
    for j=1:100
        fni=(x1.^2+0.1*randn(1,n1)).'; % interpolation
        fne=(x2.^2+0.1*randn(1,n2)).'; % extrapolation

        ani=pinv(phi_i)*fni; fnai=phi_i*ani;
        Eni(j,jj)=norm(ftrain-fnai)/norm(ftrain);

        fnae=phi_e*ani;  % use loadings from x in [0,4]
        Ene(j,jj)=norm(ftest-fnae)/norm(ftest);
    end
end

subplot(3,2,3), boxplot(Eni), axis([0.5 30.5 0 0.7]), set(gca,'Xlim',[0.5 30.5],'Xtick',1:30,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30'})
subplot(3,2,4), boxplot(Eni), axis([0.5 30.5 0 0.02]), set(gca,'Xlim',[0.5 30.5],'Xtick',1:30,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30'})
subplot(3,2,5), boxplot(Ene), set(gca,'Xlim',[0.5 30.5],'Xtick',1:30,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30'})
subplot(3,2,6), boxplot(log(Ene+1)), axis([0.5 30.5 0 30]), set(gca,'Xtick',1:30,'Xticklabel',{'1','','','','5','','','','','10','','','','','15','','','','','20','','','','','25','','','','','30'})

author: dynamicslab


text:
function c = differ_solve ( n, stencil, order )

%*****************************************************************************80
%
%% DIFFER_SOLVE solves for finite difference coefficients.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of stencil points.
%
%    Input, real STENCIL(N), the stencil vector.
%    The entries in this vector must be distinct.
%    No entry of STENCIL may be 0.
%
%    Input, integer ORDER, the order of the derivative to be approximated.
%    1 <= ORDER <= N.
%
%    Output, real C(N,1), the coefficients to be used
%    to multiply U(STENCIL(I))-U(0), so that the sum forms an
%    approximation to the derivative of order ORDER, with error
%    of order H^(N+1-ORDER).
%
  a = differ_matrix ( n, stencil );

  b = zeros ( n, 1 );
  b(order) = 1.0;

  c = a \ b;

  return
end


text:
%                           runBenchmark1DGaussianMN.m
%
% This example benchmarks algorithms based on their ability to reconstruct
% a synthetic signal (random Gaussian) using synthetic measurements
% (random Gaussian).  The benchmark shows how the
% different methods behave as the runtime increases.
%
% This script does the following:
%
% 1. Set up parameters and create a list of algorithm structs.
%
% 2. Invoke the general benchmark function benchmarkPR. A graph of errors
% (under specified error metrics) of different algorithms at each level
% of allowed runtime will be shown.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

%% -----------------------------START----------------------------------

clc
clear
close all

%% 1.Set up parameters
% Choose x label (values shown on the x axis of the benchmark plot) and
% y label (values shown on the y-axis). The value on the x axis is the
% runtime. The value on the y axis is 'reconerror', which is the relative
% 2-norm difference between the true and recovered signal.
xitem = 'time';
xvalues = [.1 1]; % The time in seconds
yitem = 'reconerror';

% Choose Dataset and set up dataSet '1DGaussian' specific parameters
dataSet = '1DGaussian';

% Set up general parameters
params.verbose = false;
params.numTrials = 2;        % run several random trials for each scenario, and report average results
params.n = 500;              % num of unknown elements
params.m = 4*params.n;       % number of measurements
params.isComplex = true;     % use complex matrices? or just stick to real?
params.policy = 'median';

% Create a list of algorithms structs
wf = struct('initMethod','spectral','algorithm','wirtflow');
twf = struct('algorithm','twf');
rwf = struct('algorithm','rwf');
ampflow = struct('algorithm','amplitudeflow');
taf = struct('initMethod','orthogonal','algorithm','taf');
raf = struct('initMethod','weighted','algorithm','raf');
fienup = struct('algorithm','fienup');
gs = struct('algorithm','gerchbergsaxton');
cd = struct('algorithm','coordinatedescent','maxIters',300*2*params.n);
kac = struct('algorithm','kaczmarz','maxIters',1000);
pmax = struct( 'algorithm','phasemax', 'maxIters',1000);
plamp = struct('algorithm', 'phaselamp');
scgm = struct('algorithm','sketchycgm');
plift = struct('algorithm','phaselift','maxIters',1000);

% Grab your pick of algorithms.
algorithms = {raf,fienup,ampflow,plift,pmax,plamp};

% Run benchmark
benchmarkSynthetic(xitem, xvalues, yitem, algorithms, dataSet, params);

author: tomgoldstein

repo: phasepack-matlab


text:
function estimation_results = notREKF_SLAM(data, NumberOfSteps)
% R-EKF SLAM

% load pre-given data: odometry and observations
if nargin < 1
    load('./data.mat');
end

data_matrix = data.state;

odom_sigma = data.odom_sigma;
obsv_sigma = data.obsv_sigma;

% odoCov = data.odom_cov;   % constant variable
% obsCov = data.obse_cov;   % constant variable

%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% In every step, all elements of Estimation_X will be changed %%%%%%%%%%%%
Estimation_X.orientation = data.poses.orientation(1:3,1:3);
Estimation_X.position    = data.poses.position(:,1);
Estimation_X.cov         = sparse(6,6);
Estimation_X.landmarks   = [];       % the landmarks observed until this step (included), 4*N format, the 4-th row is the index
Estimation_X0.IndexObservedNew=[];
Estimation_X0.IndexObservedAlreadyThis=[];
%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%

% Initialize
if nargin < 2
    NumberOfSteps = max(data_matrix(:,4));  % step instead of pose,  hence, it does not include pose 0
elseif NumberOfSteps > max(data_matrix(:,4))
    NumberOfSteps = max(data_matrix(:,4));
end
estimation_results = cell(1, NumberOfSteps+1);
estimation_results{1} = Estimation_X;
row_idx = (data_matrix(:, end) <= NumberOfSteps+1);
data_matrix = data_matrix(row_idx, :);

for i = 0:NumberOfSteps
    IndexOfCurrentStepInDataMatrix = find(data_matrix(:,4) == i);
    m = size(IndexOfCurrentStepInDataMatrix, 1);
    if ( mod(i, 50) == 0 )
        disp(['Processing pose ', int2str(i)]);
    end
    % det(Estimation_X.cov)
    if i==NumberOfSteps-1
    a=1;
    end

    if i ~= NumberOfSteps
        OdometryFromThis2Next = data_matrix(IndexOfCurrentStepInDataMatrix(m-5):IndexOfCurrentStepInDataMatrix(m),1);
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 1 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 3 )];
           [Estimation_X] = notREKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end

        estimation_results{i+1} = Estimation_X;

%        propagation using odometry info
        [Estimation_X] = notREKF_propagate(Estimation_X, OdometryFromThis2Next, odom_sigma );

    else
        a=2;
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 1 ) , ...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 3 )];
            [Estimation_X] = notREKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end
        estimation_results{i+1} = Estimation_X;
    end
end
clearvars -except estimation_results

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function res = imGeodesicPropagation(img, varargin)
%IMGEODESICPROPAGATION Compute geodesic propagation for each foreground pixel
%
%   RES = imGeodesicPropagation(IMG);
%   IMG is a binary image. For each foreground pixel, the geodesic
%   progagation is defined as the maximum geodesic distance to another
%   pixel of the foreground. If the foreground is not connected this
%   distance equals infinity.
%
%   RES = imGeodesicPropagation(IMG, WEIGHTS);
%   use different weights for the computation of distances. See
%   imChamferDistance for further details.
%
%   Note:
%   * As the algorithm propagates geodesic distances from each foreground
%       pixel, the computation time may be expensive.
%   * the function is defined for both 2D and 3D images
%
%
%   Example
%     % Compute geodesic propagation in a L-shape
%     img = zeros(20, 20);
%     img(4:16, 4:9) = 1;
%     img(11:16, 4:16) = 1;
%     prop = imGeodesicPropagation(img);
%     imagesc(prop);
%     colormap([1 1 1 ; jet]);
%
%     % Compute geodesic propagation in a set of particles
%     prop = zeros(size(img));
%     lbl = bwlabel(img);
%     for i = 1:max(lbl(:))
%         prop = max(prop, imGeodesicPropagation(lbl==i));
%     end
%     imagesc(prop);
%     colormap([1 1 1 ; jet]);
%
%   See also
%   imGeodesics, imGeodesicDistanceMap, imGeodesicRadius, imGeodesicExtremities

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2009-05-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

img = img > 0;
res = zeros(size(img));

dim = size(img);
if length(dim) == 2
    for i = 1:dim(1)
        for j = 1:dim(2)
            if ~img(i,j)
                continue;
            end

            marker = false(size(img));
            marker(i, j) = true;

            dist = imGeodesicDistanceMap(img, marker, varargin{:});
            res(i, j) = max(dist(img));
        end
    end

elseif length(dim) == 3
    for k = 1:dim(3)
        for j = 1:dim(2)
            for i = 1:dim(1)
                if ~img(i,j,k)
                    continue;
                end

                marker = false(size(img));
                marker(i, j, k) = true;

                dist = imGeodesicDistanceMap3d(img, marker, varargin{:});
                res(i, j, k) = max(dist(img));
            end
        end
    end

else
    error('Requires a 2D or a 3D image');
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [x,M]=semidefQuadProgEqConst(B,a,A,b)
%%SEMIDEFQUADPROGEQCONST Solve a quadratic programming problem involving a
%       semidefinite matrix (or a positive defintie matrix) where there are
%       only equality constraints, no inequality constraints. Specifically,
%       this algorithm solves the optimization problem:
%        minimize_x x'*B*x+2*a'*x
%        subject to A*x=b
%       If B is positive semidefinite, then if a solution exists, an
%       infinite number actually exist. This function returns one solution
%       and a matrix that can be used to find the other solutions.
%
%INPUTS: B An nXn real positive definite or positive semidefinite,
%          symmetric matrix.
%        a An nX1 real vector.
%        A A numConstXn real matrix. This can be an omitted or empty matrix
%          can be passed if there are no constraints. numConst<=n. There
%          shouldn't be any redundant constraints.
%        b A numConstX1 real vector.
%
%OUTPUTS: x The nX1 solution.
%         M An nXn matrix. All other solutions are of the form
%           xOther=x+M*y, where y is a an nX1 real vector.
%
%This function implements the equations of Chapter IV of [1]. Note that
%based on the constraints, it is possible for no solution to exist. This
%function is run assuming that a solution exists; it does not check for
%soluion non-existence.
%
%EXAMPLE:
%This is example 4.1 in [1]. x will be [9;18;35]/20.
% B=[1, 2,-1;
%    2, 4,-2;
%   -1,-2,1];
% a=[2;4;1];
% A=[1,2,1];
% b=4;
% [x,M]=semidefQuadProgEqConst(B,a,A,b)
%One sees that the third row and column of M is all zeros. Also,
%M(:,2)=-M(:,1)/2, so the matrix is singular. Thus, the second column of M
%can be used as a basis. So, to get any new  point, we can use x+M(:,2)*v
%(for a scalar v) or we could use x+5*M(:,2)*v, which is the same as the
%solution given in [1].
%
%REFERENCES:
%[1] D. L. Nelson, "Quadratic programming techniques using matrix
%    pseudoinverses," Ph.D. dissertation, Texas Tech University, Lubbock,
%    TX, Dec. 1969.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(A))
    %If there are no constraints, call the unconstrained algorithm.
    [x,M]=semidefQuadProgUnconst(B,a);
    return;
end

n=length(a);
I=eye(n,n);
pinvA=pinv(A);
diffA=(I-pinvA*A);
C=diffA*B*diffA;
pinvC=pinv(C);

%A combination of Equation 4.3 and 4.5.
x=pinvA*b+diffA*pinv(C)*diffA*(a-B*pinvA*b);
M=diffA*(I-pinvC*C);

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 yFiltM = applyRecursGaussFilter(xM,coeffS,dim)
% function yFiltM = applyRecursGaussFilter(xM,coeffS,dim)
%
% This function applies recursive filter of order 4 to xM along the passed dim.
% coeffS is struture with filter coefficients.
% For example for LoG filter:
%     scan3M: 3D scanArray (for example, planC{indexS.scan}.scanArray - ctOffset).
%     pixelSizeV = [1,1,3]; % voxel size
%     sigmaVal = 3; % sigma
%     dim = 1; % filter along rows
%     derivativeOrder = 'zero';
%     coeffS.sigmad = sigmaVal / pixelSizeV(dim);
%     coeffS = setGaussOrder(coeffS,derivativeOrder);
%     filt3M = applyRecursGaussFilter(scan3M,coeffS,dim);
%
% APA, 6/12/2018

ykPlusM = zeros(size(xM));
ykMinusM = zeros(size(xM));

switch dim
    case 1

        % Causal
        ykPlusM(1,:,:) = coeffS.N0 * xM(1,:,:) + coeffS.N1 * xM(1,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(2,:,:) = coeffS.N0 * xM(2,:,:) + coeffS.N1 * xM(1,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(3,:,:) = coeffS.N0 * xM(3,:,:) + coeffS.N1 * xM(2,:,:) + coeffS.N2 * xM(1,:,:)...
            + coeffS.N3 * xM(1,:,:);
        ykPlusM(4,:,:) = coeffS.N0 * xM(4,:,:) + coeffS.N1 * xM(3,:,:) + coeffS.N2 * xM(2,:,:)...
            + coeffS.N3 * xM(1,:,:);

        ykPlusM(1,:,:) = ykPlusM(1,:,:) - coeffS.BN1 * xM(1,:,:) - coeffS.BN2 * xM(1,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(2,:,:) = ykPlusM(2,:,:) + -coeffS.D1 * ykPlusM(1,:,:) - coeffS.BN2 * xM(1,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(3,:,:) = ykPlusM(3,:,:)  + -coeffS.D1 * ykPlusM(1,:,:) -+ coeffS.D2 * ykPlusM(2,:,:) - coeffS.BN3 * xM(1,:,:)...
            - coeffS.BN4 * xM(1,:,:);
        ykPlusM(4,:,:) = ykPlusM(4,:,:) + -coeffS.D1 * ykPlusM(1,:,:) -+ coeffS.D2 * ykPlusM(2,:,:) -+ coeffS.D3 * ykPlusM(3,:,:)...
            - coeffS.BN4 * xM(1,:,:);

        for ind = 5:size(xM,dim)
            ykPlusM(ind,:,:) = coeffS.N0 * xM(ind,:,:) + coeffS.N1 * xM(ind-1,:,:) + coeffS.N2 * xM(ind-2,:,:)...
                + coeffS.N3 * xM(ind-3,:,:) -+ coeffS.D1 * ykPlusM(ind-1,:,:) -+ coeffS.D2 * ykPlusM(ind-2,:,:) ...
                -+ coeffS.D3 * ykPlusM(ind-3,:,:) -+ coeffS.D4 * ykPlusM(ind-4,:,:);
        end

        % Anticausal
        ykMinusM(end,:,:) = coeffS.M1 * xM(end,:,:) + coeffS.M2 * xM(end,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-1,:,:) = coeffS.M1 * xM(end,:,:) + coeffS.M2 * xM(end,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-2,:,:) = coeffS.M1 * xM(end-1,:,:) + coeffS.M2 * xM(end-1,:,:) + coeffS.M3 * xM(end,:,:)...
            + coeffS.M4 * xM(end,:,:);
        ykMinusM(end-3,:,:) = coeffS.M1 * xM(end-2,:,:) + coeffS.M2 * xM(end-2,:,:) + coeffS.M3 * xM(end-1,:,:)...
            + coeffS.M4 * xM(end,:,:);

        ykMinusM(end,:,:) = ykMinusM(end,:,:) - coeffS.BM1 * xM(end,:,:) - coeffS.BM2 * xM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-1,:,:) = ykMinusM(end-1,:,:) + -coeffS.D1 * ykMinusM(end,:,:) - coeffS.BM2 * xM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-2,:,:) = ykMinusM(end-2,:,:) + -coeffS.D1 * ykMinusM(end-1,:,:) +- coeffS.D2 * ykMinusM(end,:,:) - coeffS.BM3 * xM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);
        ykMinusM(end-3,:,:) = ykMinusM(end-3,:,:) + -coeffS.D1 * ykMinusM(end-2,:,:) +- coeffS.D2 * ykMinusM(end-1,:,:) +- coeffS.D3 * ykMinusM(end,:,:)...
            - coeffS.BM4 * xM(end,:,:);

        for ind = size(xM,1)-4:-1:1
            ykMinusM(ind,:,:) = coeffS.M1 * xM(ind+1,:,:) + coeffS.M2 * xM(ind+2,:,:) + coeffS.M3 * xM(ind+3,:,:)...
                + coeffS.M4 * xM(ind+4,:,:) - coeffS.D1 * ykMinusM(ind+1,:,:) - coeffS.D2 * ykMinusM(ind+2,:,:) ...
                - coeffS.D3 * ykMinusM(ind+3,:,:) - coeffS.D4 * ykMinusM(ind+4,:,:);
        end

    case 2

        % Causal
        ykPlusM(:,1,:) = coeffS.N0 * xM(:,1,:) + coeffS.N1 * xM(:,1,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,2,:) = coeffS.N0 * xM(:,2,:) + coeffS.N1 * xM(:,1,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,3,:) = coeffS.N0 * xM(:,3,:) + coeffS.N1 * xM(:,2,:) + coeffS.N2 * xM(:,1,:)...
            + coeffS.N3 * xM(:,1,:);
        ykPlusM(:,4,:) = coeffS.N0 * xM(:,4,:) + coeffS.N1 * xM(:,3,:) + coeffS.N2 * xM(:,2,:)...
            + coeffS.N3 * xM(:,1,:);

        ykPlusM(:,1,:) = ykPlusM(:,1,:) - coeffS.BN1 * xM(:,1,:) - coeffS.BN2 * xM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,2,:) = ykPlusM(:,2,:) + -coeffS.D1 * ykPlusM(:,1,:) - coeffS.BN2 * xM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,3,:) = ykPlusM(:,3,:) + -coeffS.D1 * ykPlusM(:,2,:) -+ coeffS.D2 * ykPlusM(:,1,:) - coeffS.BN3 * xM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);
        ykPlusM(:,4,:) = ykPlusM(:,4,:) + -coeffS.D1 * ykPlusM(:,3,:) -+ coeffS.D2 * ykPlusM(:,2,:) -+ coeffS.D3 * ykPlusM(:,1,:)...
            - coeffS.BN4 * xM(:,1,:);

        for ind = 5:size(xM,dim)
            ykPlusM(:,ind,:) = coeffS.N0 * xM(:,ind,:) + coeffS.N1 * xM(:,ind-1,:) + coeffS.N2 * xM(:,ind-2,:)...
                + coeffS.N3 * xM(:,ind-3,:) -+ coeffS.D1 * ykPlusM(:,ind-1,:) -+ coeffS.D2 * ykPlusM(:,ind-2,:) ...
                -+ coeffS.D3 * ykPlusM(:,ind-3,:) -+ coeffS.D4 * ykPlusM(:,ind-4,:);
        end

        % Anticausal
        ykMinusM(:,end,:) = coeffS.M1 * xM(:,end,:) + coeffS.M2 * xM(:,end,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-1,:) = coeffS.M1 * xM(:,end,:) + coeffS.M2 * xM(:,end,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-2,:) = coeffS.M1 * xM(:,end-1,:) + coeffS.M2 * xM(:,end-1,:) + coeffS.M3 * xM(:,end,:)...
            + coeffS.M4 * xM(:,end,:);
        ykMinusM(:,end-3,:) = coeffS.M1 * xM(:,end-2,:) + coeffS.M2 * xM(:,end-2,:) + coeffS.M3 * xM(:,end-1,:)...
            + coeffS.M4 * xM(:,end,:);

        ykMinusM(:,end,:) = ykMinusM(:,end,:) - coeffS.BM1 * xM(:,end,:) - coeffS.BM2 * xM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-1,:) = ykMinusM(:,end-1,:) + -coeffS.D1 * ykMinusM(:,end,:) - coeffS.BM2 * xM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-2,:) = ykMinusM(:,end-2,:) + -coeffS.D1 * ykMinusM(:,end-1,:) +- coeffS.D2 * ykMinusM(:,end,:) - coeffS.BM3 * xM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);
        ykMinusM(:,end-3,:) = ykMinusM(:,end-3,:) + -coeffS.D1 * ykMinusM(:,end-2,:) +- coeffS.D2 * ykMinusM(:,end-1,:) +- coeffS.D3 * ykMinusM(:,end,:)...
            - coeffS.BM4 * xM(:,end,:);

        for ind = size(xM,dim)-4:-1:1
            ykMinusM(:,ind,:) = coeffS.M1 * xM(:,ind+1,:) + coeffS.M2 * xM(:,ind+2,:) + coeffS.M3 * xM(:,ind+3,:)...
                + coeffS.M4 * xM(:,ind+4,:) - coeffS.D1 * ykMinusM(:,ind+1,:) - coeffS.D2 * ykMinusM(:,ind+2,:) ...
                - coeffS.D3 * ykMinusM(:,ind+3,:) - coeffS.D4 * ykMinusM(:,ind+4,:);
        end

    case 3

        % Causal
        ykPlusM(:,:,1) = coeffS.N0 * xM(:,:,1) + coeffS.N1 * xM(:,:,1) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,2) = coeffS.N0 * xM(:,:,2) + coeffS.N1 * xM(:,:,1) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,3) = coeffS.N0 * xM(:,:,3) + coeffS.N1 * xM(:,:,2) + coeffS.N2 * xM(:,:,1)...
            + coeffS.N3 * xM(:,:,1);
        ykPlusM(:,:,4) = coeffS.N0 * xM(:,:,4) + coeffS.N1 * xM(:,:,3) + coeffS.N2 * xM(:,:,2)...
            + coeffS.N3 * xM(:,:,1);

        ykPlusM(:,:,1) = ykPlusM(:,:,1) - coeffS.BN1 * xM(:,:,1) - coeffS.BN2 * xM(:,:,1) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,2) = ykPlusM(:,:,2) + -coeffS.D1 * ykPlusM(:,:,1) - coeffS.BN2 * xM(:,:,1) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,3) = ykPlusM(:,:,3) + -coeffS.D1 * ykPlusM(:,:,1) -+ coeffS.D2 * ykPlusM(:,:,2) - coeffS.BN3 * xM(:,:,1)...
            - coeffS.BN4 * xM(:,:,1);
        ykPlusM(:,:,4) = ykPlusM(:,:,4) + -coeffS.D1 * ykPlusM(:,:,1) -+ coeffS.D2 * ykPlusM(:,:,2) -+ coeffS.D3 * ykPlusM(:,:,3)...
            - coeffS.BN4 * xM(:,:,1);

        for ind = 5:size(xM,dim)
            ykPlusM(:,:,ind) = coeffS.N0 * xM(:,:,ind) + coeffS.N1 * xM(:,:,ind-1) + coeffS.N2 * xM(:,:,ind-2)...
                + coeffS.N3 * xM(:,:,ind-3) -+ coeffS.D1 * ykPlusM(:,:,ind-1) -+ coeffS.D2 * ykPlusM(:,:,ind-2) ...
                -+ coeffS.D3 * ykPlusM(:,:,ind-3) -+ coeffS.D4 * ykPlusM(:,:,ind-4);
        end

        % Anticausal
        ykMinusM(:,:,end) = coeffS.M1 * xM(:,:,end) + coeffS.M2 * xM(:,:,end) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-1) = coeffS.M1 * xM(:,:,end) + coeffS.M2 * xM(:,:,end) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-2) = coeffS.M1 * xM(:,:,end-1) + coeffS.M2 * xM(:,:,end-1) + coeffS.M3 * xM(:,:,end)...
            + coeffS.M4 * xM(:,:,end);
        ykMinusM(:,:,end-3) = coeffS.M1 * xM(:,:,end-2) + coeffS.M2 * xM(:,:,end-2) + coeffS.M3 * xM(:,:,end-1)...
            + coeffS.M4 * xM(:,:,end);

        ykMinusM(:,:,end) = ykMinusM(:,:,end) - coeffS.BM1 * xM(:,:,end) - coeffS.BM2 * xM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-1) = ykMinusM(:,:,end-1) + -coeffS.D1 * ykMinusM(:,:,end) - coeffS.BM2 * xM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-2) = ykMinusM(:,:,end-2) + -coeffS.D1 * ykMinusM(:,:,end-1) +- coeffS.D2 * ykMinusM(:,:,end) - coeffS.BM3 * xM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);
        ykMinusM(:,:,end-3) = ykMinusM(:,:,end-3) + -coeffS.D1 * ykMinusM(:,:,end-2) +- coeffS.D2 * ykMinusM(:,:,end-1) +- coeffS.D3 * ykMinusM(:,:,end)...
            - coeffS.BM4 * xM(:,:,end);

        for ind = size(xM,dim)-4:-1:1
            ykMinusM(:,:,ind) = coeffS.M1 * xM(:,:,ind+1) + coeffS.M2 * xM(:,:,ind+2) + coeffS.M3 * xM(:,:,ind+3)...
                + coeffS.M4 * xM(:,:,ind+4) - coeffS.D1 * ykMinusM(:,:,ind+1) - coeffS.D2 * ykMinusM(:,:,ind+2) ...
                - coeffS.D3 * ykMinusM(:,:,ind+3) - coeffS.D4 * ykMinusM(:,:,ind+4);
        end

end

yFiltM = ykPlusM + ykMinusM;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%                       plotRecoveredVSOriginal.m
%
% This function plots the real part of the recovered signal against
% the real part of the original signal.
% It is used in all the test*.m scripts.
%
% Inputs:
%       x:  a n x 1 vector. Recovered signal.
%       xt: a n x 1 vector. Original signal.

%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function plotRecoverdVSOriginal(x,xt)
    figure;
    scatter(real(x),real(xt)); hold on;
    plot([-3 3], [-3 3]);
    title('Visual Correlation of Recovered signal with True Signal')
    xlabel('Recovered Signal')
    ylabel('True Signal')
end


text:
%COREG_INV invert coreg translations for coefficients
%
%    Andy Hooper, Aug 2005
%
%   ======================================================================
%   04/2008 AH: Updated for compatibility with matlab 2008a
%   11/2008 AH: Deltaline/pixel added for Doris 3.96 compatibility
%   ======================================================================

cpmname=dir('CPM_Data.*');

load coreg_parms
Lrange=coreg_parms(1)-1
Prange=coreg_parms(2)-1
N=coreg_parms(3)^2
%osf=coreg_parms(4)
osf=1
n_ifg=coreg_parms(5)

G=sparse(0,n_ifg*12);
d=zeros(0,1);

for i=1:length(cpmname);

    thisname=cpmname(i).name;
    CPM_Data=load(thisname);

    posL1=((CPM_Data(:,2)-1)*4/Lrange)-2;
    posP1=((CPM_Data(:,3)-1)*4/Prange)-2;
    offL=CPM_Data(:,4);
    offP=CPM_Data(:,5);
    posL2=((CPM_Data(:,2)+offL-1)*4/Lrange)-2;
    posP2=((CPM_Data(:,3)+offP-1)*4/Prange)-2;
    n_pos=size(posL2,1);
    corrf=CPM_Data(:,6);
    stdev=sqrt(3/2/N).*sqrt(1-corrf.^2)/pi./corrf*osf^(3/2);
    weighting=1./stdev;
    ifgs=sscanf(thisname,'CPM_Data.%d.%d');
    ifg1=ifgs(1);
    ifg2=ifgs(2);

    Gblock1=[ones(size(posL1)),posL1,posP1,posL1.^2,posL1.*posP1,posP1.^2];
    Gblock1=Gblock1.*repmat(weighting,1,6);
    Gblock2=[ones(size(posL2)),posL2,posP2,posL2.^2,posL2.*posP2,posP2.^2];
    Gblock2=Gblock2.*repmat(weighting,1,6);
    Gnew=sparse(n_pos*2,n_ifg*12);
    if ifg1~=0
       Gnew(1:n_pos,(ifg1-1)*12+1:(ifg1-1)*12+6)=Gblock1;
       Gnew(n_pos+1:n_pos*2,(ifg1-1)*12+7:ifg1*12)=Gblock1;
    end
    Gnew(1:n_pos,(ifg2-1)*12+1:(ifg2-1)*12+6)=-Gblock2;
    Gnew(n_pos+1:n_pos*2,(ifg2-1)*12+7:ifg2*12)=-Gblock2;

    G=[G;Gnew];
    d=[d;weighting.*offL;weighting.*offP];
    %d=[d;offL;offP];
end

% coeff_s gives mapping of slave to master w.r.t. slave position
coeff_s=G\d;

load slave_corners.txt
deltaL=zeros(4,size(slave_corners,1));
deltaP=zeros(4,size(slave_corners,1));
for i = 1:size(slave_corners,1)
    posL1=((slave_corners(i,1:2))*4/Lrange)-2;
    posL1=[posL1(1);posL1(1);posL1(2);posL1(2)];
    posP1=((slave_corners(i,3:4))*4/Prange)-2;
    posP1=[posP1(1);posP1(2);posP1(1);posP1(2)];
    G=[ones(4,1),posL1,posP1,posL1.^2,posL1.*posP1,posP1.^2];
    deltaL(:,i)=G*coeff_s((i-1)*12+1:(i-1)*12+6);
    deltaP(:,i)=G*coeff_s((i-1)*12+7:(i)*12);
end
corner_offsets=[deltaL(:)';deltaP(:)'];
corner_offsets=corner_offsets(:);
save('corner_offsets.txt','-ascii','corner_offsets')

% we want mapping of slave to master w.r.t. master position
l=[-2:0.05:2];
p=l;
[Ls,Ps]=meshgrid(l,p);
Ls=Ls(:);
Ps=Ps(:);
nsynth=length(Ls);

Gblock=[ones(nsynth,1),Ls,Ps,Ls.^2,Ls.*Ps,Ps.^2];
Gsynth=sparse(2*nsynth,12);
Gsynth(1:nsynth,1:6)=Gblock;
Gsynth(nsynth+1:end,7:12)=Gblock;
coeff_m=zeros(size(coeff_s));

for i=1:n_ifg
    ifg_coeff_s=coeff_s((i-1)*12+1:i*12);
    dsynth=Gsynth*ifg_coeff_s;
    Lm=Ls+dsynth(1:nsynth)*4/Lrange;
    Pm=Ps+dsynth(nsynth+1:end)*4/Prange;
    Gblock=[ones(nsynth,1),Lm,Pm,Lm.^2,Lm.*Pm,Pm.^2];
    Gsynth(1:nsynth,1:6)=Gblock;
    Gsynth(nsynth+1:end,7:12)=Gblock;
    coeff_m((i-1)*12+1:i*12)=Gsynth\-dsynth;
end

save('coreg_coeffs.txt','-ascii','coeff_m')

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% EEG_DECODECHAN - given an input EEG dataset structure, output a new EEG data structure
%                retaining and/or excluding specified time/latency, data point, channel,
%                and/or epoch range(s).
% Usage:
%   >> [chaninds chanlist] = eeg_decodechan(chanlocs, chanlist);
%
% Inputs:
%   chanlocs  - channel location structure
%   chanlist  - list of channels, numerical indices [1 2 3 ...] or string
%               'cz pz fz' or cell array { 'cz' 'pz' 'fz' }. Can also be
%               a list of channel types (see below)
%   field     - ['labels'|'type'] channel field to match. Default is
%               'labels'
%   ignoremissing - [true|false] ignore channel in channel list that would
%               be missing in the dataset.
%
% Outputs:
%   chaninds  - integer array with the list of channel indices
%   chanlist  - cell array with a list of channel labels
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, 2009-
%
% see also: EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ chaninds, chanlist ] = eeg_decodechan(chanlocs, chanstr, field, ignoremissing)

if nargin < 2
    help eeg_decodechan;
    return;
end
if nargin < 3
    field = 'labels';
end
if nargin < 4
    ignoremissing = false;
end
if isstruct(chanlocs) && isfield(chanlocs, 'chanlocs')
    chanlocs = chanlocs.chanlocs;
end

if isempty(chanlocs) && ischar(chanstr)
    chaninds = str2num(chanstr);
    chanlist = chaninds;
    return;
end

if ischar(chanstr)
    % convert chanstr
    % ---------------
    chanstr(find(chanstr == ']')) = [];
    chanstr(find(chanstr == '[')) = [];
    chanlistnum = [];
    chanstr  = [ ' ' chanstr ' ' ];
    try
        chanlist = eval( [ '{' chanstr '}' ] );
        chanlistnum = cellfun(@str2double, chanlist);
        if isnumeric(chanlist{1}) chanlist = [ chanlist{:} ]; end
    catch
        chanlist = {};
        sp = find(chanstr == ' ');
        for i = 1:length(sp)-1
            c = chanstr(sp(i)+1:sp(i+1)-1);
            if ~isempty(c)
                chanlist{end+1} = c;
                if isnan(str2double(chanlocs(1).(field))) % channel labels are not numerical
                    if ~isnan(str2double(c))
                        chanlistnum(end+1) = str2double(c);
                    end
                end
            end
        end
        if length(chanlistnum) == length(chanlist)
            chanlist = chanlistnum;
        end
    end
else
    chanlist = chanstr;
end

% convert to values
% -----------------
chanval = 0;
if isnumeric(chanlist)
    chanval = chanlist;
end
% chanval  = [];
% if iscell(chanlist)
%     for ind = 1:length(chanlist)
%
%         valtmp = str2double(chanlist{ind});
%         if ~isnan(valtmp)
%              chanval(end+1) = valtmp;
%         else chanval(end+1) = 0;
%         end
%     end
% else
%     chanval = chanlist;
% end

% convert to numerical
% --------------------
if all(chanval) > 0
    chaninds = chanval;
    chanlist = chanval;
else
    chaninds = [];
    alllabs  = lower({ chanlocs.(field) });
    chanlist = lower(chanlist);
    for ind = 1:length(chanlist)
        indmatch = find(strcmp(alllabs,chanlist{ind})); %#ok<STCI>
        if ~isempty(indmatch)
            for tmpi = 1:length(indmatch)
                chaninds(end+1) = indmatch(tmpi);
            end
        else
            try
                eval([ 'chantmp = ' chanlist{ind} ';' ]);
                if isempty(chantmp)
                     error([ 'Channel ''' chanlist{ind} ''' not found' ]);
                else
                    if chantmp <= length(chanlist)
                        chaninds(end+1) = chantmp;
                    elseif ~ignoremissing
                        error([ 'Channel ''' chanlist{ind} ''' not found' ]);
                    end
                end
            catch
                if ~ignoremissing
                    error([ 'Channel ''' chanlist{ind} ''' not found' ]);
                end
            end
        end
    end
end
chaninds = sort(chaninds);
if ~isempty(chanlocs)
    if any(chaninds > length(chanlocs)) || any(chaninds <= 0)
        error('Channel index out of range')
    end
    chanlist = { chanlocs(chaninds).(field) };
else
    chanlist = {};
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function line_num = sphere_llt_grid_line_num ( lat_num, long_num )

%*****************************************************************************80
%
%% SPHERE_LLT_GRID_LINE_NUM counts lines for an LLT grid.
%
%  Discussion:
%
%    An LLT grid is a grid of triangles bounded by latitude and longitude
%    lines over the surface of a sphere in 3D.
%
%    The number returned is the number of pairs of points to be connected.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LAT_NUM, LONG_NUM, the number of latitude and
%    longitude lines to draw.  The latitudes do not include the North and South
%    poles, which will be included automatically, so LAT_NUM = 5, for instance,
%    will result in points along 7 lines of latitude.
%
%    Output, integer LINE_NUM, the number of grid lines.
%
  line_num = long_num * ( lat_num + 1 ) ...
           + long_num *   lat_num ...
           + long_num * ( lat_num - 1 );

  return
end


text:
function c = comp_fwt(f,h,a,J,ext)
%COMP_FWT Compute DWT using FWT
%   Usage:  c=comp_fwt(f,h,J,a,Lc,ext);
%
%   Input parameters:
%         f     : Input data - L*W array.
%         h     : Analysis Wavelet filters - cell-array of length *filtNo*.
%         J     : Number of filterbank iterations.
%         a     : Subsampling factors - array of length *filtNo*.
%         ext   : 'per','zero','even','odd' Type of the forward transform boundary handling.
%
%   Output parameters:
%         c     : Cell array of length M. Each element is Lc(m)*W array.
%

% This could be removed with some effort. The question is, are there such
% wavelet filters? If your filterbank has different subsampling factors following the first two filters, please send a feature request.
assert(a(1)==a(2),'First two elements of *a* are not equal. Such wavelet filterbank is not suported.');

% Time-reversed, complex conjugate impulse responses.
filtNo = length(h);
%hCell = cellfun(@(hEl) conj(flipud(hEl.h(:))),h,'UniformOutput',0);
hCell = cellfun(@(hEl) hEl.h(:),h,'UniformOutput',0);

if(strcmp(ext,'per'))
   % Initial shift of the filter to compensate for it's delay.
   % "Zero" delay transform is produced
   % offset = cellfun(@(hEl) 1-numel(hEl.h)-hEl.offset,h);
   offset = cellfun(@(hEl) hEl.offset,h);
elseif strcmp(ext,'valid')
   offset = -cellfun(@(hEl) numel(hEl.h)-1,h);
else
   % No compensation for the filter delay (filters are made causal with respect to the output sequences).
   % This creates relative shift between levels of coefficients.
   % Initial shift determines type of subsampling.
   % This is even subsampling. e.g. subs. [1,2,3,4,5,6] by a factor 3 becomes [3,6]
   % The number of output coefficients depends on it.
   offset = -(a-1);
   % For odd subsampling skip = 0; but it requires slight touches
   % elsewhere.
end

M = (filtNo-1)*J+1;
c = cell(M,1);
runPtr = M-filtNo+2;
ctmp = f;
for jj=1:J
    % Run filterbank
    ctmp = comp_filterbank_td(ctmp,hCell,a,offset,ext);
    % Bookkeeping
    c(runPtr:runPtr+filtNo-2) = ctmp(2:end);
    ctmp = ctmp{1};
    runPtr = runPtr - (filtNo - 1);
end
% Save final approximation coefficients
c{1} = ctmp;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function sos = sosprogram(vartable,decvartable)
% SOSPROGRAM --- Initialize a new sum of squares program.
%
% SOSP = sosprogram(VARTABLE,DECVARTABLE)
%
% SOSP is a new sum of squares program.
% VARTABLE is a vector of independent variables.
% DECVARTABLE is a vector of decision variables (optional).
%
% Both VARTABLE and DECVARTABLE are either symbolic or polynomial objects.
%

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%

% 12/24/01 - SP
% 01/07/02 - SP
% 02/21/02 - SP -- Symbolic polynomial
% 10/10/02 - SP -- Path checking

if ~exist('sedumi') & ~exist('sqlp') & ~exist('csdp') & ~exist('sdpnal')
    error('No SDP solvers found.');
end;

if ~exist('getpolysym')
    dd = which('sosprogram');
    dd = strrep(dd,'/sosprogram.m','/');
    dd = strrep(dd,'\sosprogram.m','\');
    pp = genpath(dd);
    addpath(pp);
end;

if isa(vartable,'sym')

	sos.var.num = 0;
	sos.var.type = {};
	sos.var.Z = {};
	sos.var.ZZ = {};
	sos.var.T = {};

	sos.expr.num = 0;
	sos.expr.type = {};
	sos.expr.At = {};
	sos.expr.b = {};
	sos.expr.Z = {};

	sos.extravar.num = 0;
	sos.extravar.Z = {};
	sos.extravar.ZZ = {};
	sos.extravar.T = {};
	sos.extravar.idx = {};

	sos.objective = [];        % 01/07/02

	sos.solinfo.x = [];
	sos.solinfo.y = [];
	sos.solinfo.RRx = [];
	sos.solinfo.RRy = [];
	sos.solinfo.info = [];

	sos.vartable = sym2chartable(vartable);    % 02/21/02
	if size(vartable,2) > 1
        vartable = vartable.';
	end;
	sos.symvartable = vartable;

    sos.varmat.vartable = '[]';                 %JA&GV 06/06/13
    sos.varmat.symvartable = [];
    sos.varmat.count = 0;

	if (nargin == 2 & ~isempty(decvartable))
        sos.objective = sparse(length(decvartable),1);
        sos.decvartable = sym2chartable(decvartable);     % 03/01/02
        setofCommaBrackets = [1 strfind(sos.decvartable,',') strfind(sos.decvartable,']')];%26/04/13

        if size(decvartable,2) > 1
            decvartable = decvartable.';%the transpose of a matrix of decision varibles
                                        %it is put under this form in
                                        %sos.symdecvartable (next line)
        end;
        sos.symdecvartable = decvartable;
        sos.var.idx{1} = length(find(sos.decvartable==','))+2;

	else
        sos.decvartable = '[]';
        sos.symdecvartable = [];
        sos.var.idx{1} = 1;

	end;

else

	sos.var.num = 0;
	sos.var.type = {};
	sos.var.Z = {};
	sos.var.ZZ = {};
	sos.var.T = {};

	sos.expr.num = 0;
	sos.expr.type = {};
	sos.expr.At = {};
	sos.expr.b = {};
	sos.expr.Z = {};

	sos.extravar.num = 0;
	sos.extravar.Z = {};
	sos.extravar.ZZ = {};
	sos.extravar.T = {};
	sos.extravar.idx = {};

	sos.objective = [];        % 01/07/02

	sos.solinfo.x = [];
	sos.solinfo.y = [];
	sos.solinfo.RRx = [];
	sos.solinfo.RRy = [];
	sos.solinfo.info = [];

    sos.vartable = sort(vartable.varname);

    sos.varmat.vartable = [];                 % PJS 9/9/2013
    sos.varmat.symvartable = [];
    sos.varmat.count = 0;

	if (nargin == 2 & ~isempty(decvartable))
        sos.objective = sparse(length(decvartable),1);
        sos.decvartable = sort(decvartable.varname);

        sos.var.idx{1} = length(decvartable)+1;
	else
        sos.decvartable = {};
        sos.var.idx{1} = 1;
	end;

end;

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [samples, energies, diagn] = hmc2(f, x, opt, gradf, varargin)
%HMC2   Hybrid Monte Carlo sampling.
%
%       Description
%       SAMPLES = HMC2(F, X, OPTIONS, GRADF) uses a  hybrid Monte Carlo
%       algorithm to sample from the distribution P ~ EXP(-F), where F is the
%       first argument to HMC2. The Markov chain starts at the point X, and
%       the function GRADF is the gradient of the `energy' function F.
%
%       HMC2(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to
%       be passed to F() and GRADF().
%
%       [SAMPLES, ENERGIES, DIAGN] = HMC2(F, X, OPTIONS, GRADF) also returns a
%       log of the energy values (i.e. negative log probabilities) for the
%       samples in ENERGIES and DIAGN, a structure containing diagnostic
%       information (position, momentum and acceptance threshold) for each
%       step of the chain in DIAGN.POS, DIAGN.MOM and DIAGN.ACC respectively.
%       All candidate states (including rejected ones) are stored in
%       DIAGN.POS. The DIAGN structure contains fields:
%
%       pos
%        the position vectors of the dynamic process
%       mom
%        the momentum vectors of the dynamic process
%       acc
%        the acceptance thresholds
%       rej
%        the number of rejections
%       stp
%        the step size vectors
%
%       S = HMC2('STATE') returns a state structure that contains
%       the used random stream and its state and the momentum of
%       the dynamic process. These are contained in fields stream,
%       streamstate and mom respectively. The momentum state is
%       only used for a persistent momentum update.
%
%       HMC2('STATE', S) resets the state to S. If S is an integer,
%       then it is used as a seed for the random stream and the
%       momentum variable is randomised. If S is a structure
%       returned by HMC2('STATE') then it resets the random stream
%       to exactly the same state.
%
%       See HMC2_OPT for the optional parameters in the OPTIONS structure.
%
%       See also
%       HMC2_OPT, METROP
%

%       Copyright (c) 1996-1998 Christopher M Bishop, Ian T Nabney
%       Copyright (c) 1998-2000,2010 Aki Vehtari

%  The portion of the HMC algorithm involving "windows" is derived
%  from the C code for this function included in the Software for
%  Flexible Bayesian Modeling written by Radford Neal
%  <http://www.cs.toronto.edu/~radford/fbm.software.html>.

% Global variable to store state of momentum variables: set by set_state
% Used to initialize variable if set

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

global HMC_MOM

%Return state or set state to given state.
if nargin <= 2
  if ~strcmp(f, 'state')
    error('Unknown argument to hmc2');
  end
  switch nargin
   case 1
    samples = get_state(f);
    return;
   case 2
    set_state(f, x);
    return;
  end
end

% Set empty options to default values
opt=hmc2_opt(opt);
% Reference to structures is much slower, so...
opt_nsamples=opt.nsamples;
opt_window=opt.window;
opt_steps=opt.steps;
opt_display = opt.display;
opt_persistence = opt.persistence;

if opt_persistence
  alpha = opt.decay;
  salpha = sqrt(1-alpha*alpha);
end

% Stepsizes, varargin gives the opt.stepsf arguments (net, x ,y)
% where x is input data and y is a target data.
if ~isempty(opt.stepsf)
  epsilon = feval(opt.stepsf,varargin{:}).*opt.stepadj;
else
  epsilon = opt.stepadj;
end

% Force x to be a row vector
x = x(:)';
nparams = length(x);

%Set up strings for evaluating potential function and its gradient.
f = fcnchk(f, length(varargin));
gradf = fcnchk(gradf, length(varargin));

% Check the gradient evaluation.
if (opt.checkgrad)
  % Check gradients
  feval('gradcheck', x, f, gradf, varargin{:});
end

% Matrix of returned samples.
samples = zeros(opt_nsamples, nparams);
% Return energies?
if nargout >= 2
  en_save = 1;
  energies = zeros(opt_nsamples, 1);
else
  en_save = 0;
end
% Return diagnostics?
if nargout >= 3
  diagnostics = 1;
  diagn_pos = zeros(opt_nsamples, nparams);
  diagn_mom = zeros(opt_nsamples, nparams);
  diagn_acc = zeros(opt_nsamples, 1);
else
  diagnostics = 0;
end

if (~opt_persistence | isempty(HMC_MOM) | nparams ~= length(HMC_MOM))
  % Initialise momenta at random
  p = randn(1, nparams);
else
  % Initialise momenta from stored state
  p = HMC_MOM;
end

% Evaluate starting energy.
E = feval(f, x, varargin{:});

% Main loop.
nreject = 0;         %number of rejected samples
window_offset=0;     %window offset initialised to zero
k = - opt.nomit + 1;       %nomit samples are omitted, so we store
while k <= opt_nsamples    %samples from k>0

  % Store starting position and momenta
  xold = x;
  pold = p;
  % Recalculate Hamiltonian as momenta have changed
  Eold = E;
  Hold = E + 0.5*(p*p');

  % Decide on window offset, if windowed HMC is used
  if opt_window>1
    window_offset=fix(opt_window*rand(1));
  end

  have_rej = 0;
  have_acc = 0;
  n = window_offset;
  dir = -1;           %the default value for dir (=direction)
                      %assumes that windowing is used
  while (dir==-1 | n~=opt_steps)

    %if windowing is not used or we have allready taken
    %window_offset steps backwards...
    if (dir==-1 & n==0)
      % Restore, next state should be original start state.
      if window_offset > 0
	x = xold;
	p = pold;
	n = window_offset;
      end
      %set dir for forward steps
      E = Eold;
      H = Hold;
      dir = 1;
      stps = dir;
    else
      if (n*dir+1<opt_window | n>(opt_steps-opt_window))
	% State in the accept and/or reject window.
	stps = dir;
      else
	% State not in the accept and/or reject window.
	stps = opt_steps-2*(opt_window-1);
      end
      % First half-step of leapfrog.
      p = p - dir*0.5*epsilon.*feval(gradf, x, varargin{:});
      x = x + dir*epsilon.*p;
      % Full leapfrog steps.
      for m = 1:(abs(stps)-1)
	p = p - dir*epsilon.*feval(gradf, x, varargin{:});
	x = x + dir*epsilon.*p;
      end
      % Final half-step of leapfrog.
      p = p - dir*0.5*epsilon.*feval(gradf, x, varargin{:});

      E = feval(f, x, varargin{:});
      H = E + 0.5*(p*p');

      n=n+stps;
    end

    if (opt_window~=opt_steps+1 & n<opt_window)
      % Account for state in reject window.  Reject window can be
      % ignored if windows consist of the entire trajectory.
      if ~have_rej
	rej_free_energy = H;
      else
	rej_free_energy = -addlogs(-rej_free_energy, -H);
      end
      if (~have_rej | rand(1) < exp(rej_free_energy-H));
	E_rej=E;
	x_rej=x;
	p_rej=p;
        have_rej = 1;
      end
    end

    if (n>(opt_steps-opt_window))
      % Account for state in the accept window.
      if ~have_acc
	acc_free_energy = H;
      else
	acc_free_energy = -addlogs(-acc_free_energy, -H);
      end
      if (~have_acc | rand(1) < exp(acc_free_energy-H))
	E_acc=E;
	x_acc=x;
	p_acc=p;
	have_acc = 1;
      end
    end
  end

  % Acceptance threshold.
  a = exp(rej_free_energy - acc_free_energy);

  if (diagnostics & k > 0)
    diagn_pos(k,:) = x_acc;
    diagn_mom(k,:) = p_acc;
    diagn_acc(k,:) = a;
  end
  if (opt_display > 1)
    fprintf(1, 'New position is\n');
    disp(x);
  end

  % Take new state from the appropriate window.
  if a > rand(1)
    % Accept
    E=E_acc;
    x=x_acc;
    p=-p_acc; % Reverse momenta
    if (opt_display > 0)
      fprintf(1, 'Finished step %4d  Threshold: %g\n', k, a);
    end
  else
    % Reject
    if k > 0
      nreject = nreject + 1;
    end
    E=E_rej;
    x=x_rej;
    p=p_rej;
    if (opt_display > 0)
      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', k, a);
    end
  end

  if k > 0
    % Store sample
    samples(k,:) = x;
    if en_save
      % Store energy
      energies(k) = E;
    end
  end

  % Set momenta for next iteration
  if opt_persistence
    % Reverse momenta
    p = -p;
    % Adjust momenta by a small random amount
    p = alpha.*p + salpha.*randn(1, nparams);
  else
    % Replace all momenta
    p = randn(1, nparams);
  end

  k = k + 1;
end

if (opt_display > 0)
  fprintf(1, '\nFraction of samples rejected:  %g\n', ...
    nreject/(opt_nsamples));
end

% Store diagnostics
if diagnostics
  diagn.pos = diagn_pos;   %positions matrix
  diagn.mom = diagn_mom;   %momentum matrix
  diagn.acc = diagn_acc;   %acceptance treshold matrix
  diagn.rej = nreject/(opt_nsamples);   %rejection rate
  diagn.stps = epsilon;    %stepsize vector
end

% Store final momentum value in global so that it can be retrieved later
if opt_persistence
  HMC_MOM = p;
else
  HMC_MOM = [];
end
return

function state = get_state(f)
%GET_STATE           Return complete state of sampler
%                    (including momentum)

global HMC_MOM
state.stream=setrandstream();
state.streamstate = state.stream.State;
state.mom = HMC_MOM;
return

function set_state(f, x)
%SET_STATE    Set complete state of sample
%
%          Description
%          Set complete state of sampler (including momentum)
%          or just set randn and rand with integer argument.
global HMC_MOM
if isnumeric(x)
  setrandstream(x);
else
  if ~isstruct(x)
    error('Second argument to hmc must be number or state structure');
  end
  if (~isfield(x, 'stream') | ~isfield(x, 'streamstate') ...
      | ~isfield(x, 'mom'))
    error('Second argument to hmc must contain correct fields')
  end
  setrandstream(x.stream);
  x.State=x.streamstate;
  HMC_MOM = x.mom;
end
return

function c=addlogs(a,b)
%ADDLOGS(A,B)  Add numbers represented by their logarithms.
%
%        Description
%        Add numbers represented by their logarithms.
%        Computes log(exp(a)+exp(b)) in such a fashion that it
%        works even when a and b have large magnitude.
if a>b
  c = a + log(1+exp(b-a));
else
  c = b + log(1+exp(a-b));
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Comparison of different learning methods of Hawkes processes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

options.N = 200; % the number of sequences
options.Nmax = 100; % the maximum number of events per sequence
options.Tmax = 100; % the maximum size of time window
options.tstep = 0.1;
options.dt = 0.1;
options.M = 250;
options.GenerationNum = 5;
D = 3; % the dimension of Hawkes processes
nTest = 1;
nSeg = 5;
nNum = options.N/nSeg;

disp('Approximate simulation of Hawkes processes via branching process')
disp('Simple exponential kernel')
para1.kernel = 'exp';
para1.w = 1;
para1.landmark = 0;
L = length(para1.landmark);
para1.mu = rand(D,1)/D;
para1.A = zeros(D, D, L);
for l = 1:L
    para1.A(:,:,l) = (0.5^l)*(0.5+rand(D));
end
para1.A = 0.5*para1.A./max(abs(eig(sum(para1.A,3))));
para1.A = reshape(para1.A, [D, L, D]);
Seqs1 = Simulation_Branch_HP(para1, options);
%Seqs1 = SimulationFast_Thinning_ExpHP(para1, options);

%%
disp('Learning Hawkes processes via different methods')
alg.LowRank = 0;
alg.Sparse = 1;
alg.alphaS = 1;
alg.GroupSparse = 0;
alg.outer = 5;
alg.rho = 0.1;
alg.inner = 8;
alg.thres = 1e-5;
alg.Tmax = [];
alg.storeErr = 0;
alg.storeLL = 0;

Err = zeros(nTest, nSeg);

for n = 1:nTest
    for i = nSeg

        [A, Phi] = ImpactFunc( para1, options );

        disp('Maximum likelihood estimation and basis representation')
        model1 = Initialization_Basis(Seqs1);
        model1 = Learning_MLE_Basis( Seqs1(1:i*nNum), model1, alg );
        [A1, Phi1] = ImpactFunc( model1, options );

        disp('Least squares and discretization')
        model2.D = D;
        model2.h = 1;
        model2.k = floor(options.M * options.dt/model2.h);
        model2 = Initialization_Discrete(Seqs1);
        model2 = Learning_LS_Discrete( Seqs1(1:i*nNum), model2 );

        figure
        title('Simple exponetial kernels')
        for u = 1:D
            for v = 1:D
                subplot(D,D,D*(u-1)+v)
                hold on
                plot(options.dt*(0:(size(Phi,2)-1)), Phi(v,:,u), 'k-')
                plot(options.dt*(0:(size(Phi1,2)-1)), Phi1(v,:,u), 'r-')
                plot(model2.h*(0:(size(model2.A,2)-1)), model2.A(v,:,u), 'b-')
                hold off
                axis tight
                legend('Real', 'MLE', 'LS')%, 'LS2')
                xlabel('Time interval between events')
                ylabel(['\phi', sprintf('%d%d', u, v)])
            end
        end

    end
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function pco = ge(p,q)
% function pco = ge(p,q)
%
% DESCRIPTION
%   Creates a polynomial constraint object representing the
%   constraint p>=q.  sosopt interprets this constraint as
%           p-q is sum of squares
%
% INPUTS
%   p: polynomial
%   q: polynomial
%
% OUTPUT
%   pco: polynomial constraint object
%
% SYNTAX
%   p>=q
%     Creates a polynomial constraint object.  If p is a scalar and
%     q is a vector then p will be expanded to have the same dimension as
%     q.  Similarly, if q is a scalar and p is a vector then q will be
%     expanded.
%

% 10/22/2010:   PJS  Initial Coding

% Create polynomial constraint object
pco = polyconstr(p,q,'>=');


text:
function s = setfield(s,field,v)
%SETFIELD Set mutable structure field contents.
% setfield(s,'field',v) sets the contents of the specified field to the
% value V.  This is equivalent to the syntax S.field = V.
%
% See also GETFIELD, RMFIELD, ISFIELD, FIELDNAMES.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

jv = s.obj;
if ~jv.containsKey(field)
  % add a new field
  jv.get('_fields').addElement(field);
end
jv.put(field,asJava(v));

author: andrewssobral

repo: lrslibrary


text:
classdef LvdDataPoint < AbstractGeometricPoint
    %LvdDataPoint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        allTimes cell = {};
        allRVVect cell = {};
        allFrames(1,:) AbstractReferenceFrame
        allGIs cell = {};

        name(1,:) char = 'New Point';

        inputLvdData LvdData
        lvdData LvdData

        %marker
        markerColor(1,1) ColorSpecEnum = ColorSpecEnum.Red;
        markerShape(1,1) MarkerStyleEnum = MarkerStyleEnum.RightTriangle;

        %track line
        plotTrkLine(1,1) logical = true;
        trkLineColor(1,1) ColorSpecEnum = ColorSpecEnum.Black;
        trkLineSpec(1,1) LineSpecEnum = LineSpecEnum.DottedLine;
    end

    methods
        function obj = LvdDataPoint(inputLvdData, lvdData, name)
            arguments
                inputLvdData(1,1) LvdData %the lvd data that contains the trajectory to be loaded
                lvdData(1,1) LvdData %the lvd data for the current case/scenario/mission
                name(1,:) char %the name of the point
            end

            obj.name = name;
            obj.lvdData = lvdData;

            obj.loadLvdData(inputLvdData, lvdData);
        end

        function tf = hasValidLvdData(obj)
            tf = true;

            if(isempty(obj.inputLvdData) || ...
               obj.inputLvdData == obj.lvdData)
                tf = false;

                return;
            end
        end

        function cartElems = getPositionAtTime(obj, inputTimes, ~, inFrame)
            cartElems = obj.getCartElemForTime(inputTimes, inFrame);
        end

        function loadLvdData(obj, inputLvdData, lvdData)
            obj.inputLvdData = inputLvdData;

            stateLog = inputLvdData.script.executeScript(false, inputLvdData.script.getEventForInd(1), true, false, false, false);

            obj.allTimes = {};
            obj.allRVVect = {};
            obj.allFrames = AbstractReferenceFrame.empty(1,0);
            obj.allGIs = {};

            for(i=1:inputLvdData.script.getTotalNumOfEvents()) %#ok<*NO4LP>
                evt = inputLvdData.script.getEventForInd(i);
                entries = stateLog.getAllStateLogEntriesForEvent(evt);

                [~,ia,~] = unique([entries.time], 'sorted');
                entries = entries(ia);

                if(numel(entries) > 2)
                    useFrame = entries(1).centralBody.getBodyCenteredInertialFrame();

                    times = [];
                    rvVects = [];
                    for(j=1:length(entries))
                        entry = entries(j);
                        maEntry = entry.getMAFormattedStateLogMatrix(true);

                        time = maEntry(1);
                        rVect = maEntry(2:4)';
                        vVect = maEntry(5:7)';

                        bodyInfo = lvdData.celBodyData.getBodyInfoById(maEntry(8));
                        frame = bodyInfo.getBodyCenteredInertialFrame();

                        ce = CartesianElementSet(time, rVect(:), vVect(:), frame);
                        ce = ce.convertToFrame(useFrame, true);

                        rvVect = [ce.rVect; ce.vVect];

                        times(j) = time; %#ok<AGROW>
                        rvVects(:,j) = rvVect; %#ok<AGROW>
                    end

                    obj.allTimes{end+1} = times;
                    obj.allRVVect{end+1} = rvVects;
                    obj.allFrames(end+1) = useFrame;
                    obj.allGIs{end+1} = griddedInterpolant(times,[rvVects;times]', "makima", "nearest");
                end
            end
        end

        function name = getName(obj)
            name = obj.name;
        end

        function setName(obj, name)
            obj.name = name;
        end

        function listboxStr = getListboxStr(obj)
            listboxStr = sprintf('%s (LVD Trajectory)', obj.getName());
        end

        function useTf = openEditDialog(obj, ~)
            output = AppDesignerGUIOutput({false});
            lvd_EditLvdTrajectoryPointGUI_App(obj, output, obj.lvdData);
            useTf = output.output{1};
        end

        function tf = isVehDependent(obj) %#ok<MANU>
            tf = false;
        end

        function tf = canBePlotted(obj) %#ok<MANU>
            tf = true;
        end

        function bodyInfo = getOriginBody(obj)
            bodyInfo = obj.bodyInfo;
        end

        function tf = usesGroundObj(~, ~)
            tf = false;
        end

        function tf = usesGeometricPoint(~, ~)
            tf = false;
        end

        function tf = usesGeometricVector(~, ~)
            tf = false;
        end

        function tf = usesGeometricCoordSys(~, ~)
            tf = false;
        end

        function tf = usesGeometricRefFrame(~, ~)
            tf = false;
        end

        function tf = usesGeometricAngle(~, ~)
            tf = false;
        end

        function tf = usesGeometricPlane(~, ~)
            tf = false;
        end

        function tf = isInUse(obj, lvdData)
            tf = lvdData.usesGeometricPoint(obj);
        end
    end

    methods(Access=private)
        function cartElem = getCartElemForTime(obj, qTimes, inFrame)
            fillerCe = CartesianElementSet(0, [0;0;0], [0;0;0], inFrame);
            cartElem = repmat(fillerCe, [1, numel(qTimes)]);

            for(i=1:length(obj.allTimes))
                times = obj.allTimes{i};
                minTime = min(times);
                maxTime = max(times);

                bool = qTimes >= minTime & qTimes <= maxTime;
                if(any(bool))
                    subQTimes = qTimes(bool);
                    gi = obj.allGIs{i};
                    frame = obj.allFrames(i);

                    rvVects = gi(subQTimes)';
                    subCe = CartesianElementSet(subQTimes, rvVects(1:3,:), rvVects(4:6,:), frame);
                    subCe = convertToFrame(subCe, inFrame, true);
                    cartElem(bool) = subCe;
                end
            end

            bool = cartElem == fillerCe;
            if(any(bool))
                outQTimes = qTimes(bool);

                for(i=1:length(obj.allTimes))
                    times = obj.allTimes{i};
                    minTime = min(times);
                    maxTime = max(times);

                    distToLb = abs(outQTimes - minTime);
                    distToUb = abs(outQTimes - maxTime);
                    distToBnd(:,i) = min([distToLb(:), distToUb(:)], [], 2); %#ok<AGROW>
                end

                boolInds = find(bool);

                [~,I] = min(distToBnd, [], 2);
                for(i=1:length(I))
                    qTime = outQTimes(i);
                    Ii = I(i);

                    gi = obj.allGIs{Ii};
                    frame = obj.allFrames(Ii);

                    rvVects = gi(qTime)';
                    rVects = rvVects(1:3,:);
                    vVects = rvVects(4:6,:);
                    nearestTimes = rvVects(7,:);
                    subCe = CartesianElementSet(nearestTimes, rVects, vVects, frame);
                    subCe = convertToFrame(subCe, inFrame, true);
                    cartElem(boolInds(i)) = subCe;
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function p24_title ( )

%*****************************************************************************80
%
% P24_TITLE prints a title for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 24\n' );
  fprintf ( 1, '  Name:       (|4X-2|+C)/(1+C)\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = prod ( (|4*X(i)-2|+C(i)) / (1+C(i)) )\n' );
  fprintf ( 1, '  Parameters:\n' );
  fprintf ( 1, '              C(1:DIM_NUM) defaults to 0.0\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function unicycle_test ( )

%*****************************************************************************80
%
%% UNICYCLE_TEST tests the UNICYCLE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNICYCLE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the UNICYCLE library.\n' );

  unicycle_test01 ( );
  unicycle_test02 ( );
  unicycle_test03 ( );
  unicycle_test04 ( );
  unicycle_test05 ( );
  unicycle_test06 ( );
  unicycle_test07 ( );
  unicycle_test08 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNICYCLE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% std_rmalldatafields - remove all data fields from STUDY (before saving
%                       it for instance.
%
% Usage:
%   STUDY = std_rmalldatafields(STUDY, type);
%
% Input:
%   STUDY   - EEGLAB study structure
%
% Optional input:
%   type    - ['chan'|'clust'|'both'] remove from changrp channel location
%             structure, cluster structure or both. Default is 'both'.
%
% Output:
%   STUDY   - updated EEGLAB study structure
%
% Author: Arnaud Delorme, CERCO/CNRS, UCSD, 2010-

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function STUDY = std_rmalldatafields(STUDY, chanorcomp);

if nargin < 1
    help std_rmalldatafields;
    return;
end
if nargin < 2
    chanorcomp = 'both';
end

fields = { 'erpdata'  'erptimes'  'erpdatatrials'  'erptrialinfo' ...
               'specdata' 'specfreqs' 'specdatatrials' 'spectrialinfo' ...
               'erspdata' 'erspbase'  'erspfreqs'      'ersptimes' 'erspdatatrials' 'ersptrialinfo'  ...
               'topo' 'topox' 'topoy' 'topoall' 'topopol' ...
               'itcdata'  'itcfreqs'  'itctimes'       'itcdatatrials' 'itctrialinfo'  ...
               'erpimdata' 'erpimtrials' 'erpimevents' 'erpimtimes' 'erspsubjinds' 'itcsubjinds' 'dipoles' ...
               'data' 'datatimes' 'datasortvals' 'datacontinds' 'centroid' };
for ff = 1:length(fields)
    if strcmpi(chanorcomp, 'data') || strcmpi(chanorcomp, 'both')
        if isfield(STUDY.changrp, fields{ff})
            STUDY.changrp = rmfield(STUDY.changrp, fields{ff} );
        end
    end
    if strcmpi(chanorcomp, 'clust') || strcmpi(chanorcomp, 'both')
        if isfield(STUDY.cluster, fields{ff})
            STUDY.cluster = rmfield(STUDY.cluster, fields{ff} );
        end
    end
end


text:
function y = upsample2(x,K)
for i = 1:size(x,3)
    y(:,:,i) = upsample(upsample(x(:,:,i),K)',K)';
end
end

author: cszn

repo: IRCNN

save_path: github-repos/MATLAB/cszn-IRCNN


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%
Sigma=(1./m).*X'*X;
[U, S, V]=svd(Sigma);

% =========================================================================

end

author: yhyap

repo: machine-learning-coursera


text:
%% Define Feedforward Convolutional Code
SNR_dB = 3
LM = 1600         % Mesg Length excluding pre-determined bits for starting & ending trellis at state 0
TREL_TYPE = 'Feedback'  % {'Feedback', 'Feedforward'}

%% Matlab Generator Polynomial convention
% Build a binary number representation by placing a 1 in each spot where a connection line from the shift register feeds into the adder,
% and a 0 elsewhere. The leftmost spot in the binary number represents the current input, while the rightmost spot represents the
% oldest input that still remains in the shift register.

if strcmp(TREL_TYPE, 'Feedback')
%     CL = 2		% Rate 1/2 Feedback encoder with 2 states
%     GenPoly0_1by2 = 3 % in octal
%     GenPoly1_1by2 = 2
%     FeedBackCoef = [1,1];  % For computing i/p bits needed to terminate trellis at state =0.
%     TREL = poly2trellis(CL, [GenPoly0_1by2, GenPoly1_1by2], GenPoly0_1by2)

%% Rate 1/2 Feedback encoder with 8-states used in 3GPP cellular 3G/4G standard
    CL = 4
    GenPoly0_1by2 = 13 % in octal
    GenPoly1_1by2 = 15
    FeedBackCoef = [1,0,1,1];% For computing i/p bits needed to terminate trellis at state =0.
    TREL = poly2trellis(CL, [GenPoly0_1by2, GenPoly1_1by2], GenPoly0_1by2)

else
    %% Rate 1/3 Feedforward encoder with 4 states
%     CL = 3; % constraint length
%     GenPoly0_third = 4; % in octal
%     GenPoly1_third = 5;
%     GenPoly2_third = 7;
%     TREL = poly2trellis( CL, [GenPoly0_third, GenPoly1_third, GenPoly2_third])

%     CL = 4 % constraint length
%     GenPoly0_1by4 = 13 % in octal
%     GenPoly1_1by4 = 15
%     GenPoly2_1by4 = 15
%     GenPoly3_1by4 = 17
%     TREL = poly2trellis( CL, [GenPoly0_1by4, GenPoly1_1by4, GenPoly2_1by4, GenPoly3_1by4])

end
LM = LM + 2*(CL-1); % space for start & tail bits

%% Verify trellis-structure is OK
[isok, status] = istrellis(TREL)

randn('state', sum(100*clock)); % initialize to random state

%% The encoder is assumed to have both started and ended at the all-zeros state
%% Ensure msg is s.t. trellis starts at "0" state and ends at "0" state.
%%  Generate Random binary stream & encode
if strcmp(TREL_TYPE, 'Feedback')
    msg1(1 : CL-1) = zeros(1, CL-1); % 1st (CL-1) bits must be 0
    msg1(CL : LM -CL+1) = randint(LM -2*CL +2, 1, 2)' ; % Random data
% 	msg1(CL : LM -CL+1) = [0 1];

    % Encode first part of msg, recording final state for later use.
    [cenc_o1, final_state1] = convenc(msg1, TREL);

    % Rest of msg depends on final_state1; it makes trellis terminate at final_state=0.
    bvec_fs = bitget(final_state1, (CL-1) : -1 : 1); % All possible binary-vectors of length EncoderN
    for idx = 1 : (CL-1)
        msg2(idx) = rem( [0, bvec_fs]*FeedBackCoef', 2);
        bvec_fs = [0, bvec_fs(1: (end-1))];
    end
%     msg2(1 : CL -1) = bvec_fs; % Last (CL-1) bits depend on state at time=(LM-CL+1)
    [cenc_o2, final_state] = convenc(msg2, TREL, final_state1);

    msg = [msg1, msg2]; clear msg1 msg2
    [cenc_o] = [cenc_o1, cenc_o2]; clear cenc_o1 cenc_o2
    final_state
else % Feedforward
    msg(1 : CL-1) = zeros(1,CL-1); % 1st (CL-1) bits must be 0
    msg(CL : LM -CL+1) = randint(LM -2*CL +2, 1, 2)' ; % Random data
    msg(LM -CL + 2 : LM) = zeros(CL-1, 1) % Last (CL-1) bits must be 0
    [cenc_o, final_state] = convenc(msg, TREL);
end

if final_state ~= 0
    disp('trellis not terminated properly: check last (CL-1) bits of mesg')
    return
end

%%	Map to BPSK constellation: bit0 -> 1, bit1 -> -1
signal_power = 1;
chan_in = (1 - 2*cenc_o)*sqrt(signal_power);

%%	Generate and add Gaussian-noise mean=0, variance= noise_power
noise_power = signal_power / 10^(SNR_dB/10);
noise = randn(size(chan_in))*sqrt(noise_power);
chan_o = chan_in + noise;
% chan_o = chan_in;

Es = log2(TREL.numOutputSymbols)*signal_power;
No = noise_power;
Lc = 4*Es/No

%%	Soft-Decision decoding, map to decision values
soft_in = chan_o;
% decodeds = vitdec_htm(TREL, soft_in');

[LLR, Alpha, Beta] = LogMAPdecode_htm(TREL, chan_o, Lc);

decd_msg = (1 - sign(LLR))/2;
err_vec = abs(decd_msg - msg);

sprintf('# of Bit-Errors = %d out of %d info-bits ', sum(abs(err_vec)), LM -CL -2)

author: Sable


text:
function[varargout]=printall(varargin)
%PRINTALL  Print and close all open figures.
%
%   PRINTALL prints all open figures to files in the current directory,
%   and closes the figures.
%
%   Files are named with the time stamp YYYY-MM-DD+HH|MM|SS followed
%   by their window number.
%
%   For image file output (png, jpeg, or tiff), the files are also cropped
%   using CROP by A. Bliss.
%
%   PRINTALL(NAME) prepends the string NAME to the file names.
%
%   PRINTALL(NAME,STR) passes the string STR as the printer device, for
%   example, PRINTALL('-djpg -r200').  The default choice is STR='-dpng'.
%
%   PRINTALL(H,...) only prints those figures whose handles are given in H.
%
%   Usage: printall samplefigures
%          printall samplefigures -depsc
%          printall([],' -depsc')
%          printall(h,'-djpeg')
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2020 J.M. Lilly --- type 'help jlab_license' for details

str='-dpng';
name=[];
h=[];
if ~isempty(varargin)
    if ~ischar(varargin{1})
        h=varargin{1};
        varargin=varargin(2:end);
    end
end
if ~isempty(varargin)
    name=varargin{1};
    varargin=varargin(2:end);
end
if ~isempty(varargin)
    str=varargin{1};
    varargin=varargin(2:end);
end

if ~isempty(name)
    name=[name '@'];
end

ext=[];
if contains(str,'jpeg')
    ext='.jpg';
elseif contains(str,'png')
    ext='.png';
elseif contains(str,'tiff')
    ext='.tif';
end

if isempty(h)
    a=findobj('Type','Figure');
    for i=length(a):-1:1
        figure(a(i).Number)
        filename= datestr(now,31);filename(11)='+';filename(14:3:end)='|';
        namestr=[name filename '-' int2str(a(i).Number)];
        eval(['print ' str ' ' namestr])
        if ~isempty(ext)
            eval(['crop ' namestr ext])
        end
        close(a(i).Number);
    end
else
    for i=1:length(h)
        figure(h(i))
        filename= datestr(now,31);filename(11)='+';filename(14:3:end)='|';
        namestr=[name filename '-' int2str(h(i))];
        eval(['print ' str ' ' namestr])
        if ~isempty(ext)
            eval(['crop ' namestr ext])
        end
        close(h(i));
    end
end
%close all

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%%*****************************************************************
%% min sum_k bk*yk
%% s.t. sum yk*Hk  <= 0
%%      y1 = 1
%% Hk = -hankel(ek) if   1 <= k <= n
%%    = -hankel(0,e(k-n+1)) if n+1 <= k <=2*n-1
%%
%% [blk,At,C,b] = sdphankel(n);
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [blk,At,C,b] = sdphankel(n);

   randn('seed',0);
   tmp = randn(n,n);
   tmp = tmp+tmp';
   X{1} = tmp + norm(tmp,'fro')*speye(n,n);
%%
   for k = 1:n
      ek = zeros(n,1); ek(k) = -1;
      AA{k} = sparse(hankel(ek));
   end
   zz = zeros(n,1);
   for k = n+1:2*n-1
      ek = zeros(n,1); ek(k-n+1) = -1;
      AA{k} = sparse(hankel(zz,ek));
   end
   blk{1,1} = 's'; blk{1,2} = n;
   At = svec(blk,AA,1);
   C{1} = spconvert([n n 0]);
   b = AXfun(blk,At,[],X);
%%
   blk{2,1} = 'u'; blk{2,2} = 1;
   ee = zeros(1,2*n-1); ee(1) = 1;
   At{2,1} = ee;
   C{2,1} = 1;
%%***************************************************