text:
function [xUpdate,SUpdate,innov,Szz,W]=sqrtEKFUpdateWithPred(z,SR,zPred,otherInfo,innovTrans,stateTrans)
%%SQRTEKFUPDATEWITHPRED Given the output of the measurement prediction step
%           from sqrtEKFMeasPred and a measurement, complete the measurement
%           update step of the first-order square root extended Kalman
%           filter. Separating the measurement prediction step from the
%           rest of the update step can make the creation of multiple
%           measurement association hypotheses from a single target
%           prediction more efficient. The full measurement update function
%           is sqrtEKFUpdate.
%
%INPUTS: z The zDimX1 measurement vector.
%       SR The zDimXzDim lower-triangular square root of the measurement
%          covariance matrix in the native coordinate system of the
%          measurement.
%    zPred The zDimXnumComp measurement prediction from the filter.
% otherInfo The intermediate results returned in the otherInfo output of
%          the sqrtEKFMeasPred function.
% innovTrans An optional function handle that computes and optionally
%          transforms the value of the difference between the observation
%          and any predicted points. This is called as innovTrans(a,b) and
%          the default if omitted or an empty matrix is passed is
%          @(a,b)bsxfun(@minus,a,b). This must be able to handle sets of
%          values. For a zDimX1 measurement, either of the inputs could be
%          zDimXN in size while one of the inputs could be zDimX1 in size.
%          This only needs to be supplied when a measurement difference
%          must be restricted to a certain range. For example, the
%          innovation between two angles will be 2*pi if one angle is zero
%          and the other 2*pi, even though they are the same direction. In
%          such an instance, a function handle to the
%          wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%          appropriate parameters should be passed for innovTrans.
% stateTrans An optional function that takes a state estimate and
%          transforms it. This is useful if one wishes the elements of the
%          state to be bound to a certain domain. For example, if an
%          element of the state is an angle, one should generally want to
%          bind it to the region +/-pi.
%
%OUTPUTS: xUpdate The xDimXnumComp updated state vectors.
%         SUpdate The updated xDimXxDimXnumComp lower-triangular square-
%                 root state covariance matrices.
%      innov, Szz The zDimXnumComp innovations and the zDimXzDimXnumComp
%                 square-root innovation covariance matrices are returned
%                 in case one wishes to analyze the consistency of the
%                 estimator or use those values in gating or likelihood
%                 evaluation.
%               W The xDimXzDimXnumComp gains used in the update.
%
%See the comments to the function sqrtEKFMeasPred for an example of usage
%of this function. See the comments to sqrtEKFUpdate for more information
%on the algorithm.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(stateTrans))
    stateTrans=@(x)x;
end

if(nargin<5||isempty(innovTrans))
    innovTrans=@(a,b)bsxfun(@minus,a,b);
end

H=otherInfo.H;
Pxz=otherInfo.Pxz;
xPred=otherInfo.xPred;
SPred=otherInfo.SPred;

xDim=size(xPred,1);
numComp=size(xPred,2);
zDim=size(z,1);

xUpdate=zeros(xDim,numComp);
SUpdate=zeros(xDim,xDim,numComp);
innov=zeros(zDim,numComp);
Szz=zeros(zDim,zDim,numComp);
W=zeros(xDim,zDim,numComp);
for k=1:numComp
    Szz(:,:,k)=tria([H(:,:,k)*SPred(:,:,k),SR]);
    W(:,:,k)=(Pxz(:,:,k)/Szz(:,:,k)')/Szz(:,:,k);

    innov(:,k)=innovTrans(z,zPred(:,k));
    xUpdate(:,k)=stateTrans(xPred(:,:,k)+W(:,:,k)*innov(:,k));

    temp=W(:,:,k)*H(:,:,k);
    SUpdate=tria([(eye(size(temp))-temp)*SPred(:,:,k),W(:,:,k)*SR]);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%  Figure 5.45  Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami

clf
numG=[1 0];
denG=[1 1 4];
K1=0:.1:8;   % K = 2 at breakin point
K2=100;
K=[K1 K2];
axis('square')
r=rlocus(numG,denG,K);    % roots vs. K
plot(r,'-'),nicegrid;
axis([-6 2 -4 4])
hold on
p=roots(denG);             % find and plot OL poles
z=roots(numG);             % find and plot OL zeros
plot(z(1),0,'o')
plot(real(p(1)),imag(p(1)),'x',real(p(2)),imag(p(2)),'x')
title('Fig. 5.45  Root Locus vs. K_T')
xlabel('Re(s)')
ylabel('Im(s)')
hold off
axis('normal')

author: Sable

repo: mcbench-benchmarks


text:
function jac = p16_jac ( option, nvar, x )

%*****************************************************************************80
%
%% P16_JAC computes the jacobian for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, X(NVAR), the point where the jacobian is evaluated.
%
%    Output, real JAC(NVAR-1,NVAR), the jacobian matrix evaluated
%    at X.  The NVAR-th row is not set by this routine.
%
  jac = zeros ( nvar, nvar );

  beta = 0.4;
  gamma = 20.0;

  lambda = x(nvar);
  m = nvar - 1;
  h = 1.0 / ( m - 1  );

  for i = 1 : m

    s = h * ( i - 1 );
    trapezoid = 0.0;

    for j = 1 : m

      t = h * ( j - 1 );

      arg = beta * gamma * ( 1.0 - x(j) ) / ( 1.0 + beta * ( 1.0 - x(j) ) );

      if ( j == 1 )
        factor = 0.5;
      elseif ( j < m - 1 )
        factor = 1.0;
      elseif ( j == m )
        factor = 0.5;
      end

      trapezoid = trapezoid + h * factor * x(j) * exp ( arg ) * ...
        ( max ( s, t ) - 1.0 );

      dg = - beta * gamma / ( 1.0 + beta * ( 1.0 - x(j) ) )^2;

      jac(i,j) = jac(i,j) - lambda * h * factor * exp ( arg ) * ...
        ( 1.0 + x(j) * dg ) * ( max ( s, t ) - 1.0 );

    end

    jac(i,i) = jac(i,i) + 1.0;

    jac(i,nvar) = - trapezoid;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestCalcBlurriness
    %TestCalcBlurriness

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
            f = cv.calcBlurriness(img);
            validateattributes(f, {'numeric'}, {'scalar', 'real'});
        end

        function test_error_argnum
            try
                cv.calcBlurriness();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function i4vec_aminz_index_test ( )

%*****************************************************************************80
%
%% I4VEC_AMINZ_INDEX_TEST tests I4VEC_AMINZ_INDEX;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_AMINZ_INDEX_TEST\n' );
  fprintf ( 1, '  For an integer vector:\n' );
  fprintf ( 1, '  I4VEC_AMINZ_INDEX: index of minimum nonzero absolute entry;\n' );

  seed = 123456789;

  b = -n;
  c = n;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Input vector:' );

  fprintf ( 1, '\n' );

  ival = i4vec_aminz_index ( n, a );

  fprintf ( 1, '  Minimum abs nonzero index: %d\n', ival );

  return
end


text:
function [sFile, ChannelMat] = in_fopen_ant(DataFile)
% IN_FOPEN_ANT: Open an ANT EEProbe .cnt/.avr file (continuous recordings).
%
% USAGE:  [sFile, ChannelMat] = in_fopen_ant(DataFile)

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

%% ===== READ HEADER =====
% Read a small block of data, to get all the extra information
hdr = eepv4_read_info(DataFile);
% Copy some fields for backward compatibility with previous versions of the library
hdr.nsample = hdr.sample_count;

% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder  = 'l';
sFile.filename   = DataFile;
sFile.format     = 'EEG-ANT-CNT';
sFile.prop.sfreq = double(hdr.sample_rate);
sFile.device     = 'ANT';
sFile.header     = hdr;
% Comment: short filename
[fPath, fBase, fExt] = bst_fileparts(DataFile);
sFile.comment = fBase;
% Time and samples indices
sFile.prop.times = [0, hdr.sample_count - 1] ./ sFile.prop.sfreq;
sFile.prop.nAvg  = 1;
% Get bad channels
sFile.channelflag = ones(hdr.channel_count, 1);

%% ===== EVENTS =====
if isfield(hdr, 'triggers') && ~isempty(hdr.triggers)
    % Get list of events
    allNames = {hdr.triggers.label};
    [uniqueEvt, iUnique] = unique(allNames);
    uniqueEvt = allNames(sort(iUnique));
    % Initialize list of events
    events = repmat(db_template('event'), 1, length(uniqueEvt));
    % Format list
    for iEvt = 1:length(uniqueEvt)
        % Find list of occurences of this event
        iOcc = find(strcmpi(allNames, uniqueEvt{iEvt}));
        % Fill event structures
        events(iEvt).label      = uniqueEvt{iEvt};
        events(iEvt).color      = [];
        events(iEvt).times      = round([hdr.triggers(iOcc).seconds_in_file] .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
        events(iEvt).epochs     = ones(1, length(events(iEvt).times));  % Epoch: set as 1 for all the occurrences
        events(iEvt).reactTimes = [];
        events(iEvt).select     = 1;
        events(iEvt).channels   = [];
        events(iEvt).notes      = [];
    end
    % Import this list
    sFile = import_events(sFile, [], events);

%% ===== EXTERNAL EVENT FILE =====
else
    % If a .trg file exists with the same name: load it
    [fPath, fBase, fExt] = bst_fileparts(DataFile);
    TrgFile = bst_fullfile(fPath, [fBase '.trg']);
    % If file exists
    if file_exist(TrgFile)
        sFile = import_events(sFile, [], TrgFile, 'ANT');
    end
end

%% ===== CREATE DEFAULT CHANNEL FILE =====
% Create channel structure
Channel = repmat(db_template('channeldesc'), [1 hdr.channel_count]);
for i = 1:hdr.channel_count
    Channel(i).Name    = hdr.channels(i).label;
    Channel(i).Type    = 'EEG';
    Channel(i).Orient  = [];
    Channel(i).Weight  = 1;
    Channel(i).Comment = [];
    Channel(i).Loc = [0; 0; 0];
end
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'ANT standard position';
ChannelMat.Channel = Channel;

author: brainstorm-tools


text:
domain = [-1,1,-1,1,-1,1];
% bm = 1; bp = 100;
% rx = 1; ry = 0.075; rz = 3;
% pde = elli3DorthocircIntf(bm,bp,rx,ry,rz);
%intf = @(x,y,z) x.^2+y.^2+z.^2-1;epsm = 8.85*10^(-3)*2;
epsp = 8.85*10^(-3);
sigm = 1;
sigp = 0.01;
mum = (4*pi)*3;
mup = (4*pi);
x0=0; y0=0; z0=-0.3; r1=0.2; r2=pi/5; omega = 1; a = omega*sqrt(epsp*mup); b= 150; intPt = -1;
pde = TorusTimeInitial3(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2,a,b,intPt);

fimplicit3(pde.intf,domain)
material metal

colormap(jet)

%%%%%%%%%%%%%%%%%%%%%%%
nx = 10;  h=(domain(2) - domain(1))/nx;
ny = nx;
nz = nx;

mesh = genMesh3D(domain, nx, ny, nz);
mesh = enrichMesh3D(mesh,2); % Mesh detail level = 1 (for IFE).
mesh = genIntfMesh3D(mesh,pde.intf);
meshI = genIVmesh(mesh);

mdpt = (mesh.p(mesh.t(:,1),:) + mesh.p(mesh.t(:,2),:) +...
    mesh.p(mesh.t(:,3),:) + mesh.p(mesh.t(:,4),:))/4;
tidshow = find(mdpt(:,2)>0);
tetramesh(mesh.t(tidshow,:),mesh.p,'FaceAlpha',0)
hold on
trisurf(meshI.iface,meshI.node(:,1),meshI.node(:,2),meshI.node(:,3)) %'edgecolor','none'
box on
%axis off
axis equal
set(gca,'XTick',[],'YTick',[],'ZTick',[])

view(50,30)

repo: ifem


text:
function i = binary_to_i4 ( s )

%*****************************************************************************80
%
%% BINARY_TO_I4 converts a binary representation into an integer value.
%
%  Example:
%
%        S        I
%
%      '101'      5
%    '-1000'     -8
%        '1'      1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the binary representation.
%
%    Output, integer I, the integer whose representation was input.
%
  nchar = s_len_trim ( s );

  i = 0;
  ichr = 1;
  istate = 0;
  isgn = 1;

  while ( ichr <= nchar )

    c = s(ichr);
%
%  Blank.
%
    if ( c == ' ' )

      if ( istate == 2 )
        istate = 3;
      end
%
%  Sign, + or -.
%
    elseif ( c == '-' )

      if ( istate == 0 )
        istate = 1;
        isgn = - 1;
      else
        istate = - 1;
      end

    elseif ( c == '+' )

      if ( istate == 0 )
        istate = 1;
      else
        istate = - 1;
      end
%
%  Digit, 0 or 1.
%
    elseif ( c == '1' )

      i = 2 * i;
      i = i + 1;
      istate = 2;

    elseif ( c == '0' )

      i = 2 * i;
      istate = 2;
%
%  Illegal or unknown sign.
%
    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'BINARY_TO_I4 - Serious error!\n' );
      fprintf ( 1, '  Illegal digit = "%c"\n', c );
      fprintf ( 1, '  Conversion halted prematurely!\n' );
      error ( 'BINARY_TO_I4 - Serious error!' );

    end

    if ( istate == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'BINARY_TO_I4 - Serious error!\n' );
      fprintf ( 1, '  Unable to decipher input!\n' );
      error ( 'BINARY_TO_I4 - Serious error!' );
    end

    if ( 3 <= istate )
      break;
    end

    ichr = ichr + 1;

  end
%
%  Apply the sign.
%
  i = isgn * i;

  return
end


text:
function [anybad] = check_spm_matfiles(P)
% Check a list of SPM-style images to see if they all have the same
% space, dims.
%
% :Usage:
% ::
%
%     [anybad] = check_spm_matfiles(P)
%

V=spm_vol(P);
for i=2:length(V),
    chk = V(1).mat - V(i).mat; chk = chk(:);
    chk = chk(1:end-1);             % eliminate SPM scale factor
    chk = any(chk);
    notok(i) = chk;
end

if any(notok)
    wh = find(notok);

    disp(['The following images'' mat files differend from the first:'])
    disp(num2str(wh));

    disp(['First mat:']);
    disp(V(1).mat);

    disp(['Bad mats:']);
    for i = wh
        disp(V(i).fname);
        disp(V(i).mat);
    end
end

anybad = any(notok);
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function M = rmCompareModelsGUI_selectVoxels(M, whichVoxels, plotFlag);
% Sub-select voxels in the ROI of rmCompareModelsGUI.
%
% M = rmCompareModelsGUI_selectVoxels([M=get from cur figure], [whichVoxels=dialog], [plotFlag=2]);
%
%
% INPUTS:
%	M: rmCompareModelsGUI structure. See rmCompareModelsGUI_getData. Can also be the handle
%	to an rmCompareModelsGUI figure.
%
%	whichVoxels: index of voxels within M.coords to select.
%	[default: open a dialog allowing the user to select voxels based on
%	basic pRF parameter ranges (e.g., varexp, eccentricity, polar angle,
%	sigma)].
%
%	plotFlag: flag indicating where to show the new analysis on the
%	restricted set of voxels. Values are:
%		0 -- just return the M structure and don't make/update a GUI.
%		1 -- replace the existing GUI with the new data (throw out the
%			 eliminated voxels).
%		2 -- open a new GUI for the restricted subset.
%		[default 2, open a new figure]
%
% Returns a modified version of M.
%
% ras, 04/2009.
if notDefined('M'),			M = get(gcf, 'UserData');		end
if ishandle(M),				M = get(M, 'UserData');			end
if notDefined('plotFlag'),	plotFlag = 2;					end
if notDefined('whichVoxels'),
	whichVoxels = selectVoxelsDialog(M);
end

% is M the proper structure?
if ~checkfields(M, 'modelList') | ~checkfields(M, 'tSeries')
	error('Need an rmCompareModelsGUI structure. See rmCompareModelsGUI.')
end

% are the coordinates empty? If so, warn and do nothing, but don't error.
if isempty(whichVoxels)
	warning('No voxels selected.')
	return
end

% sub-select the voxels
M.roi.coords = M.roi.coords(:,whichVoxels);
for m = 1:M.nModels
	for f = {'tSeries' 'x0' 'y0' 'sigma' 'pol' 'ecc' 'varexp'}
		M.(f{1}){m} = M.(f{1}){m}(:,whichVoxels);
	end
	M.beta{m} = M.beta{m}(whichVoxels,:);
end
M.roi.name = ['Subset of ' M.roi.name];
M.nVoxels = length(whichVoxels);

% update the GUI as desired
switch plotFlag
	case 0, % do nothing

	case 1, % replace data in existing GUI
		set(M.fig, 'UserData', M);

	case 2,
		M = rmCompareModelsGUI_openFig(M);

	otherwise, error('Invalid plot flag.')
end

return
% /-----------------------------------------------------------/ %

% /-----------------------------------------------------------/ %
function whichVoxels = selectVoxelsDialog(M);
%% put up a dialog to get pRF criteria for the voxels to keep, then find
%% the voxels which satisfy these criteria.

% build the dialog
dlg(1).fieldName = 'modelNum';
dlg(end).style = 'popup';
dlg(end).value = 1;
dlg(end).list  = M.dtList;
dlg(end).string = 'Select based on which model''s parameters?';

dlg(end+1).fieldName = 'notes';
dlg(end).style = 'text';
dlg(end).value = '';
dlg(end).string = 'Leave a field empty to omit thresholding:';

dlg(end+1).fieldName = 'varexp';
dlg(end).style = 'number';
dlg(end).value = 0.10;
dlg(end).string = 'Variance explained threshold:';

dlg(end+1).fieldName = 'ecc';
dlg(end).style = 'number';
dlg(end).value = [0 30];
dlg(end).string = 'Eccentricity range (blank for all values):';

dlg(end+1).fieldName = 'pol';
dlg(end).style = 'number';
dlg(end).value = [0 360];
dlg(end).string = 'Polar Angle range (deg CW from up):';

dlg(end+1).fieldName = 'sigma';
dlg(end).style = 'number';
dlg(end).value = [0 30];
dlg(end).string = 'pRF Diameter (major axis) range:';

% put up the dialog
resp = generalDialog(dlg, 'rmCompareModelsGUI: Select Voxels');

% get the model number
m = cellfind(M.dtList, resp.modelNum);
m = m(1); % in case 2 models from the same data type

% find voxels which pass the crtieria
ok = ones(1, size(M.roi.coords, 2));

if ~isempty(resp.varexp)
	ok = ok & (M.varexp{m} >= resp.varexp);
end

if length(resp.pol) > 1
	ok = ok & (M.pol{m} >= resp.pol(1) & M.pol{m} <= resp.pol(2));
end

if length(resp.ecc) > 1
	ok = ok & (M.ecc{m} >= resp.ecc(1) & M.ecc{m} <= resp.ecc(2));
end

if length(resp.sigma) > 1
	ok = ok & (M.sigma{m} >= resp.sigma(1) & M.sigma{m} <= resp.sigma(2));
end

whichVoxels = find(ok);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Q = cost_allstim(V,t,tc,Run)
%
% Least-squares cost function for the IL model
%
% INPUT:
% Run = stick function
% tc = time course
% t = vector of time points
% V = parameters
%
% OUTPUT:
% Q = cost
%
% By Martin Lindquist and Tor Wager
% Edited 12/12/06
% Further edited by Christian Waugh 2/15/08 to include multiple trialtypes

numstim = length(Run);
len = length(Run{1});
h = zeros(length(t),numstim);
yhatt =zeros(len,numstim);

for k = 1:numstim
    h(:,k) = Get_Logit(V(k*7-6:k*7),t);            % Get IL model corresponding to parameters V
    yhat(:,k) = conv(Run{k}, h(:,k));              % Convolve IL model with stick function
    yhatt(:,k) = yhat(1:len,k);
end

yhat2 = sum(yhatt,2); %Sum models together to get overall estimate

Q = sum((yhat2-tc).^2);              % Calculate cost function

return


text:
function out = BER(hidden, retrieved)
%BER Bit Error Rate to measure error rate

y = getBits(hidden); x = getBits(retrieved);
len = min(length(x), length(y));

ber = 0;
for i=1:len
    err = (x(i)~= y(i));
    ber = ber + err;
end
out = 100*(ber/len);
end

author: ktekeli

repo: audio-steganography-algorithms


text:
function [y1]=rs(y1,T,f2)
% [y2] = rs(y1, T) resamples y1 to the target sampling rate y2 using T
% [y2] = rs(y1, f1, f2) resamples y1 with f1 to the target sampling rate f2
%
% RS does not require overlap data.
%
% see also: SOPEN, SREAD, SCLOSE, MAT2SEL, SAVE2TXT, SAVE2BKR
%
% Reference(s):

%	$Revision: 1.2 $
%	$Id: rs.m 2202 2009-10-27 12:06:45Z schloegl $
%	Copyright (C) 1997-2004 by Alois Schloegl
%	a.schloegl@ieee.org
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% Version 2 of the License, or (at your option) any later version.
%
% This library 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
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.

if nargin==3,
        f1=T;
        if f1==f2
                return;
        elseif f1>f2
                D=f1/f2;
                [yr,yc]=size(y1);
                LEN=yr/D;
		y2=zeros(yr*f2/f1,yc);
                for k=0:LEN-1
                        y2(k+1,:)=sum(y1(k*D+(1:D),:),1)/D;
                end;
		y1=y2;
        else %f1<f2
		y1=y1(ceil((1:size(y1,1)*f2/f1)/f2*f1),:);
        end;

elseif nargin==2,
        [f1,f2]=size(T);
        if f1==f2,
                return;
        end;
        [yr,yc]=size(y1);
        LEN=yr/f1;
	y2=zeros(yr*f2/f1,yc);
        for k=0:LEN-1
                y2(k*f2+(1:f2),:)=T'*y1(k*f1+(1:f1),:);
        end;
	y1=y2;
end;

author: PatternRecognition

repo: OpenBMI


text:
function glottal_pulse_GUI25
% Modifiable runGUI file
clc;
clear all;
fileName = 'glottal_pulse.mat';    %USER - ENTER FILENAME
fileData=load(fileName);
temp=fileData(1).temp;

f = figure('Visible','on',...
'Units','normalized',...
'Position',[0,0,1,1],...
'MenuBar','none',...
'NumberTitle','off');

Callbacks_glottal_pulse_GUI25(f,temp);    %USER - ENTER PROPER CALLBACK FILE
%panelAndButtonEdit(f, temp);       % Easy access to Edit Mode

% Note comment PanelandBUttonCallbacks(f,temp) if panelAndButtonEdit is to
% be uncommented and used
end

% glottal_pulse_rev1_gui25 design
% 2 Panels
%   #1 - parameters
%   #2 - graphics
% 2 Graphic Panels
%   #1 - glottal pulse impulse response
%   #2 - glottal pulse frequency response
% 1 TitleBox
% 8 Buttons
%   #1 - editable button - alpha1: rising glottal pulse cycle
%   #2 - editable button - alpha2: falling glottal pulse cycle
%   #3 - editable button - period for periodic excitation
%   #4 - popupmenu - list of vowel sounds for synthesis
%   #5 - pushbutton - Play impulse train
%   #6 - pushbutton - Play glottal pulses
%   #7 - pushbutton - Play vowel sound
%   #8 - pushbutton - Close GUI


text:
function [max_err, mean_err] = tperror(lpv, S, U, domain, n)
%TPERROR Calculate the error of the discretized TP model at random points
%	[max_err, mean_err] = TPERROR(lpv, S, U, domain, n)
%
%	lpv      - LPV model
%	S        - core tensor of the TP model
%	U        - weight function data of the TP model
%	domain   - parameter domain of the LPV model
%	n        - number of random test points
%
%	max_err  - maximum L2 error over the n test points
%	mean_err - mean L2 error over the n test points

% TODO: use queryw instead of queryw1

% n random points in the given parameter space
P = size(domain,1);
x = zeros(n, P);
for i = 1:P
	x(:,i) = rand(n,1) * (domain(i,2) - domain(i,1)) + domain(i,1);
end

err = zeros(n, 1);
for i = 1:n
	% S from original model
	S_lpv = querylpv(lpv, x(i,:));

	% S from tp model
	W = queryw1(U, domain, x(i,:));
	S_tp = squeeze(tprod(S, W));

	% L2 error
	err(i) = norm(S_lpv - S_tp, 2);
	if err(i) > 100
		disp('Huge error')
		disp(x(i,:))
		S_lpv
		S_tp
	end
end

max_err = max(err);
mean_err = mean(err);


text:
function min_f = minf(datafe,dataor)
T=length(dataor);
sum=0;
for i=1:T
   sum=sum+ abs((datafe(i)-dataor(i))*100/dataor(i));
end

min_f=sum;

author: sfvsfv

repo: Mathematical-modeling


text:
function [n]=nnz(a)
% Returns TT-ranks of the given tensor in TT-format
%   [N]=NNZ(A) Returns number of non zero elements in the TT representation of A
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

n=nnz(a.core);
return
end

author: oseledets

repo: TT-Toolbox


text:
function varargout = fillPolygon3d(varargin)
%FILLPOLYGON3D Fill a 3D polygon specified by a list of vertex coords.
%
%   fillPolygon3d(COORD, COLOR)
%   packs coordinates in a single [N*3] array.
%   COORD can also be a cell array of polygon, in this case each polygon is
%   drawn using the same color.
%
%   fillPolygon3d(PX, PY, PZ, COLOR)
%   specifies coordinates in separate numeric vectors (either row or
%   columns)
%
%   fillPolygon3d(..., PARAM, VALUE)
%   allows to specify some drawing parameter/value pairs as for the plot
%   function.
%
%   H = fillPolygon3d(...)
%   also returns a handle to the list of created patch objects.
%
%   Example
%     t = linspace(0, 2*pi, 100)';
%     xt = 10 * cos(t);
%     yt = 5 * sin(t);
%     zt = zeros(1,100);
%     figure; fillPolygon3d(xt, yt, zt, 'c');
%
%   See Also:
%   polygons3d, drawPolygon3d, drawPolyline3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-01-05
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% check case we want to draw several curves, stored in a cell array
var1 = varargin{1};
if iscell(var1)
    hold on;
    h = [];
    for i = 1:length(var1(:))
        h = [h; fillPolygon3d(var1{i}, varargin{2:end})]; %#ok<AGROW>
    end
    if nargout>0
        varargout{1}=h;
    end
    return;
end

% extract vertex coordinates
if min(size(var1)) == 1
    % if first argument is a vector (either row or column), then assumes
    % first argument contains x coords, second argument contains y coords
    % and third one the z coords
    px = var1;
    if length(varargin) < 3
        error('geom3d:fillPolygon3d:Wrong number of arguments in fillPolygon3d');
    end
    py = varargin{2};
    pz = varargin{3};
    varargin = varargin(4:end);
else
    % first argument contains all three coordinates
    px = var1(:, 1);
    py = var1(:, 2);
    pz = var1(:, 3);
    varargin = varargin(2:end);
end

% extract color information
if isempty(varargin)
    color = 'c';
else
    color = varargin{1};
    varargin = varargin(2:end);
end

% fill the polygon
h = fill3(px, py, pz, color, varargin{:});

if nargout>0
    varargout{1}=h;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Function to derive short term power contour as is done in Ishi et al 2008
%
% Description
% Function to derive short term power contour as is done in Ishi et al 2008
%
%
% Inputs
%  x        : [samples] [Nx1]  Speech signal
%  fs       : [Hz]      [1x1]  Sampling frequency
%
% Outputs
%  pow      : [dB]      [Mx1]  Power contour
%  pow_std  : [samples] [Px1]  Standard deviation of power contour
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Ishi, C., Sakakibara, K-I, Ishiguro, H., (2008) `A method for
%       automatic detection of vocal fry', IEEE TASLP, 16(1), 47-56.
%  [2] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [3] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [4] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  This program is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function [pow,pow_std,pow_std_inter] = get_short_pow(x,fs)

% Get very short term power contour
veryShort_len = 4*(fs/1000); % 4ms frame length for "very short-term" analysis
veryShort_shift = 2*(fs/1000); % 2ms shift for "very short-term" analysis
veryShort_powCont = zeros(1,ceil((length(x)-veryShort_len)/veryShort_shift));
start=1;
finish=start+veryShort_len-1;

n=1;
x2 = x.^2;
while finish <= length(x)
    veryShort_powCont(n) = mean(x2(start:finish));
    start = start + veryShort_shift;
    finish=start+veryShort_len-1;
    n=n+1;
end
clear x2;

pow = 20*log10(veryShort_powCont);
inf_idx=isinf(pow);
pow(inf_idx)=min(pow(~inf_idx));

pow_std=zeros(1,length(pow));
std_len=16;

for n=std_len+1:length(pow)-std_len

    pow_std(n) = std(pow(n-std_len:n+std_len));
end
pow_std=medfilt1(pow_std,13);

pow_std_inter=interp1(linspace(1,length(x),length(pow)),pow_std,1:length(x));

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%ISAFFINE Test affine mapping
%
%    I = ISAFFINE(W)
%    ISAFFINE(W)
%
% I is true if W is an affine mapping.
% If called without an output argument ISAFFINE generates an error
% if W is not an affine mapping.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function hout=suptitle(str, fs)
%SUPTITLE Puts a title above all subplots.
%	SUPTITLE('text') adds text to the top of the figure
%	above all subplots (a "super title"). Use this function
%	after all subplot commands.

% Drea Thomas 6/15/95 drea@mathworks.com

% Warning: If the figure or axis units are non-default, this
% will break.

% Parameters used to position the supertitle.

% Amount of the figure window devoted to subplots
plotregion = .92;

% Y position of title in normalized coordinates
titleypos  = .95;

% Fontsize for supertitle
if nargin < 2
  fs = get(gcf,'defaultaxesfontsize')+4;
end

% Fudge factor to adjust y spacing between subplots
fudge=1;

haold = gca;
figunits = get(gcf,'units');

% Get the (approximate) difference between full height (plot + title
% + xlabel) and bounding rectangle.

	if (~strcmp(figunits,'pixels')),
		set(gcf,'units','pixels');
		pos = get(gcf,'position');
		set(gcf,'units',figunits);
	else,
		pos = get(gcf,'position');
	end
	ff = (fs-4)*1.27*5/pos(4)*fudge;

        % The 5 here reflects about 3 characters of height below
        % an axis and 2 above. 1.27 is pixels per point.

% Determine the bounding rectange for all the plots

% h = findobj('Type','axes');

% findobj is a 4.2 thing.. if you don't have 4.2 comment out
% the next line and uncomment the following block.

h = findobj(gcf,'Type','axes');  % Change suggested by Stacy J. Hills

% If you don't have 4.2, use this code instead
%ch = get(gcf,'children');
%h=[];
%for i=1:length(ch),
%  if strcmp(get(ch(i),'type'),'axes'),
%    h=[h,ch(i)];
%  end
%end

max_y=0;
min_y=1;

oldtitle =0;
for i=1:length(h),
	if (~strcmp(get(h(i),'Tag'),'suptitle')),
		pos=get(h(i),'pos');
		if (pos(2) < min_y), min_y=pos(2)-ff/5*3;end;
		if (pos(4)+pos(2) > max_y), max_y=pos(4)+pos(2)+ff/5*2;end;
	else,
		oldtitle = h(i);
	end
end

if max_y > plotregion,
	scale = (plotregion-min_y)/(max_y-min_y);
	for i=1:length(h),
		pos = get(h(i),'position');
		pos(2) = (pos(2)-min_y)*scale+min_y;
		pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
		set(h(i),'position',pos);
	end
end

np = get(gcf,'nextplot');
set(gcf,'nextplot','add');
if (oldtitle),
	delete(oldtitle);
end
ha=axes('pos',[0 1 1 1],'visible','off','Tag','suptitle');
ht=text(.5,titleypos-1,str);set(ht,'horizontalalignment','center','fontsize',fs);
set(gcf,'nextplot',np);
axes(haold);
if nargout,
	hout=ht;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = subsref( F, ref )
%SUBSREF   DISKFUNV subsref.
%   F(X,Y) or F(X, Y, 'cart') returns values of the DISKFUNV F evaluated on
%   the array (X, Y) in Cartesian coordinates.
%
%   F(T, R, 'polar') returns the values of the DISKFUNV F evaluated on the
%   array (T,R) in polar coordinates.
%
%   F(k) returns the kth component of F.
%
%   F.PROP returns the property PROP of F as defined by GET(F, 'PROP').
%
% See also DISKFUN/GET, DISKFUNV/SET

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

% Empty check:
if ( isempty( F ) )
   varargout = { [] };
   return
end

% Get indexing reference:
indx = ref(1).subs;

switch ( ref(1).type )

    case '.'
        if ( numel( ref ) == 1 )

            % This is a get call to get a property.
            varargout = { get(F, indx) };

        else

            % Probably .^ or maybe .*
            %t2 = index(2).type;
            t2 = ref(2).type;
            if ( strcmp(t2, '.') )
                out = get(F, indx, ref(2).subs{:});
            else
                out = get(F, indx);
                out = out( ref(2).subs{:} );
            end
            if ( numel(ref) > 2 )
                varargout = {subsref(out, ref(3:end))};
            else
                varargout = { out };
            end

        end

    case '()'
        % Either evaluation or the user wants the kth component of F:
        if ( length(indx) > 1 )
            % EVALUATION, call FEVAL:
            x = indx{1};
            y = indx{2}; % where to evaluate
            vals = feval(F, x, y, ref(1).subs{:});
            varargout = { vals };

        else
            % COMPONENTS, extract the component:
            if ( isa(indx{1},'double') )
                if all( indx{1} == 1  )
                    % F(1):
                    varargout = F.components(1);
                elseif ( all( indx{1} == 2 ) )
                    % F(2):
                    varargout = F.components(2);
                else
                    error('CHEBFUN:DISKFUNV:subsref:index', ...
                        'DISKFUNV only contains two components');
                end
            end
        end

    otherwise
        error('CHEBFUN:DISKFUNV:subsref:unexpectedType', ...
            ['??? Unexpected index.type of ' index(1).type]);

end

% Recurse down the subsref tree, if appropriate:
if ( numel( ref ) > 1 )
   ref(1) = [];
   varargout = { subsref( varargout{ : }, ref ) };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Z = variable_replace(X,Y,W)

% Check so that Y is a simple unit variable
Ybase = getbase(Y);
Yvariables = getvariablesSORTED(Y);
Xbase = getbase(X);
Xvariables = getvariables(X);
[i,j,k] = find(Ybase);
if ~isequal(sort(i),1:length(i))
end
if ~isequal(sort(j),2:(length(i)+1))
end
if ~all(k == 1)
end

[mt,variabletype] = yalmip('monomtable');
% Linear, or at least linear in Y
if all(variabletype(Xvariables) == 0) %| all(sum(mt(any(mt(getvariables(X),getvariables(Y)),2),:),2)==1)
    % Simple linear replacement
    v = 1;
    v1 = [];
    v2 = [];
    i1 = [];
    i2 = [];
    for i = 1:length(Xvariables)
        XisinY = find(Xvariables(i) == Yvariables);
        if ~isempty(XisinY)
         %   v = [v;W(XisinY)];
            v1 = [v1 XisinY];
            i1 = [i1 i];
        else
         %   v = [v;recover(Xvariables(i))];
            v2 = [v2 Xvariables(i)];
            i2 = [i2 i];
        end
    end
    v = sparse(i1,ones(length(i1),1),W(v1),length(Xvariables),1);
    v = v + sparse(i2,ones(length(i2),1),recover(v2),length(Xvariables),1);

    Z = Xbase*[1;v];
    %Z = Xbase*[v];
    Z = reshape(Z,size(X,1),size(X,2));
else
    if nnz(mt(getvariables(X),getvariables(Y)))==0
        Z = X;
    else
        Z = nonlinearreplace(X,Y,W);
    end
    return
%    error('Nonlinear replacement not supported yet')
end

% This has not been tested (copied from variable_replace) so it is placed
% in a catch to be safe.
try
    Xvariables = getvariables(Z);
    extvar = yalmip('extvariables');
    Xext = find(ismember(Xvariables,extvar));
    if ~isempty(Xext)
        %We must dig down in extended operators to see if they use the replaced
        %set of variables
        for i = 1:length(Xext)
            extstruct = yalmip('extstruct',Xvariables(Xext(i)));
            anychanged = 0;
            for j = 1:length(extstruct.arg)
                if isa(extstruct.arg{j},'sdpvar')
                    XinY = find(ismembc(getvariables(extstruct.arg{j}),Yvariables));
                    if ~isempty(XinY)
                        anychanged = 1;
                        extstruct.arg{j} = replace(extstruct.arg{j},Y,W);
                    else
                    end
                end
            end
            if anychanged
                Zi = yalmip('define',extstruct.fcn,extstruct.arg{:});
                Xvariables(Xext(i)) = getvariables(Zi);
            end
        end
        % And now recover this sucker
        Z = struct(Z);
        Z.lmi_variables = Xvariables;
        % Messed up order (lmi_variables should be sorted)
        if any(diff(Z.lmi_variables)<0)
            [i,j]=sort(Z.lmi_variables);
            Z.basis = [Z.basis(:,1) Z.basis(:,j+1)];
            Z.lmi_variables = Z.lmi_variables(j);
        end
        Z = sdpvar(Z.dim(1),Z.dim(2),[],Z.lmi_variables,Z.basis);
    end
catch
end

function Yvariables = getvariablesSORTED(Y);
Y = Y(:);
for i = 1:length(Y)
    Yvariables(i) = getvariables(Y(i));
end

author: shamilmamedov


text:
classdef LinkSphere < frost.Animator.Sphere
    % Sphere depicting center of mass of a rigid link for
    % frost.Animator.Display
    %
    % @author omar @date 2017-06-01
    %
    % Copyright (c) 2017, UMICH Biped Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/BSD-3-Clause

    methods
        function obj = LinkSphere(ax, model, linkFrame, name, varargin)
            obj = obj@frost.Animator.Sphere(ax, model, linkFrame, name, varargin{:});

            obj.radius = 0.03;
            obj.surface.FaceColor = 'blue';
        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function c=ref_dft(f)
%REF_DFT  Reference Discrete Fourier Transform
%   Usage:  c=ref_dft(f);
%
%   REF_DFT(f) computes the unitary discrete Fourier transform of f.
%
%   AUTHOR: Jordy van Velthoven

L=size(f,1);
W=size(f,2);
c= zeros(L,W);

for w=1:W
for k=0:L-1
  for l=0:L-1
    c(k+1,w) = c(k+1,w) + f(l+1,w) * exp(-2*pi*i*k*l/L);
  end;
end;
end

c = c./sqrt(L);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


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

    enumeration
        AngleBetweenVectors('Angle Between Vectors');
        AngleBetweenVectorPlane('Angle Between Vector and Plane');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = GeometricAngleEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('GeometricAngleEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('GeometricAngleEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

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


text:
function mtrPDBScoreHist(filename,chColor,lowThresh,dt6Filename)

% lowThresh is value between 0 and 1, i.e. 0.1 will cut off the bottom 10%
% from the histogram
% chColor can be usual 'r', 'g', 'b', ...

if ieNotDefined('dt6Filename')
  outPathName = pwd;
  [f,p] = uigetfile('*.mat', 'Select dt6 file...', outPathName);
  if(isnumeric(f)) error('User cancelled.'); end
  dt6Filename = fullfile(p,f);
end

if ieNotDefined('lowThresh')
    lowThresh = 0.1;
end

% Load pathways
dt6 = load(dt6Filename,'xformToAcPc');
disp('Loading pathways ...');
fg = mtrImportFibers(filename, dt6.xformToAcPc);

% Find length param
scoreVec = mtrGetFGScoreVec(fg);
scoreVec = sort(scoreVec,'descend');
[n,xout] = hist(scoreVec(1:end*(1-lowThresh)),linspace(-50,150,20));
%bar(xout,n,chColor);
nWithNoZeros = n;
nWithNoZeros(n==0) = 1;
plot(xout,log(nWithNoZeros).*double(n>0),chColor);
xlabel('ln(score)');
ylabel('ln(count)');


text:
% pop_extractlag() - If less than
%           4 arguments are given, a window pops up
%           to ask for the value of the additional
%           parameters.
%
% Usage:
%    >>    [OUTEEG, COM, ORDER, LAGS_MS, EVENT_TYPE] = pop_extractlag( INEEG, USE_ICA, CHANNEL_ID, TIME_WIN, OPTIONS );
%
% Inputs:
%    INEEG      - input EEG dataset
%    use_ica    - type of processing. 1 process the raw data and 0 the ICA components.
%    channel_id - channel or component id to consider.
%    time_win   - time window in ms on which lag extraction should be performed
%    options    : structure containing options with default values
%        Some important options :
%           options.use_maximum = true (resp. false) : lag pass through a maximum (resp. lag pass through a minimum)
%           options.sigma = 0.05 : controls the similarity measures between time series
%           options.alpha = 0.01 (example of value) : controls regularity of lag function
%
% Outputs:
%    OUTEEG    - output dataset
%
% See also:
%    EEGLAB

% (c) Copyright 2008 Alexandre Gramfort INRIA. All Rights Reserved.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.    See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA    02111-1307    USA

% $Id: pop_extractlag.m 4 2009-08-15 21:10:35Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-08-15 17:10:35 -0400 (Sam, 15 aoû 2009) $
% $Revision: 4 $

function [EEG, com, order, lags_ms, event_type, E_lags] = pop_extractlag( EEG, use_ica, channel_id, time_win, options )

% the command output is a hidden output that does not have to
% be described in the header

com = ''; % this initialization ensure that the function will return something
    % if the user press the cancel button

% display help if not enough arguments
% ------------------------------------
if nargin < 1
    help pop_extractlag;
    return;
end;

% Set default options

if nargin<5
    options.null = 0;
end

if ~isfield(options, 'NORMALIZE_DATA')
    options.NORMALIZE_DATA = 1;
end
NORMALIZE_DATA = options.NORMALIZE_DATA;

if ~isfield(options, 'CENTER_DATA')
    options.CENTER_DATA = 0;
end
CENTER_DATA = options.CENTER_DATA;

if ~isfield(options, 'USE_ADAPTIVE_SIGMA')
    options.USE_ADAPTIVE_SIGMA = 0;
end
USE_ADAPTIVE_SIGMA = options.USE_ADAPTIVE_SIGMA;

if ~isfield(options, 'DISPLAY_SIGNAL')
    options.DISPLAY_SIGNAL = 0;
end
DISPLAY_SIGNAL = options.DISPLAY_SIGNAL;

if ~isfield(options, 'sigma')
    sigma_txt = '0.01:0.01:0.2';
    options.sigma = eval([ '[' sigma_txt ']' ]);
end
sigma = options.sigma;

if ~isfield(options, 'distance')
    options.distance = 'corr2';
end
distance = options.distance;

if ~isfield(options, 'CONNECT_EMBEDDED_POINTS')
    options.CONNECT_EMBEDDED_POINTS = 0;
end
CONNECT_EMBEDDED_POINTS = options.CONNECT_EMBEDDED_POINTS;

if ~isfield(options, 'alpha')
    alpha_txt = '0.001,0.01,0.1,1';
    options.alpha = eval([ '[' alpha_txt ']' ]);
end
alpha = options.alpha;

if ~isfield(options, 'exponent')
    options.exponent = 1;
end
exponent = options.exponent;

if ~isfield(options, 'use_maximum')
    options.use_maximum = true;
end
use_maximum = options.use_maximum;

if ~isfield(options, 'show_pca')
    options.show_pca = false;
end
show_pca = options.show_pca;

if ~isfield(options, 'show_trial_number')
    options.show_trial_number = false;
end
show_trial_number = options.show_trial_number;

if ~isfield(options, 'verbose')
    options.verbose = true;
end
verbose = options.verbose;

if ~isfield(options, 'disp_log')
    options.disp_log = true;
end
disp_log = options.disp_log;

% pop up window
% -------------

if nargin < 4

    cb_plotcomps = [ 'pop_topoplot(EEG);' ];
    cb_erpimage = [ 'pop_erpimage(EEG,1);' ];

    listui = { { 'style' 'text'       'string' 'Use ICA components' } ...
               { 'style' 'checkbox'   'string' '' } ...
               { 'style' 'text'       'string' '(check=yes)' } ...
               { 'style' 'text'       'string' 'Channel or Component' } ...
               { 'style' 'edit'       'string' [ num2str(1) ] } ...
               { 'style' 'pushbutton' 'string' 'Plot topo.' 'callback' cb_plotcomps }  ...
               { 'style' 'text'       'string' 'Time window (ms)' } ...
               { 'style' 'edit'       'string' [ num2str(EEG.times(1)) ':' num2str(EEG.times(end)) ] } ...
               { 'style' 'pushbutton' 'string' 'Erpimage' 'callback' cb_erpimage } ...
               { 'style' 'text'       'string' 'Use minimum' } ...
               { 'style' 'checkbox'   'string' '' } ...
               { 'style' 'text'       'string' '(check=yes)' } ...
               { 'style' 'text'       'string' 'Sigma' } ...
               { 'style' 'edit'       'string' sigma_txt } ...
               { 'style' 'text'       'string' '' } ...
               { 'style' 'text'       'string' 'Alpha' } ...
               { 'style' 'edit'       'string' alpha_txt } ...
               { 'style' 'text'       'string' '' } ...
               { 'style' 'text'       'string' 'Verbose' } ...
               { 'style' 'checkbox'   'string' '' } ...
               { 'style' 'text'       'string' '(check=yes)' } ...
             };
    geom = { [1 1 0.5] [0.8 0.2 1] [1 1 1] [1 1 1] [1 1 1] [1 1 1] [1 1 1] };

    results = inputgui('geometry', geom, 'uilist', listui, 'helpcom', 'pophelp(''pop_extractlag'')', ...
        'title', 'Extract lags');

    if isempty(results), return; end;

    use_ica             = eval( [ '[' num2str(results{1}) ']' ] );
    channel_id          = eval( [ '[' results{2} ']' ] );
    time_win            = eval( [ '[' results{3} ']' ] );
    time_win = [time_win(1),time_win(end)];

    options.use_maximum = not(results{4});
    options.sigma       = eval( [ '[' results{5} ']' ] );
    options.alpha      = eval( [ '[' results{6} ']' ] );
    options.verbose     = results{7};

end;

tmin = find(EEG.times > time_win(1),1,'first');
tmax = find(EEG.times < time_win(2),1,'last');

% call function
% ---------------------------------------------------
if size(EEG.data,3) == 1
    error('EEG dataset is not epoched (EEG.data should of dimensions : nb channels x nb times x nb epochs)')
end

if use_ica == 0
    data = squeeze(EEG.data(channel_id,:,:))';
    points = squeeze(EEG.data(channel_id,tmin:tmax,:))';
else
    if ~isempty( EEG.icadata )
        data = squeeze(EEG.icadata(channel_id,:,:))';
        points = squeeze(EEG.icadata(channel_id,tmin:tmax,:))';
    else
        error('You must run ICA first');
    end;
end;
npoints = size(data,1);

%% Cross Validation
alphas = options.alpha;

if length(alphas) > 1 % Use Cross validation error if multiple alphas
    best_CVerr = -Inf;

    K = 10;
    disp(['--- Running K Cross Validation (K = ',num2str(K),')']);

    block_idx = fix(linspace(1,npoints,K+1)); % K cross validation
    for jj=1:length(alphas)
        options.alpha = alphas(jj);

        CVerr = 0;
        for kk=1:K
            bidx = block_idx(jj):block_idx(jj+1);
            idx = 1:npoints;
            idx(bidx) = [];

            data_k = data(idx,:);
            points_k = points(idx,:);
            [order,lags,coords,E_lags] = extractlag(points_k,options);

            data_reordered = data_k(order,:);
            lags = lags+tmin;
            [data_aligned,data_aligned_times] = perform_realign( data_reordered, EEG.times, lags );
            ep_evoked = mean(data_aligned);
            ep_evoked(isnan(ep_evoked)) = 0;
            ep_evoked(ep_evoked==Inf) = 0;

            ep_evoked = ep_evoked ./ norm(ep_evoked);

            data_k = data(bidx,:);
            data_k = normalize_rows(data_k);

            ep_raw = mean(data);
            ep_raw = ep_raw ./ norm(ep_raw);
            for pp=1:length(bidx)
                c = xcorr(ep_raw,data_k(pp,:));
                % max(c(:))
                c = xcorr(ep_evoked,data_k(pp,:));
                % max(c(:))
                [mc,mci] = max(c(:));
                % keyboard

                CVerr = CVerr + max(c(:));
            end
        end

        CVerr = CVerr/npoints;

        if CVerr > best_CVerr
            best_CVerr = CVerr;
            best_alpha = alphas(jj);
        end
    end

    options.alpha = best_alpha;
end

if use_maximum
    [order,lags,coords,E_lags] = extractlag( points, options );
else
    [order,lags,coords,E_lags] = extractlag( -points, options );
end

disp(['---------- Using alpha = ',num2str(options.alpha)]);

data_reordered = data(order,:);
lags = lags+tmin;
[data_aligned,data_aligned_times] = perform_realign( data_reordered, EEG.times, lags );

order_inv = perminv(order);
lags_no_order = lags(order_inv);
lags_ms = EEG.times(lags_no_order);

event_type = [num2str(time_win(1)),' - ',num2str(time_win(2)),' ms'];

if isfield(EEG,'chanlocs') % Real dataset (not synthetic)
    event_type = [EEG.chanlocs(channel_id).labels,' (',num2str(channel_id),') : ',event_type];
end

if use_maximum
    event_type = [event_type ' (max)'];
else
    event_type = [event_type ' (min)'];
end

% try

    clear events
    events = zeros(length(lags),3);

    for ii=1:length(lags)

        clear event
        event.type = event_type;
        event.latency = (ii-1)*length(EEG.times) + lags_no_order(ii);

        if ~isfield(EEG,'urevent')
            nurevent = 0;
            EEG.urevent(1) = event;
        else
            nurevent = length(EEG.urevent);
            EEG.urevent(end+1) = event;
        end

        event.code = 'Realignement';
        event.duration = 1;
        % event.channel = channel_id; % or 0 for all channels
        event.channel = 0; % or 0 for all channels
        event.time = [];

        % event.urevent = nurevent+1;
        event.epoch = ii;

        if isfield(EEG,'event')
            evidx = length(EEG.event)+1;
        else
            evidx = 1;
        end
        EEG.event(evidx).type = event.type;
        EEG.event(evidx).latency = event.latency;
        EEG.event(evidx).epoch = event.epoch;

        if ~isfield(EEG,'event')
            EEG.event(1) = event;
        else
            if isfield(EEG.event,'code'), EEG.event(evidx).code = event.code; end
            if isfield(EEG.event,'duration'), EEG.event(evidx).duration = event.duration; end
            if isfield(EEG.event,'channel'), EEG.event(evidx).channel = event.channel; end
            if isfield(EEG.event,'time'), EEG.event(evidx).time = event.time; end
            if isfield(EEG.event,'epoch'), EEG.event(evidx).epoch = event.epoch; end
        end
    end

    EEG = eeg_checkset(EEG);

% catch
%     warning('Problem with EEG set structure : Markers could not be correctly added to the dataset !')
% end

% return the string command
% -------------------------
com = sprintf('pop_extractlag( %s, %d, [%s], [%s] );', inputname(1), use_ica, int2str(channel_id), num2str(time_win));

%% view embedding results

t0 = find(EEG.times > 0,1,'first');

if options.verbose
    smart_figure('data');
    % plot(EEG.times,data(fix(linspace(1,size(data,1),4)),:));
    plot(EEG.times,data(fix(linspace(1,size(data,1),2)),:));
    axis tight
    xlabel('Time (ms)')
    line([time_win(1) time_win(1)],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([time_win(2) time_win(2)],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([0 0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');

    smart_figure('erp_data');
    imagesc(data,[min(data(:)),max(data(:))]);
    % colorbar
    xlabel('Time (ms)')
    ylabel('Trial')
    title('Data')
    set(gca,'XTickLabel',fix(EEG.times(get(gca,'XTick'))))
    line([tmin tmin],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([tmax tmax],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([t0 t0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');

    smart_figure('erp_data_ordered');
    imagesc(data_reordered,[min(data(:)),max(data(:))]);
    % colorbar
    xlabel('Time (ms)')
    ylabel('Trial')
    title('Data Ordered')
    set(gca,'XTickLabel',fix(EEG.times(get(gca,'XTick'))))
    % set(gcf,'WindowButtonDownFcn','callbackRaster(EEG)'); % set the callback
    line([tmin tmin],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([tmax tmax],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([t0 t0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');

    smart_figure('2D_embedding');
    plot(coords(:,1),coords(:,2),'+');
    xlabel('y_1')
    ylabel('y_2')
    if show_trial_number
        dxtext = std(coords(:,1))*0.1;
        dytext = std(coords(:,2))*0.1;
        for i=1:size(coords,1)
            text(coords(i,1)+dxtext,coords(i,2)+dytext,num2str(i));
        end
    end
    set(get(gca,'Title'),'String','2D Embedding');
end

if options.show_pca
    % [coeff, SCORE, LATENT] = princomp(points(order,:));
    [coeff, SCORE, LATENT] = princomp(zscore(points(order,:)));

    smart_figure('PCA embedding'); clf
    scatter(SCORE(:,1),SCORE(:,2),25,jet(npoints),'filled');
    % scatter3(SCORE(:,1),SCORE(:,2),SCORE(:,3),25,jet(npoints),'filled'); cameramenu
    axis square
end

%% view lag extraction results

if 1 & options.verbose
    smart_figure('erp_data_ordered_lags');
    imagesc(data_reordered,[min(data(:)),max(data(:))]);
    % colorbar
    xlabel('Time (ms)')
    ylabel('Trial')
    title('Data Ordered')
    set(gca,'XTickLabel',fix(EEG.times(get(gca,'XTick'))))
    % set(gcf,'WindowButtonDownFcn','callbackRaster(EEG)'); % set the callback
    line([tmin tmin],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([tmax tmax],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
    line([t0 t0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');
    hold on % super impose the lags and the smoothed maximums
    plot(lags,[1:npoints],'k+')
    plot(lags,[1:npoints],'k')
    line([mean(lags) mean(lags)],get(gca,'ylim'),'Color','b','LineWidth',1,'LineStyle','--')
    hold off
end

if 1 & options.verbose
    smart_figure('erp_data_aligned');
    imagesc(data_aligned,[min(data_aligned(:)),max(data_aligned(:))])
    xlabel('Time (ms)')
    ylabel('Trial')
    title('Data Ordered and aligned')
    set(gca,'XTickLabel',fix(EEG.times(get(gca,'XTick'))))
    line([mean(lags) mean(lags)],get(gca,'ylim'),'Color','b','LineWidth',1,'LineStyle','--')
    line([t0 t0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');
end

%% view evoked potential after realigning

ep_orig = mean(data_reordered); % evoked potential init
ep_aligned = mean(data_aligned); % evoked potential after realignment

smart_figure('evoked_potentials');
plot(EEG.times,ep_orig,'r');
hold on
plot(EEG.times,ep_aligned,'b');
hold off
xlabel('Time (ms)')
ylim = get(gca,'ylim');
ylim(1) = ylim(1) - 0.1 * abs(ylim(2)-ylim(1));
ylim(2) = ylim(2) + 0.1 * abs(ylim(2)-ylim(1));
set(gca,'ylim',ylim);
title('Naive averaging vs realigned averaging')
line([time_win(1) time_win(1)],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
line([time_win(2) time_win(2)],get(gca,'ylim'),'Color','k','LineWidth',1,'LineStyle','--');
line([0 0],get(gca,'ylim'),'Color','k','LineWidth',2,'LineStyle','-');
return;


text:
function out=downsample(inp, nsam, ave)
[m,n]=size(inp);
dim=2;
if (m>n)
	dim=1;
end

mx_a=cumsum(double(inp), dim);
if (dim==1)
    out=diff(mx_a(1:nsam:end,:),1,dim)/ave;
else
    out=diff(mx_a(:,1:nsam:end),1,dim)/ave;
end

repo: ekfmonoslam


text:
function logEmissionProb = emission(poseData,G,P,N,K)
logEmissionProb = zeros(N,K);
  for i = 1:N
	  data = reshape(poseData(i,:,:),10,3);
	  for j = 1:10
		  if G(j,1) == 1
			  parent = data(G(j,2),:);
			  for k = 1:K
				  theta = P.clg(j).theta(k,:);
				  mu_y = sum(theta(1:4).*[1,parent]);
				  mu_x = sum(theta(5:8).*[1,parent]);
				  mu_a = sum(theta(9:12).*[1,parent]);
				  logEmissionProb(i,k) += lognormpdf(data(j,1),mu_y,P.clg(j).sigma_y(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,2),mu_x,P.clg(j).sigma_x(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,3),mu_a,P.clg(j).sigma_angle(k));
			  end
		  else
			  for k = 1:K
				  logEmissionProb(i,k) += lognormpdf(data(j,1),P.clg(j).mu_y(k),P.clg(j).sigma_y(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,2),P.clg(j).mu_x(k),P.clg(j).sigma_x(k));
				  logEmissionProb(i,k) += lognormpdf(data(j,3),P.clg(j).mu_angle(k),P.clg(j).sigma_angle(k));
			  end
		  end
	  end
  end
  end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function mbasis = basis_matrix_b_uni ( )

%*****************************************************************************80
%
%% BASIS_MATRIX_B_UNI sets up the uniform B spline basis matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Foley, van Dam, Feiner, Hughes,
%    Computer Graphics: Principles and Practice,
%    page 493.
%
%  Parameters:
%
%    Output, real MBASIS(4,4), the basis matrix.
%
  mbasis = [
    -1.0 / 6.0,   3.0 / 6.0, -3.0 / 6.0, 1.0 / 6.0; ...
     3.0 / 6.0, - 6.0 / 6.0,  3.0 / 6.0, 0.0; ...
   - 3.0 / 6.0,   0.0,        3.0 / 6.0, 0.0; ...
     1.0 / 6.0,   4.0 / 6.0,  1.0 / 6.0, 0.0 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Fig. 9.8  Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
 n=[1 1];
 d = [1 0 0];
 rlocus(n,d)
 axis([-6 2 -3 3])
 hold on
 r=roots([1 1 1]);
 plot(r,'*')
  z=0:.1:.9;
 wn= 1:6;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function [pos2, tlight] = geocen (pos1, pe)

% this function moves the origin of coordinates from the
% barycenter of the solar system to the observer (or the
% geocenter).  i.e., this function accounts for parallax
% (annual + geocentric or just annual).

% input

%  pos1 = position vector of star or planet, with respect to
%         origin at solar system barycenter, components in au

%  pe   = position vector of observer (or the geocenter),
%         with respect to origin at solar system barycenter,
%         components in au

% output

%   pos2   = position vector of star or planet, with respect to
%            origin at observer (or the geocenter), components in au

%   tlight = light-time from star or planet to observer (or the
%            geocenter) in days

% ported from NOVAS 3.1

%%%%%%%%%%%%%%%%%%%%%%%

% speed of light in au/day

c = astcon('c(au/day)', 1.0);

for j = 1:1:3

    pos2(j) = pos1(j) - pe(j);

end

tlight = sqrt(pos2(1)^2 + pos2(2)^2 + pos2(3)^2) / c;

author: goGPS-Project


text:
function [ MRI_signal, bvalue ] = dt_from_smoldyn_pt2( positions, voxel_size, step_time, delta, Delta, gradient_directions, gradient_strength, bvalue )
%DT_FROM_SMOLDYN_PT2 Calculates the MRI signal
%
% PART 2 OF THE CALCULATION OF THE DIFFUSION TENSOR FROM THE SMOLDYN MODEL DATA
% This part calculates the diffusion MRI signal from the positions data.
% Note that providing each set of positions data contains the same number
% of molecules, the MRI signal from different sets of positions data can be
% combined by just taking the mean.
%
% Inputs are:
%   positions: matrix of size no of molecules x 3 x no of timesteps generated
%       from the Smoldyn program and cut down in part 1
%   voxel_size: scalar. the total (isotropic) size of the volume simulated (um)
%   step_time: scalar. the time of each step as used to generate the Smoldyn diffusion (us)
%   delta: scalar. the gradient pulse time (us)
%   Delta: scalar. the time between gradient pulses (us)
%   gradient_directions: matrix of size no of directions x 3 which gives
%       the x,y,z components of the gradient directions
%   gradient_strength: IF UNKNOWN, SET TO 0 & IT WILL BE CALCULATED FROM
%       BVALUE. scalar. (T/um)
%   bvalue: IF UNKNOWN, SET TO 0 & IT WILL BE CALCULATED FROM
%       GRADIENT_STRENGTH. scalar. (um^2/us)
%
% Outputs are:
%   MRI_signal: relative signal for each of the gradient directions
%   bvalue: as input, included as an output as needed later

% Author: Jo Bates <jobates81@gmail.com>
% Copyright © 2014 University of Oxford
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

no_of_directions = size(gradient_directions,1);
no_of_molecules = size(positions,1);
% number of jumps
no_of_timesteps = size(positions, 3);

% define the gyromagnetic ratio of water (rad/us/T)
gyro_ratio = 267.5;

% calculate the bvalue or gradient_strength if required
if bvalue == 0;
    bvalue = gradient_strength^2 * gyro_ratio^2 * delta^2 *(Delta - (delta/3));
end
if gradient_strength == 0;
    gradient_strength = sqrt(bvalue/(gyro_ratio^2 * delta^2 *(Delta - (delta/3))));
end

% resolution & range of the co-ordinates
xres = 1;
yres = 1;
zres = 1;
xmax = (voxel_size-1)/2;
ymax = (voxel_size-1)/2;
zmax = (voxel_size-1)/2;

% parameters independent of gradient direction or voxel
pulse_1_end_step = (delta/step_time) + 1;
pulse_2_start_step = (Delta/step_time)+1;
pulse_2_end_step = ((Delta+delta)/step_time) + 1;

% calc parameters which are independent of voxel
start_corner = zeros(no_of_directions, 3);
for g = 1:no_of_directions
    % work out the starting corner of the image, where the MRI pulse with be 0
    sign_x = sign(gradient_directions(g,1));
    sign_y = sign(gradient_directions(g,2));
    sign_z = sign(gradient_directions(g,3));
    if sign_x < 0
        x_1 = xmax*2+1;
    else
        x_1 = 1;
    end
    if sign_y < 0
        y_1 = ymax*2+1;
    else
        y_1 = 1;
    end
    if sign_z < 0
        z_1 = zmax*2+1;
    else
        z_1 = 1;
    end
    start_corner (g,:) = [x_1, y_1, z_1];
end

MRI_signal = zeros(no_of_directions,1);
% for every gradient direction
for g = 1:no_of_directions
    x_1 = start_corner (g,1);
    y_1 = start_corner (g,2);
    z_1 = start_corner (g,3);
        % for each molecule, for each position, the position in gradient
        % direction is the dot product of the position times the gradient
        % direction. since A.B = a1*b1 + a2*b2 + a3*b3, calculate like this as
        % need to include the potential for different resolutions in each
        % direction and can't seem to do dot product of correct parts for these
        % 2 matrices.

        x_pos = (positions(:,1,:)-x_1).*xres.*gradient_directions(g,1);
        y_pos = (positions(:,2,:)-y_1).*yres.*gradient_directions(g,2);
        z_pos = (positions(:,3,:)-z_1).*zres.*gradient_directions(g,3);
        position_grad_dirn = x_pos + y_pos + z_pos;

        % calculate phase shift for each molecule
        mol_pulse = zeros(size(position_grad_dirn));
        mol_pulse(:,:,1:pulse_1_end_step) = position_grad_dirn(:,:,1:pulse_1_end_step).*gradient_strength;
        mol_pulse(:,:,pulse_2_start_step:pulse_2_end_step) = -position_grad_dirn(:,:,pulse_2_start_step:pulse_2_end_step).*gradient_strength;
        molpulsesum = sum(mol_pulse,3);

        % calc MRI signal
        phasech = gyro_ratio*molpulsesum*step_time;
        cosphase = cos(phasech);
        MRI_signal(g) = mean(cosphase);
end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function save_mat(key,outfilename)
% Saves the Key to a mat file.
% Inputs:
%   key: An object of type Key.
%   outfilename: The name for the output mat file.

assert(nargin==2)
assert(isa(key,'Key'))
assert(isa(outfilename,'char'))
assert(key.validate())

warning('off','MATLAB:structOnObject')
key = struct(key);
warning('on','MATLAB:structOnObject')
save(outfilename,'key');

author: nesl


text:
addpath('../core');
addpath('../helpers');

run ../vlfeat/toolbox/vl_setup ;
run ../matconvnet/matlab/vl_setupnn ;

if(~exist('data', 'dir'))
  mkdir('data');
end

if(~exist('data/results', 'dir'))
  mkdir('data/results');
end

if(~exist('networks/imagenet-caffe-ref.mat', 'file'))
  cmd='wget http://www.robots.ox.ac.uk/~aravindh/imagenet-caffe-ref.mat -O networks/imagenet-caffe-ref.mat';
  system(cmd);
  clear cmd;
end

addpath('networks')

author: aravindhm

repo: deep-goggle


text:
function test_bug1652

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_datatype_segmentation

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1652'));

load seg1.mat
% this contains TPMs
seg1_i  = ft_datatype_segmentation(seg1, 'segmentationstyle', 'indexed');
seg1_p  = ft_datatype_segmentation(seg1, 'segmentationstyle', 'probabilistic');
% seg1_ib = ft_datatype_segmentation(seg1, 'hasbrain', 'yes', 'segmentationstyle', 'indexed')  % this one is not supported due to the combination of the TPMs with the brain
seg1_pb = ft_datatype_segmentation(seg1, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic');
clear seg1*

load seg2.mat
% this one has three BEM tissues
% note that the skull and scalp are not consistent (one is cumulative, the other exclusive)
seg2_i  = ft_datatype_segmentation(seg2, 'segmentationstyle', 'indexed');
seg2_p  = ft_datatype_segmentation(seg2, 'segmentationstyle', 'probabilistic');
seg2_ib = ft_datatype_segmentation(seg2, 'hasbrain', 'yes', 'segmentationstyle', 'indexed');
seg2_pb = ft_datatype_segmentation(seg2, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic');
clear seg2*

load seg3.mat
% this contains three cumulative or overlapping BEM tissues
seg3_i  = ft_datatype_segmentation(seg3, 'segmentationstyle', 'indexed');
seg3_p  = ft_datatype_segmentation(seg3, 'segmentationstyle', 'probabilistic');
seg3_ib = ft_datatype_segmentation(seg3, 'hasbrain', 'yes', 'segmentationstyle', 'indexed');
seg3_pb = ft_datatype_segmentation(seg3, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic');
clear seg3*

load seg4.mat
% this contains three exclusive BEM tissues
seg4_i  = ft_datatype_segmentation(seg4, 'segmentationstyle', 'indexed');
seg4_p  = ft_datatype_segmentation(seg4, 'segmentationstyle', 'probabilistic');
seg4_ib = ft_datatype_segmentation(seg4, 'hasbrain', 'yes', 'segmentationstyle', 'indexed');
seg4_pb = ft_datatype_segmentation(seg4, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic');
clear seg4*

load seg5.mat
% this one is an indexed representation of the BEM tissues, without labels
seg5_i  = ft_datatype_segmentation(seg5, 'segmentationstyle', 'indexed');
seg5_p  = ft_datatype_segmentation(seg5, 'segmentationstyle', 'probabilistic');

seg5.seglabel = {'scalp', 'skull', 'brain'};
seg5_ib = ft_datatype_segmentation(seg5, 'hasbrain', 'yes', 'segmentationstyle', 'indexed')         % this fails without the labels;
seg5_pb = ft_datatype_segmentation(seg5, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic')   % this fails without the labels;
clear seg5*

load seg6.mat
% this only contains a single boolean brain
seg6_i  = ft_datatype_segmentation(seg6, 'segmentationstyle', 'indexed');
seg6_p  = ft_datatype_segmentation(seg6, 'segmentationstyle', 'probabilistic');
seg6_ib = ft_datatype_segmentation(seg6, 'hasbrain', 'yes', 'segmentationstyle', 'indexed');
seg6_pb = ft_datatype_segmentation(seg6, 'hasbrain', 'yes', 'segmentationstyle', 'probabilistic');
clear seg6*

% try it with the AFNI atlas
atlas = ft_read_atlas('TTatlas+tlrc.BRIK');
atlas_i  = ft_datatype_segmentation(atlas, 'segmentationstyle', 'indexed');
atlas_p  = ft_datatype_segmentation(atlas, 'segmentationstyle', 'probabilistic');


text:
function [p,priors]=create_parameters12(ncp_shocks,start_at_mode)
if nargin<2
    start_at_mode=false;
    if nargin<1
        ncp_shocks=false;
    end
end

p=struct();

p.beta = 0.99;
p.mu = 0.34;
p.gam = 2;
p.phi_d = 0.001;
p.delta = 0.025;
p.alpha = 0.36;
p.theta = 1.05;
p.omega = 0.10;

if ncp_shocks
    p.mss_H = 0.9997;
    p.mss_F = 0.9997;
    p.m_hf_ss=1;
    p.zss_H = 1.0045;
else
    p.mss_H = 0.9;
    p.mss_F = 0.9;
    p.zss_H = 1.0042;
end

p.zss_F = p.zss_H;
p.vss_H = 1.0004;
p.vss_F = p.vss_H;
p.z_hf_ss = 1;
p.v_hf_ss = 1;
p.gss_H = 0.20;
p.gss_F = 0.20;
p.rho_m_HF = 0;
p.rho_m_FH = 0;
p.rho_z_HF = 0;
p.rho_z_FH = 0;
p.rho_g_HF = 0;
p.rho_g_FH = 0;
p.rho_v_HF = 0;
p.rho_v_FH = 0;

priors=struct();

if start_at_mode
    priors.theta = {1.5709,0,5};

    priors.philtr = {17.113/100,0,3};
    priors.phiktr = {2.971/100,0,3};

    if ncp_shocks
        priors.kappamhtr = {0.0077294,0,3};
        priors.kappa_m_F = {8.0501e-06,0,3};
    end
    priors.kappazhtr = {0.002383,0,3};
    priors.kappa_z_F = {1.6159e-05,0,3};

    priors.kappavhtr = {0.017416,0,3};
    priors.kappa_v_F = {sqrt(eps),0,3};

    priors.rho_m_HH = {0.9671,0,1};
    priors.rho_m_FF = {0.9924,0,1};

    priors.rho_z_HH = {0.1752,0,1};
    priors.rho_z_FF = {0.35979,0,1};

    priors.rho_v_HH = {0.15793,0,1};
    priors.rho_v_FF = {0.9834,0,1};

    priors.rho_g_HH = {0.9695,0,1};
    priors.rho_g_FF = {0.9414,0,1};

    priors.sig_m_H = {0.017084,0,2};
    priors.sig_m_F = {0.0063409,0,2};

    priors.sig_z_H = {0.01176,0,2};
    priors.sig_z_F = {0.0085657,0,2};

    priors.sig_v_H = {0.01436,0,2};
    priors.sig_v_F = {0.00060192,0,2};

    priors.sig_g_H = {0.018465,0,2};
    priors.sig_g_F = {0.0098922,0,2};
else
    thetatr = 1.5709;
    priors.theta = {abs(thetatr),0,5};

    priors.philtr = {17.1310/100,0,3};
    priors.phiktr = {2.7056/100,0,3};

    if ncp_shocks
        kappamftr = 0.0010;
        priors.kappamhtr = {0.0077,0,3};
        priors.kappa_m_F = {abs(kappamftr),0,3};
    end

    priors.kappazhtr = {0.0018,0,3};
    kappazftr = 0.0010;
    priors.kappa_z_F = {abs(kappazftr),0,3};

    priors.kappavhtr = {0.0089,0,3};
    kappavftr = 0.0010;
    priors.kappa_v_F = {abs(kappavftr),0,3};

    rhomhhtr = sqrt(0.9671/(1-0.9671));
    priors.rho_m_HH = {rhomhhtr^2/(1+rhomhhtr^2),0,1};
    rhomfftr = sqrt(0.9924/(1-0.9924));
    priors.rho_m_FF = {rhomfftr^2/(1+rhomfftr^2),0,1};

    rhozhhtr = sqrt(0.1752/(1-0.1752));
    priors.rho_z_HH = {rhozhhtr^2/(1+rhozhhtr^2),0,1};
    rhozfftr = sqrt(0.3598/(1-0.3598));
    priors.rho_z_FF = {rhozfftr^2/(1+rhozfftr^2),0,1};

    rhovhhtr = sqrt(0.1581/(1-0.1581));
    priors.rho_v_HH = {rhovhhtr^2/(1+rhovhhtr^2),0,1};
    rhovfftr = sqrt(0.9834/(1-0.9834));
    priors.rho_v_FF = {rhovfftr^2/(1+rhovfftr^2),0,1};

    rhoghhtr = sqrt(0.9695/(1-0.9695));
    priors.rho_g_HH = {rhoghhtr^2/(1+rhoghhtr^2),0,1};
    rhogfftr = sqrt(0.9414/(1-0.9414));
    priors.rho_g_FF = {rhogfftr^2/(1+rhogfftr^2),0,1};

    sigmhtr = 0.01;
    priors.sig_m_H = {abs(sigmhtr),0,2};
    sigmftr = 0.01;
    priors.sig_m_F = {abs(sigmftr),0,2};

    sigzhtr = 0.01;
    priors.sig_z_H = {abs(sigzhtr),0,2};
    sigzftr = 0.01;
    priors.sig_z_F = {abs(sigzftr),0,2};

    sigvhtr = 0.01;
    priors.sig_v_H = {abs(sigvhtr),0,2};
    sigvftr = 0.01;
    priors.sig_v_F = {abs(sigvftr),0,2};

    sigghtr = 0.01;
    priors.sig_g_H = {abs(sigghtr),0,2};
    siggftr = 0.01;
    priors.sig_g_F = {abs(siggftr),0,2};
end

% add the initial conditions of the priors to the parameters
%------------------------------------------------------------
fields=fieldnames(priors);
for ip=1:numel(fields)
    name=fields{ip};
    p.(name)=priors.(name){1};
end

end

author: jmaih


text:
classdef MimRemoveDataset < MimGuiPlugin
    % MimRemoveDataset. Gui Plugin for removing a dataset.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimRemoveDataset is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will delete all results files from the current
    %     dataset results cache folder. Certain internal cache files will not be
    %     removed.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Delete dataset'
        SelectedText = 'Delete dataset'
        ToolTip = 'Remove this dataset from the Toolkit'
        Category = 'Dataset'
        Visibility = 'Dataset'
        Mode = 'Toolbar'
        Icon = 'bin.png'
        Location = 3

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            % Delete files from the disk cache
            gui_app.DeleteThisImageInfo;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = false;
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
% KM_DEMO_KPCA_U Kernel principal component analysis (KPCA) on a U-shaped
% two-dimensional data set.
%
% This program implements the example shown in Figure 2.4 of "Kernel
% Methods for Nonlinear Identification, Equalization and Separation of
% Signals", Ph.D. dissertation by S. Van Vaerenbergh.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2010.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

close all
clear

%% PARAMETERS

N_split = [30, 30, 60];		% number of data points in 2 straight parts and curve of "U"
corners = [3 4 2 4];	% corners: x1 x2 y1 y2
nvar = 0.05;	% noise variance

kernel.type = 'gauss';
kernel.par =  .5;
numeig = 4;	% number of eigenvalues to plot

Ntest = [50,50];	% number of test grid divisions, in each dimension
border = [0 6 0 6];	% test grid border

%% PROGRAM
tic

%% generate data
N = sum(N_split);
N1 = N_split(1); N2 = N_split(2); N3 = N_split(3);

X1 = [linspace(corners(1),corners(2),N1)' repmat(corners(3),N1,1)];	% lower straight part
X2 = [linspace(corners(1),corners(2),N2)' repmat(corners(4),N2,1)];	% upper straight part

angles = linspace(pi/2,3*pi/2,N3)';
rad = (corners(4)-corners(3))/2;
m3 = [corners(1) corners(3)+rad];
X3 = repmat(m3,N3,1) + rad*[cos(angles) sin(angles)];

n = nvar*randn(N,2);
X = [X1; X2; X3] + n;

%% generate test grid data
Nt1 = Ntest(1); Nt2 = Ntest(2);
Xtest = zeros(Nt1*Nt2,2);
absc = linspace(border(1),border(2),Nt1);
ordi = linspace(border(3),border(4),Nt2);
for i=1:Nt1,
	for j=1:Nt2,
		Xtest((i-1)*Nt2+j,:) = [absc(i) ordi(j)];
	end
end

%% calculate kernel principal components and projections of Xtest
[E,v] = km_kpca(X,numeig,kernel.type,kernel.par);

Kt = km_kernel(X,Xtest,kernel.type,kernel.par);	% kernels of test data set
Xtestp = E'*Kt;	 % projections of test data set on the principal directions

Y = cell(numeig,1);
for i=1:numeig,
 	Y{i} = reshape(Xtestp(i,:),Nt2,Nt1);	% shape into 2D grid
end

toc
%% OUTPUT

figure;
plot(X(:,1),X(:,2),'.')
axis(border)

% fireworks!
for i=1:numeig,
	figure; hold on
	[C,h] = contourf(-interp2(Y{i},2),25);
	plot(X(:,1)/border(2)*(Nt1*4-1)+1,X(:,2)/border(4)*(Nt2*4-1)+1,'o',...
		'MarkerFaceColor','White','MarkerEdgeColor','Black')
	set(gca, 'XTick', [])
	set(gca, 'YTick', [])
end

author: epfl-lasa


text:
O = 3;
Q = 2;

% "true" parameters
prior0 = normalise(rand(Q,1));
transmat0 = mk_stochastic(rand(Q,Q));
obsmat0 = mk_stochastic(rand(Q,O));

% training data
T = 5;
nex = 10;
data = dhmm_sample(prior0, transmat0, obsmat0, T, nex);

% initial guess of parameters
prior1 = normalise(rand(Q,1));
transmat1 = mk_stochastic(rand(Q,Q));
obsmat1 = mk_stochastic(rand(Q,O));

% improve guess of parameters using EM
[LL, prior2, transmat2, obsmat2] = dhmm_em(data, prior1, transmat1, obsmat1, 'max_iter', 5);
LL

% use model to compute log likelihood
loglik = dhmm_logprob(data, prior2, transmat2, obsmat2)
% log lik is slightly different than LL(end), since it is computed after the final M step


text:
report_this_filefun(mfilename('fullpath'));

str = [];
out = reshape(bv2,33759,1,1);
out = out*100;

[newmatfile, newpath] = uiputfile([ hodo '*.dat'], 'Save As'); %Syntax change Matlab Version 7, no window positioning on macs

s = [33  ;33 ; 31  ;out(:,1)   ];
s = s'
fid = fopen([newpath newmatfile],'w') ;;
fwrite(fid,s,'integer*8');
fclose(fid);
clear s
return

author: CelsoReyes


text:
function [n,dndr,T,P]=SinclairAtmos(h,phl0,Rh,P0,T0,wl,ht)
%SINCLAIRATMOS   Compute atmospheric parameters for the Sinclair
%                atmospheric model. This simple model is often used for
%                computing atmospheric refraction when viewing stars. The
%                index of refraction in the model depends on the height
%                above mean sea level. In this implementation, mean sea
%                level is approximated by the surface of the WGS-84
%                reference ellipsoid. The model is not sufficiently precise
%                that geoid undulations should matter and it also shouldn't
%                matter if a geocentric latitude is substituted for the
%                geodetic latitude.
%
%INPUTS: h  The vector or matrix of heights above sea level where the
%           indices of refraction and other parameters are desired.
%      phl0 The location of an observer in WGS-84 ellipsoidal coordinates
%           [latitude;longitude;ellipsoidal height] in the troposphere
%           measuring the relative humidity, ambient air pressure and the
%           temperature. Only the latitude in radians and the ellipsoidal
%           height in meters matter; the longitude is ignored.
%       Rh  The relative humidity at the observer (between 0 and 1). If
%           this parameter is omitted or an empty matrix is passed, then
%           Constants.standardRelHumid is used.
%       P0  The atmospheric pressure at the observer in Pascals (N/m^2). If
%           this parameter is omitted or an empty matrix is passed, then
%           Constants.standardAtmosphericPressure is used.
%       T0  The air temperature at the observer in degrees Kelvin. If this
%           parameter is omitted or an empty matrix is passed, then
%           Constants.standardTemp is used.
%       wl  The wavelength at which the observation is made in units of
%           meters. If this parameter is omitted or an empty matrix is
%           passed, then a wavelength of 0.574 micrometers is used, which
%           is in the visible spectrum (a rather yellow color). This
%           parameter only matters for determining the narrowband index of
%           refraction at a given .
%       ht  The assumed height of the troposphere in meters. This parameter
%           specifies when a stratospheric model begins being used. If this
%           parameter is omitted, then the default value of 11000m is used.
%
%OUTPUTS: n A matrix providing the index of refraction at all of the
%           distances given in r.
%      dndr A matrix providing the derivative of the index of refraction
%           with respect to r evaluated at all of the values given in r.
%        T  A matrix giving the temperature in degrees Kelvin in the
%           atmospheric model at the distances given in r.
%        P  A matrix giving the barometric pressures in Pascals  in the
%           atmospheric model at the distances given in r.
%
%The Sinclair atmospheric model is described in Chapter 7.2 of [1] and in
%[2]. The original source is cited in [2] as being [3]. Reference 3 was not
%consulted in implementing this algorithm.
%
%In the references, the atmosphere is defined in terms of radial
%distances from the center of the Earth. However, the radials distances
%only appear as differences, so heights above mean sea level (MSL) can be
%substituted. MSL in this implementations is just taken to be the surface
%of the WGS-84 reference ellipsoid. The model makes a number of local
%spherical Earth approximations.
%
%REFERENCES:
%[1] S. E. Urban and K. P.Seidelmann, Eds.,Explanatory Supplement to the
%    Astronomical Almanac, 3rd ed. Mill Valley, CA: University Science
%    Books, 2013.
%[2] C. Y. Hohenkerk and A. T. Sinclair, "The computation of angular
%    atmospheric refraction at large zenith angles," United Kingdom
%    Hydrographic Office, HM Nautical Almanac, Tech. Rep. 63, Apr. 1985.
%    http://astro.ukho.gov.uk/data/tn/naotn63.pdf
%[3] A. T. Sinclair, "The effect of atmospheric refraction on laser ranging
%    data," United Kingdom Hydrographic Office, HM Nautical Almanac, Tech.
%    Rep. 59, 1982.
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<7||isempty(ht))
       ht=11000;%Assumed top of the troposphere in meters.
    end

    if(nargin<6||isempty(wl))
       wl=0.574e-6;
    end

    if(nargin<5||isempty(T0))
       T0=Constants.standardTemp;
    end

    if(nargin<4||isempty(P0))
        P0=Constants.standardAtmosphericPressure;
    end

    if(nargin<3||isempty(Rh))
        Rh=Constants.standardRelHumid;
    end

    %The geodetic latitude in radians.
    phi=phl0(1);
    %The height above the reference ellipsoid in meters is used in
    %place of the height above the geoid.
    h0=phl0(3);
    phi=ellips2Sphere(phi);

    %%SET ATMOSPHERIC MODEL PARAMETERS.
    %Convert the pressure from Pascals to millibars.
    P0=P0*0.01;
    %Convert the wavelength from meters to micrometers.
    lambda=wl*1e6;
    %The set of constants is from Equation 7.82 in [1].
    R=1000*Constants.molarGasConstant;%(R) J/(kilomol K)
    Md=28.966;%Assumed molecular mass (im amu) of dry air.
    %The molecular mass (in amu) of water.
    Mw=2*Constants.elementAMU(1)+Constants.elementAMU(8);
    %Exponent of the temperature dependence of water vapor pressure.
    %This is used in an approximate conversion to get the partial
    %pressure of water vapor from the relative humidity.
    delta=18.36;
    alpha=0.0065;%Tropospheric lapse rate of temperature in K/m.

    %The following set of Equations is from Equation 7.83 in [1].
    %The partial pressure of water vapor at the observer in millibars.
    Pw0=Rh*(T0/247.1)^delta;
    gBar=9.784*(1-0.0026*cos(2*phi)-0.00000028*h0);
    A=(287.604+1.6288/lambda^2+0.0136/lambda^4)*(273.15/1013.25)*1e-6;
    C2=gBar*Md/R;
    gamma=C2/alpha;
    C5=Pw0*(1-Mw/Md)*gamma/(delta-gamma);
    C6=A*(P0+C5)/T0;
    C7=(A*C5+11.2684e-6*Pw0)/T0;
    C8=alpha*(gamma-1)*C6/T0;
    C9=alpha*(delta-1)*C7/T0;

    %Allocate space for the return variables.
    n=zeros(size(h));
    dndr=zeros(size(h));
    T=zeros(size(h));
    P=zeros(size(h));

    %First, deal with all of the points that are in the stratosphere.
    sel=h>ht;

    %First, compute the refraction and temperature at the top of the
    %troposphere using Equation 7.85 in [1].
    Tt=T0-alpha*(ht-h0);
    T(sel)=Tt;
    TRat=Tt/T0;
    nt=1+(C6*TRat^(gamma-2)-C7*TRat^(delta-2))*TRat;
    %Approximate partial vapor pressure of water at the tropopause,
    %taken from page 4 of [2].
    Pwt=Pw0*TRat^delta;
    %Air pressure at the tropopause taken from page 4 of [2].
    Pt=(P0+C5)*TRat^(gamma)-Pwt*(1-Mw/Md)*gamma/(delta-gamma);

    %The modelled refraction in the stratosphere.
    %Equation 7.86 in [1].
    n(sel)=1+(nt-1).*exp(-C2*(h(sel)-ht)./Tt);
    dndr(sel)=-(C2./Tt).*(nt-1).*exp(-C2*(h(sel)-ht)./Tt);

    %Stratospheric air pressure, taken from page 5 of [2]. The 0.01
    %term converts the pressure from millibars to Pascals.
    P(sel)=Pt.*exp(-C2*(h(sel)-ht)./Tt)/0.01;

    %Next, deal with the points in the troposphere.
    sel=~sel;
    %The modelled refraction in the troposphere.
    %Equation 7.85 in [1].
    T(sel)=T0-alpha*(h(sel)-h0);
    TRat=T(sel)/T0;
    n(sel)=1+(C6*TRat.^(gamma-2)-C7*TRat.^(delta-2)).*TRat;
    dndr(sel)=-C8*TRat.^(gamma-2)+C9*TRat.^(delta-2);

    %Approximate partial vapor pressure of water at altitude, taken
    %from page 4 of [2].
    Pw=Pw0*TRat.^delta;
    %Air pressure at altitude taken from page 4 of [2]. The 0.01 term
    %converts the pressure from millibars to Pascals.
    P(sel)=((P0+C5).*TRat.^(gamma)-Pw*(1-Mw/Md)*gamma/(delta-gamma))/0.01;
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 y = gsp_smooth_downstep(x, a, o)
%GSP_SMOOTH_DOWNSTEP Smooth downstep function from 1 to 0
%   Usage: y = gsp_smooth_downstep(x);
%          y = gsp_smooth_downstep(x, a);
%          y = gsp_smooth_downstep(x, a, o);
%
%   Input parameters:
%          x        : input value
%          a        : smoothing parameter (default 1)
%          o        : offset (default 1)
%   Ouput parameters:
%          y        : output value of the function
%
%   This function is a smooth downstep from 1 to 0 arround *o* using the
%   following function:
%
%   .. f(x) = s(o-x)
%
%   .. math:: f(x) = s(o-x)
%
%   where
%
%   ..          /   0                                      if x < -1
%   ..  s(x) = | exp(-a/x) / ( exp(-a/x) + exp(-a/(1-x)) ) if x in [-1, 1]
%   ..          \   1                                      if x > 1
%
%   .. math:: s(x)=\begin{cases} 0 & \mbox{if }x<-1 \\ \frac{e^{-\frac{a}{x}}}{e^{-\frac{a}{x}}+e^{-\frac{a}{1-x}}} & \mbox{if }x\in[-1,1]\\ 1 & \mbox{if }x>1 \end{cases}
%

% Author: Nathanael Perraudin
% Date  : 5 mai 2016

if nargin<3
    o = 1;
end

if nargin<3
    a = 1;
end

y = gsp_smooth_step(o-x, a);

end

repo: gspbox


text:
function row_sum = r8row_sum ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_SUM returns the sums of an R8ROW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the R8ROW
%
%    Output, real ROW_SUM(M), the sum of the entries of
%    each row.
%
  for i = 1 : m
    row_sum(i) = sum ( a(i,1:n) );
  end

  return
end


text:
function log_L = Kalman_Estimation(y, psi, matur, dt, a0, P0, N, nobs, locked_parameters)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extracting initial parameter values from initial psi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = psi(1,1);
sigmax = psi(2,1);
lambdax = psi(3,1);
mu = psi(4,1);
sigmae = psi(5,1);
rnmu = psi(6,1);
pxe = psi(7,1);

if sum(locked_parameters) == 0
    k = psi(1,1);
    sigmax = psi(2,1);
    lambdax = psi(3,1);
    mu = psi(4,1);
    sigmae = psi(5,1);
    rnmu = psi(6,1);
    pxe = psi(7,1);

    s = zeros(1, size(psi,1)-7);
    for i = 1:size(s,2)
        s(1, i) = psi(i+7,1);
    end
end

if sum(locked_parameters) ~= 0
    s = zeros(1, size(psi,1)-7+size(locked_parameters,1));
    j = 1;
    for i = 1:size(s,2)
        if all(abs(i-(locked_parameters))) == 1
             s(1, i) = psi(7+j,1);
             j = j+1;
        end
    end
end

% m = Number of state variables (number of rows in a0)
m = size(a0,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE TRANSITION EQUATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S&S NOTATION: x(t)=c+G*x(t-1)+w(t)        w~N(0,W)    Equation (14)
% NEW NOTATION: a(t)=c+T*a(t-1)+R(t)*n(t)   n~N(0,Q)

% c is a {m x 1} Vector
% T is a {m x m} Matrix
c=[0;mu*dt];
T=[exp(-k*dt),0;0,1];

% Defining Q = var[n(t)] and R
xx=(1-exp(-2*k*dt))*(sigmax)^2/(2*k);
xy=(1-exp(-k*dt))*pxe*sigmax*sigmae/k;
yx=(1-exp(-k*dt))*pxe*sigmax*sigmae/k;
yy=(sigmae)^2*dt;
Q=[xx,xy;yx,yy];
R=eye(size(Q,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE MEASUREMENT EQUATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S&S NOTATION: y(t)=d(t)+F(t)'x(t)+v(t)    v~N(0,V) Equation (15)
% NEW NOTATION: y(t)=d(t)+Z(t)a(t)+e(t)     e~N(0,H)

% d is a {N x 1} Vector
% Z is a {N x m} Matrix
    for i=1:N
        p1=(1-exp(-2*k*matur(i)))*(sigmax)^2/(2*k);
        p2=(sigmae)^2*matur(i);
        p3=2*(1-exp(-k*matur(i)))*pxe*sigmax*sigmae/k;
        d(i,1)=rnmu*matur(i)-(1-exp(-k*matur(i)))*lambdax/k+.5*(p1+p2+p3);
        Z(i,1)=exp(-k*matur(i));
        Z(i,2)=1;
    end

% Measurment errors Var-Cov Matrix: Cov[e(t)]=H
H=diag(s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUNNING THE KALMAN FILTER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating placeholder vectors/matrices for variables to be stored in
global save_vt save_att save_dFtt_1 save_vFv save_vtt save_Ptt_1 save_Ftt_1 save_Ptt

save_ytt_1 = zeros(nobs,N);
save_vtt = zeros(nobs,N);
save_vt    = zeros(nobs,N);
save_att_1 = zeros(nobs,m);
save_att   = zeros(nobs,m);
save_Ptt_1 = zeros(nobs,m*m);
save_Ptt   = zeros(nobs,m*m);
save_Ftt_1 = zeros(nobs,N*N);
save_dFtt_1 = zeros(nobs,1);
save_vFv    = zeros(nobs,1);
%save_log_Lt   = zeros(nobs,1);

Ptt = P0;
att = a0;

% Running the kalman filter for t = 1,...,nobs
    for t = 1:nobs
        Ptt_1   = T*Ptt*T'+R*Q*R';
        Ftt_1   = Z*Ptt_1*Z'+H;
        dFtt_1  = det(Ftt_1);

        %Ptt_1_test = [Ptt_1(1,1) 0; 0 Ptt_1(2,2)];
        %Ftt_1_test   = Z*Ptt_1_test*Z'+H;
        %dFtt_1_test  = det(Ftt_1_test);

        att_1   = T*att + c;
        yt      = y(t,:)';
        ytt_1   = Z*att_1+d;
        vt      = yt-ytt_1;

        att = att_1 + Ptt_1*Z'*inv(Ftt_1)*(vt);
        Ptt = Ptt_1 - Ptt_1*Z'*inv(Ftt_1)*Z*Ptt_1;

        ytt = Z*att+d;
        vtt  = yt-ytt;

        % save_ytt_1(t,:) = ytt_1';
        save_vtt(t,:) = vtt';
        save_vt(t,:)    = (vt)';
        % save_att_1(t,:) = att_1';
        save_att(t,:)   = att';
        save_Ptt_1(t,:) = [Ptt_1(1,1), Ptt_1(1,2), Ptt_1(2,1), Ptt_1(2,2)];
        save_Ptt(t,:)   = [Ptt(1,1), Ptt(1,2), Ptt(2,1), Ptt(2,2)];
        % save_Ftt_1(t,:) = [Ftt_1(1,1), Ftt_1(1,2), Ftt_1(1,3), Ftt_1(1,4), Ftt_1(1,5), Ftt_1(2,1), Ftt_1(2,2), Ftt_1(2,3), Ftt_1(2,4), Ftt_1(2,5), Ftt_1(3,1), Ftt_1(3,2), Ftt_1(3,3), Ftt_1(3,4), Ftt_1(3,5), Ftt_1(4,1), Ftt_1(4,2), Ftt_1(4,3), Ftt_1(4,4), Ftt_1(5,5), Ftt_1(5,1), Ftt_1(5,2), Ftt_1(5,3), Ftt_1(5,4), Ftt_1(5,5)];

        %save_dFtt_1(t,:)= dFtt_1_test;
        %save_vFv(t,:)   = vt'*inv(Ftt_1_test)*vt;
        save_dFtt_1(t,:)= dFtt_1;
        save_vFv(t,:)   = vt'*inv(Ftt_1)*vt;

    end

logL = -(N*nobs/2)*log(2*pi)-0.5*sum(log(save_dFtt_1))-0.5*sum(save_vFv);
% logL = -(N*nobs/2)*log(2*pi)-0.5*sum(save_vFv);
% logL = sum(diag(save_vt'*save_vt));
log_L = -logL;


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

    enumeration
        ConicalSensor('Conical Sensor');
        RectangularSensor('Rectangular Sensor');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = SensorEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('SensorEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('SensorEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

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

author: Arrowstar


text:
function X = full(X)
%FULL Convert to a (dense) tensor.
%
%   FULL(X) returns X, i.e., does nothing.
%
%   See also SPTENSOR/FULL, KTENSOR/FULL, TTENSOR/FULL.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

return;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [w, state] = adadelta(w, state, grad, opts, ~)
%ADADELTA
%   Example AdaDelta solver, for use with CNN_TRAIN and CNN_TRAIN_DAG.
%
%   AdaDelta sets its own learning rate, so any learning rate set in the
%   options of CNN_TRAIN and CNN_TRAIN_DAG will be ignored.
%
%   If called without any input argument, returns the default options
%   structure.
%
%   Solver options: (opts.train.solverOpts)
%
%   `epsilon`:: 1e-6
%      Small additive constant to regularize variance estimate.
%
%   `rho`:: 0.9
%      Moving average window for variance update, between 0 and 1 (larger
%      values result in slower/more stable updating).

% Copyright (C) 2016 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin == 0  % Return the default solver options
  w = struct('epsilon', 1e-6, 'rho', 0.9) ;
  return ;
end

if isequal(state, 0)  % First iteration, initialize state struct
  state = struct('g_sqr', 0, 'delta_sqr', 0) ;
end

rho = opts.rho ;

state.g_sqr = state.g_sqr * rho + grad.^2 * (1 - rho) ;
new_delta = -sqrt((state.delta_sqr + opts.epsilon) ./ ...
                  (state.g_sqr + opts.epsilon)) .* grad ;
state.delta_sqr = state.delta_sqr * rho + new_delta.^2 * (1 - rho) ;

w = w + new_delta ;

repo: LapSRN