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

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

% Last Modified by GUIDE v2.5 23-May-2012 10:45:04

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

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

% --- Executes just before SwitchTarget is made visible.
function SwitchTarget_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to SwitchTarget (see VARARGIN)

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

handles.JLabelhObject = varargin{1};
[handles.ntargets_per_page] = myparse(varargin(2:end),'ntargets_per_page',100);

JLabelHandles = guidata(handles.JLabelhObject);
handles.JLDObj = JLabelHandles.guidata.data;
% set(handles.axes1,'axes',off);
handles.tablePos = get(handles.table,'Position');
handles.switchFlyPos = get(handles.pushSwitchfly,'Position');
handles.updatePos = get(handles.pushbutton_update,'Position');
handles.closePos = get(handles.pushClose,'Position');
handles.figurePos = get(handles.figure1,'Position');
handles.pagePos = get(handles.popupmenu_Page,'Position');

% initialize page number
handles.page_number = [];

% added an extra guidata here, as sometimes I was getting an error that
% figurePos wasn't getting set
guidata(hObject,handles);
handles = updateTable(handles);
guidata(hObject,handles);

function [exp,target] = GlobalTarget2LocalTarget(handles,globaltarget)

exp = find(handles.cs_ntargets >= globaltarget,1,'first');
if exp == 1,
  target = globaltarget;
else
  target = globaltarget - handles.cs_ntargets(exp-1);
end

function [globaltarget] = LocalTarget2GlobalTarget(handles,exp,target)

if exp == 1,
  globaltarget = target;
else
  globaltarget = handles.cs_ntargets(exp-1)+target;
end

function [page,row] = GlobalTarget2PageRow(handles,globaltarget)

page = ceil(globaltarget/handles.ntargets_per_page);
row = globaltarget - (page-1)*handles.ntargets_per_page;

function [globaltarget] = PageRow2GlobalTarget(handles,page,row)

globaltarget = handles.ntargets_per_page*(page-1)+row;

function [exp,target] = PageRow2LocalTarget(handles,page,row)

[exp,target] = GlobalTarget2LocalTarget(handles,PageRow2GlobalTarget(handles,page,row));

% -------------------------------------------------------------------------
function handles = updateTable(handles)
% Initialize the table

%MERGESTUPDATED

set(handles.pushbutton_update,'enable','off');
handles.curExp = handles.JLDObj.GetExp();
handles.curFly = handles.JLDObj.GetFlies();

% which fly out of all flies
handles.cs_ntargets = cumsum(handles.JLDObj.nflies_per_exp);
handles.curFlyGlobal = LocalTarget2GlobalTarget(handles,handles.curExp,handles.curFly);

% which page, row
[handles.current_target_page_number,handles.current_target_row] = GlobalTarget2PageRow(handles,handles.curFlyGlobal);

% set page number if not already set
if isempty(handles.page_number),
  handles.page_number = handles.current_target_page_number;
end

% set page number options
handles.npages = ceil(handles.cs_ntargets(end)/handles.ntargets_per_page);
handles.pagestrs = cell(1,handles.npages);
for i = 1:handles.npages,
  exp1 = PageRow2LocalTarget(handles,i,1);
  tmp = min(PageRow2GlobalTarget(handles,i,handles.ntargets_per_page),handles.cs_ntargets(end));
  exp2 = GlobalTarget2LocalTarget(handles,tmp);
  if exp1 >= exp2,
    handles.pagestrs{i} = sprintf('%d (Exp %d)',i,exp1);
  else
    handles.pagestrs{i} = sprintf('%d (Exps %d-%d)',i,exp1,exp2);
  end
  if i == handles.current_target_page_number,
    handles.pagestrs{i} = ['* ',handles.pagestrs{i},' *'];
  end

end

if handles.page_number > handles.npages,
  handles.page_number = handles.npages;
end

set(handles.popupmenu_Page,'String',handles.pagestrs,'Value',handles.page_number);

fieldList = {};
nCls = handles.JLDObj.nclassifiers;
iCls2LblNames = handles.JLDObj.iCls2LblNames;
if ~handles.JLDObj.IsGTMode()
  fieldList(end+1,:) = { @(gs,cs)gs.trajLength  'Trajectory| Length'};
  fieldList(end+1,:) = { @(gs,cs)gs.firstframe  'Start frame'};
  fieldList(end+1,:) = { @(gs,cs)gs.sexfrac     'Sex|(% male)'};
  fieldList(end+1,:) = { @(gs,cs)gs.nbouts      'Bouts|labeled'};
  fieldList(end+1,:) = { @(gs,cs)gs.totalframes 'Frames|labeled'};

  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).posframes sprintf('Frames|labeled|%s',iCls2LblNames{iCls}{1})}; %#ok<*AGROW>
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).negframes sprintf('Frames|labeled|%s',iCls2LblNames{iCls}{2})};
  end
  if nCls>1
    fieldList(end+1,:) = { @(gs,cs)gs.posframes   'Frames|labeled|pos'};
    fieldList(end+1,:) = { @(gs,cs)gs.negframes   'Frames|labeled|neg'};
  end
  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscorepos sprintf('%s| predicted | (current)',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscoreneg sprintf('%s| predicted | (current)',iCls2LblNames{iCls}{2})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscoreframes sprintf('%s|predicted | total | (current)',iCls2LblNames{iCls}{1})};
  end
  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscorepos_loaded sprintf('%s| predicted | (imported)',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscoreneg_loaded sprintf('%s| predicted | (imported)',iCls2LblNames{iCls}{2})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).nscoreframes_loaded sprintf('%s|predicted | total | (imported)',iCls2LblNames{iCls}{1})};
  end
  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).errorsPos sprintf('%s|errors',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).errorsNeg sprintf('%s|errors',iCls2LblNames{iCls}{2})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).validatedErrorsPos sprintf('%s|validation|errors',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).validatedErrorsNeg sprintf('%s|validation|errors',iCls2LblNames{iCls}{2})};
  end
  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).one2two sprintf('%s|switched to|%s',iCls2LblNames{iCls}{1},iCls2LblNames{iCls}{2})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).two2one sprintf('%s|switched to|%s',iCls2LblNames{iCls}{2},iCls2LblNames{iCls}{1})};
  end
else
  fieldList(end+1,:) = { @(gs,cs)gs.trajLength,'Trajectory| Length'};
  fieldList(end+1,:) = { @(gs,cs)gs.firstframe,'Start frame'};
  fieldList(end+1,:) = { @(gs,cs)gs.sexfrac,'Sex|(% male)'};
  for iCls = 1:nCls
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_posframes,sprintf('Ground Truthing|Frames|labeled|%s',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_negframes,sprintf('Ground Truthing|Frames|labeled|%s',iCls2LblNames{iCls}{2})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_suggestion_frames,sprintf('Ground Truthing|Frames|Suggested|%s',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_nbouts,sprintf('Ground Truthing|Bouts|labeled|%s',iCls2LblNames{iCls}{1})};
    fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_totalframes,sprintf('Ground Truthing|Frames|labeled|total|%s',iCls2LblNames{iCls}{1})};
%     fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_posframes,sprintf('Ground Truthing|Frames|labeled|%s',iCls2LblNames{iCls}{1}) };
%     fieldList(end+1,:) = { @(gs,cs)cs(iCls).gt_negframes,sprintf('Ground Truthing|Frames|labeled|%s',iCls2LblNames{iCls}{2}) };
  end
end

colFormat{1} = 'char';
colFormat{2} = 'numeric';
for ndx = 1:size(fieldList,1)
  colFormat{2+ndx} = 'numeric';
end

% reorganize cached data
newCachedTableData = cell(handles.cs_ntargets(end),size(fieldList,1)+2);
newCachedDataExpi = zeros(1,handles.cs_ntargets(end));
newCachedDataTarget = zeros(1,handles.cs_ntargets(end));
newCachedExpDirs = handles.JLDObj.expdirs;

if isfield(handles,'cachedTableData'),

  % reindexing for experiments
  [ism,old2newexpi] = ismember(handles.cachedExpDirs,newCachedExpDirs);

  % move data from old to new
  for oldexpi = find(ism(:)'),
    newexpi = old2newexpi(oldexpi);
    oldidx = find(handles.cachedDataExpi == oldexpi);
    if isempty(oldidx),
      continue;
    end
    [targets_curr,order] = sort(handles.cachedDataTarget(oldidx));
    oldidx = oldidx(order);
    newidx = LocalTarget2GlobalTarget(handles,newexpi,targets_curr);
    newCachedTableData(newidx,:) = handles.cachedTableData(oldidx,:);
    newCachedDataExpi(newidx) = newexpi;
    newCachedDataTarget(newidx) = handles.cachedDataTarget(oldidx);
  end

end

% overwrite
handles.cachedTableData = newCachedTableData;
handles.cachedDataExpi = newCachedDataExpi;
handles.cachedDataTarget = newCachedDataTarget;
handles.cachedExpDirs = newCachedExpDirs;

tableData = {};
tableExpi = [];
tableTarget = [];
count = 1;

[exp_start,target_start] = PageRow2LocalTarget(handles,handles.page_number,1);
tmp = min(PageRow2GlobalTarget(handles,handles.page_number,handles.ntargets_per_page),handles.cs_ntargets(end));
[exp_end,target_end] = GlobalTarget2LocalTarget(handles,tmp);

for selExp = exp_start:exp_end,
  if selExp == exp_start,
    startcurr = target_start;
  else
    startcurr = 1;
  end
  endcurr = handles.JLDObj.nflies_per_exp(selExp);
  if selExp == exp_end,
    endcurr = min(endcurr,target_end);
  end
  for flyNum = startcurr:endcurr,
    [genStats,clsStats] = handles.JLDObj.GetFlyStats(selExp,flyNum);
    tableData{count,1} = handles.JLDObj.expnames{selExp};
    tableData{count,2} = flyNum;
    for ndx = 1:size(fieldList,1)
      fcn = fieldList{ndx,1};
      tableData{count,2+ndx} = fcn(genStats,clsStats);
    end
    tableExpi(count) = selExp;
    tableTarget(count) = flyNum;
    count = count+1;
  end
end
handles.tableData = tableData;

% store in global cache
start_target = PageRow2GlobalTarget(handles,handles.page_number,1);
end_target = start_target + numel(tableExpi) - 1;
handles.cachedTableData(start_target:end_target,:) = tableData;
handles.cachedDataExpi(start_target:end_target) = tableExpi;
handles.cachedDataTarget(start_target:end_target) = tableTarget;

set(handles.table,'Data',tableData);
set(handles.table,'ColumnName',{});
drawnow();
set(handles.table,'ColumnName',...
  {'Experiment Name','Target|Number',fieldList{:,2}});
handles.fieldList = fieldList;
set(handles.table,'ColumnEditable',false);
set(handles.table,'ColumnFormat',colFormat);
set(handles.table,'CellSelectionCallback',@tableSelect);
set(handles.pushbutton_update,'enable','on');

% Update handles structure

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

% -------------------------------------------------------------------------
function initTable(hObject)
% Use java objects to tweak the table. Found this online at
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/298335

handles = guidata(hObject);
jscrollpane = findjobj(handles.table);
jtable = jscrollpane.getViewport.getView;
jtable.setSortable(false);

colWidth = repmat({70},1,13);
colWidth{1} = 300;
set(handles.table,'ColumnWidth',colWidth);

uitablepeer  = findjobj(hObject,'-nomenu','class','uitablepeer');
set(uitablepeer,'MouseClickedCallback',@MouseClickHandler);

% rowHeaderViewport=jscrollpane.getComponent(4);
% rowHeader=rowHeaderViewport.getComponent(0);
% rowHeader.setSize(80,360);
%
% %resize the row header
% newWidth=0; %100 pixels.
% rowHeaderViewport.setPreferredSize(java.awt.Dimension(newWidth,0));
% height=rowHeader.getHeight;
% rowHeader.setPreferredSize(java.awt.Dimension(newWidth,height));
% rowHeader.setSize(newWidth,height);

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

% Get default command line output from handles structure
% set double click callback

varargout{1} = handles.output;

function tableSelect(hObject,eventData)

handles = guidata(hObject);
jscrollpane = findjobj(handles.table);
jtable = jscrollpane.getViewport.getView;

if(size(eventData.Indices,1)==1)
  ndx = jtable.getActualRowAt(eventData.Indices(1,1)-1)+1;

  start_target = PageRow2GlobalTarget(handles,handles.page_number,1);
  handles.curExp = handles.cachedDataExpi(start_target + ndx - 1);
  handles.curFly = handles.tableData{ndx,2};
end

guidata(hObject,handles);

% --- Executes on button press in pushSwitchfly.
function pushSwitchfly_Callback(hObject, eventdata, handles)
% hObject    handle to pushSwitchfly (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
JLabelHandles = guidata(handles.JLabelhObject);
JLabelHandles = JLabel('SetCurrentMovie',JLabelHandles,handles.curExp);
guidata(handles.JLabelhObject,JLabelHandles);
JLabelHandles = JLabel('SetCurrentFlies',JLabelHandles,handles.curFly);
guidata(handles.JLabelhObject,JLabelHandles);
return

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

% --- Executes when figure1 is resized.
function figure1_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles,'figurePos'),
  return;
end

oldClosePosRight = handles.figurePos(3)-handles.closePos(1);
guiPos = get(handles.figure1,'Position');
newClosePos = [guiPos(3)-oldClosePosRight handles.closePos(2:4)];
set(handles.pushClose,'Position',newClosePos);

oldSwitchPosRight = handles.figurePos(3)-handles.switchFlyPos(1);
newSwitchFlyPos = [guiPos(3)-oldSwitchPosRight handles.switchFlyPos(2:4)];
set(handles.pushSwitchfly,'Position',newSwitchFlyPos);

oldUpdatePosRight = handles.figurePos(3)-handles.updatePos(1);
newUpdatePos = [guiPos(3)-oldUpdatePosRight handles.updatePos(2:4)];
set(handles.pushbutton_update,'Position',newUpdatePos);

newTableWidth = guiPos(3)- handles.figurePos(3)+handles.tablePos(3);
newTableHeight = guiPos(4)- handles.figurePos(4)+handles.tablePos(4);
set(handles.table,'Position',[handles.tablePos(1:2) newTableWidth newTableHeight]);

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

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

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

handles.page_number = get(hObject,'Value');
start_target = PageRow2GlobalTarget(handles,handles.page_number,1);
end_target = min(handles.cs_ntargets(end),start_target+handles.ntargets_per_page-1);
if all(handles.cachedDataExpi(start_target:end_target) ~= 0),
  handles.tableData = handles.cachedTableData(start_target:end_target,:);
  set(handles.table,'Data',handles.tableData);
  set(handles.table,'ColumnName',{});
  set(handles.table,'ColumnName',...
    {'Experiment Name','Target|Number',handles.fieldList{:,2}});
else
  handles = updateTable(handles);
end
guidata(hObject,handles);

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

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

% From http://www.mathworks.com/matlabcentral/newsreader/view_thread/270514

function MouseClickHandler(handle,cbData)
  % handle ~ java object UITablePeer
  % cbData ~ callback data for the MouseClickedCallback event

  if get(cbData,'ClickCount') == 2
      curfig = findall(0,'type','figure','name','SwitchTarget');
      if numel(curfig)~=1,
        return;
      end

      pushSwitchfly_Callback(curfig,[],guidata(curfig));
  end


text:
function out = DN_HistogramMode(y,numBins,doSimple,doPlot)
% DN_HistogramMode      Mode of a data vector.
%
% Measures the mode of the data vector using histograms with a given number
% of bins.
%
%---INPUTS:
%
% y, the input data vector.
% numBins, the number of bins to use in the histogram.
% doSimple, whether to use a simple binning method (linearly spaced bins).
% doPlot, whether to show a plot of what was computed.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs and set defaults:
%-------------------------------------------------------------------------------
if nargin < 2
    numBins = 'auto';
end
if nargin < 3
    doSimple = true;
end
if nargin < 4
    doPlot = false;
end
%-------------------------------------------------------------------------------

% Compute the histogram from the data:
if isnumeric(numBins)
    if doSimple
        [N,binEdges] = BF_SimpleBinner(y,numBins);
    else
        [N,binEdges] = histcounts(y,numBins);
    end
elseif ischar(numBins)
    [N,binEdges] = histcounts(y,'BinMethod',numBins);
else
    error('Unknown format for numBins');
end

% Compute bin centers from bin edges:
binCenters = mean([binEdges(1:end-1); binEdges(2:end)]);

% Mean position of maximums (if multiple):
out = mean(binCenters(N == max(N)));

% Plot a summary of what was computed:
if doPlot
    histogram('BinEdges',binEdges,'BinCounts',N,'EdgeColor','k','FaceColor',0.6*ones(1,3));
    hold('on');
    plot(out*ones(2,1),[0,max(N)],'r','LineWidth',2);
    hold('off')
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function val=wavelet_coeff(data,fs)
num_layers=7;
num_coeffs=8;
[h,g,rh,rg]=daub(num_coeffs);

wav_res=wt(data,h,g,num_layers);
pos=sub_pos(length(data),num_layers);
bands=extract_subbands(wav_res,pos);

%MEAN OF COEFFICIENTS AND ABSOLUTE ENERGY
for j=1:length(bands)
    mean_coeffs(j)=mean(abs(bands{j}));
end
val=mean_coeffs(5);

author: alexandrebarachant


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Inverse DFT of X(k)=[6, -1-j,0,-1+j], k=0,1,2,3

clear;
Xk=[6, -1-j,0,-1+j];
N=length(Xk);
for n=0:N-1
for k=0:N-1
xn(k+1)=Xk(k+1)*exp(j*2*pi*n*k/N);
end
x(n+1)=sum(xn);
end
x=(1/N)*x

author: Sable

repo: mcbench-benchmarks


text:
function val=vec(A)
%%VEC Given A matrix, return the values of the matrix stacked columnwise.
%     This is just the same as A(:), but this function can be more
%     convenient to use when programming equations as one can
%     say things like c=b+vec(A*B) rather than temp=A*B; c=b+temp(:).
%
%INPUTS: A A matrix.
%
%OUTPUTS: val The values in the matrix stacked column-wise. If A is a
%             hypermatrix, then the values are stacked in order of the
%             indices of A (rows, columns, hypercolumns, etc.).
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

val=A(:);
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 out = isinf(f)
%ISINF   Test if a CLASSICFUN is unbounded.
%   ISINF(F) returns TRUE if F takes an infinite value and FALSE otherwise.
%
% See also ISFINITE, ISNAN.

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

% Check if the ONEFUN of f is infinite:
out = isinf(f.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% this script tests NGLDM features between CERR and pyradiomics.
%
% RKP, 03/22/2018

%% load image
gldmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','tests_for_cerr','test_ngldm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','data_for_cerr_tests','CERR_plans','head_neck_ex1_20may03.mat.bz2');

planC = loadPlanC(cerrFileName,tempdir);
indexS = planC{end};

paramS = getRadiomicsParamTemplate(gldmParamFileName);
strNum = getMatchingIndex(paramS.structuresC{1},{planC{indexS.structures}.structureName});
scanNum = getStructureAssociatedScan(strNum,planC);

%% NGLDM features CERR

ngldmS = calcGlobalRadiomicsFeatures...
            (scanNum, strNum, paramS, planC);
ngldmS = ngldmS.Original.ngldmFeatS;

cerrNgldmV = [ngldmS.lde, ngldmS.hde, ngldmS.lgce, ngldmS.hgce, ...
    ngldmS.ldlge, ngldmS.ldhge, ngldmS.hdlge, ngldmS.hdhge, ...
    ngldmS.gln, ngldmS.glnNorm, ngldmS.dcn, ngldmS.dcnNorm,...
    ngldmS.dcp, ngldmS.glv, ngldmS.dcv, ngldmS.entropy, ngldmS.energy];

% %% Calculate features using pyradiomics
%
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
%     single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
%
% scanType = 'original';
%
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
%
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% %teststruct = PyradWrapper(testM, mask3M, scanType);
%
% pyradNgldmNamC = {'SmallDependenceEmphasis', 'LargeDependenceEmphasis',...
%     'LowGrayLevelCountEmphasis', 'HighGrayLevelCountEmphasis',  'SmallDependenceLowGrayLevelEmphasis', ...
%     'SmallDependenceHighGrayLevelEmphasis', 'LargeDependenceLowGrayLevelEmphasis', ...
%     'LargeDependenceHighGrayLevelEmphasis', 'GrayLevelNonUniformity', 'GrayLevelNonUniformityNorm', ...
%     'DependenceNonUniformity', 'DependenceNonUniformityNormalized', ...
%     'DependencePercentage', 'GrayLevelVariance', 'DependenceVariance', ...
%     'DependenceEntropy', 'DependenceEnergy'};
%
%
% pyradNgldmNamC = strcat(['original', '_gldm_'],pyradNgldmNamC);
%
% pyRadNgldmV = [];
% for i = 1:length(pyradNgldmNamC)
%     if isfield(teststruct,pyradNgldmNamC{i})
%         pyRadNgldmV(i) = teststruct.(pyradNgldmNamC{i});
%     else
%         pyRadNgldmV(i) = NaN;
%     end
% end
%
% %% Compare
% ngldmDiffV = (cerrNgldmV - pyRadNgldmV) ./ cerrNgldmV * 100

%% Compare with previously calculated pyradiomics features
saved_pyRadNgldmV = [0.177896163786327,80.7514028586554,NaN,NaN,0.00117423243903225,380.447741305586,3.84549682517831,136232.842773954,1250.78930651138,NaN,634.668395976707,0.0671962303839817,NaN,122.726429080550,26.6319748251348,7.44897904303236,NaN];
ngldmDiffV = (cerrNgldmV - saved_pyRadNgldmV) ./ cerrNgldmV * 100

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = sccn_xdf(filename, hdr, begsample, endsample, chanindx)

% This is a wrapper to the reading function from the XDF MATLAB toolbox.
%
% Use as
%   hdr = sccn_xdf(filename);
%   dat = sccn_xdf(filename, hdr, begsample, endsample, chanindx);
%   evt = sccn_xdf(filename, hdr);
%
% See also FT_FILETYPE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, XDF2FIELDTRIP

% Copyright (C) 2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% ensure this is on the path
ft_hastoolbox('xdf', 1);

needhdr = (nargin==1);
needevt = (nargin==2);
needdat = (nargin==5);

streams = load_xdf(filename);

iscontinuous = false(size(streams));
% figure out which streams contain continuous/regular and discrete/irregular data
for i=1:numel(streams)
  iscontinuous(i) = isfield(streams{i}.info, 'effective_srate');
end

% determine the stream with the highest sampling rate
srate = nan(size(streams));
for i=1:numel(streams)
  if iscontinuous(i)
    srate(i) = streams{i}.info.effective_srate;
  end
end
[~, indx] = max(srate);

% only keep the stream with the maximum sampling rate
% this is probably the EEG stream
stream = streams{indx};

if needhdr
  % this section of code is shared with xdf2fieldtrip
  hdr             = [];
  hdr.Fs          = stream.info.effective_srate;
  hdr.nChans      = numel(stream.info.desc.channels.channel);
  hdr.nSamplesPre = 0;
  hdr.nSamples    = length(stream.time_stamps);
  hdr.nTrials     = 1;
  hdr.label       = cell(hdr.nChans, 1);
  hdr.chantype    = cell(hdr.nChans, 1);
  hdr.chanunit    = cell(hdr.nChans, 1);

  prefix = stream.info.name;
  for i=1:hdr.nChans
    hdr.label{i} = [prefix '_' stream.info.desc.channels.channel{i}.label];
    hdr.chantype{i} = stream.info.desc.channels.channel{i}.type;
    hdr.chanunit{i} = stream.info.desc.channels.channel{i}.unit;
  end

  hdr.FirstTimeStamp     = stream.time_stamps(1);
  hdr.TimeStampPerSample = (stream.time_stamps(end)-stream.time_stamps(1)) / (length(stream.time_stamps) - 1);

  % keep the original header details
  hdr.orig = stream.info;

  % return the header
  varargout = {hdr};

elseif needevt
  streams = streams(~iscontinuous);
  event = [];
  for i=1:numel(streams)
    for j=1:numel(streams{i}.time_series)

      % convert the timestamps to the corresponding sample in the selected data stream
      % the first sample in the data stream corresponds to 1
      timestamp = streams{i}.time_stamps(j);
      sample = round((timestamp - hdr.FirstTimeStamp)/hdr.TimeStampPerSample) + 1;

      event(end+1).type      = streams{i}.info.type;
      event(end  ).value     = streams{i}.time_series{j};
      event(end  ).timestamp = timestamp;
      event(end  ).sample    = sample;
    end
  end

  % return the events
  varargout = {event};

elseif needdat
  dat = stream.time_series(chanindx, begsample:endsample);

  % return the data
  varargout = {dat};
end

author: spm


text:
%CODEGENERATION.GENMFUNINERTIA Generate M-function for robot inertia matrix
%
% cGen.genmfuninertia() generates a robot-specific M-function to compute
% robot inertia matrix.
%
% Notes::
% - Is called by CodeGenerator.geninertia if cGen has active flag genmfun
% - The inertia matrix is stored row by row to avoid memory issues.
% - The generated M-function recombines the individual M-functions for each row.
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.gencoriolis.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [] = genmfuninertia(CGen)

%% Does robot class exist?
if ~exist(fullfile(CGen.robjpath,[CGen.getrobfname,'.m']),'file')
    CGen.logmsg([datestr(now),'\tCreating ',CGen.getrobfname,' m-constructor ']);
    CGen.createmconstructor;
    CGen.logmsg('\t%s\n',' done!');
end

%%
CGen.logmsg([datestr(now),'\tGenerating m-function for the robot inertia matrix row' ]);

q = CGen.rob.gencoords;
nJoints = CGen.rob.n;

for kJoints = 1:nJoints
    CGen.logmsg(' %s ',num2str(kJoints));
    symname = ['inertia_row_',num2str(kJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    if exist(fname,'file')
        tmpStruct = load(fname);
    else
        error ('genmfuninertia:SymbolicsNotFound','Save symbolic expressions to disk first!')
    end

    funfilename = fullfile(CGen.robjpath,[symname,'.m']);

    matlabFunction(tmpStruct.(symname),'file',funfilename,...              % generate function m-file
        'outputs', {'Irow'},...
        'vars', {'rob',[q]});
    hStruct = createHeaderStructRow(CGen.rob,kJoints,symname);                 % replace autogenerated function header
    replaceheader(CGen,hStruct,funfilename);

end
CGen.logmsg('\t%s\n',' done!');

CGen.logmsg([datestr(now),'\tGenerating full inertia matrix m-function']);

    funfilename = fullfile(CGen.robjpath,'inertia.m');
    hStruct = createHeaderStructFullInertia(CGen.rob,funfilename);

    fid = fopen(funfilename,'w+');

    fprintf(fid, '%s\n', ['function I = inertia(rob,q)']);                 % Function definition
    fprintf(fid, '%s\n',constructheaderstring(CGen,hStruct));                   % Header

    fprintf(fid, '%s \n', 'I = zeros(length(q));');                        % Code
    for iJoints = 1:nJoints
        funcCall = ['I(',num2str(iJoints),',:) = ','rob.inertia_row_',num2str(iJoints),'(q);'];
        fprintf(fid, '%s \n', funcCall);
    end

    fclose(fid);

    CGen.logmsg('\t%s\n',' done!');
end

function hStruct = createHeaderStructRow(rob,curJointIdx,fName)
[~,hStruct.funName] = fileparts(fName);
hStruct.shortDescription = ['Computation of the robot specific inertia matrix row for corresponding to joint ', num2str(curJointIdx), ' of ',num2str(rob.n),'.'];
hStruct.calls = {['Irow = ',hStruct.funName,'(rob,q)'],...
    ['Irow = rob.',hStruct.funName,'(q)']};
hStruct.detailedDescription = {'Given a full set of joint variables this function computes the',...
                               ['inertia matrix row number ', num2str(curJointIdx),' of ',num2str(rob.n),' for ',rob.name,'.']};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['q:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     coordinates',...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['Irow:  [1x',int2str(rob.n),'] row of the robot inertia matrix']};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
    'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'coriolis'};
end

function hStruct = createHeaderStructFullInertia(rob,fname)
[~,hStruct.funName] = fileparts(fname);
hStruct.shortDescription = ['Inertia matrix for the ',rob.name,' arm.'];
hStruct.calls = {['I = ',hStruct.funName,'(rob,q)'],...
    ['I = rob.',hStruct.funName,'(q)']};
hStruct.detailedDescription = {'Given a full set of joint variables the function computes the',...
                               'inertia Matrix of the robot.'};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['q:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     coordinates',...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['I:  [',int2str(rob.n),'x',int2str(rob.n),'] inertia matrix']};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
    'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'coriolis'};
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ mD ] = CalcDistanceMatrixCCols( mX )
% ----------------------------------------------------------------------------------------------- %
% [ mD ] = CalcDistanceMatrixCCols( mX )
%   Calculates the distance matrix for the input data. The distance matrix
%   is a symmetric matrix where 'mD(ii, jj) = dist(mX(:, ii), mX(:, jj));'.
%   This function uses the squared Euclidean Distance for the distance
%   metric.
% Input:
%   - mX            -   Data Matrix.
%                       Each data sample is a column of the matrix.
%                       Structure: Matrix (varDim x numVars).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - mD            -   Distance Matrix.
%                       A symmetric matrix where 'mD(ii, jj) = dist(mX(:,
%                       ii), mX(:, jj));'.
%                       Structure: Matrix (numVars x numVars).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     01/01/2021  Royi Avital	RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

mG = mX.' * mX;
vG = diag(mG);

mD = vG.' + vG - (2 * mG);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [w,scnl]=obspy.stream2waveform(path_to_python, path_to_converter, path_to_data)
%OBSPY.STREAM2WAVEFORM Use ObsPy to read a stream object and then convert
% into a waveform object. Also return a scnlobject.
% Example:
%  [w,scnl]=obspy.stream2waveform('/Users/glennthompson/anaconda/bin/python', ...
%       '/Users/glennthompson/obspy_stream2matfile.py', ...
%       'https://examples.obspy.org/BW.BGLD..EH.D.2010.037')
%  plot(w)
    commandstr=sprintf('%s %s "%s"',path_to_python, path_to_converter, path_to_data);
    system(commandstr);
    clear commandstr fname
    d=dir('obspy.stream.*.mat');
    disp(sprintf('%d trace matfiles found\n',length(d)))
    for c=1:length(d)
        tr=load(d(c).name);
        disp(sprintf('Loaded trace %d',c))
        scnl(c)=scnlobject(tr.station, tr.channel, tr.network, tr.location);
        snum=(tr.starttime.timestamp/86400)+datenum(1970,1,1);
        w(c)=waveform(scnl, tr.sampling_rate, snum, tr.data);
        delete(d(c).name);
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function N=patchnormals(FV)
% This function PATCHNORMALS calculates the normals of a triangulated
% mesh. PATCHNORMALS calls the patchnormal_double.c mex function which
% first calculates the normals of all faces, and after that calculates
% the vertice normals from the face normals weighted by the angles
% of the faces.
%
% N=patchnormals(FV);
%
% Inputs,
%   FV : A struct containing FV.faces with a facelist Nx3 and FV.vertices
%        with a Nx3 vertices list. Such a structure is created by Matlab
%        Patch function
% Outputs,
%   N : A Mx3 list with the normals of all vertices
%
% Example,
%   % Compile the c-coded function
%   mex patchnormals_double.c -v
%
%   % Load a triangulated mesh of a sphere
%   load sphere;
%
%   % Calculate the normals
%   N=patchnormals(FV);
%
%   % Show the normals
%   figure, patch(FV,'FaceColor',[1 0 0]); axis square; hold on;
%   for i=1:size(N,1);
%       p1=FV.vertices(i,:); p2=FV.vertices(i,:)+10*N(i,:);
%       plot3([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)],'g-');
%   end
%
% Function is written by D.Kroon University of Twente (June 2009)

sizev=size(FV.vertices);
% Check size of vertice array
if((sizev(2)~=3)||(length(sizev)~=2))
    error('patchnormals:inputs','The vertice list is not a m x 3 array')
end

sizef=size(FV.faces);
% Check size of vertice array
if((sizef(2)~=3)||(length(sizef)~=2))
    error('patchnormals:inputs','The vertice list is not a m x 3 array')
end

% Check if vertice indices exist
if(max(FV.faces(:))>size(FV.vertices,1))
    error('patchnormals:inputs','The face list contains an undefined vertex index')
end

% Check if vertice indices exist
if(min(FV.faces(:))<1)
    error('patchnormals:inputs','The face list contains an vertex index smaller then 1')
end

[Nx,Ny,Nz]=patchnormals_double(double(FV.faces(:,1)),double(FV.faces(:,2)),double(FV.faces(:,3)),double(FV.vertices(:,1)),double(FV.vertices(:,2)),double(FV.vertices(:,3)));

N=zeros(length(Nx),3);
N(:,1)=Nx;
N(:,2)=Ny;
N(:,3)=Nz;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef VademecumComputerAndPlotter < handle

    properties (GetAccess = public, SetAccess = private)
        vademecumData
    end

    properties (Access = private)
        fileName
        outPutPath
        vademecumPlotter
        vadVariables
    end

    methods (Access = public)

        function obj = VademecumComputerAndPlotter(d)
            obj.init(d)
        end

        function compute(obj)
            obj.calculateVademecumPlotterData();
            obj.plotVademecum();
            obj.printData();
        end

        function fI = getFeasibleIndex(obj)
            vp = obj.vademecumPlotter;
            fI = vp.feasibleIndex;
        end
    end

    methods (Access = private)

        function init(obj,d)
            obj.fileName   = d.fileName;
            obj.outPutPath = d.outPutPath;
        end

        function calculateVademecumPlotterData(obj)
            obj.calculateVademecumData();
            obj.transformVademecumCellDataInMatrix();
            d = obj.vademecumData;
            d.fileName   = obj.fileName;
            d.outPutPath = obj.outPutPath;
            obj.vademecumData = d;
        end

        function calculateVademecumData(obj)
            d.fileName = obj.fileName;
            d.outPutPath = obj.outPutPath;
            vc = VademecumPinvComputer(d);
            vc.compute();
            obj.vadVariables = vc.vademecumData;
        end

        function transformVademecumCellDataInMatrix(obj)
            mxV = obj.vadVariables.domVariables.mxV;
            myV = obj.vadVariables.domVariables.myV;
            for imx = 1:length(mxV)
                for imy = 1:length(myV)
                    vad = obj.vadVariables.variables{imx,imy};
                    variables.C(:,:,imx,imy)     = vad.Ctensor;
                    variables.invP(:,:,imx,imy)  = vad.Pinvtensor;
                    variables.Ptensor(:,:,imx,imy) = vad.Ptensor;
                    variables.volume(imx,imy)    = vad.volume;
                end
            end
            obj.vademecumData = variables;
            obj.vademecumData.mxV = mxV;
            obj.vademecumData.myV = myV;
        end

        function plotVademecum(obj)
            d = obj.vademecumData;
            obj.vademecumPlotter = VademecumPlotters(d);
            obj.vademecumPlotter.plot();
        end

        function printData(obj)
            d.outPutPath = [obj.outPutPath,obj.fileName];
            d.Ctensor  = obj.vademecumData.C;
            d.Ptensor  = obj.vademecumData.Ptensor;
            d.volume   = obj.vademecumData.volume;
            p = VademecumDataPrinter(d);
           % p.print();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [L,S,G,error,time]=GreBackground(X,rank,tau,power,isize,tol,k)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Background Modeling by Greedy Semi-Soft GoDec (GreBsmo)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%isize: image size of each frame, for example, 1024x768
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Tianyi Zhou, 2013, All rights reserved.

[L,S,error,time]=GreGoDec(X,rank,tau,tol,power,k);
G=X-L-S;

figure;
for i=1:min([200,size(X,2)])
    subplot(2,2,1);imagesc(reshape(X(:,i),isize));colormap(gray);axis image;axis off;title('X(Sample)');
    subplot(2,2,2);imagesc(reshape(L(:,i),isize));colormap(gray);axis image;axis off;title('L(Low-rank)');
    %subplot(1,4,3);imagesc(reshape(X(:,i).*(~~S(:,i)),isize));colormap(gray);axis image;axis off;title('S(Sparse)');
    subplot(2,2,4);imagesc(reshape(S(:,i),isize));colormap(gray);axis image;axis off;title('S(Sparse)');
    subplot(2,2,3);imagesc(reshape(G(:,i),isize));colormap(gray);axis image;axis off;title('G(Noise)');
    pause(0.1);
    if i<10
        print('-djpeg', '-r250', sprintf('%s%s','00',num2str(i)));
    elseif i<100
        print('-djpeg', '-r250', sprintf('%s%s','0',num2str(i)));
    else
        print('-djpeg', '-r250', sprintf('%s',num2str(i)));
    end
end

author: andrewssobral

repo: lrslibrary


text:
function [thetamat, dthetamat] ...
         = ForwardDynamicsTrajectory(thetalist, dthetalist, taumat, g, ...
                                     Ftipmat, Mlist, Glist, Slist, dt, ...
                                     intRes)
% *** CHAPTER 8: DYNAMICS OF OPEN CHAINS ***
% Takes thetalist: n-vector of initial joint variables,
%       dthetalist: n-vector of initial joint rates,
%       taumat: An N x n matrix of joint forces/torques, where each row is
%               the joint effort at any time step,
%       g: Gravity vector g,
%       Ftipmat: An N x 6 matrix of spatial forces applied by the
%                end-effector (If there are no tip forces, the user should
%                input a zero and a zero matrix will be used),
%       Mlist: List of link frames {i} relative to {i-1} at the home
%              position,
%       Glist: Spatial inertia matrices Gi of the links,
%       Slist: Screw axes Si of the joints in a space frame, in the format
%              of a matrix with the screw axes as the columns,
%       dt: The timestep between consecutive joint forces/torques,
%       intRes: Integration resolution is the number of times integration
%               (Euler) takes places between each time step. Must be an
%               integer value greater than or equal to 1.
% Returns thetamat: The N x n matrix of robot joint angles resulting from
%                   the specified joint forces/torques,
%         dthetamat: The N x n matrix of robot joint velocities.
% This function simulates the motion of a serial chain given an open-loop
% history of joint forces/torques. It calls a numerical integration
% procedure that uses ForwardDynamics.
% Example Inputs (3 Link Robot):
%
% clc; clear;
% thetalist = [0.1; 0.1; 0.1];
% dthetalist = [0.1; 0.2; 0.3];
% taumat = [[3.63, -6.58, -5.57]; [3.74, -5.55, -5.5]; ...
%         [4.31, -0.68, -5.19]; [5.18, 5.63, -4.31]; ...
%         [5.85, 8.17, -2.59]; [5.78, 2.79, -1.7]; ...
%         [4.99, -5.3, -1.19]; [4.08, -9.41, 0.07]; ...
%         [3.56, -10.1, 0.97]; [3.49, -9.41, 1.23]];
% %Initialise robot description (Example with 3 links)
% g = [0; 0; -9.8];
% Ftipmat = ones(size(taumat, 1), 6);
% M01 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.089159]; [0, 0, 0, 1]];
% M12 = [[0, 0, 1, 0.28]; [0, 1, 0, 0.13585]; [-1, 0 ,0, 0]; [0, 0, 0, 1]];
% M23 = [[1, 0, 0, 0]; [0, 1, 0, -0.1197]; [0, 0, 1, 0.395]; [0, 0, 0, 1]];
% M34 = [[1, 0, 0, 0]; [0, 1, 0, 0]; [0, 0, 1, 0.14225]; [0, 0, 0, 1]];
% G1 = diag([0.010267, 0.010267, 0.00666, 3.7, 3.7, 3.7]);
% G2 = diag([0.22689, 0.22689, 0.0151074, 8.393, 8.393, 8.393]);
% G3 = diag([0.0494433, 0.0494433, 0.004095, 2.275, 2.275, 2.275]);
% Glist = cat(3, G1, G2, G3);
% Mlist = cat(3, M01, M12, M23, M34);
% Slist = [[1; 0; 1;      0; 1;     0], ...
%        [0; 1; 0; -0.089; 0;     0], ...
%        [0; 1; 0; -0.089; 0; 0.425]];
% dt = 0.1;
% intRes = 8;
% [thetamat, dthetamat] ...
% = ForwardDynamicsTrajectory(thetalist, dthetalist, taumat, g, ...
%                           Ftipmat, Mlist, Glist, Slist, dt, intRes);
% %Output using matplotlib to plot the joint forces/torques
% Tf = size(taumat, 1);
% time=0: (Tf / size(thetamat, 1)): (Tf - (Tf / size(thetamat, 1)));
% plot(time,thetamat(:, 1),'b')
% hold on
% plot(time,thetamat(:, 2), 'g')
% plot(time,thetamat(:, 3), 'r')
% plot(time,dthetamat(:, 1), 'c')
% plot(time,dthetamat(:, 2), 'm')
% plot(time,dthetamat(:, 3), 'y')
% title('Plot of Joint Angles and Joint Velocities')
% xlabel('Time')
% ylabel('Joint Angles/Velocities')
% legend('Theta1', 'Theta2', 'Theta3', 'DTheta1', 'DTheta2', 'DTheta3')
%

taumat = taumat';
Ftipmat = Ftipmat';
thetamat = taumat;
thetamat(:, 1) = thetalist;
dthetamat = taumat;
dthetamat(:, 1) = dthetalist;
for i = 1: size(taumat, 2) - 1
    for j = 1: intRes
       ddthetalist ...
       = ForwardDynamics(thetalist, dthetalist, taumat(:,i), g, ...
                         Ftipmat(:, i), Mlist, Glist, Slist);
       [thetalist, dthetalist] = EulerStep(thetalist, dthetalist, ...
                                           ddthetalist, dt / intRes);
    end
    thetamat(:, i + 1) = thetalist;
    dthetamat(:, i + 1) = dthetalist;
end
thetamat = thetamat';
dthetamat = dthetamat';
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function assertAlmostEqual(A, B, reltol, message)
%assertEqual Assert that inputs are equal within relative tolerance
%   assertEqual(A, B, RELTOL) throws an exception of any of the values in A and
%   B are not equal within the specified tolerance.  NaN values are considered
%   to be equal.  A and B have to have the same class and sparsity to be
%   considered equal.
%
%   assertEqual(A, B) uses the following relative tolerance value:
%
%       100 * eps(class(A))
%
%   assertEqual(A, B, RELTOL, MESSAGE) uses the specified message string when
%   throwing the exception.  With this syntax, use RELTOL = [] to specify the
%   default relative tolerance.
%
%   Note that if either A or B are not floating-point arrays, then A and B are
%   compared using ISEQUALWITHEQUALNANS and the relative tolerance value is not
%   used.
%
%   Examples
%   --------
%   % This call returns silently.
%   assertAlmostEqual(1.0, 1.0 + eps);
%
%   % This call throws an error.
%   assertAlmostEqual(1.0, 1.1);
%
%   See also assertEqual, mtest.utils.isAlmostEqual

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

if ~(issparse(A) == issparse(B))
   throw(MException('assertAlmostEqual:sparsityNotEqual', message));
end

if ~strcmp(class(A), class(B))
   throw(MException('assertAlmostEqual:classNotEqual', message));
end

if nargin < 3 || isempty(reltol)
    reltol = 100 * eps(class(A));
end

if nargin < 4
    message = sprintf('Inputs are not equal within relative tolerance: %g', ...
        reltol);
end

if ~mtest.utils.isAlmostEqual(A, B, reltol)
   throw(MException('assertAlmostEqual:tolExceeded', message));
end

author: ricedsp


text:
function value = index3_col ( i_min, i, i_max, j_min, j, j_max, ...
  k_min, k, k_max, index_min )

%*****************************************************************************80
%
%% INDEX3_COL indexes a 3D array by columns.
%
%  Discussion:
%
%    Entries of the array are indexed starting at entry (I_MIN,J_MIN,K_MIN),
%    and increasing the row index first, then the column index.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I_MIN, I, I_MAX, for row indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer J_MIN, J, J_MAX, for column indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer K_MIN, K, K_MAX, for plane indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer INDEX_MIN, the index of (I_MIN,J_MIN,K_MIN).
%    Typically, this is 0 or 1.
%
%    Output, integer VALUE, the index of element (I,J,K).
%
  value = index_min ...
             + ( i - i_min ) ...
             + ( j - j_min ) * ( i_max + 1 - i_min ) ...
             + ( k - k_min ) * ( j_max + 1 - j_min ) * ( i_max + 1 - i_min );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% compute the source zones for landers

report_this_filefun(mfilename('fullpath'));

% compute the source zones for landers

% define Input parameters
if ~exist('nb', 'var'); nb = 10; end
if ~exist('tmin1', 'var'); tmin1 = 0.05; end
if ~exist('tpre', 'var'); tpre = 30; end
if ~exist('Tobs', 'var'); Tobs = 30; end
if ~exist('prol', 'var'); prol = 0.5; end
if isempty(storedcat); storedcat=a; end

def = {num2str(nb),num2str(tmin1),num2str(tpre),num2str(Tobs), num2str(prol)};

tit ='Hazard analysis input parameters';
prompt={ 'Number of source zones',...
    'Minimum Time (tmin1) [days]',...
    'Prediction Time (tpre) [days]',...
    'Data period considered (Tobs) [days] ',...
    'Probability level (prol)',...
    };

ni2 = inputdlg(prompt,tit,1,def);

l = ni2{1}; nb= str2double(l);
l = ni2{2}; tmin1= str2double(l);
l = ni2{3}; tpre= str2double(l);
l = ni2{4}; Tobs = str2double(l);
l = ni2{5}; prol= str2double(l);

tlen = Tobs;
l=storedcat.Date < min(maepi.Date) + days(Tobs);
a=storedcat.subset(l);
da = []; anz = [];
B = [];
dt = 1;
update(mainmap())

cd /home2/stefan/ZMAP/aspar
fid2 = fopen('sourczones.txt','w');

l1 = [];

for kk = 1:nb-1
    l1 = [l1 prctile2(a.Latitude,kk*100/nb)];
end

bo = [34.35 l1 34.88 ];
for ii = 1:length(bo)-1
    l = a.Latitude >= bo(ii) & a.Latitude < bo(ii+1);
    figure_w_normalized_uicontrolunits(map)
    hold on
    plot(a(l,1),a(l,2),'k+')

    b1 = prctile2(a(l,1),25);
    b2 = prctile2(a(l,1),75);
    plot(b1,bo(ii),'+r','Markersize',12)
    plot(b2,bo(ii),'+r','Markersize',12)
    plot(b1,bo(ii+1),'+r','Markersize',12)
    plot(b2,bo(ii+1),'+r','Markersize',12)
    drawnow

    newt2 = a.subset(l);
    calcp
    B = [B ; b1 bo(ii)  b2 bo(ii)];

    anz = [];
    for m2 = 4.25:0.5:6.25
        M2 = maepi(1,6) - m2;
        t0 = days(max(a.Date) - min(maepi));
        pla = 0; pla2 = 0;

        for t = t0:dt:t0+tpre
            pla = pla + (10^(A + bv*(M2)) * (t + c)^(-p))  *dt;
            pla2 = pla2 + (10^(A + bv*(M2-0.5)) * (t + c)^(-p))  *dt;
        end

        %fac = length(a(l,1))/a.Count
        anz = [anz ;  m2+0.25  (pla-pla2)/tpre];

    end

    % write info to file
    s = ['0    1.     -1          zn035.00']; s = s';
    fprintf(fid2,'%s\n',s);
    s = ['2 1 1']; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(-b1,5) ' ' num2str(bo(ii)) ' ' num2str(-b2,5) ' ' num2str(bo(ii))]; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(-b1,5) ' ' num2str(bo(ii+1)) ' ' num2str(-b2,5) ' ' num2str(bo(ii+1))]; s = s';
    fprintf(fid2,'%s\n',s);
    fprintf(fid2,'%7.6f    ',anz(:,2));
    fprintf(fid2,'\n');
    fprintf(fid2,'\n');
    fprintf(fid2,'%3.2f  ',anz(:,1));
    fprintf(fid2,'\n');
end

fclose(fid2)

%return
cd /home2/stefan/ZMAP/aspar
do = [' ! cat head5.txt | sed -e"s/sub1/' num2str(prol) ' 1 ' num2str(tpre) '/" > head2.txt ' ]; eval(do)
do = [' ! cat head2.txt sourczones.txt tail.txt > /home2/stefan/srisk/myrisk.inp']; eval(do)
cd /home2/stefan/srisk/

do = [ '! /nfs/alaska/home2/stefan/srisk/seis4b.exe  myrisk.inp myrisk.out f2 f3' ]; eval(do)
do = [' !cat myrisk.out | grep -e "LAT   "  -e "' num2str(tpre) ' YE"  > tmp2 ']; eval(do)
try
    condata2
catch ME
    error_handler(ME, @do_nothing);
end

[X2,Y2] = meshgrid((-116.8:0.02:-115.7),(33.92:0.01:35.3));

Z = griddata(-da(:,1),da(:,2),da(:,3),X2,Y2,'linear');
figure

axes('pos',[0.05 0.1 0.4 0.7])
pcolor(X2,Y2,Z);
set(gca,'TickDir','out');
colorbar; shading flat
hold on
shading interp
overlay
[ca1, ca2] = caxis;

cd /home2/stefan/ZMAP

mati = decyear(min(maepi.Date));
l = obs(:,3) > mati + tlen/365 &  obs(:,3) <  mati+ tlen/365 + tpre/365;
b = obs(l,:);
l = b(:,6) > 3;
b = b(l,:);

ve = [];dx = 0.02; dy = 0.02

for x = -116.8:dx:-115.7
    for y = 33.92:dy:35.3
        ve =    [ ve ; x y ];
    end
end

le = length(ve);
Y0 = zeros(le,1);

for i = 1:length(b)
    di2 = deg2km((distance(ve(:,2),ve(:,1),repmat(b(i,2),le,1),repmat(b(i,1),le,1))));
    R = di2;
    r = sqrt(R.^2 + 5.57^2);
    M = b(i,6);
    Y = -0.136 + 0.229*(M-6) - 0.778 * log10(r) + 0.161*1 ;
    Y = 10.^Y;
    c = [Y , Y0];
    mapga = max(c');
    Y0 = mapga';
end

mapga = mapga';

l1 = length(-116.8:dx:-115.7);
l2 = length(33.92:dy:35.3);

axes('pos',[0.55 0.1 0.4 0.7])
re = reshape(mapga,l2,l1);
rey = reshape(ve(:,2),l2,l1);
rex = reshape(ve(:,1),l2,l1);
pcolor(rex,rey,re);
set(gca,'YTicklabels',[],'TickDir','out');
caxis([ca1 ca2]);
shading interp

colorbar;
overlay
hold on

author: CelsoReyes


text:
function varargout=bss_decomp_tvfilt(varargin)

% decompose an estimated source into target/interference/noise/artefacts components, assuming the admissible distortion is a time-varying filter.
%
% Usage:
%
% [s_target,e_interf[,e_noise],e_artif]=bss_decomp_tvfilt(se,index,S[,N],tvshape,tvstep,L)
%
% Input:
%   - se: row vector of length T containing the estimated source,
%   - index: points which component of S se has to be compared to,
%   - S: n x T matrix containing the original sources,
%   - N: m x T matrix containing the noise on the obseravtions (if any).
%   - tvshape : row vector of length V at most T containing the shape of the elementary
%     allowed time variations of the filter coefficients
%   - tvstep  : hop size (in number of samples) between two consecutive
%     variations of the filter coefficients
%   - L: the number of lags
%
% Output:
%   - s_target: row vector of length T containing the target source(s)
%   contribution,
%   - e_interf: row vector of length T containing the interferences
%   contribution,
%   - e_noise: row vector of length T containing the noise contribution (if
%   any),
%   - e_artif: row vector of length T containing the artifacts
%   contribution.
%
% Developers:  - Cedric Fevotte (cf269@cam.ac.uk) - Emmanuel Vincent
% (vincent@ircam.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

se=varargin{1}; index=varargin{2}; S=varargin{3};

switch nargin
    case 5
        N=[]; tvshape = varargin{4}; tvstep = varargin{5}; L = varargin{6};
    case 6
        N=varargin{4}; tvshape = varargin{5}; tvstep = varargin{6}; L = varargin{7};
    otherwise
        disp('Wrong number of arguments.')
end

[ne,Te]=size(se);
[n,T]=size(S);

%%%%%%%%%% WARNINGS %%%%%%%%%%%%%
switch isempty(N)
    case 1
        if n>T | ne>Te, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: se and S have different lengths.'), return; end
    case 0
        [m,Tm]=size(N);
        if n>T | ne>Te | m>Tm, disp('Watch out: signals must be in rows.'), return; end
        if ne~=1, disp('Watch out: se must contain only one row.'), return; end
        if T~=Te, disp('Watch out: S and Se have different lengths.'), return; end
        if T~=Tm, disp('Watch out: N, S and Se have different lengths.'), return; end
end

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

% Create the space of target source(s)
target_space = bss_make_lags(S(index,:),L);
% Create the space of sources
sources_space= bss_make_lags(S,L);
% Create the noise space
noise_space  = bss_make_lags(N,L);

s_target=zeros(1,T);
e_interf=zeros(1,T);
e_artif=zeros(1,T);
if isempty(noise_space)==0, e_noise=zeros(1,T); end

%%% Target source(s) contribution %%%
s_target = bss_tvproj(se,target_space,tvshape,tvstep);

%%% Interferences contribution %%%
P_S_se = bss_tvproj(se,[sources_space],tvshape,tvstep);
e_interf = P_S_se - s_target;

switch isempty(noise_space)
    case 1 % No noise
        %%% Artifacts contribution %%%
        e_artif= se - P_S_se;

        %%% Output %%%
        varargout{1}=s_target;
        varargout{2}=e_interf;
        varargout{3}=e_artif;

    case 0 % Noise
        %%% Noise contribution %%%
        P_SN_se= bss_tvproj(se,[sources_space;noise_space],tvshape,tvstep);
        e_noise=P_SN_se-P_S_se;

        %%% Artifacts contribution %%%
        e_artif=se-P_SN_se;

        %%% Output %%%
        varargout{1}=s_target;
        varargout{2}=e_interf;
        varargout{3}=e_noise;
        varargout{4}=e_artif;
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function [Xpoints1, Ypoints1, Xpoints2, Ypoints2, tms, yc] = ProcessVideo(videoFilePath)
if nargin < 1
    videoFilePath = fullfile(pwd, 'video/video.avi');
end
% clc; clear all; close all;
% videoFilePath = fullfile(pwd, 'video/video.avi');
time_start = cputime;
[pathstr, name, ext] = fileparts(videoFilePath);
foldername = fullfile(pwd, sprintf('%s_images', name));
T = 1;
P = 5;
W1 = [75 95];
L1 = [360 17];
W2 = [55 55];
L2 = [35 1565];
Xpoints1 = [];
Ypoints1 = [];
Xpoints2 = [];
Ypoints2 = [];
Xpointst1 = [];
Ypointst1 = [];
Xpointst2 = [];
Ypointst2 = [];
figure('Position', get(0, 'ScreenSize'));
hg1 = subplot(1, 2, 1);
hg2 = subplot(1, 2, 2);
for frame = 1:146
    filename = fullfile(foldername, sprintf('%04d.jpg', frame));
    R = imread(filename);
    Imi = R;
    xc1 = 0;
    yc1 = 0;
    xc2 = 0;
    yc2 = 0;
    if frame > 72
        I = rgb2hsv(Imi);
        I = I(:,:,1);
        I = roicolor(I, 0.1, 0.17);
        MeanConverging1 = 1;
        while MeanConverging1
            M00 = 0.0;
            for i = L1(1)-P : (L1(1)+W1(1)+P)
                for j = L1(2)-P : (L1(2)+W1(2)+P)
                    if i > size(I,1) || j > size(I,2) || i < 1 || j < 1
                        continue;
                    end
                    M00 = M00 + double(I(i,j));
                end
            end
            M10 = 0.0;
            for i = L1(1)-P : (L1(1)+W1(1)+P)
                for j = L1(2)-P : (L1(2)+W1(2)+P)
                    if i > size(I,1) || j > size(I,2) || i < 1 || j < 1
                        continue;
                    end
                    M10 = M10 + i * double(I(i,j));
                end
            end
            M01 = 0.0;
            for i = L1(1)-P : (L1(1)+W1(1)+P)
                for j = L1(2)-P : (L1(2)+W1(2)+P)
                    if i > size(I,1) || j > size(I,2)|| i < 1 || j < 1
                        continue;
                    end
                    M01 = M01 + j * double(I(i,j));
                end
            end
            xc1 = round(M10 / M00);
            yc1 = round(M01 / M00);
            oldL = L1;
            L1 = [floor(xc1 - (W1(1)/2)) floor(yc1 - (W1(2)/2))];
            if abs(oldL(1)-L1(1)) < T || abs(oldL(2)-L1(2)) < T
                MeanConverging1 = 0;
            end
        end
        s = round(1.1 * sqrt(M00));
        W1 = [ s floor(1.2*s) ];
        L1 = [floor(xc1 - (W1(1)/2)) floor(yc1 - (W1(2)/2))];
        Xpoints1 = [Xpoints1 xc1];
        Ypoints1 = [Ypoints1 yc1];
        yc1t = yc1+randi(2,1,1)*25;
        xc1t = xc1+randi(2,1,1)*25;
        Xpointst1 = [Xpointst1 xc1t];
        Ypointst1 = [Ypointst1 yc1t];
    else
        Xpoints1 = [Xpoints1 NaN];
        Ypoints1 = [Ypoints1 NaN];
        Xpointst1 = [Xpointst1 NaN];
        Ypointst1 = [Ypointst1 NaN];
    end
    if frame > 90 && frame < 146
        R = Imi;
        I = rgb2ycbcr(R);
        I = I(:,:,1);
        I = mat2gray(I);
        I = roicolor(I, 0.05, 0.3);
        MeanConverging2 = 1;
        while MeanConverging2
            M00 = 0.0;
            M00 = 0.0;
            for i = L2(1)-P : (L2(1)+W2(1)+P),
                for j = L2(2)-P : (L2(2)+W2(2)+P),
                    if i > size(I,1) || j > size(I,2) || i < 1 || j < 1
                        continue;
                    end
                    M00 = M00 + double(I(i,j));
                end
            end
            M10 = 0.0;
            for i = L2(1)-P : (L2(1)+W2(1)+P),
                for j = L2(2)-P : (L2(2)+W2(2)+P),
                    if i > size(I,1) || j > size(I,2) || i < 1 || j < 1
                        continue;
                    end
                    M10 = M10 + i * double(I(i,j));
                end
            end
            M01 = 0.0;
            for i = L2(1)-P : (L2(1)+W2(1)+P),
                for j = L2(2)-P : (L2(2)+W2(2)+P),
                    if i > size(I,1) || j > size(I,2)|| i < 1 || j < 1
                        continue;
                    end
                    M01 = M01 + j * double(I(i,j));
                end
            end
            xc2 = round(M10 / M00);
            yc2 = round(M01 / M00);
            oldL = L2;
            L2 = [floor(xc2 - (W2(1)/2)) floor(yc2 - (W2(2)/2))];
            if abs(oldL(1)-L2(1)) < T || abs(oldL(2)-L2(2)) < T
                MeanConverging2 = 0;
            end
        end
        s = round(1.1 * sqrt(M00));
        W2 = [ s      floor(1.2*s) ];
        L2 = [floor(xc2 - (W2(1)/2)) floor(yc2 - (W2(2)/2))];
        Xpoints2 = [Xpoints2 xc2];
        Ypoints2 = [Ypoints2 yc2];
        yc2t = yc2+randi(2,1,1)*25;
        xc2t = xc2+randi(2,1,1)*25;
        Xpointst2 = [Xpointst2 xc2t];
        Ypointst2 = [Ypointst2 yc2t];
    else
        Xpoints2 = [Xpoints2 NaN];
        Ypoints2 = [Ypoints2 NaN];
        Xpointst2 = [Xpointst2 NaN];
        Ypointst2 = [Ypointst2 NaN];
    end
    axes(hg1); cla;
    imshow(Imi, []); hold on;
    if xc1 > 0 && yc1 > 0
        plot(yc1, xc1, 'go', 'MarkerFaceColor', 'g');
        plot(yc1t, xc1t, 'g+', 'MarkerFaceColor', 'g');
    end
    if xc2 > 0 && yc2 > 0
        plot(yc2, xc2, 'bo', 'MarkerFaceColor', 'b');
        plot(yc2t, xc2t, 'b+', 'MarkerFaceColor', 'b');
    end
    hold off; title(sprintf('%04d֡', frame));
    bg = true(size(Imi,1), size(Imi,2));
    axes(hg2); cla; imshow(bg);
    hold on; box on;
    plot(Ypoints1, Xpoints1, 'go-', 'MarkerFaceColor', 'g');
    plot(Ypoints2, Xpoints2, 'bo-', 'MarkerFaceColor', 'b');
    hold off; title(sprintf('%04d֡', frame));
    pause(0.001);
end
time_end = cputime;
tms = time_end - time_start;
yc.Xpointst1 = Xpointst1;
yc.Ypointst1 = Ypointst1;
yc.Xpointst2 = Xpointst2;
yc.Ypointst2 = Ypointst2;

author: decouples


text:
function spm_preproc_write(p,opts)
% Write out VBM preprocessed data
% FORMAT spm_preproc_write(p,opts)
% p    - results from spm_prep2sn
% opts - writing options.  A struct containing these fields:
%        biascor - write bias corrected image
%        cleanup - level of brain segmentation cleanup
%        GM      - flags for which images should be written
%        WM      - similar to GM
%        CSF     - similar to GM
%__________________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_preproc_write.m 6894 2016-09-30 16:48:46Z spm $

if ischar(p), p = cellstr(p); end

if nargin==1
    opts = spm_get_defaults('old.preproc.output');
end

for i=1:numel(p)
    if iscellstr(p), q = load(p{i}); else q = p(i); end
    if i==1, b0 = spm_load_priors(q.VG); end
    preproc_apply(q,opts,b0);
end

%==========================================================================
% function preproc_apply(p,opts,b0)
%==========================================================================
function preproc_apply(p,opts,b0)

nclasses = size(fieldnames(opts),1) - 2 ;
switch nclasses
    case 3
        sopts = [opts.GM ; opts.WM ; opts.CSF];
    case 4
        sopts = [opts.GM ; opts.WM ; opts.CSF ; opts.EXTRA1];
    case 5
        sopts = [opts.GM ; opts.WM ; opts.CSF ; opts.EXTRA1 ; opts.EXTRA2];
    otherwise
        error('Unsupported number of classes.')
end

T    = p.flags.Twarp;
bsol = p.flags.Tbias;
d2   = [size(T) 1];
d    = p.VF.dim(1:3);

[x1,x2,o] = ndgrid(1:d(1),1:d(2),1);
x3  = 1:d(3);
d3  = [size(bsol) 1];
B1  = spm_dctmtx(d(1),d2(1));
B2  = spm_dctmtx(d(2),d2(2));
B3  = spm_dctmtx(d(3),d2(3));
bB3 = spm_dctmtx(d(3),d3(3),x3);
bB2 = spm_dctmtx(d(2),d3(2),x2(1,:)');
bB1 = spm_dctmtx(d(1),d3(1),x1(:,1));

mg  = p.flags.mg;
mn  = p.flags.mn;
vr  = p.flags.vr;
K   = length(p.flags.mg);
Kb  = length(p.flags.ngaus);

for k1=1:size(sopts,1)
    dat{k1}    = zeros(d(1:3),'uint8');
    if sopts(k1,3)
        Vt     = struct('fname',   spm_file(p.VF.fname,'prefix',['c', num2str(k1)]),...
                        'dim',     p.VF.dim,...
                        'dt',      [spm_type('uint8') spm_platform('bigend')],...
                        'pinfo',   [1/255 0 0]',...
                        'mat',     p.VF.mat,...
                        'n',       [1 1],...
                        'descrip', ['Tissue class ' num2str(k1)]);
        Vt     = spm_create_vol(Vt);
        VO(k1) = Vt;
    end
end
if opts.biascor
    VB = struct('fname',   spm_file(p.VF.fname,'prefix','m'),...
                'dim',     p.VF.dim(1:3),...
                'dt',      [spm_type('float32') spm_platform('bigend')],...
                'pinfo',   [1 0 0]',...
                'mat',     p.VF.mat,...
                'n',       [1 1],...
                'descrip', 'Bias Corrected');
    VB = spm_create_vol(VB);
end

lkp = []; for k=1:Kb, lkp = [lkp ones(1,p.flags.ngaus(k))*k]; end

spm_progress_bar('init',length(x3),['Working on ' spm_file(p.VF.fname,'basename')],'Planes completed');
M = p.VG(1).mat\p.flags.Affine*p.VF.mat;

for z=1:length(x3)

    % Bias corrected image
    f          = spm_sample_vol(p.VF,x1,x2,o*x3(z),0);
    cr         = exp(transf(bB1,bB2,bB3(z,:),bsol)).*f;
    if opts.biascor
        % Write a plane of bias corrected data
        VB = spm_write_plane(VB,cr,z);
    end

    if any(sopts(:))
        msk        = (f==0) | ~isfinite(f);
        [t1,t2,t3] = defs(T,z,B1,B2,B3,x1,x2,x3,M);
        q          = zeros([d(1:2) Kb]);
        bt         = zeros([d(1:2) Kb]);
        for k1=1:Kb
            bt(:,:,k1) = spm_sample_priors(b0{k1},t1,t2,t3,k1==Kb);
        end
        b = zeros([d(1:2) K]);
        for k=1:K
            b(:,:,k) = bt(:,:,lkp(k))*mg(k);
        end
        s = sum(b,3);
        for k=1:K
            p1            = exp((cr-mn(k)).^2/(-2*vr(k)))/sqrt(2*pi*vr(k)+eps);
            q(:,:,lkp(k)) = q(:,:,lkp(k)) + p1.*b(:,:,k)./s;
        end
        sq = sum(q,3)+eps;
        sw = warning('off','MATLAB:divideByZero');
        for k1=1:size(sopts,1)
            tmp            = q(:,:,k1);
            tmp(msk)       = 0;
            tmp            = tmp./sq;
            dat{k1}(:,:,z) = uint8(round(255 * tmp));
        end
        warning(sw);
    end
    spm_progress_bar('set',z);
end
spm_progress_bar('clear');

if opts.cleanup > 0
    [dat{1},dat{2},dat{3}] = clean_gwc(dat{1},dat{2},dat{3}, opts.cleanup);
end
if any(sopts(:,3))
    for z=1:length(x3)
        for k1=1:size(sopts,1)
            if sopts(k1,3)
                tmp = double(dat{k1}(:,:,z))/255;
                spm_write_plane(VO(k1),tmp,z);
            end
        end
    end
end

for k1=1:size(sopts,1)
    if any(sopts(k1,1:2))
        so      = struct('wrap',[0 0 0],...
                         'interp',1,...
                         'vox',[NaN NaN NaN],...
                         'bb',ones(2,3)*NaN,...
                         'preserve',0);
        ovx     = abs(det(p.VG(1).mat(1:3,1:3)))^(1/3);
        fwhm    = max(ovx./sqrt(sum(p.VF.mat(1:3,1:3).^2))-1,0.1);
        dat{k1} = decimate(dat{k1},fwhm);
        dim     = [size(dat{k1}) 1];
        VT      = struct('fname', spm_file(p.VF.fname,'prefix',['c', num2str(k1)]),...
                         'dim',   dim(1:3),...
                         'dt',    [spm_type('uint8') spm_platform('bigend')],...
                         'pinfo', [1/255 0]',...
                         'mat',   p.VF.mat,...
                         'dat',   dat{k1});
        if sopts(k1,2)
            spm_write_sn(VT,p,so);
        end
        so.preserve = 1;
        if sopts(k1,1)
           VN       = spm_write_sn(VT,p,so);
           VN.fname = spm_file(p.VF.fname,'prefix',['mwc', num2str(k1)]);
           spm_write_vol(VN,VN.dat);
        end
    end
end

%==========================================================================
% function [x1,y1,z1] = defs(sol,z,B1,B2,B3,x0,y0,z0,M)
%==========================================================================
function [x1,y1,z1] = defs(sol,z,B1,B2,B3,x0,y0,z0,M)
x1a = x0    + transf(B1,B2,B3(z,:),sol(:,:,:,1));
y1a = y0    + transf(B1,B2,B3(z,:),sol(:,:,:,2));
z1a = z0(z) + transf(B1,B2,B3(z,:),sol(:,:,:,3));
x1  = M(1,1)*x1a + M(1,2)*y1a + M(1,3)*z1a + M(1,4);
y1  = M(2,1)*x1a + M(2,2)*y1a + M(2,3)*z1a + M(2,4);
z1  = M(3,1)*x1a + M(3,2)*y1a + M(3,3)*z1a + M(3,4);

%==========================================================================
% function t = transf(B1,B2,B3,T)
%==========================================================================
function t = transf(B1,B2,B3,T)
if ~isempty(T)
    d2 = [size(T) 1];
    t1 = reshape(reshape(T, d2(1)*d2(2),d2(3))*B3', d2(1), d2(2));
    t  = B1*t1*B2';
else
    t  = zeros(size(B1,1),size(B2,1),size(B3,1));
end

%==========================================================================
% function dat = decimate(dat,fwhm)
%==========================================================================
function dat = decimate(dat,fwhm)
% Convolve the volume in memory (fwhm in voxels).
lim = ceil(2*fwhm);
x  = -lim(1):lim(1); x = spm_smoothkern(fwhm(1),x); x  = x/sum(x);
y  = -lim(2):lim(2); y = spm_smoothkern(fwhm(2),y); y  = y/sum(y);
z  = -lim(3):lim(3); z = spm_smoothkern(fwhm(3),z); z  = z/sum(z);
i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;
spm_conv_vol(dat,dat,x,y,z,-[i j k]);

%==========================================================================
% function [g,w,c] = clean_gwc(g,w,c,level)
%==========================================================================
function [g,w,c] = clean_gwc(g,w,c,level)
if nargin<4, level = 1; end

b    = w;
b(1) = w(1);

% Build a 3x3x3 seperable smoothing kernel
%--------------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

th1 = 0.15;
if level==2, th1 = 0.2; end
% Erosions and conditional dilations
%--------------------------------------------------------------------------
niter = 32;
spm_progress_bar('Init',niter,'Extracting Brain','Iterations completed');
for j=1:niter
    if j>2, th=th1; else th=0.6; end  % Dilate after two its of erosion
    for i=1:size(b,3)
        gp       = double(g(:,:,i));
        wp       = double(w(:,:,i));
        bp       = double(b(:,:,i))/255;
        bp       = (bp>th).*(wp+gp);
        b(:,:,i) = uint8(round(bp));
    end
    spm_conv_vol(b,b,kx,ky,kz,-[1 1 1]);
    spm_progress_bar('Set',j);
end
th = 0.05;
for i=1:size(b,3)
    gp       = double(g(:,:,i))/255;
    wp       = double(w(:,:,i))/255;
    cp       = double(c(:,:,i))/255;
    bp       = double(b(:,:,i))/255;
    bp       = ((bp>th).*(wp+gp))>th;
    g(:,:,i) = uint8(round(255*gp.*bp./(gp+wp+cp+eps)));
    w(:,:,i) = uint8(round(255*wp.*bp./(gp+wp+cp+eps)));
    c(:,:,i) = uint8(round(255*(cp.*bp./(gp+wp+cp+eps)+cp.*(1-bp))));
end
spm_progress_bar('Clear');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [objfiles, timestamp_out] = cs_make (f)
%CS_MAKE compiles CSparse for use in MATLAB.
%   Usage:
%       cs_make
%       [objfiles, timestamp] = cs_make (f)
%
%   With no input arguments, or with f=0, only those files needing to be
%   compiled are compiled (like the Unix/Linux/GNU "make" command, but not
%   requiring "make").  If f is a nonzero number, all files are compiled.
%   If f is a string, only that mexFunction is compiled.  For example,
%   cs_make ('cs_add') just compiles the cs_add mexFunction.  This option is
%   useful when developing a single new mexFunction.  This function can only be
%   used if the current directory is CSparse/MATLAB/CSparse.  Returns a list of
%   the object files in CSparse, and the latest modification time of any source
%   codes.
%
%   To add a new function and its MATLAB mexFunction to CSparse:
%
%       (1) Create a source code file CSparse/Source/cs_mynewfunc.c.
%       (2) Create a help file, CSparse/MATLAB/CSparse/cs_mynewfunc.m.
%           This is very useful, but not strictly required.
%       (3) Add the prototype of cs_mynewfunc to CSparse/Include/cs.h.
%       (4) Create its MATLAB mexFunction, CSparse/MATLAB/cs_mynewfunc_mex.c.
%       (5) Edit cs_make.m, and add 'cs_mynewfunc' to the 'cs' and 'csm' lists.
%       (6) Type 'cs_make' in the CSparse/MATLAB/CSparse directory.
%           If all goes well, your new function is ready for use in MATLAB.
%
%       (7) Optionally add 'cs_mynewfunc' to CSparse/Source/Makefile
%           and CSparse/MATLAB/CSparse/Makefile, if you want to use the
%           Unix/Linux/GNU make command instead of cs_make.m.  See where
%           'cs_add' and 'cs_add_mex' appear in those files, and add
%           'cs_mynewfunc' accordingly.
%       (8) Optionally add 'cs_mynewfunc' to Tcov/Makefile, and add additional
%           test code to cs_test.c, and add MATLAB test code to MATLAB/Test/*.
%
%   Example:
%       cs_make                  % compile everything
%       cs_make ('cs_chol') ;    % just compile cs_chol mexFunction
%
%   See also MEX.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

fprintf ('Compiling CSparse\n') ;

% CSparse source files, in ../../Source, such as ../../Source/cs_add.c.
% Note that not all CSparse source files have their own mexFunction.
cs = { 'cs_add', 'cs_amd', 'cs_chol', 'cs_cholsol', 'cs_counts', ...
    'cs_cumsum', 'cs_dfs', 'cs_dmperm', 'cs_droptol', 'cs_dropzeros', ...
    'cs_dupl', 'cs_entry', 'cs_etree', 'cs_fkeep', 'cs_gaxpy', 'cs_happly', ...
    'cs_house', 'cs_ipvec', 'cs_load', 'cs_lsolve', 'cs_ltsolve', 'cs_lu', ...
    'cs_lusol', 'cs_malloc', 'cs_maxtrans', 'cs_multiply', 'cs_norm', ...
    'cs_permute', 'cs_pinv', 'cs_post', 'cs_print', 'cs_pvec', 'cs_qr', ...
    'cs_qrsol', 'cs_scatter', 'cs_scc', 'cs_schol', 'cs_sqr', 'cs_symperm', ...
    'cs_tdfs', 'cs_transpose', 'cs_compress', 'cs_updown', 'cs_usolve', ...
    'cs_utsolve', 'cs_util', 'cs_reach', 'cs_spsolve', 'cs_ereach', ...
    'cs_leaf', 'cs_randperm' } ;
    % add cs_mynewfunc to the above list

details = 1 ;
kk = 0 ;
csm = { } ;
if (nargin == 0)
    force = 0 ;
elseif (ischar (f))
    fprintf ('cs_make: compiling ../../Source files and %s_mex.c\n', f) ;
    force = 0 ;
    csm = {f} ;
else
    force = f ;
    details = details | (force > 1) ;                                       %#ok
    if (force & details)                                                    %#ok
        fprintf ('cs_make: re-compiling everything\n') ;
    end
end
if (force)
    fprintf ('Compiling CSparse\n') ;
end

if (isempty (csm))
    % mexFunctions, of the form cs_add_mex.c, etc, in this directory
    csm = { 'cs_add', 'cs_amd', 'cs_chol', 'cs_cholsol', 'cs_counts', ...
        'cs_dmperm', 'cs_droptol', 'cs_etree', 'cs_gaxpy', 'cs_lsolve', ...
        'cs_ltsolve', 'cs_lu', 'cs_lusol', 'cs_multiply', 'cs_permute', ...
        'cs_print', 'cs_qr', 'cs_qrsol', 'cs_scc', 'cs_symperm', 'cs_thumb', ...
        'cs_transpose', 'cs_sparse', 'cs_updown', 'cs_usolve', ...
        'cs_utsolve', 'cs_randperm', 'cs_sqr' } ;
        % add cs_mynewfunc to the above list
end

try
    % ispc does not appear in MATLAB 5.3
    pc = ispc ;
catch
    % if ispc fails, assume we are on a Windows PC if it's not unix
    pc = ~isunix ;
end

if (pc)
    obj = '.obj' ;
else
    obj = '.o' ;
end

srcdir = '../../Source/' ;
hfile = '../../Include/cs.h' ;

% compile each CSparse source file
[anysrc timestamp kk] = compile_source ('', 'cs_mex', obj, hfile, force, ...
    kk, details) ;
CS = ['cs_mex' obj] ;
if (nargout > 0)
    objfiles = ['..' filesep 'CSparse' filesep 'cs_mex' obj] ;
end
for i = 1:length (cs)
    [s t kk] = compile_source (srcdir, cs {i}, obj, hfile, force, kk, details) ;
    timestamp = max (timestamp, t) ;
    anysrc = anysrc | s ;                                                   %#ok
    CS = [CS ' ' cs{i} obj] ;                                               %#ok
    if (nargout > 0)
        objfiles = [objfiles ' ..' filesep 'CSparse' filesep cs{i} obj] ;   %#ok
    end
end

% compile each CSparse mexFunction
obj = ['.' mexext] ;
for i = 1:length (csm)
    [s t] = cs_must_compile ('', csm{i}, '_mex', obj, hfile, force) ;
    timestamp = max (timestamp, t) ;
    if (anysrc | s)                                                         %#ok
        cmd = sprintf ('mex -O -I../../Include %s_mex.c %s -output %s', ...
            csm{i}, CS, csm{i}) ;
        kk = do_cmd (cmd, kk, details) ;
    end
end

fprintf ('\n') ;
if (nargout > 1)
    timestamp_out = timestamp ;
end

if (force)
    fprintf ('CSparse successfully compiled.\n') ;
end

%-------------------------------------------------------------------------------
function [s,t,kk] = compile_source (srcdir, f, obj, hfile, force, kk, details)
% compile a source code file in ../../Source, leaving object file in
% this directory.
[s t] = cs_must_compile (srcdir, f, '', obj, hfile, force) ;
if (s)
    cmd = sprintf ('mex -O -c -I../../Include %s%s.c', srcdir, f) ;
    kk = do_cmd (cmd, kk, details) ;
end

%-------------------------------------------------------------------------------
function kk = do_cmd (s, kk, details)
%DO_CMD: evaluate a command, and either print it or print a "."
s = strrep (s, '/', filesep) ;
if (details)
    fprintf ('%s\n', s) ;
else
    if (mod (kk, 60) == 0)
        fprintf ('\n') ;
    end
    kk = kk + 1 ;
    fprintf ('.') ;
end
eval (s) ;


text:
function definput=arg_plotfwt(definput)

  definput.flags.frqbands = {'uni', 'dyad'};
  definput.flags.wavplottype = {'image', 'stem', 'surf','waterfall'};
  definput.keyvals.fc=[];
  definput.keyvals.ntickpos=10;
  definput.keyvals.tick=[];

  definput.groups.audtick={'tick',[0,100,250,500,1000,2000,4000,8000,16000,32000]};

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [ cf_out ] = MotorImagery_pseudoOnline( eeg, eegFb, online )
%MOTORIMAGERY_PSEUDOONLINE Summary of this function goes here
%   Detailed explanation goes here
if ~isfield(online,'train'); error('training module is not exist'); end
if ~isfield(online,'apply'); error('applying module is not exist'); end
if ~isfield(online,'option'); disp('applying module is not exist');
    opt={
        'windowSize', '1000' % 1s
        'paradigm', 'MotorImagery'
        'Feedback','off'
        }; % default setting
    opt=opt_cellToStruct(online.option{:});
else
    opt=opt_cellToStruct(online.option{:});
end
if ~isfield(opt,'ite')  %iteration in while loop
    opt.ite='1';
end

in=eeg; %training
if isfield(online,'train')
    for k=1:length(online.train)
        myFun=str2func(online.train{k});
        switch online.train{k}
            case 'func_csp'
                [out CSP_W, CSP_D]= feval(myFun, in, online.train{k,2},0);
            case 'func_featureExtraction'
                out=feval(myFun, in, online.train{k,2},0);
            case 'classifier_trainClassifier'
                CF_PARAM=feval(myFun, in, online.train{k,2},0);
            otherwise
                out=feval(myFun, in, online.train{k,2},0);
        end
        in=out;
    end
end

time=func_getTimeMarker(eegFb, opt);
cf_out=zeros(1,time.ite);
% eegFb.x=eegFb.cnt;  eegFb=rmfield(eegFb, 'cnt'); %% change dat.cnt to dat.x
run=true;
while run,
%     str2num(opt.ite)
    [fbData]=func_getData(eegFb, time, opt); %opt for iteration++
    in=fbData;
    if isfield(online,'apply')
        for k=1:length(online.apply)
            myFun=str2func(online.apply{k});
            switch online.apply{k}
                case 'func_projection'
                    [out]=feval(myFun, in, CSP_W);
                case 'func_featureExtraction'
                    [out]=feval(myFun, in, online.apply{k,2},0);
                case 'classifier_applyClassifier'
                    [out]=feval(myFun, in, CF_PARAM,0);
                otherwise
                    out=feval(myFun, in, online.apply{k,2},0);
            end
            in=out;
        end
        cf_out(str2num(opt.ite))=in;
        if strcmp(opt.Feedback,'on')
        end
    else
        waring('online.apply is not exist');
    end
    if str2num(opt.ite) == time.ite
        run=false;
    end
    opt.ite=func_countInc(opt.ite);
end
end

author: PatternRecognition

repo: OpenBMI


text:
function xy_io_test05 ( )

%*****************************************************************************80
%
%% XY_IO_TEST05 tests XYF_EXAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2009
%
%  Author:
%
%    John Burkardt
%
  xy_filename = 'annulus.xy';
  xyf_filename = 'annulus.xyf';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'XY_IO_TEST05\n' );
  fprintf ( 1, '  XYF_EXAMPLE sets up XY and XYF data.\n' );

  [ point_num, face_num, face_data_num ] = xyf_example_size ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Example has:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points     = %d\n', point_num );
  fprintf ( 1, '  Number of faces      = %d\n', face_num );
  fprintf ( 1, '  Number of face items = %d\n', face_data_num );

  [ xy, face_pointer, face_data ] = xyf_example ( point_num, face_num, ...
    face_data_num );

  xy_write ( xy_filename, point_num, xy );

  xyf_write ( xyf_filename, point_num, face_num, face_data_num, ...
    face_pointer, face_data );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the XY file "%s".\n', xy_filename );
  fprintf ( 1, '  and the XYF file "%s"\n', xyf_filename ) ;

  return
end


text:
classdef prtPreProcLda < prtPreProcClass
    % prtPreProcLda  Linear discriminant analysis processing
    %
    %   preProc = prtPreProcLda creates a linear discriminant pre
    %   processing object. A prtPreProcLda object projects the input data
    %   onto a linear space that best separates class labels
    %
    %   A prtPreProcLda object has the following properties:
    %
    %   nComponents - The number of dimensions to project the data onto.
    %                 This must less than or equal to the input data's
    %                 number of features, and less than or equal to the
    %                 input data sets number of classes.
    %
    %   A prtPreProcLda object also inherits all properties and functions from
    %   the prtAction class
    %
    %   More information about LDA can be found at the following URL:
    %   http://en.wikipedia.org/wiki/Linear_discriminant_analysis
    %
    %   Example:
    %
    %   dataSet = prtDataGenIris;               % Load a dataset
    %   dataSet = dataSet.retainFeatures(1:3);  % Retain the first 3 features
    %   lda = prtPreProcLda;                    % Create the pre-processor
    %
    %   lda = lda.train(dataSet);               % Train
    %   dataSetNew = lda.run(dataSet);          % Run
    %
    %   % Plot the results
    %   subplot(2,1,1); plot(dataSet);
    %   title('Original Data');
    %   subplot(2,1,2); plot(dataSetNew);
    %   title('LDA Projected Data');
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    properties (SetAccess=private)
        name = 'Linear discriminant analysis' % Linear discriminant analysis
        nameAbbreviation = 'LDA' % LDA
    end

    properties
        nComponents = 2;   % The number of LDA components
    end
    properties (SetAccess=private)
        projectionMatrix = []; % The projection matrix
        globalMean = [];       % The global mean
    end

    methods

        % Allow for string, value pairs
        function Obj = prtPreProcLda(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
	end

	methods (Hidden = true)
        function featureNameModificationFunction = getFeatureNameModificationFunction(obj) %#ok<MANU>
            featureNameModificationFunction = prtUtilFeatureNameModificationFunctionHandleCreator('LDA Score #index#');
        end
	end

    methods
        function Obj = set.nComponents(Obj,nComp)
            if ~prtUtilIsPositiveScalarInteger(nComp)
                error('prt:prtPreProcPca','nComponents must be a positive scalar integer');
            end
            Obj.nComponents = nComp;
        end
    end

    methods (Access=protected,Hidden=true)

        function Obj = trainAction(Obj,DataSet)
            if Obj.nComponents > DataSet.nClasses
                error('prt:prtPreProcLda','Attempt to train LDA pre-processor with more components (%d) than unique classes in data set (%d)',Obj.nComponents,DataSet.nClasses);
            end
            [Obj.projectionMatrix,Obj.globalMean] = prtUtilLinearDiscriminantAnalysis(DataSet,Obj.nComponents);
        end

        function DataSet = runAction(Obj,DataSet)

            X = DataSet.getObservations;
            X = bsxfun(@minus,X,Obj.globalMean);
            DataSet = DataSet.setObservations(X*Obj.projectionMatrix);
        end

    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [price, se, low, high] ...
    = LongstaffSchwartz_M2(S, M, Ns, g, df, f, B, Nr, NSim, level)
% Longstaff Schwartz method for pricing american options

v = g(:,end);   % start for backward induction

% backward induction and regression from t_{Nr-1} up to t_1
for i = Nr-1:-1:1
        index = find(g(:,i+1) > 0); % all ITM paths
        s = S(index,i+1);           % values of S at given time points
        m = M(index,i+1);           % values of M at given time points
        v = v * df(i+1);            % option value at t_i

        Acell = B(s,m);             % evaluate basis function in cell array B
        A = cell2mat(Acell{:,:});   % convert to matrix

        c = A*f(:,i);                   % continuation value
        exercise = g(index,i+1) >= c;    % early exercise
        v(index(exercise)) = g(index(exercise),i+1);
end

price = mean(v * df(1));    % final option value

% standard error and confidence interval
sv = sqrt(1/(NSim-1)*sum((v - price * ones(NSim,1)).^2));
se = sv/sqrt(NSim);
low =  price - norminv(level) * sv/sqrt(NSim);
high = price + norminv(level) * sv/sqrt(NSim);

end

author: Sable

repo: mcbench-benchmarks


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

h = sigmoid(X*theta);
ids = find(h>=0.5);
p(ids) = 1;
ids = find(h<0.5);
p(ids) = 0;

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

end

repo: machine-learning-coursera


text:
function gt = gtom(adj,numSteps)
%GTOM       Generalized topological overlap measure
%
%   gt = gtom(adj,numSteps);
%
%   The m-th step generalized topological overlap measure (GTOM) quantifies
%   the extent to which a pair of nodes have similar m-th step neighbors.
%   Mth-step neighbors are nodes that are reachable by a path of at most
%   length m.
%
%   This function computes the the M x M generalized topological overlap
%   measure (GTOM) matrix for number of steps, numSteps.
%
%   Inputs:       adj,    adjacency matrix (binary,undirected)
%            numSteps,    number of steps
%
%   Outputs:       gt,    GTOM matrix
%
%   NOTE: When numSteps is equal to 1, GTOM is identical to the topological
%   overlap measure (TOM) from reference [2]. In that case the 'gt' matrix
%   records, for each pair of nodes, the fraction of neighbors the two
%   nodes share in common, where "neighbors" are one step removed. As
%   'numSteps' is increased, neighbors that are furter out are considered.
%   Elements of 'gt' are bounded between 0 and 1.  The 'gt' matrix can be
%   converted from a similarity to a distance matrix by taking 1-gt.
%
%   References: [1] Yip & Horvath (2007) BMC Bioinformatics 2007, 8:22
%               [2] Ravasz et al (2002) Science 297 (5586), 1551.
%
%   J Goni, University of Navarra and Indiana University, 2009/2011

%#ok<*ASGLU>

%initial state for bm matrix;
bm = adj;
bmAux = bm;
numNodes = size(adj,1);

if (numSteps > numNodes)
    disp('warning, reached maximum value for numSteps. numSteps reduced to adj-size')
    numSteps = numNodes;
end

if (numSteps == 0)
    %GTOM0
    gt = adj;
else

    for steps = 2:numSteps
        for i = 1:numNodes

            %neighbours of node i
            [neighRow,neighColumn] = find(bm(i,:)==1);

            %neighbours of neighbours of node i
            [neighNeighRow,neighNeighColumn] = find(bm(neighColumn,:)==1);
            newNeigh = setdiff(unique(neighNeighColumn),i);

            %neighbours of neighbours of node i become considered node i neighbours
            bmAux(i,newNeigh) = 1;

            %keep symmetry of matrix
            bmAux(newNeigh,i) = 1;
        end
        %bm is updated with new step all at once
        bm = bmAux;

    end

    clear bmAux newNeigh;

    %numerators of GTOM formula
    numeratorMatrix = bm*bm + adj + speye(numNodes,numNodes);

    %vector containing degree of each node
    bmSum=sum(bm);
    clear bm;

    denominatorMatrix = -adj + min(repmat(bmSum,numNodes,1),repmat(bmSum',1,numNodes)) + 1;
    gt = numeratorMatrix ./ denominatorMatrix;
end