text:
function out_projector_fif(FifFile, ChannelNames, Projectors)
% OUT_PROJECTOR_FIF: Saves a list of SSP projectors in Brainstorm format in a FIF file.

% @=============================================================================
% 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

% Get fiff constants
global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

% Convert projectors to FIF format
projs = repmat(struct('kind',[], 'active',[], 'desc',[], 'data', []), 0);
% Loop on all the projector categories
for iCat = 1:length(Projectors)
    % Skip projectors that are not selected
    if (Projectors(iCat).Status == 0)
        continue;
    end
    % Get the selected components for this projector category
    iSelComp = find(Projectors(iCat).CompMask);
    % Save each component as a projector in the FIF file
    for i = 1:length(iSelComp)
        iComp = iSelComp(i);
        iProj = length(projs) + 1;
        projData = Projectors(iCat).Components(:,iComp)';
        iRows    = find(projData);
        projs(iProj).kind   = FIFF.FIFFV_PROJ_ITEM_FIELD;
        projs(iProj).active = (Projectors(iCat).Status == 2);
        projs(iProj).desc   = Projectors(iCat).Comment;
        projs(iProj).data.nrow      = 1;
        projs(iProj).data.ncol      = length(iRows);
        projs(iProj).data.row_names = [];
        projs(iProj).data.col_names = ChannelNames(iRows);
        projs(iProj).data.data      = projData(iRows);
    end
end

% Write FIF file
fid = fiff_start_file(FifFile);
fiff_write_proj(fid, projs);
fiff_end_file(fid);

author: brainstorm-tools


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This files simplifies the generation of FAIR data
% Input
%    example        the name of the game (for figures etc)
%    fileT        name of a image file for the template  image
%    fileR        name of a image file for the reference image
%    varargin    list of additional variables such as omega, m, LM
%==============================================================================

function expfile = jpgs2data(example,fileT,fileR,varargin)

if nargin == 0,
  help(mfilename)
  runMinimalExample
  return;
end

omega     = [];
m         = [];
LM        = [];
viewPara  = {'viewImage','viewImage2D','colormap','gray(256)'};
imgPara   = {'imgModel','linearInter'};
overwrite = 0;

for k=1:2:length(varargin), % overwrite defaults
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

if isempty(example),
  caller  = dbstack;
  expfile = caller(min(2,length(caller))).name;
  example = expfile(6:end-5);
else
  expfile = sprintf('setup%sData',example);
end
expfile = fullfile(FAIRpath,'temp',expfile);

FAIRmessage(expfile)
if exist([expfile,'.mat'],'file') & ~overwrite,
  fprintf('[%s] %s exists\n',mfilename,expfile);
  return
end;

% do whatever needed to be done to get your data here
image = @(str) double(flipud(imread(str))'); % reads and converts

% load the original data, set domain, initial discretization, and grid
dataT = image(fileT);
dataR = image(fileR);

if isempty(omega), omega = [0,size(dataT,1),0,size(dataT,2)];  end;
if isempty(m),     m     = size(dataR);                        end;

viewImage('reset',viewPara{:});
imgModel('reset',imgPara{:});

omegaI = @(i) omega(min(i,size(omega,1)),:);
xc     = @(k) getCellCenteredGrid(omegaI(k),m);

viewData  = @(I,k) viewImage(imgModel(I,omegaI(k),xc(k)),omegaI(k),m);

FAIRfigure(1,'figname',expfile); clf;
subplot(1,2,1); viewData(dataT,1); title('template');
subplot(1,2,2); viewData(dataR,2); title('reference');

ML = getMultilevel({dataT,dataR},omega,m,'fig',2);

% save to outfile
save(expfile,'example','dataT','dataR','omega','m','ML','viewPara','imgPara');
if ~isempty(LM),
  save(expfile,'-append','LM','example');
end;

%------------------------------------------------------------------------------

function runMinimalExample
setup2DhandData;

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


text:
%
%   simple implementation of N-Dimensional PDF Transfer
%
%   [DR] = pdf_transferND(D0, D1, rotations);
%
%     D0, D1 = NxM matrix containing N-dimensional features
%     rotations = { {R_1}, ... , {R_n} } with R_i PxN
%
%     note that we can use more than N projection axes. In this case P > N
%     and the inverse transformation is done by least mean square.
%     Using more than N axes leads to a more stable (but also slower)
%     convergence.
%
%  (c) F. Pitie 2007
%
%  see reference:
%  Automated colour grading using colour distribution transfer. (2007)
%  Computer Vision and Image Understanding.
%
function [DR] = pdf_transfer(D0, D1, Rotations, varargin)

nb_iterations = length(Rotations);

numvarargs = length(varargin);
if numvarargs > 1
    error('pdf_transfer:TooManyInputs', ...
        'requires at most 1 optional input');
end

optargs = {1};
optargs(1:numvarargs) = varargin;
[relaxation] = optargs{:};

prompt = '';

for it=1:nb_iterations
    fprintf(repmat('\b',[1, length(prompt)]))
    prompt = sprintf('IDT iteration %02d / %02d', it, nb_iterations);
    fprintf(prompt);

    R = Rotations{it};
    nb_projs = size(R,1);

    % apply rotation

    D0R = R * D0;
    D1R = R * D1;
    D0R_ = zeros(size(D0));

    % get the marginals, match them, and apply transformation
    for i=1:nb_projs
        % get the data range
        datamin = min([D0R(i,:) D1R(i,:)])-eps;
        datamax = max([D0R(i,:) D1R(i,:)])+eps;
        u = (0:(300-1))/(300-1)*(datamax - datamin) + datamin;

        % get the projections
        p0R = hist(D0R(i,:), u);
        p1R = hist(D1R(i,:), u);

        % get the transport map
        f = pdf_transfer1D(p0R, p1R);

        % apply the mapping
        D0R_(i,:) = (interp1(u, f', D0R(i,:))-1)/(300-1)*(datamax-datamin) + datamin;
    end

    D0 = relaxation * (R \ (D0R_ - D0R)) + D0;
end

fprintf(repmat('\b',[1, length(prompt)]))

DR = D0;

end

%
% 1D - PDF Transfer
%
function f = pdf_transfer1D(pX,pY)
    nbins = max(size(pX));

    eps = 1e-6; % small damping term that faciliates the inversion

    PX = cumsum(pX + eps);
    PX = PX/PX(end);

    PY = cumsum(pY + eps);
    PY = PY/PY(end);

    % inversion

    f = interp1(PY, 0:nbins-1, PX, 'linear');
    f(PX<=PY(1)) = 0;
    f(PX>=PY(end)) = nbins-1;
    if sum(isnan(f))>0
        error('colour_transfer:pdf_transfer:NaN', ...
              'pdf_transfer has generated NaN values');
    end
end

author: frcs

repo: colour-transfer


text:
function ROI(img, nroi)

% 1) Goal: Draw & process multiple ROIs interactively within an image.
%
% 2) Usage: ROI(img, nroi), where 'img' is your image, and 'nroi' is a
% total number to ROIs to be processed.  The opened image will be processed
% BY DEFAULT.  A prefered procedure is as follows: a) img = imread(...); b)
% imh = imshow(img) and c) ROI(img, nroi).
% Alternatively, if there is no image in your WorkSpace, you MUST use
% square brackets to occupy the argument space for img, such as, ROI([],5),
% will let you open a new image and process with 5 ROIs.
%
% 3) Since getline('closed') is used to get the polygon interactively,
% please click left mouse button to select, and right button to finish up a
% ROI (Backspace to delete the lastest click). You may repeat this process
% till all ROIs were processed.  For more infos, please may see help
% getline.
%
% 4) Results: ROI statistics are displayed on screen or output to a text
% file (optional).
%
% 5) The colors of line/text label are generated by 'jet' colormap,
% therefore, certain color may be too close to tell, especially when you
% select too many ROIs. In that case, you may need to edit the color after
% ROI processing. For more infos, please see help jet.
%
% 6) For your convenience, you may use the following m-scripts to clear
% your enviroments: 'dt', 'dl' and 'df' to delete all text characters,
% lines and figures, respectively.
%
% 7) This package includes 5 m-files: 'ROI.m', 'MultiROI', 'df.m', 'dl.m' and 'dt.m'.
%    MultiROI is similar to ROI. The only difference is that it uses spline
%    interprate to smooth the ROI edge.
%
% Shanrong Zhang
% Department of Radiology
% University of Washington
% 02/09/2004
%
% email: zhangs@u.washington.edu

if nargin == 2

    imh = findobj(0, 'Type', 'Image');

    % if no image opened, open a new one
    if isempty(imh)
        if length(img) == 0
            [infn inpn] = uigetfile('*.*','Please select an image file');
            if infn ~= 0
                img = imread([inpn infn]);
                imh = imshow(img);
                axis image; axis off;
            else
                disp('Cancel by user!')
                return
            end
        else
            imh = imshow(img);
        end
    end

    [nrows, ncols, ncolors] = size(img);

    % Save ROIs to a file (optional)
    SaveIt = questdlg('Do you want to save ROI outputs ?');
    switch SaveIt
        case 'Yes'
            [outfn, outpn] = uiputfile('*', 'Select an output file');
            if outfn == 0
                disp('Cancel by user !');
            else
                fid = fopen([outpn outfn], 'w+');
                fprintf(fid, '%20s\t %-50s\n', 'Date\time = ', datestr(now));
            end
        otherwise
            outfn = 0;
    end

    % generate a jet colormap according to nroi
    cmap = jet(nroi);
    rndp = randperm(nroi);

    hold on;

    croi = 1;

    while croi <= nroi
        [x,y] =getline('closed');

        XData = get(imh, 'XData');
        YData = get(imh, 'YData');
        xmingrid = max( XData(1), floor(min(x)) );
        xmaxgrid = min( XData(2),  ceil(max(x)) );
        ymingrid = max( YData(1), floor(min(y)) );
        ymaxgrid = min( YData(2),  ceil(max(y)) );
        xgrid = xmingrid : xmaxgrid;
        ygrid = ymingrid : ymaxgrid;
        [X, Y] = meshgrid(xgrid, ygrid);
        k_inside = inpolygon(X, Y, x, y);
        Xin = X(k_inside);
        Yin = Y(k_inside);

        cdata = get(imh, 'CData');
        smallcdata = double(cdata(ygrid, xgrid, :));

        roi.index = croi;
        roi.area = polyarea(x,y);
        roi.center =  [mean(Xin(:)), mean(Yin(:))];

        for i=1:ncolors
            roicidata     = smallcdata(:, :, i);
            roi.mean(i)   =   mean(roicidata(k_inside));
            roi.std(i)    =    std(roicidata(k_inside));
            roi.min(i)    =    min(roicidata(k_inside));
            roi.max(i)    =    max(roicidata(k_inside));
            roi.median(i) = median(roicidata(k_inside));
        end;

        plot(x,y,'Color',cmap(rndp(croi), :));
        text(roi.center(1), roi.center(2), num2str(croi), 'Color', cmap(rndp(croi), :), 'FontWeight','Bold');

        % write ROI statistics into file if necessary
        if outfn ~= 0
            fprintf(fid, '\n');
            fprintf(fid, '%20s\t %10.0f\n', 'ROI Index = ', roi.index);

            fprintf(fid, '%20s\t ', 'area = ');
            fprintf(fid, '%10.2f\t', roi.area);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'mean = ');
            fprintf(fid, '%10.2f\t', roi.mean);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'std = ');
            fprintf(fid, '%10.2f\t', roi.std);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'min = ');
            fprintf(fid, '%10.2f\t', roi.min);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'max = ');
            fprintf(fid, '%10.2f\t', roi.max);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'median = ');
            fprintf(fid, '%10.2f\t', roi.median);
            fprintf(fid, '\n');

            fprintf(fid, '%20s\t ', 'roicenter = ');
            fprintf(fid, '%10.2f\t', roi.center);
            fprintf(fid, '\n');
        end

        disp(' ');
        disp(sprintf('%20s\t %10.0f', 'ROI index = ', roi.index));
        disp(sprintf('%20s\t %10.2f', 'area = ', roi.area));
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'mean = ', roi.mean) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'std = ',  roi.std) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'min = ', roi.min) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'max = ', roi.max) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t %10.2f', 'median = ', roi.median) );
        disp(sprintf('%20s\t %10.2f\t %10.2f\t', 'roicenter [x,y] = ', roi.center));
        disp(' ');

        croi = croi + 1;

    end

    if outfn == 0
        disp('Done, but ROI statistics were not saved !!!');
        disp(' ' );
    else
        disp(['Done, ROI statistics been output to ', outfn]);
        disp('But the image with ROI lines/labels was not saved yet !!!');
        fclose(fid);
    end

else
    disp(' ')
    disp('  Number of arguments is incorrect !!!')
    disp(' ')
    help ROI
end

% end of code

author: Sable

repo: mcbench-benchmarks


text:
function world = thickCloister(x,X,y,Y,h,n)

% THICKCLOISTER  Generates features in a 3D cloister shape.
%   THICKCLOISTER(XMIN,XMAX,YMIN,YMAX,H,N) generates a 3D cloister in the
%   limits indicated as parameters, with height H and with N points per
%   side.
%
%   See also CLOISTER.

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

switch nargin
    case {1, 2, 3, 4}
        n = 9;
        h = 1;
    case 5
        n = 9;
end

plane = cloister(x,X,y,Y,n);
nlm = length(plane);
low = zeros(1,nlm);
high  = low + h;
world = [plane plane;low high];

% nid = 2*nlm;

% ids = 1:nid;

% id = 0;
% ids = zeros(1,nid);
% for i = 1:nid
%     id = id+floor(10*rand+1);
%     ids(i) = id;
% end

% world = [ids;world];

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
close all;
%findBestAndWorstMCMCRunByLogPr( 1515157:1515161, 10 );

% plotLogPrVsTime( [1515157 1515159 1515160], 1:10, {'SM+DD one', 'Prior u5', 'SM+DD u5'} )
% ylim( [-5.8 -5.33]*1e4 );
% set( gcf, 'Name', 'LogPrCompare_ValidSamplers' );
% set(gca,'FontSize', 20);
%
% plotLogPrVsTime( [1515157 1515158], 1:10, {'SM+DD one', 'SMnoq+DD one'} )
% ylim( [-5.8 -5.33]*1e4 );
% set( gcf, 'Name', 'LogPrCompare_HastingsFactor_ONE' );
% set(gca,'FontSize', 20);
%
% plotLogPrVsTime( [1515160 1515161], 1:10, {'SM+DD unique5', 'SMnoq+DD unique5'} )
% ylim( [-5.8 -5.33]*1e4 );
% set( gcf, 'Name', 'LogPrCompare_HastingsFactor_Unique5' );
% set(gca,'FontSize', 20);

plotStateSeq( 1515157, 'best' );
set(gcf,'Name', 'StateSeq_SMvalid_ONE_Best' );

plotStateSeq( 1515158, 'best' );
set(gcf,'Name', 'StateSeq_SMnoq_ONE_Best' );

plotStateSeq( 1515160, 'best' );
set(gcf,'Name', 'StateSeq_SMvalid_UNIQUE5_Best' );

plotStateSeq( 1515161, 'best' );
set(gcf,'Name', 'StateSeq_SMnoq_UNIQUE5_Best' );

author: michaelchughes

repo: NPBayesHMM


text:
function a = c8vec_unity ( n )

%*****************************************************************************80
%
%% C8VEC_UNITY returns the N roots of unity.
%
%  Discussion:
%
%    X(1:N) = exp ( 2 * PI * (0:N-1) / N )
%
%    X(1:N)**N = ( (1,0), (1,0), ..., (1,0) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Output, complex A(N), the N roots of unity.
%

%
%  Yes, I realize that in MATLAB "i" already means sqrt(-1).
%  But to me, "i" means an array index, and nothing else.
%
  imaginary = sqrt ( -1.0 );

  for i = 1 : n
    theta = 2 * pi * ( i - 1 ) / n;
    a(i) = cos ( theta ) + imaginary * sin ( theta );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ rexp, sexp ] = poly_q9 ( )

%*****************************************************************************80
%
%% POLY_Q9 returns the monomials associated with a 9 node quadrilateral.
%
%  Reference Element Q9:
%
%    |
%    1  4--7--3
%    |  |     |
%    |  |     |
%    S  8  9  6
%    |  |     |
%    |  |     |
%    0  1--5--2
%    |
%    +--0--R--1-->
%
%  Formula:
%
%    Given coefficients A(I), the polynomial interpolant at (R,S) is
%
%      P(R,S) = sum ( 1 <= I <= N ) A(I) * R**REXP(I) * S**SEXP(I)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer REXP(9), SEXP(9), the powers of R and S associated
%    with each monomial.
%
  rexp(1:9) = [ 0, 0, 1, 0, 1, 2, 1, 2, 2 ];
  sexp(1:9) = [ 0, 1, 0, 2, 1, 0, 2, 1, 2 ];

  return
end


text:
function window_sz = get_search_window( target_sz, im_sz)
% GET_SEARCH_WINDOW

% if(target_sz(1)/target_sz(2) > 2)
%     % For objects with large height, we restrict the search window with padding.height
%     window_sz = floor(target_sz.*[1+padding.height, 1+padding.generic]);
%
% elseif(prod(target_sz)/prod(im_sz(1:2)) > 0.05)
%     % For objects with large height and width and accounting for at least 10 percent of the whole image,
%     % we only search 2x height and width
%     window_sz=floor(target_sz*(1+padding.large));
%
% else
%     %otherwise, we use the padding configuration
%     window_sz = floor(target_sz * (1 + padding.generic));

ratio=target_sz(1)/target_sz(2);
if ratio>1
    window_sz=round(target_sz.*[5,5*ratio]);
else
    window_sz=round(target_sz.*[5/ratio,5]);
end

%window_sz=round(target_sz.*[5,9]);
window_sz=window_sz-mod(window_sz,2)+1;

end

repo: CREST-Release


text:
function neqn = p01_neqn ( )

%*****************************************************************************80
%
%% P01_NEQN returns the number of equations for problem P01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 1;

  return
end


text:
function varargout = chebyball(F)
%CHEBYBALL Computes Chebyshev ball of a constraint object
%
% If two outputs are requested, the numerical data is returned
% [xc,R] = chebyball(F)
%
% If three outputs are requrested, the symbolic model (x-xc)'(x-xc)<R^2 is
% appended
% [xc,R,C] = chebyball(F)
%
% If only one output is requested, only the symbolic constraint is returned
% C = chebyball(F)

switch nargout
    case 0
        f = chebyball(lmi(F))
    case 1
        [f] = chebyball(lmi(F));
        varargout{1} = f;
    case 2
        [xc,R] = chebyball(lmi(F));
        varargout{1} = xc;
        varargout{2} = R;
    case 3
        [xc,R,f] = chebyball(lmi(F));
        varargout{1} = xc;
        varargout{2} = R;
        varargout{3} = f;
    otherwise
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%----------------------------------------------------------------
% Depth Map Super-Resolution by Deep Multi-Scale Guidance, ECCV16
%                     written by T.-W. HUI
%----------------------------------------------------------------
function Z = normalizeIm(Z, range)

if nargin < 2
    range = [0.01 1];
end

validMap = Z > 0;
Z = Z - min(Z(validMap));
Z = Z/max(Z(validMap));
Z(~validMap) = 0;
Z(validMap) = Z(validMap) + range(1);
Z = range(2)*Z/max(max(Z));

end

author: twhui

repo: MSG-Net

save_path: github-repos/MATLAB/twhui-MSG-Net


text:
function dms=rad2dms(rad)
% RAD2DMS  Converts radians to degrees-minutes-seconds. Vectorized.
% Version: 12 Mar 00
% Useage:  dms=rad2dms(rad)
% Input:   rad - vector of angles in radians
% Output:  dms - [d m s] array of angles in deg-min-sec, where
%                d = vector of degrees
%                m = vector of minutes
%                s = vector of seconds

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

d=abs(rad).*180/pi;
id=floor(d);
rm=(d-id).*60;
im=floor(rm);
s=(rm-im).*60;

%if rad<0
%  if id==0
%    if im==0
%      s = -s;
%    else
%      im = -im;
%    end
%  else
%    id = -id;
%  end
%end

ind=(rad<0 & id~=0);
id(ind)=-id(ind);

ind=(rad<0 & id==0 & im~=0);
im(ind)=-im(ind);

ind=(rad<0 & id==0 & im==0);
s(ind)=-s(ind);

dms=[id im s];


text:
function test_suite = test_hierprior

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_HIERPRIOR
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

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

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testDerivatesRegression(testCase)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{1}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01, 'AbsTol', 0.4)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{2}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01, 'AbsTol', 0.4)
end

function testDerivatesLaplace(testCase)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{3}, ...
    'AbsTol', 1e-5)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{4}, ...
    'AbsTol', 1e-5)
end

function testDerivatesEP(testCase)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{5}, ...
    'AbsTol', 1e-5)
  verifyVarsEqual(testCase, getName(), {'derivs'}, @(x)x{6}, ...
    'AbsTol', 1e-5)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [angle_SAM,map] = SAM(ref,tar)
%--------------------------------------------------------------------------
% Spectral angle mapper (SAM)
%
% USAGE
%   [angle_SAM,map] = SAM(msoriginal,msfused)
%
% INPUT
%   ref : reference HS data (rows,cols,bands)
%   tar : target HS data (rows,cols,bands)
%
% OUTPUT
%   angle_SAM : average SAM in degree (scalar)
%   map       : 2-D map of SAM
%
%--------------------------------------------------------------------------
[rows,cols,bands] = size(tar);
prod_scal = dot(ref,tar,3);
norm_orig = dot(ref,ref,3);
norm_fusa = dot(tar,tar,3);
prod_norm = sqrt(norm_orig.*norm_fusa);
prod_map = prod_norm;
prod_map(prod_map==0)=eps;
map = acos(prod_scal./prod_map);
prod_scal = reshape(prod_scal,rows*cols,1);
prod_norm = reshape(prod_norm, rows*cols,1);
z=find(prod_norm==0);
prod_scal(z)=[];prod_norm(z)=[];
angolo = sum(sum(acos(prod_scal./prod_norm)))/(size(prod_norm,1));
angle_SAM = real(angolo)*180/pi;

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function plaza = create_plaza(B, L)
global plazalength;
topgap = 5;
bottomgap = 1;
plaza = zeros(plazalength,B+2);
if mod(B-L,2)==0
for row = 1:plazalength
plaza(row,1) = -888;
plaza(row,2+B) = -888;
end
for col = 2:B/2 - L/2 + 1
for row = 1:(plazalength-1)/2 - topgap * (col-1)
plaza(row,col) = -888;
plaza(row,B+3-col) = -888;
end
end
for col = 2:B/2 - L/2 + 1
for row = (plazalength+3)/2 + bottomgap*(col-1):plazalength
plaza(row,col) = -888;
plaza(row,B+3-col) = -888;
end
end
else
for row = 1:plazalength
plaza(row,1) = -888;
plaza(row,3+B) = -888;
end
for col = 2:(B+1)/2 - L/2 + 1
for row = 1:(plazalength-1)/2 - topgap * (col-1)
plaza(row,col) = -888;
plaza(row,B+4-col) = -888;
end
end
for col = 2:(B+1)/2 - L/2 + 1
for row = (plazalength+3)/2 + bottomgap*(col-1):plazalength
plaza(row,col) = -888;
plaza(row,B+4-col) = -888;
end
end
for row = 1:plazalength
plaza(row,2+B) = -888;
end
end

author: Eurus-Holmes

path: github-repos/MATLAB/Eurus-Holmes-Mathematical_Modeling/Mathematical_Modeling-cb9dd53af84ffbd455ec62ab89886351e0ec98d9/Cellular automaton/元胞自动机/1月13日课件（元胞自动机）/程序/The Booth Tolls for Thee/create_plaza.m


text:
function tagS = ROI_contour_module_tags
%"ROI_contour_module_tags"
%   Return the tags used to represent a ROI contour specified by
%   C.8.8.6 in PS3.3 of 2006 DICOM specification. Tags are returned in a
%   struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array
%          containing the tags of child fields.
%
%JRA 06/06/06
%
%Usage:
%   tagS = ROI_contour_module_tags
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
template = tagS;

%Add tags based on PS3.3 attribute lists.

%ROI Contour Sequence
tagS(end+1) = struct('tag', '30060039', 'tagdec', 805699641, 'type', '1', 'children', []);
child_1     = template;

    %Referenced ROI Number
    child_1(end+1) = struct('tag', '30060084', 'tagdec', 805699716, 'type', '1', 'children', []);

    %ROI Display Color
    child_1(end+1) = struct('tag', '3006002A', 'tagdec', 805699626, 'type', '3', 'children', []);

    %Contour Sequence
    child_1(end+1) = struct('tag', '30060040', 'tagdec', 805699648, 'type', '3', 'children', []);
    child_2        = template;

        %Contour Number
        child_2(end+1) = struct('tag', '30060048', 'tagdec', 805699656, 'type', '3', 'children', []);

        %Attached Contours
        child_2(end+1) = struct('tag', '30060049', 'tagdec', 805699657, 'type', '3', 'children', []);

        %Contour Image Sequence
        child_2(end+1) = struct('tag', '30060016', 'tagdec', 805699606, 'type', '3', 'children', []);
        child_3        = template;

            %Image SOP Instance Reference Macro
            child_3 = image_SOP_instance_reference_macro_tags;
            child_2(end).children = child_3;

        %Contour Geometric Type
        child_2(end+1) = struct('tag', '30060042', 'tagdec', 805699650, 'type', '1C', 'children', []);

        %Contour Slab Thickness
        child_2(end+1) = struct('tag', '30060044', 'tagdec', 805699652, 'type', '3', 'children', []);

        %Contour Offset Vector
        child_2(end+1) = struct('tag', '30060045', 'tagdec', 805699653, 'type', '3', 'children', []);

        %Number of Contour Points
        child_2(end+1) = struct('tag', '30060046', 'tagdec', 805699654, 'type', '1C', 'children', []);

        %Contour Data
        child_2(end+1) = struct('tag', '30060050', 'tagdec', 805699664, 'type', '1C', 'children', []);

        child_1(end).children = child_2;

   tagS(end).children = child_1;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = acosh(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('positive','callback');
        operator.convexity = @convexity;
        operator.bounds = @bounds;
        operator.domain = [1 inf];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

function [L,U] = bounds(xL,xU)
if xU<=-1
    L = acosh(xU);
    U = acosh(xL);
elseif xL>=1
    L = acosh(xL);
    U = acosh(xU);
elseif xU <= 1
    L = 0;
    U = acosh(xL);
else
    L = 0;
    U = acosh(xU);
end

function vexity = convexity(xL,xU)
if xL >= 1
    vexity = 'concave';
elseif xU <= -1
    vexity = 'concave';
else
    vexity = 'none';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function t = readsqdinfo(t,fileinfo,acq_offset,raw_offset)
% t = readsqdinfo(t,fileinfo,acq_offset,raw_offset); Read info from sqd-file
% and write it to ACQPARAM object

SEEK_SET=-1;    % Start of the file for fseek: -1 => at the beginning

if ischar(fileinfo) % If filename is passed
    % open file as binary and little endian
    fid = fopen(fileinfo,'rb','l');
    if fid==-1,
        disp('File not found');
        [f,p] = uigetfile('*.sqd','MEG160 files');
        if ((~isstr(f)) | ~min(size(f))), return; end;
        fname = [p f];
        fid = fopen(fname,'rb','l');
    end;
elseif isnumeric(fileinfo)% If filepointer is passed
    if fileinfo>0
        fid = fileinfo;
    else
        error('Invalid fid');
    end;
end;

% Initialize the offsets
switch nargin
case 2
    % Get offset of acquisition parameters
    fseek( fid, 128, SEEK_SET );
    acq_offset = fread(fid,1,'long');
    % Get Acquisition type
    fseek(fid,acq_offset,SEEK_SET);
    t.AcquisitionType       = fread(fid,1,'long');
    % Get Raw Offset
    fseek(fid,t.Private.Offset2Raw(t.AcquisitionType),SEEK_SET);
    raw_offset  = fread(fid,1,'long');
case 3
    % Get Acquisition type
    fseek(fid,acq_offset,SEEK_SET);
    t.AcquisitionType       = fread(fid,1,'long');
    % Get Raw Offset
    fseek(fid,t.Private.Offset2Raw(t.AcquisitionType),SEEK_SET);
    raw_offset  = fread(fid,1,'long');
end;
%Save Offsets
t.RawOffset         = raw_offset;

% Read acquisition type
fseek( fid, acq_offset, SEEK_SET );
t.AcquisitionType       = fread(fid,1,'long');
t.SampleRate            = fread(fid,1,'double');
t.Datatype              = t.Private.Datatypes{t.AcquisitionType};

switch t.AcquisitionType
    case 1
        % Read acquisition parameters
        t.SampleInfo.SamplesAcquired    = fread(fid,1,'long');
        t.SampleInfo.ActSamplesAcquired = fread(fid,1,'long');
        t.SampleInfo.SamplesAvailable   = t.SampleInfo.ActSamplesAcquired;
    case {2,3}
        % Read acquisition parameters
        t.SampleInfo.FrameLength        = fread(fid,1,'long');
        t.SampleInfo.SamplesAvailable   = t.SampleInfo.FrameLength;
        t.SampleInfo.PreTrigger         = fread(fid,1,'long');
        t.SampleInfo.AverageCount       = fread(fid,1,'long');
        t.SampleInfo.ActAverageCount    = fread(fid,1,'long');
end;
if ischar(fileinfo);fclose(fid);end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Demo for DAVB Van Der Pol oscillator.
% This demo inverts a model of nonlinear Van Dr Pol oscillator, which is
% observed through a nonlinear sigmoid function.

clear variables
close all

% Choose basic settings for simulations
n_t = 2e2;
deltat = 1e-1;
f_fname = @f_vanDerPol;
g_fname = @g_sigmoid;
alpha   = 1e2;
sigma   = 1e1;
theta   = [1e0];
phi     = [];
u       = [];

% Build options structure for temporal integration of SDE
inG.scale = 50;
inG.slope = 5;
inF.deltat = deltat;
options.inF     = inF;
options.inG     = inG;
options.backwardLag = 5;

% Build priors for model inversion
priors.muX0 = 1e-1*ones(2,1);
priors.SigmaX0 = 1e0*eye(2);
priors.muTheta = 0.*ones(1,1);
priors.SigmaTheta = 1e0*eye(1);
priors.a_alpha = 1e0;
priors.b_alpha = 1e0;
priors.a_sigma = 1e0;
priors.b_sigma = 1e0;

% Build options and dim structures for model inversion
options.priors      = priors;
dim.n_theta         = 1;
dim.n_phi           = 0;
dim.n               = 2;

% Build time series of hidden states and observations
[y,x,x0,eta,e] = VBA_simulate (n_t,f_fname,g_fname,theta,phi,u,alpha,sigma,options);

% display time series of hidden states and observations
displaySimulations(y,x,eta,e);
% disp('--paused--')
% pause

% [posterior,out] = VBA_onlineWrapper(y,u,f_fname,g_fname,dim,options);

% Call inversion routine
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% Display results
displayResults(posterior,out,y,x,x0,theta,phi,alpha,sigma);

% Make predictions
try
    options = out.options;
    [xs,ys,xhat,vx,yhat,vy] = VBA_comparePredictions(...
        n_t,theta,phi,u,alpha,sigma,options,posterior,dim);
catch
    disp('------!!Unable to form predictions!!------')
end

author: MBB-team

repo: VBA-toolbox


text:
function D = spm_mesh_distmtx(M,order)
% Compute the distance matrix of a triangle mesh
% FORMAT D = spm_mesh_distmtx(M,order)
% M        - patch structure
% order    - 0: adjacency matrix
%            1: first order distance matrix [default]
%            2: second order distance matrix
%
% D        - distance matrix
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_distmtx.m 4035 2010-08-05 18:54:32Z guillaume $

if nargin < 2, order = 1; end

if order > 2, error('High order distance matrix not handled.'); end

%-Adjacency matrix
%--------------------------------------------------------------------------
if order == 0
    D = spm_mesh_adjacency(M);
    return;
end

%-First order distance matrix
%--------------------------------------------------------------------------
if isstruct(M) && ~isa(M.faces,'double')
    M.faces = double(M.faces);
elseif isa(M,'gifti')
    M = export(M,'patch');
    M.faces = double(M.faces);
end
d = M.vertices(M.faces(:,[1 2 3]),:) - M.vertices(M.faces(:,[2 3 1]),:);

D = sparse([M.faces(:,1); M.faces(:,2); M.faces(:,3)], ...
           [M.faces(:,2); M.faces(:,3); M.faces(:,1)], sqrt(sum(d.^2, 2)));

D = (D + D')/2;

if order == 1, return; end

%-Second order distance matrix
%--------------------------------------------------------------------------
D2 = D;
for i=1:size(M.vertices,1)
    a      = find(D2(i,:));
    [b,c]  = find(D2(a,:));
    D(i,c) = D(i,a(b)) + diag(D(a(b),c))';
    D(c,i) = D(i,c)';
end

author: spm


text:
function p = PlotMean(X,Y,L,U,style,color)

%PlotMean - Plot mean and confidence intervals.
%
%  USAGE
%
%    p = PlotMean(X,Y,L,U,style,color)
%
%    X              abscissae
%    Y              ordinates
%    L              lower error level
%    U              upper error level
%    style          optional style (':' or '-')
%    color          optional color specification (see plot)
%
%  NOTES
%
%    Notice that L and U specifiy actual error levels rather than (relative)
%    ranges, contrary to the built-in 'errorbar' where the error is in [Y-L Y+U].
%
%    If X is empty, it is set to 1:N (where N is the length of Y). If Y is
%    empty, only error limits are plotted. Conversely, if L and U are empty, only
%    the mean is plotted.

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% 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.

if nargin < 4,
  error('Incorrect number of parameters (type ''help PlotMean'' for details).');
end

if nargin < 6,
	color = 'b';
end
if nargin < 5,
	style = ':';
end

nX = length(X);
nY = length(Y);
nL = length(L);
nU = length(U);

% Make sure all non-empty inputs have the same length
n = [nX nY nL nU];
n = n(n~=0);
if isempty(n),
	error('At least one input must contain data (type ''help PlotMean'' for details).');
end
m = n - n(1);
if any(m),
	error('All non-empty inputs must have the same length (type ''help PlotMean'' for details).');
end

if nX == 0,
	X = (1:n(1))';
end

hold on;
% Plot mean
if nY ~= 0,
	p = plot(X,Y);
	if strcmp(style,':'),
		set(p,'LineStyle','-','color',color);
	else
		set(p,'LineWidth',3,'color',color);
	end
end
% Plot lower error level
if nL ~= 0,
	p = plot(X,L);
	if strcmp(style,':'),
		set(p,'LineStyle',':','color',color);
	else
		set(p,'LineWidth',1,'color',color);
	end
end
% Plot upper error level
if nU ~= 0,
	p = plot(X,U);
	if strcmp(style,':'),
		set(p,'LineStyle',':','color',color);
	else
		set(p,'LineWidth',1,'color',color);
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ m, p, t ] = naca4_mpt ( code )

%*****************************************************************************80
%
%% NACA4_MPT returns the parameters stored in a NACA 4 digit airfoil code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer CODE, the NACA4 code.
%    0 <= CODE <= 9999.
%
%    Output, real M, the maximum camber, as a percentage of the chord length.
%    0 <= M <= 1.0.
%
%    Output, real P, the relative distance of the occurrence of the maximum
%    camber from the beginning of the chord.
%    0 <= P <= 1.0.
%
%    Output, real T, the maximum thickness relative to the chord length.
%    0 <= T <= 1.0.
%
  if ( code < 0 || 9999 < code )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NACA4_MPT - Fatal error!\n' );
    fprintf ( 1, '  CODE should be an integer between 0 and 9999.\n' );
    error ( 'NACA4_MPT - Fatal error!' );
  end

  m = floor ( code / 1000 );
  code = code - m * 1000;
  m = m / 100.0;

  p = floor ( code / 100 );
  code = code - p * 100;
  p = p / 10.0;

  t = code / 100.0;

  return
end


text:
function [out_profile,out_IMU_bias_est,out_KF_SD,out_R_matrix,out_Q_matrix,corrections] =...
    Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, LC_KF_config, est_IMU_bias)
%Loosely_coupled_INS_GNSS - Simulates inertial navigation using ECEF
% navigation equations and kinematic model, GNSS using a least-squares
% positioning algorithm, and loosely-coupled INS/GNSS integration.
%
% Software for use with "Principles of GNSS, Inertial, and Multisensor
% Integrated Navigation Systems," Second Edition.
%
% This function created 12/4/2012 by Paul Groves
%
% Significant changes made to data-formatting and time-handling by Adam
% Werries, 2015/2016.
%
% Inputs:
%   LC_KF_config
%     .init_att_unc           Initial attitude uncertainty per axis (rad)
%     .init_vel_unc           Initial velocity uncertainty per axis (m/s)
%     .init_pos_unc           Initial position uncertainty per axis (m)
%     .init_b_a_unc           Initial accel. bias uncertainty (m/s^2)
%     .init_b_g_unc           Initial gyro. bias uncertainty (rad/s)
%     .gyro_noise_PSD         Gyro noise PSD (rad^2/s)
%     .accel_noise_PSD        Accelerometer noise PSD (m^2 s^-3)
%     .accel_bias_PSD         Accelerometer bias random walk PSD (m^2 s^-5)
%     .gyro_bias_PSD          Gyro bias random walk PSD (rad^2 s^-3)
%     .pos_meas_SD            Position measurement noise SD per axis (m)
%     .vel_meas_SD            Velocity measurement noise SD per axis (m/s)
%     .init_biases            IMU bias initialization
%
% Outputs:
%   out_profile        Navigation solution as a motion profile array
%   out_IMU_bias_est   Kalman filter IMU bias estimate array
%   out_clock          GNSS Receiver clock estimate array
%   out_KF_SD          Output Kalman filter state uncertainties
%
% Format of motion profiles:
%  Column 1: time (sec)
%  Column 2: latitude (rad)
%  Column 3: longitude (rad)
%  Column 4: height (m)
%  Column 5: north velocity (m/s)
%  Column 6: east velocity (m/s)
%  Column 7: down velocity (m/s)
%  Column 8: roll angle of body w.r.t NED (rad)
%  Column 9: pitch angle of body w.r.t NED (rad)
%  Column 10: yaw angle of body w.r.t NED (rad)
%
% Format of output IMU biases array:
%  Column 1: time (sec)
%  Column 2: estimated X accelerometer bias (m/s^2)
%  Column 3: estimated Y accelerometer bias (m/s^2)
%  Column 4: estimated Z accelerometer bias (m/s^2)
%  Column 5: estimated X gyro bias (rad/s)
%  Column 6: estimated Y gyro bias (rad/s)
%  Column 7: estimated Z gyro bias (rad/s)
%
% Format of KF state uncertainties array:
%  Column 1: time (sec)
%  Column 2: X attitude error uncertainty (rad)
%  Column 3: Y attitude error uncertainty (rad)
%  Column 4: Z attitude error uncertainty (rad)
%  Column 5: X velocity error uncertainty (m/s)
%  Column 6: Y velocity error uncertainty (m/s)
%  Column 7: Z velocity error uncertainty (m/s)
%  Column 8: X position error uncertainty (m)
%  Column 9: Y position error uncertainty (m)
%  Column 10: Z position error uncertainty (m)
%  Column 11: X accelerometer bias uncertainty (m/s^2)
%  Column 12: Y accelerometer bias uncertainty (m/s^2)
%  Column 13: Z accelerometer bias uncertainty (m/s^2)
%  Column 14: X gyro bias uncertainty (rad/s)
%  Column 15: Y gyro bias uncertainty (rad/s)
%  Column 16: Z gyro bias uncertainty (rad/s)

old_time = filter_time(1);

% Initialize INS from GPS
old_est_r_eb_e = init_cond(1:3)';
old_est_v_eb_e = init_cond(4:6)';
[old_est_L_b,old_est_lambda_b,old_est_h_b,old_est_v_eb_n] = pv_ECEF_to_NED(old_est_r_eb_e,old_est_v_eb_e);
est_L_b = old_est_L_b;

% Initialize estimated attitude solution
% old_est_C_b_n = Initialize_NED_attitude(true_C_b_n,initialization_errors);
old_est_C_b_n = Euler_to_CTM(init_cond(7:9));
[~,~,old_est_C_b_e] = NED_to_ECEF(old_est_L_b,old_est_lambda_b,old_est_h_b,old_est_v_eb_n,old_est_C_b_n);
old_est_r_eb_e = old_est_r_eb_e + old_est_C_b_e*LC_KF_config.lever_arm + old_est_C_b_e*LC_KF_config.gps_correction;

% Initialize output profile record and errors record
out_profile = zeros(length(filter_time),10);

% Generate output profile record
out_profile(1,1) = old_time;
out_profile(1,2:4) = old_est_r_eb_e;
out_profile(1,5:7) = old_est_v_eb_e';
out_profile(1,8:10) = CTM_to_Euler(old_est_C_b_n')';

% Initialize Kalman filter P matrix and IMU bias states
P_matrix = Initialize_LC_P_matrix(LC_KF_config);

% Generate IMU bias and clock output records
out_IMU_bias_est(1,1) = old_time;
out_IMU_bias_est(1,2:7) = est_IMU_bias';

% Generate KF uncertainty record
out_KF_SD(1,1) = old_time;
for n = 1:15
    out_KF_SD(1,n+1) = sqrt(P_matrix(n,n));
end % for i

% Initialize R (not really used but it makes me feel better
pos_variance = LC_KF_config.gps_pos_stddev.^2*gps(1,16).^2.*ones(1,3);
vel_variance = LC_KF_config.gps_vel_stddev.^2*gps(1,16).^2.*ones(1,3);
R_matrix(1:3,1:3) = diag(max(LC_KF_config.pos_sd_min, ...
                         min(LC_KF_config.pos_sd_max,...
                             pos_variance)));
R_matrix(1:3,4:6) = zeros(3);
R_matrix(4:6,1:3) = zeros(3);
R_matrix(4:6,4:6) = diag(max(LC_KF_config.vel_sd_min,...
                         min(LC_KF_config.vel_sd_max,...
                         vel_variance)));

out_R_matrix = zeros(size(gps,1), 6);
for n = 1:6
    out_R_matrix(1,n) = R_matrix(n,n);
end % for i
out_Q_matrix = zeros(size(gps,1), 15);
Q_matrix = zeros(15,15);
for n = 1:15
    out_Q_matrix(1,n) = Q_matrix(n,n);
end % for i
% Main loop
GNSS_epoch = 2;
last_GNSS_epoch = GNSS_epoch;
corrections = zeros(15,size(gps,1)-1);
last_imu_index = 1;
numIMU = size(imu, 1);
imu_index_windowsize = floor(20*epoch/mean(diff(imu(:,1))));
for i = 2:length(filter_time)
    time = filter_time(i);
    % find range of imu measurements to use
    endcap = min(numIMU, last_imu_index+1 + imu_index_windowsize);
    indices = find(imu(last_imu_index+1:endcap,1) < time);
    imu_range_end = last_imu_index + indices(end) - 1;
    % Apply IMU bias estimates
    meas_f_ib_b = mean(imu(last_imu_index+1:imu_range_end,2:4))' - est_IMU_bias(1:3);
    meas_omega_ib_b = mean(imu(last_imu_index+1:imu_range_end,5:7))' - est_IMU_bias(4:6);

    % Update estimated navigation solution
    [est_r_eb_e,est_v_eb_e,est_C_b_e] = Nav_equations_ECEF(epoch,...
        old_est_r_eb_e,old_est_v_eb_e,old_est_C_b_e,meas_f_ib_b,...
        meas_omega_ib_b);
%     while GNSS_epoch < size(gps,1)+1 && ~gps(GNSS_epoch,3)
%         GNSS_epoch = GNSS_epoch + 1;
%     end
    % Determine whether to update GNSS simulation and run Kalman filter
    if GNSS_epoch <= size(gps,1) && time > gps(GNSS_epoch,1)
        tor_s = gps(GNSS_epoch,1) - gps(last_GNSS_epoch,1);  % KF time interval
        GNSS_r_eb_e = gps(GNSS_epoch,9:11)';
        GNSS_v_eb_e = gps(GNSS_epoch,12:14)';
        est_L_b = lla(GNSS_epoch,1);
        % Use the GPS-reported standard deviation values, but clamping
        % min/max values according to configuration
        pos_variance = LC_KF_config.gps_pos_stddev.^2*gps(GNSS_epoch,16).^2;
        vel_variance = LC_KF_config.gps_vel_stddev.^2*gps(GNSS_epoch,16).^2;
        R_matrix(1:3,1:3) = diag(max(LC_KF_config.pos_sd_min, ...
                                     min(LC_KF_config.pos_sd_max,...
                                         pos_variance.*ones(1,3)*tor_s)));
        R_matrix(4:6,4:6) = diag(max(LC_KF_config.vel_sd_min,...
                                     min(LC_KF_config.vel_sd_max,...
                                         vel_variance.*ones(1,3)*tor_s)));
%         disp(size(R_matrix));
%         disp(diag(R_matrix));

        % Run Integration Kalman filter
        [est_C_b_e,est_v_eb_e,est_r_eb_e,est_IMU_bias,P_matrix_new,corrections(:,GNSS_epoch-1), Phi_matrix, Q_matrix] =...
            LC_KF_Epoch(GNSS_epoch, GNSS_r_eb_e,GNSS_v_eb_e,tor_s,est_C_b_e,...
            est_v_eb_e,est_r_eb_e,est_IMU_bias,P_matrix,meas_f_ib_b,...
            est_L_b,LC_KF_config,Q_matrix,R_matrix,meas_omega_ib_b);
        if any(any(isnan(P_matrix))) || any(any(isinf(P_matrix)))
            disp('Filter instability detected. Time to kick the bucket.');
            fprintf('Died on iteration %d/%d\n',i,length(filter_time));
            fprintf('and GNSS epoch %d/%d\n',GNSS_epoch,size(gps,1));
            return
        end

        % Run adaptive algorithm
        [Q_matrix] = adapt_noise_covariance(Phi_matrix, P_matrix_new, P_matrix, Q_matrix, ...
                                            LC_KF_config.n, GNSS_epoch, corrections);
        P_matrix = P_matrix_new;

        % Generate IMU bias and clock output records
        out_IMU_bias_est(GNSS_epoch,1) = time;
        out_IMU_bias_est(GNSS_epoch,2:7) = est_IMU_bias';

        % Generate KF uncertainty output record
        out_KF_SD(GNSS_epoch,1) = time;
        for n = 1:15
            out_KF_SD(GNSS_epoch,n+1) = sqrt(P_matrix(n,n));
        end % for i
        for n = 1:6
            out_R_matrix(GNSS_epoch,n) = sqrt(abs(R_matrix(n,n)));
        end % for i
        for n = 1:15
            out_Q_matrix(GNSS_epoch,n) = sqrt(abs(Q_matrix(n,n)));
        end % for i
        last_GNSS_epoch = GNSS_epoch;
        GNSS_epoch = GNSS_epoch + 1;
    end % if time

    % Convert navigation solution to NED
    [~,~,~,~,est_C_b_n] =...
        ECEF_to_NED(est_r_eb_e,est_v_eb_e,est_C_b_e);

    % Generate output profile record
    out_profile(i,1) = time;
    out_profile(i,2:4) = est_r_eb_e;
    out_profile(i,5:7) = est_v_eb_e';
    out_profile(i,8:10) = CTM_to_Euler(est_C_b_n')';

    % Reset old values
    old_est_r_eb_e = est_r_eb_e;
    old_est_v_eb_e = est_v_eb_e;
    old_est_C_b_e = est_C_b_e;
    last_imu_index = imu_range_end;
end %epoch

% Ends

author: awerries

repo: kalman-localization


text:
function [epo, ev]= proc_selectClasses(epo, cls, varargin)
%PROC_SELECTCLASSES - select only a subset of the avialble classes
%
%Synopsis:
%[epo, ev]= proc_selectClasses(epo, classes)
%
% IN   epo     - epoched data structure
%      classes - cell array of classes names or vector of class indices
%
% OUT  epo     - epoch structure containing only selected events
%      ev      - indices of select events
%
%Description:
% This function selects events from a epoch structure that belong
% to given classes. a typical application is to select from a multi-class
% experiment a two-class subproblem.
% Class names may include regular expressions, see regexpi. This also
% implies that you might have to escape some characters, see last example.
%
% Examples
%      epo_lr= proc_selectClasses(epo, {'left', 'right'});
%      epo_12= proc_selectClasses(epo, [1 2]);
%      epo_nl= proc_selectClasses(epo, {'not', 'left*'});
%
%      epo= proc_selectClasses(epo, 'r\^2 \( a , b \)';

% Benjamin Blankertz

props_selectEvents= proc_selectEpochs;

if nargin==0,
  epo= props_selectEvents; return
end

misc_checkType(epo, 'STRUCT(className y)');
epo= misc_history(epo);

opt= opt_proplistToStruct(varargin{:});

clInd= procutil_getClassIndices(epo, cls);
ev= find(any(epo.y(clInd,:)==1,1));

%% the following is done to keep the order of classes as specified
epo.y= epo.y(clInd,:);
epo.className= epo.className(clInd);

%% the following is to avoid an error in proc_selectEpochs
if length(ev)==size(epo.y,2)
  return;
end

epo= proc_selectEpochs(epo, ev, opt);

author: bbci


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function newPose = transPoseFwd(oldPose,T)
% newPose = transPoseFwd(oldPose,T)
% produce new pose according to oldPose and T
% size(oldPose) = [m,2n], each row must in the format of X1...Xn,Y1...Yn.
% T is a array containing n structures. Each one is the trans for one
% sample.
% newPose has the same size as oldPose.

m = size(oldPose,1);
if length(T)~=m,error('length(T) must be the same as size(oldPose,1)!');end;
n = size(oldPose,2);
if mod(n,2),error('size(pose,2) is odd!');end
n = n / 2;
ind_nan = find(isnan(oldPose));
oldPose(ind_nan) = 0;

if (m<100) || (matlabpool('size')==0)
    newPose = arrayfun(@(i)reshape(tformfwd(T(i),reshape(oldPose(i,:),n,2)),1,2*n),1:m,'UniformOutput',false);
    newPose = cell2mat(newPose(:));
    if m>=100
        warning('Please launch matlabpool to speed up your program!');
    end
else
    newPose = zeros(size(oldPose));
    parfor i = 1:m
        newPose(i,:) = reshape(tformfwd(T(i),reshape(oldPose(i,:),n,2)),1,2*n);
    end
end
newPose(ind_nan) = NaN;

end

author: zhusz


function G = get_payoff_G_matrix_from_ygrid_2d( y_1, y_2, S_0s, sigmas, rho, contractParams)
%UNTITLED5 Summary of this function goes here
%   Detailed explanation goes here

payoff_type = contractParams.payoff_type;

if payoff_type == 1  % G = S_1
    payoff = @(y1,y2)S_0s(1)*exp(sigmas(1)*y1);

elseif payoff_type == 2  % G = S_2
    payoff = @(y1,y2)S_0s(2)*exp(sigmas(2)*(y2 + rho*y1));

elseif payoff_type == 3  % Exchange:  G = (S_1 - S_2)^+
    payoff = @(y1,y2) max(0, S_0s(1)*exp(sigmas(1)*y1) - S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)));

elseif payoff_type == 4  % Spread:  G = (S_1 - S_2 - K)^+
    K = contractParams.K;
    payoff = @(y1,y2) max(0, S_0s(1)*exp(sigmas(1)*y1) - S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)) - K);

elseif payoff_type == 5  % Geometric Basket Call / Put:  G = (sqrt(S_1) * sqrt(S_2) - K)^+ and
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2) max(0, sqrt(S_0s(1)*exp(sigmas(1)*y1)) * sqrt(S_0s(2)*exp(sigmas(2)*(y2 + rho*y1))) - K);
    else
        payoff = @(y1,y2) max(0, K - sqrt(S_0s(1)*exp(sigmas(1)*y1)) * sqrt(S_0s(2)*exp(sigmas(2)*(y2 + rho*y1))));
    end

elseif payoff_type == 6  % Arithmetic Basket Call / Put:  G = (sqrt(S_1) * sqrt(S_2) - K)^+
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2) max(0, 0.5*S_0s(1)*exp(sigmas(1)*y1) + 0.5*S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)) - K);
    else
        payoff = @(y1,y2) max(0, K - 0.5*S_0s(1)*exp(sigmas(1)*y1) - 0.5*S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)));
    end

elseif payoff_type == 7  % Call-on-Max and Put-on-Min
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2) max(0, max(S_0s(1)*exp(sigmas(1)*y1), S_0s(2)*exp(sigmas(2)*(y2 + rho*y1))) - K);
    else
        payoff = @(y1,y2) max(0, K - min(S_0s(1)*exp(sigmas(1)*y1), S_0s(2)*exp(sigmas(2)*(y2 + rho*y1))));
    end
elseif payoff_type == 8  % Call/put on Just S_2
    K = contractParams.K;
    if contractParams.call == 1
        payoff = @(y1,y2) max(0, S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)) - K);
    else
        payoff = @(y1,y2) max(0, K - S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)));
    end
elseif payoff_type == 9  % Best-of / worst of
    if contractParams.best == 1
        payoff = @(y1,y2) max(S_0s(1)*exp(sigmas(1)*y1), S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)));
    else  % worst of
        payoff = @(y1,y2) min(S_0s(1)*exp(sigmas(1)*y1), S_0s(2)*exp(sigmas(2)*(y2 + rho*y1)));
    end
end

m_0 = length(y_1);
G = zeros(m_0, m_0);

for i=1:m_0
    for j=1:m_0
        G(i,j) = payoff(y_1(i), y_2(j));
    end
end

end


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

N = size(X,2);                  % number of features

for i=1:N,
    feature = X(:,i);           % get feature
    mu(i) = mean(feature);      % mean
    sigma(i) = std(feature);    % standard deviation

    % normalize
    X_norm(:,i) = (feature - mu(i)) / sigma(i);
end

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

end

author: rmarquis

repo: coursera-machinelearning


text:
%% Configurable Simulink Model for DC-DC Converters
%

%% DC-DC Converters
% There are three kinds of switching mode DC-DC converters, buck, boost and
% buck-boost. The buck mode is used to reduce output voltage, whilst the
% boost mode can increase the output voltage. In the buck-boost mode, the
% output voltage can be maintained either higher or lower than the source
% but in the opposite polarity. The simplest forms of these converters are
% schematically represented in Figure 1.

%%
%
%
%
%
% Figure 1. Buck (left), Boost (middle) and Buck-Boost (right) converters.

%%
% These converters consist of the same components, an inductor, $L$, a
% capacitor, $C$ and a switch, which has two states  $u=1$ and $u=0$. All
% converters connect to a DC power source with a voltage (unregulated),
% $V_\mathrm{in}$ and provide a regulated voltage, $v_\mathrm{o}$ to the
% load resistor, $R$ by controlling the state of the switch. In some
% situations, the load also could be inductive, for example a DC motor, or
% approximately, a current load, for example in a cascade configuration.
% For simplicity, here, only current and resistive loads are to be
% considered.
%

%% Principles
% The working principles of the these DC-DC converters can be explained as
% follows. In the buck mode, when the switch is on position 1, the DC
% source supplies power to the circuit which results an output voltage
% across the resistor. When the switch changes its position to 0, the
% energy stored in the inductor and capacitor will discharge through the
% resistor. Appropriately controlling the switching position can maintain
% the output voltage at a desired level lower than the source.
%
% In the boost mode, when the switch is on position 1, the circuit is
% separated into two parts: on the left, the source is charging the
% inductor, meanwhile the capacitor on the right maintains the output
% voltage using previously stored energy. When the switch changes its
% position to 0, both the DC source and energy stored in the inductor will
% supply power to the circuit on the right, hence boost the output voltage.
% Again, the output voltage can be maintain at desired level by controlling
% the switching sequence.
%
% Finally, for the buck-boost mode, switch positions 1 and 0 represents
% charging and discharging modes of the inductor. Appropriately controlling
% the switching sequence can result in output voltage higher or lower than
% the DC source. Since the inductor cannot change the direction of current,
% the output voltage is opposite to the DC source.

%% Model under ideal assumptions
% Under ideal assumptions: ideal switch, ideal capacitor and ideal
% inductor, these converters can be described using ordinary
% differentiation equations as follows:
%
%% Buck converter:
% $C{dv_{c} \over dt} = i_{L} - v_{c}/R - i_{o}$
%
% $L{di_{L} \over dt} = uv_{in} - v_{c}$
%
% where, $i_\mathrm{o}$ is the load current.
%
%% Boost converter:
% $C{dv_\mathrm{c}\over dt}=(1-u)i_\mathrm{L}-v_\mathrm{c}/R-i_\mathrm{o}$
%
% $L{di_\mathrm{L}\over dt}=v_\mathrm{in} - (1-u) v_\mathrm{c}$
%
%% Buck-boost converter:
% $C{dv_\mathrm{c}\over dt}=(1-u)i_\mathrm{L}-v_\mathrm{c}/R-i_\mathrm{o}$
%
% $L{di_\mathrm{L}\over dt}=uv_\mathrm{in} - (1-u) v_\mathrm{c}$
%
% Introduce the following state, time and load normalization:
%
% $x_1={v_\mathrm{c}\over v_\mathrm{in}}$,
%
% $x_2={i_\mathrm{L}\over v_\mathrm{in}}\sqrt{{L}\over{C}}$,
%
% $\tau = {t \over \sqrt{LC}}$,
%
% $\gamma = {\sqrt{LC}\over R}$,
%
% $d = {i_\mathrm o \over v_\mathrm{in}}\sqrt{{L}\over {C}}$
%
% Then the normalized state equations of three converters are as follows:
%
%% Normalized buck model
% $\dot{x_1} = -\gamma x_1 + x_2 - d$
%
% $\dot{x_2}= -x_1 + u $
%
% where, with an abuse of notation, `.' represents the derivation with
% respect to the normalized time, $\tau$.
%
%% Normalized boost model
% $\dot{x_1} = -\gamma x_1 + (1-u)x_2 - d$
%
% $\dot{x_2} = -(1-u)x_1 + 1$
%
%% Normalized buck-boost model
% $\dot{x_1} = -\gamma x_1 + (1-u)x_2 -d$
%
% $\dot{x_2} = -(1-u)x_1 + u$
%

%% Model with body resistors
% In more general cases, a body resistor of the inductor, $R_\mathrm{L}$
% and an equivalent series resistor (ESR) of the capacitor, $R_\mathrm{c}$
% can be added to the above models.
%
%% Buck model with $R_\mathrm{L}$ and $R_\mathrm{c}$
% Since,
%
% $C{dv_\mathrm{c}\over dt} = i_\mathrm{L} - v_\mathrm{o}/R - i_\mathrm{o}$
%
% $v_\mathrm{o} = v_\mathrm{c} + R_\mathrm{c}C{dv_\mathrm{c}\over dt}$
%
% $L{di_\mathrm{L}\over dt} = uv_\mathrm{in} - v_\mathrm{o} -
% R_\mathrm{L}i_\mathrm{L}$
%
% Inserting the second equation into the first leads to:
%
% $C{dv_\mathrm{c}\over dt} = i_\mathrm{L} - v_\mathrm{c}/R -
% {R_\mathrm{c}\over R}C{dv_\mathrm{c}\over dt}- i_\mathrm{o}$
%
% $\left(1+{R_\mathrm{c}\over R}\right)C{dv_\mathrm{c}\over
% dt}=i_\mathrm{L} - v_\mathrm{c}/R - i_\mathrm{o}$
%
% Hence,
%
% $v_\mathrm{o}={Rv_\mathrm{c}\over R+R_\mathrm{c}}+{RR_\mathrm{c}\over
% R+R_\mathrm{c}}(i_\mathrm{L}-i_\mathrm{o})$
%
% and the overall model is
%
% $C{dv_\mathrm{c}\over dt}={R\over R+R_\mathrm{c}}\left(i_\mathrm{L} - {v_\mathrm{c}\over R} - i_\mathrm{o}\right)$
%
% $L{di_\mathrm{L}\over dt} = uv_\mathrm{in} - {Rv_\mathrm{c}\over
% R+R_\mathrm{c}} - \left(R_\mathrm{L}+{RR_\mathrm{c}\over R+R_\mathrm{c}}\right) i_\mathrm{L} + {RR_\mathrm{c}i_\mathrm{o}\over R+R_\mathrm{c}}$
%
% $v_\mathrm{o}={Rv_\mathrm{c}\over R+R_\mathrm{c}}+{RR_\mathrm{c}\over
% R+R_\mathrm{c}}(i_\mathrm{L}-i_\mathrm{o})$
%
%% Boost model with $R_\mathrm{L}$ and $R_\mathrm{c}$
% $C{dv_\mathrm{c}\over dt} = (1-u)i_\mathrm{L} - v_\mathrm{o}/R - i_\mathrm{o}$
%
% $L{di_\mathrm{L}\over dt} = v_\mathrm{in} - (1-u) v_\mathrm{o} - R_\mathrm{L}i_\mathrm{L}$
%
% $v_\mathrm{o} = {Rv_\mathrm{c}\over R+R_\mathrm{c}}+{RR_\mathrm{c}\over
% R+R_\mathrm{c}}((1-u)i_\mathrm{L}-i_\mathrm{o})$
%
%% Buck-boost model with $R_\mathrm{L}$ and $R_\mathrm{c}$
% $C{dv_\mathrm{c}\over dt} = (1-u)i_\mathrm{L} - v_\mathrm{o}/R-i_\mathrm{o}$
%
% $L{di_\mathrm{L}\over dt}=uv_\mathrm{in}-(1-u) v_\mathrm{o}-R_\mathrm{L}i_\mathrm{L}$
%
% $v_\mathrm{o}={Rv_\mathrm{c}\over R+R_\mathrm{c}}+{RR_\mathrm{c}\over R+R_\mathrm{c}}((1-u)i_\mathrm{L}-i_\mathrm{o})$
%
%% Simulink Model
% These three modes of DC-DC converters have been uniformly implemented in
% the MATLAB/Simulink as show in Figure~\ref{fig:simulinkmodel}.
%
% <<DCDCModelInside.png>>
%
% Figure 2. A uniform Simulink model of DC-DC converters.
%
% The input-output connections of the model is shown in Figure 3.
%
% <<DCDCModel.png>>
%
% Figure 3. Input and output connections of the DC-DC converter model.
%
% The first input to the model is the switch signal eight 1 or 0. The
% second one defines the DC source voltage and internal resistance. The
% third input is used to define the output current. The model has two
% outputs, the output voltage and the inductor current, which are the
% states of the system.
%
% The model can be configure with a number of parameters as shown in
% Figure~\ref{fig:simulinkmodelparameters}. These parameters are: the
% capacitance, $C$, inductance, $L$, the internal resistance of the
% capacitor and the inductor, $R_C$ and $R_L$ respectively. Three converter
% modes can be selected through the pull-down menu. One can also define
% either zero or non-zero value to the initial capacitor voltage by
% selecting or de-selecting the ``zero capacitor voltage'' option. Finally,
% the option ``Positive Inductor Current'' defines whether the condition
% $i_L\ge 0$ should be enforced or not.
%
% <<DCDCModelParameters.png>>
%
% Figure 4. Parameters of the DC-DC converter model.
%

author: Sable


text:
function [outputs, out_blob_names_total] = caffe_forward_net(net, input, out_blob_names_extra)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
% c
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

if ~exist('out_blob_names_extra', 'var'), out_blob_names_extra = {}; end
assert(iscell(out_blob_names_extra));

input_size  = caffe_get_blobs_size(net, net.inputs);
assert(numel(input_size) == numel(input));
size_in     = size(input{1});
num_feats   = size_in(end);

out_blob_names       = net.outputs;
num_out_blobs        = length(out_blob_names);
num_out_blobs_extra  = length(out_blob_names_extra);
num_out_blobs_total  = num_out_blobs + num_out_blobs_extra;
out_blob_names_total = [out_blob_names(:); out_blob_names_extra(:)];

output_size  = caffe_get_blobs_size(net, out_blob_names_total);
outputs      = cell(num_out_blobs_total,1);
out_feat_dim = zeros(num_out_blobs_total,1);

batch_size  = input_size{1}(4);
num_batches = ceil(num_feats/batch_size);

for i = 1:num_out_blobs_total
    out_feat_dim(i) = prod(output_size{i}(1:(end-1)));
    if output_size{i}(end) == batch_size
        outputs{i} = zeros(out_feat_dim(i), num_feats, 'single');
    else
        outputs{i} = zeros(out_feat_dim(i), output_size{i}(end) * num_batches, 'single');
    end
end

for i = 1:num_batches
    start_idx = (i-1) * batch_size + 1;
    stop_idx  = min(i * batch_size, num_feats);
    batch_size_this = stop_idx - start_idx + 1;

    % forward propagate batch of region images
    out = net.forward(prepare_batch(input, input_size, start_idx, stop_idx));
    for j = (num_out_blobs+1):num_out_blobs_total
        out{j} = net.blobs(out_blob_names_total{j}).get_data;
    end

    for j = 1:num_out_blobs_total
        out{j} = reshape(out{j}, [out_feat_dim(j), output_size{j}(end)]);
        if output_size{j}(end) == batch_size
            if (i == num_batches)
                out{j} = out{j}(:,1:batch_size_this);
            end
            outputs{j}(:,start_idx:stop_idx) = out{j};
        else
            start_idx0 = (i-1) * output_size{j}(end) + 1;
            stop_idx0  = i * output_size{j}(end);
            outputs{j}(:,start_idx0:stop_idx0) = out{j};
        end
    end
end

end

function batch = prepare_batch(input, input_size, start_idx, stop_idx)
assert(numel(input_size) == numel(input));
batch_size_this   = stop_idx - start_idx + 1;

batch = cell(1,numel(input_size));
for i = 1:numel(input_size)
    batch{i}          = zeros(input_size{i}, 'single');
    reshape_vector    = input_size{i};
    reshape_vector(4) = batch_size_this;
    batch{i}(:,:,:,1:batch_size_this) = reshape(input{i}(:,start_idx:stop_idx), reshape_vector);
end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
%	OVERVIEW:
%       This is demo for the HRV PhysioNet Cardiovascular Signal Toolbox
%       using RR intervals with annotations.
%       Provided data are a subset from the MIT Physionet
%       NSR dataset, which contains long-term ECG recordings of subjects
%       in normal sinus rhythm.
%       It shows how to automaticly import multiple files from a folder,
%       perfrom the HRV analysis on each of them and then store the results
%       in .csv format.
%       It uses the default parameters in the configuration file using
%       'demo_NSR' option : InitializeHRVparams('demo_NSR').
%
%   OUTPUT:
%       HRV Metrics exported to .cvs files
%
%   DEPENDENCIES & LIBRARIES:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   REFERENCE:
%       Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%       Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Giulia Da Poian
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close all;

run(['..' filesep 'startup.m'])

% Remove old files generated by this demo
OldFolder = [pwd,filesep, 'OutputData', filesep, 'ResultsNSR'];
if exist(OldFolder, 'dir')
    rmdir(OldFolder, 's');
    fprintf('Old Demo Folder deleted \n \n');
end

% Initialize settings for demo
HRVparams = InitializeHRVparams('demo_NSR');
HRVparams.MSE.on = 0; % No MSE analysis for this demo
HRVparams.DFA.on = 0; % No DFA analysis for this demo
HRVparams.HRT.on = 0; % No HRT analysis for this demo
HRVparams.output.separate = 0;   % For this demo write all the results in one file

% Check for a list of files to be analyzed in current directory in .dat format
[subjectIDs,filesTBA] = GenerateListOfFilesTBA(HRVparams.ext, HRVparams.readdata,[]);

% Prepare for parallel loop by eliminating variables
clear nummatchingfiles x i filename flag match
numsub = length(subjectIDs);
notAnalyzed = 0;

% NOTE: This loop can be run in parallel by changing the loop to a parfor
% loop.
for i_patient = 1:numsub

    thisPatient = strcat(HRVparams.readdata, filesep, subjectIDs(i_patient));
    try

        % 1. Import Patient Data
        RRwindowStartIndices = [];
        tNN = [];
        NN = [];
        [samples,annotations] = read_ann(thisPatient{1},HRVparams.ext);
        rr = diff(samples)./HRVparams.Fs;
        t = cumsum(rr);

        % Demo keeps only the first 2h

        rr = rr(t<60*60*2);
        t = t(t<60*60*2);

        % 2. Perform HRV analysis on the RR intervals
        [results, resFilenameHRV] = Main_HRV_Analysis(rr,t,'RRIntervals',HRVparams,subjectIDs(i_patient),annotations);
        currentFile = [HRVparams.writedata filesep resFilenameHRV.HRV '.csv'];

        fprintf('\n');
    catch

        results = NaN;
        col_titles = {'NaN'};
        currentFile = '';
        notAnalyzed = 1;
        fprintf('Error on subject %s \n', char(subjectIDs(i_patient)));

    end

end

% 3. Compare generated output file with the reference one

referenceFile = ['ReferenceOutput' filesep 'NSR_HRV_allwindows_allpatients.csv'];
testHRV = CompareOutput(currentFile,referenceFile);

if testHRV
    fprintf('\n ** DemoAnnotatedData: TEST SUCCEEDED ** \n ')
elseif notAnalyzed == 0
    fprintf('\n ** DemoAnnotatedData: TEST FAILED ** \n')
    fprintf('Error: generated output does not match reference \n')
elseif notAnalyzed == 1
    fprintf('\n** DemoAnnotatedData: TEST FAILED ** \n')
    fprintf('Error: analysis not performed \n');
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function Handle = GCO_Create(NumSites,NumLabels)
% GCO_Create    Create a GCoptimization object.
%    Handle = GCO_Create(NumSites,NumLabels) creates a new GCoptimization
%    object and returns a 'handle' to uniquely identify it.
%    Call GCO_Delete(Handle) to delete the object and free its memory.
%    Call GCO_Delete(GCO_ListHandles) to delete all GCO objects.

GCO_LoadLib();
if (nargin < 2), error('Expected 2 arguments'); end
Handle = gco_matlab('gco_create_general',int32(NumSites),int32(NumLabels));
end

author: andrewssobral

repo: lrslibrary