text:
function res = spm_eeg_artefact_nans(S)
% Plugin for spm_eeg_artefact doing NaN detection
% S            - input structure
% fields of S:
%    S.D       - M/EEG object
%    S.chanind - vector of indices of channels that this plugin will look at
%
%    Additional parameters can be defined specific for each plugin.
%
% Output:
% res -
%    If no input is provided the plugin returns a cfg branch for itself.
%
%    If input is provided the plugin returns a matrix of size D.nchannels x D.ntrials
%    with zeros for clean channel/trials and ones for artefacts.
%__________________________________________________________________________
% Copyright (C) 2011-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_artefact_nans.m 7132 2017-07-10 16:22:58Z guillaume $

%-This part if for creating a config branch that plugs into spm_cfg_eeg_artefact
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0
    nans      = cfg_branch;
    nans.tag  = 'nans';
    nans.name = 'Detect NaNs';
    nans.val  = {};
    nans.help = {''};

    res = nans;

    return
end

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm('FigName','M/EEG NaN detection');

D = spm_eeg_load(S.D);

chanind = S.chanind;
res = zeros(D.nchannels, D.ntrials);

if isequal(S.mode, 'reject')
    res = zeros(D.nchannels, D.ntrials);

    %-Artefact detection
    %----------------------------------------------------------------------

    spm_progress_bar('Init', D.ntrials, 'Trials checked');
    if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials,100));
    else Ibar = [1:D.ntrials]; end

    for i = 1:D.ntrials
        for j = 1:length(chanind)
            if any(isnan(squeeze(D(chanind(j), :, i))))
                res(chanind(j), i) = 1;
            end
        end
        if any(Ibar == i), spm_progress_bar('Set', i); end
    end

    spm_progress_bar('Clear');

elseif isequal(S.mode, 'mark')
    if isequal(D.type, 'continuous')
        spm_progress_bar('Init', length(chanind), 'Channels checked');
        if length(chanind) > 100, Ibar = floor(linspace(1, length(chanind),100));
        else Ibar = [1:length(chanind)]; end
    else
        spm_progress_bar('Init', D.ntrials, 'Trials checked');
        if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials,100));
        else Ibar = [1:D.ntrials]; end
    end

    for i = 1:D.ntrials
        res = [];
        for j = 1:length(chanind)
            dat  = ~isnan(squeeze(D(chanind(j), :, i)));
            if  sum(dat)/length(dat)<(1-S.badchanthresh)
                res(end+1).type   = 'artefact_nan';
                res(end).value    = char(D.chanlabels(chanind(j)));
                res(end).time     = D.trialonset(i);
                res(end).duration = D.time(end) - D.time(1) + 1;
            else
                tmp  = find(dat);
                diffs = diff([0 tmp D.nsamples]);
                onsets = find(diffs>1);

                onsetsamples = [];
                if any(onsets == 1);
                    onsetsamples = 1;
                    onsets(1)    = [];
                    onsetsamples = [onsetsamples tmp(onsets-1)+1];
                    onsets       = [1 onsets];
                else
                    onsetsamples = [onsetsamples tmp(onsets-1)+1];
                end

                k = 1;
                m = 1;
                while k<=length(onsets)
                    if m <= length(onsets)
                        ind1 = onsetsamples(k);
                        ind2 = onsetsamples(m) + diffs(onsets(m))-2;
                        if ind2 > length(dat)
                            ind2 = length(dat);
                        end
                        if (sum(dat(ind1:ind2))/(ind2-ind1+1))<0.5
                            m = m+1;
                        else
                            if m>k
                                m = m-1;
                            end

                            res(end+1).type   = 'artefact_nan';
                            res(end).value    = char(D.chanlabels(chanind(j)));
                            res(end).time     = D.time(onsetsamples(k)+1) - D.time(1) + D.trialonset(i);
                            res(end).duration = (onsetsamples(m) + diffs(onsets(m))-onsetsamples(k)-1)/D.fsample;

                            k = m+1;
                            m = k;
                        end
                    else
                        ind1 = onsetsamples(k);
                        ind2 = length(dat);
                        if (sum(dat(ind1:ind2))/(ind2-ind1+1))<0.5
                            res(end+1).type   = 'artefact_nan';
                            res(end).value    = char(D.chanlabels(chanind(j)));
                            res(end).time     = D.time(onsetsamples(k)+1) - D.time(1) + D.trialonset(i);
                            res(end).duration = (length(dat)-onsetsamples(k)+2)/D.fsample;
                        else
                            m = m-1;

                            res(end+1).type   = 'artefact_nan';
                            res(end).value    = char(D.chanlabels(chanind(j)));
                            res(end).time     = D.time(onsetsamples(k)+1) - D.time(1) + D.trialonset(i);
                            res(end).duration = (onsetsamples(m) + diffs(onsets(m))-onsetsamples(k)-1)/D.fsample;
                        end
                        break;
                    end
                end
            end

            if isequal(D.type, 'continuous')
                if any(Ibar == j), spm_progress_bar('Set', j); end
            end
        end
        if ~isempty(res)
            ev = D.events(i);
            if iscell(ev)
                ev = ev{1};
            end

            if ~S.append
                ev(strmatch('artefact_nan', {ev.type})) = [];
            end

            D = events(D, i, spm_cat_struct(ev, res));
        end

        if ~isequal(D.type, 'continuous')
            if any(Ibar == i), spm_progress_bar('Set', i); end
        end
    end

    spm_progress_bar('Clear');

    res = D;
end

spm('FigName','M/EEG NaN detection: done');


text:
function val = csa_header(s, key)
fld = 'CSAImageHeaderInfo';
if isfield(s, fld) && isfield(s.(fld), key), val = s.(fld).(key); return; end
if isfield(s, key), val = s.(key); return; end % general tag: 2nd choice
try val = s.PerFrameFunctionalGroupsSequence.Item_1.(fld).Item_1.(key); return; end
fld = 'CSASeriesHeaderInfo';
if isfield(s, fld) && isfield(s.(fld), key), val = s.(fld).(key); return; end
val = [];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function connect = mocapConnections(fileName, pointNames);

% MOCAPCONNECTIONS Give a connection matrix for the motion capture data.
% FORMAT
% DESC load a connection matrix for txt file based mocap data.
% ARG fileName : the file from which to load the connectivity data.
% ARG pointNames : the names of the points in the file, with the
% ordering matching orderings in any data files to be loaded in (as
% returned by MOCAPPARSETEXT).
%
% SEEALSO : mocapLoadTextData, mocapParseText
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% DATASETS

fid = fopen(fileName);
i = 1;
rem = fgets(fid);
while(rem ~= -1)
  [token, rem] = strtok(rem, ',');
  connections{i, 1} = fliplr(deblank(fliplr(deblank(token))));
  [token, rem] = strtok(rem, ',');
  connections{i, 2} = fliplr(deblank(fliplr(deblank(token))));
  i = i + 1;
  rem = fgets(fid);
end

connect = zeros(length(pointNames));
fclose(fid);
for i = 1:size(connections, 1);
  for j = 1:length(pointNames)
    if strcmp(pointNames{j}, connections{i, 1}) | ...
          strcmp(pointNames{j}, connections{i, 2})
      for k = 1:length(pointNames)
        if k == j
          break
        end
        if strcmp(pointNames{k}, connections{i, 1}) | ...
              strcmp(pointNames{k}, connections{i, 2})
          connect(j, k) = 1;
        end
      end
    end
  end
end
connect = sparse(connect);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function problem_num = p00_problem_num ( )

%*****************************************************************************80
%
%% P00_PROBLEM_NUM returns the number of problems available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%   Output, integer PROBLEM_NUM, the number of problems available.
%
  problem_num = 12;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%RCSSM		read CSSM statistics from google group
%		RCSSM extracts CSSM's posting activity
%		since its inauguration on jan 8, 1993
%		from google's group site and
%		- returns the data in a table
%		- plots an activity surface [abs/gradient]
%
%		CSSM = MATLAB's newsgroup: comp.soft-sys.matlab
%
%SYNTAX
%		T = RCSSM;
%
%OUTPUT
% T		output table with fields
%		.y	years: 1993 - current
%		.m	month: 1:12
%		.d	activity table
%			rows = year	[1=1993]
%			cols = month	[1= jan]
%
%EXAMPLE
%		tbl=rcssm;

% created:
%	us	10-Jan-2007 us@neurol.unizh.ch
% modified:
%	us	20-Feb-2008 13:14:18

function	t=rcssm

% currently valid google web address and rex templates
% - 20feb08	us

		url='http://groups.google.com/group/comp.soft-sys.matlab/about';
		rex='(?<=month/)\d+-\d+">\d+';
		fmt='%d-%d">%d';

		t=get_html(url,rex);
		t=get_tbl(t,fmt);
	if	~isempty(t)
		t=plt_tbl(t);
	end
end
%-------------------------------------------------------------------------------
function	t=get_html(url,rex)

		import com.mathworks.mde.desk.*;
		wb=com.mathworks.mde.webbrowser.WebBrowser.createBrowser;
		wb.setCurrentLocation(url(8:end));
		pause(1);

		t={};
	while	isempty(t)
		s=char(wb.getHtmlText);
		t=regexp(s,rex,'match').';
% make sure we can bail-out with CTR-C
		pause(.1);
	end
		desk=MLDesktop.getInstance;
		desk.removeClient(wb);
end
%-------------------------------------------------------------------------------
function	tbl=get_tbl(t,fmt)

		[t,n]=cellfun(@(x) sscanf(char(x),fmt),t,'uni',false);
	if	~all([n{:}]==3)
		disp(sprintf('RCSSM> error reading table'));
		tbl=[];
		return;
	end
		t=cat(2,t{:});
		iy=t(1,:)-t(1,1)+1;
		im=t(2,:);

		tbl.y=unique(t(1,:));
		tbl.iy=unique(iy);
		tbl.m=1:12;
		tbl.d=accumarray([iy.',im.'],t(3,:),[],[],nan);
end
%-------------------------------------------------------------------------------
function	tbl=plt_tbl(tbl)

		m={
			'jan',	'feb',	'mar',...
			'apr',	'may',	'jun',...
			'jul',	'aug',	'sep',...
			'oct',	'nov',	'dec'
		};

% surfaces
		surf(tbl.d.');
		colormap(jet(256));
		shading interp;
		hold on;
		surf(gradient(tbl.d).');
		shading interp;

		mh=mesh(tbl.d.'+10);
		set(mh,...
			'facecolor','none',...
			'edgecolor',[0,0,0]);

		brighten(.6);
		view(22,25);

% labels
		y=tbl.y;
		iy=1:2:max(tbl.iy);
	if	iy(end) ~= max(tbl.iy)
		iy=[iy,max(tbl.iy)];
		y=[y(1:2:end),y(end)];
	else
		y=y(1:2:end);
	end
		y=strread(sprintf('%2.2d\n',rem(y-1900,100)),'%s');

		set(gca,...
			'xlim',[1,iy(end)],...
			'xtick',iy,...
			'xticklabel',y,...
			'ylim',[1,12],...
			'ytick',1:12,...
			'yticklabel',m,...
			'fontsize',10,...
			'fontname','arial');

		lh(1)=xlabel('year');
		lh(2)=ylabel('month');
		lh(3)=zlabel('entries   [absolute / gradient]');
		set(lh,...
			'fontweight','bold',...
			'fontsize',12,...
			'fontname','arial');

		title('CSSM activity',...
			'fontweight','bold',...
			'fontsize',16,...
			'fontname','arial');

		shg;
		rotate3d on;
end

author: Sable

repo: mcbench-benchmarks


text:
function [varargout] = kruskal_mst(A,varargin)
% KRUSKAL_MST Compute a minimum spanning with Kruskal's algorithm.
%
% The Kruskal MST algorithm computes a minimum spanning tree for a graph.
%
% This method works on weighted symmetric graphs.
% The runtime is O(E log (E)).
%
% See the mst function for calling information.  This function just calls
% mst(...,struct('algname','kruskal'));
%
% Example:
%    load graphs/clr-24-1.mat
%    kruskal_mst(A)
%
% See also MST, PRIM_MST.

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2006-04-23: Initial version
%  2008-09-24: Code cleanup
%%

algname = 'kruskal';
if ~isempty(varargin),
    options = merge_options(struct(),varargin{:});
    options.algname= algname;
else options = struct('algname',algname);
end

varargout = cell(1,max(nargout,1));

[varargout{:}] = mst(A,options);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function ma_UndoAction(handles)
%ma_UndoAction Summary of this function goes here
%   Detailed explanation goes here
    global number_state_log_entries_per_coast num_SoI_search_revs strict_SoI_search use_selective_soi_search soi_search_tol num_soi_search_attempts_per_rev;

    undo_states = getappdata(handles.ma_MainGUI,'undo_states');
    undo_pointer = getappdata(handles.ma_MainGUI,'undo_pointer');
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');

    if(undo_pointer < 1 || undo_pointer > size(undo_states,1))
        return;
    end

    if(undo_pointer == size(undo_states,1))
        undo_states{end+1,1} = maData;
        undo_states{end,2}   = '';
    end

    maData = undo_states{undo_pointer,1};
    undo_pointer = undo_pointer - 1;

    number_state_log_entries_per_coast = maData.settings.numStateLogPtsPerCoast;
    num_SoI_search_revs = maData.settings.numSoISearchRevs;
    strict_SoI_search = maData.settings.strictSoISearch;
    use_selective_soi_search = maData.settings.useSelectiveSoISearch;
    soi_search_tol = maData.settings.soiSearchTol;
    num_soi_search_attempts_per_rev = maData.settings.numSoiSearchAttemptsPerRev;

    maData.stateLog = ma_executeScript(maData.script,handles,celBodyData,findobj('Tag','scriptWorkingLbl'));
    maData.celBodyData = celBodyData;
    setappdata(handles.ma_MainGUI,'undo_pointer',undo_pointer);
    setappdata(handles.ma_MainGUI,'undo_states',undo_states);
    setappdata(handles.ma_MainGUI,'ma_data',maData);
    ma_processData(handles);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function organizeFusedTextures_STS(pathWORK,nPatient,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% function organizeFusedTextures(pathWORK,nPatient,MRIinv_cell,MRIweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function organizes FUSED texture features computed for all
% patients, for all different combinations of the following texture
% extraction parameters:
% - MRI Inv.: Inversion of MRI intensities in the PET/MRI fusion process.
% - MRI weight: Weight given to MRI wavelet sub-bands in the PET/MRI fusion
%               process.
% - WPBF ratio 'R': Ratio specifying the ratio of weight given to band-pass
%                   coefficients to the weigth given to the rest of
%                   coefficients (HHH and LLL) in the wavelet band-pass
%                   filtering process.
% - Scale: Resolution at which 'volume' is isotropically resampled.
% - Quant. Algo.: Quantization algorithm used on 'volume'.
% - Ng: Number of gray-levels on the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to read all possible combinations. This function is
% used for FUSED scans specifically. See Ref. [1] and 'prepareVolume.m'
% for more details.
%
% Textures are organized in three cells: 'textures_PET_T1' and
% 'textures_PET_T2FS', where each cell entry organizes the data for all
% patients for a specific extraction parameter combination. T2-weighted
% fat-saturated and STIR scans are combined in the same texture category
% and simply refered to as 'T2FS' scans (T2-weighted fat-suppression scans).
%
% -- > For example, size(textures_PET_T1) is equal to:
% [numel(MRIinv_cell),numel(MRIweight_mat),numel(R_mat),numel(scale_cell),numel(algo_cell),numel(Ng_mat)]
%
% In addition, the Spearman's rank correlation between each texture
% features and lung metastases is computed and saved in the cells for each
% texture feature of each entry. Results are saved in the STS WORKSPACE
% with the names given above.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the STS WORKSPACE directory.
% - nPatient: Number of patients to read.
% - MRIinv_cell: Cell vector specifying the different 'MRI Inv.' values to test.
% - MRIweight_mat: Array vector specifying the different 'MRI weight' values to test.
% - R_mat: Array vector specifying the different 'R' ratio values to test.
% - scale_cell: Cell vector specifying the different 'Scale' values to test.
% - algo_cell: Cell vector specifying the different 'Quant. Algo.' values to test.
% - Ng_mat: Array vector specifying the different 'Ng' values to test.
% -------------------------------------------------------------------------
% EXAMPLE:
% MRIinv_cell = {'NoInv','Inv'};
% MRIweight_mat = [1/4,1/3,1/2,2/3,3/4];
% R_mat = [1/2,2/3,1,3/2,2];
% scale_cell = {'pixelW',1,2,3,4,5};
% algo_cell = {'Equal','Lloyd'};
% Ng_mat = [8,16,32,64];
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathWORK), load('outcome')

% INITIALIZATION
scans = {'PET_T1','PET_T2FS'};
nScans = numel(scans);
suffix = ['_TEXT'];

fprintf('ORGANIZING TEXTURES FROM FUSED SCANS ... ')
for scan = 1:numel(scans)
    cd([pathWORK,'/TEXTURES'])
    temp = load('Patient1_PET_T1_TEXT'); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    nameType = fieldnames(temp.Experiment1); nType = numel(nameType); % All texture types are the same
    text = cell(numel(MRIinv_cell),numel(MRIweight_mat),numel(R_mat),numel(scale_cell),numel(algo_cell),numel(Ng_mat)); % Initialization of master cell
    tempText = cell(1,nPatient); % Cell used to load patient textures only once
    for p = 1:nPatient
        load(['Patient',num2str(p),'_',scans{scan},suffix]) % Variable 'textures' is now in MATLAB workspace
        tempText{p} = textures;
    end
    experiment = 0;
    for i = 1:numel(MRIinv_cell)
        for w = 1:numel(MRIweight_mat)
            for r = 1:numel(R_mat)
                for s = 1:numel(scale_cell)
                    for a = 1:numel(algo_cell)
                        for n = 1:numel(Ng_mat)
                            text{i,w,r,s,a,n} = struct;
                            experiment = experiment + 1;
                            strExperiment = ['Experiment',num2str(experiment)];
                            for t = 1:nType-1 % Last field is 'parameters'
                                nameFeature = fieldnames(temp.(strExperiment).(nameType{t})); nFeature = numel(nameFeature);
                                for f = 1:nFeature
                                    data = zeros(nPatient,1);
                                    for p = 1:nPatient
                                        data(p,1) = tempText{p}.(strExperiment).(nameType{t}).(nameFeature{f});
                                    end
                                    text{i,w,r,s,a,n}.(nameType{t}).(nameFeature{f}).Data = data;
                                    text{i,w,r,s,a,n}.(nameType{t}).(nameFeature{f}).Spearman = corr(data,outcome,'type','Spearman');
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    cd(pathWORK), save(['textures_',scans{scan}],'text')
end
fprintf('DONE\n')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function rval=rtest_1()

n = 49;
[A,b] = testmat(n,2);
x0 = [1:n]'/(n+1);
y0 = [1:n]'/(n+1);
x = repmat(x0,1,n);
y = repmat(y0',n,1);
xy = [x(:),y(:)];

rval = 0;
try
    T = mst(A);
    T = T + diag(diag(A));
    rval = 1;
catch
    lasterr
end;

try
    A(1,2)= -1;
    A(2,1)= -1;
    T = prim_mst(A);
    rval = 0;
catch
    rval = 1;
end;

function [A,b] = testmat( n,stencil )
    h = 1/(n+1);

    % initialization
    x = [1:n]'/(n+1);
    y = [1:n]'/(n+1);
    u = zeros(n,n);

    % exact solution
    u0 = repmat(x.^4,1,n) + repmat(12*y'.^2,n,1);

    % matices
    T0 = -ones(n);
    T0 = sparse( triu(tril(T0,-1),-1) + triu(tril(T0,1),1) );
    I  = speye(n);
    if stencil == 1
        T = (-8*I - T0) / 3;
        B = (I - T0) / 3;
    else
        T = (-20*I - 4*T0) / 6;
        B = (4*I - T0) / 6;
    end

    A = kron(I,T) - kron(T0,B);

    % boundary
    b = zeros(n);
    if stencil == 1
        b(1,:) = b(1,:) + 12*y'.^2 + 12*(y-h)'.^2 + 12*(y+h)'.^2;
        b(n,:) = b(n,:) + 3 + 12*y'.^2 + 12*(y-h)'.^2 + 12*(y+h)'.^2;
        b(:,1) = b(:,1) + x.^2 + (x-h).^2 + (x+h).^2;
        b(:,n) = b(:,n) + x.^2 + (x-h).^2 + (x+h).^2 + 12*3;
    else
        b(1,:) = b(1,:) + 4* 12*y'.^2 + 12*(y-h)'.^2 + 12*(y+h)'.^2;
        b(n,:) = b(n,:) + 6 + 4* 12*y'.^2 + 12*(y-h)'.^2 + 12*(y+h)'.^2;
        b(:,1) = b(:,1) + 4* x.^2 + (x-h).^2 + (x+h).^2;
        b(:,n) = b(:,n) + 4* x.^2 + (x-h).^2 + (x+h).^2 + 12*6;
    end
    % fix corners
    b(1,n) = b(1,n) - 12;
    b(n,1) = b(n,1) - 1;
    b(n,n) = b(n,n) - 13;
    % normalize
    if stencil == 1
        b = - b / 3;
    else
        b = - b / 6;
    end
    % add source
    f = 12*x.^2 + 24;
    f = repmat(f,1,n);
    b = b + f * h^2;

    b = b(:);
return


text:
function ydata = tsne_d(D, labels, no_dims, perplexity)
%TSNE_D Performs symmetric t-SNE on the pairwise Euclidean distance matrix D
%
%   mappedX = tsne_d(D, labels, no_dims, perplexity)
%   mappedX = tsne_d(D, labels, initial_solution, perplexity)
%
% The function performs symmetric t-SNE on the NxN pairwise Euclidean
% distance matrix D to construct an embedding with no_dims dimensions
% (default = 2). An initial solution obtained from an other dimensionality
% reduction technique may be specified in initial_solution.
% The perplexity of the Gaussian kernel that is employed can be specified
% through perplexity (default = 30). The labels of the data are not used
% by t-SNE itself, however, they are used to color intermediate plots.
% Please provide an empty labels matrix [] if you don't want to plot
% results during the optimization.
% The data embedding is returned in mappedX.
%
%
% (C) Laurens van der Maaten, 2010
% University of California, San Diego

    if ~exist('labels', 'var')
        labels = [];
    end
    if ~exist('no_dims', 'var') || isempty(no_dims)
        no_dims = 2;
    end
    if ~exist('perplexity', 'var') || isempty(perplexity)
        perplexity = 30;
    end

    % First check whether we already have an initial solution
    if numel(no_dims) > 1
        initial_solution = true;
        ydata = no_dims;
        no_dims = size(ydata, 2);
    else
        initial_solution = false;
    end

    % Compute joint probabilities
    D = D / max(D(:));                                                      % normalize distances
    P = d2p(D .^ 2, perplexity, 1e-5);                                      % compute affinities using fixed perplexity

    % Run t-SNE
    if initial_solution
        ydata = tsne_p(P, labels, ydata);
    else
        ydata = tsne_p(P, labels, no_dims);
    end

repo: mobile-id


text:
function sats=satid2no(satid)

global glc;
sats=0;

[prn,n]=sscanf(satid,'%d');
if n==1
%     prn=str2double(satid);
    if     glc.MINPRNGPS<=prn&&prn<=glc.MAXPRNGPS,sys=glc.SYS_GPS;
    elseif glc.MINPRNQZS<=prn&&prn<=glc.MAXPRNQZS,sys=glc.SYS_QZS;
    else  ,return ;
    end
    sats=satno(sys,prn);
end

[str,n]=sscanf(satid,'%c%d',2);
if n<2,return;end
code=char(str(1)); prn=str(2);
switch code
    case 'G',sys=glc.SYS_GPS;prn=prn+glc.MINPRNGPS-1;
    case 'R',sys=glc.SYS_GLO;prn=prn+glc.MINPRNGLO-1;
    case 'E',sys=glc.SYS_GAL;prn=prn+glc.MINPRNGAL-1;
    case 'C',sys=glc.SYS_BDS;prn=prn+glc.MINPRNBDS-1;
    case 'J'
        sys=glc.SYS_QZS;prn=prn+glc.MINPRNQZS-1;
    otherwise, return;
end
sats=satno(sys,prn);

return

repo: GINav


text:
% Test file for chebtech/cell2mat.m

function pass = test_cell2mat(pref)

if ( nargin < 2 )
    pref = chebtech.techPref();
end

for n = 1:2

    if (n == 1);
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    f = testclass.make(@(x) [sin(x) cos(x) exp(x)], [], pref);
    g = testclass.make(@(x) sin(x), [], pref);
    h = testclass.make(@(x) [cos(x) exp(x)], [], pref);

    F = cell2mat([g h]);
    pass(n, 1) = all( sum(F - f) < max(vscale(f)*eps) );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef MeshStructure
    %MeshStructure class
    % contains information about the domain and the mesh size

    properties
        dimension
        dims
        cellsize
        cellcenters
        facecenters
        corners
        edges
    end

    methods
        function meshVar = MeshStructure(dimension, dims, cellsize, ...
          cellcenters, facecenters, corners, edges)
            if nargin>0
                meshVar.dimension = dimension;
                meshVar.dims = dims;
                meshVar.cellsize = cellsize;
                meshVar.cellcenters = cellcenters;
                meshVar.facecenters = facecenters;
                meshVar.corners= corners;
                meshVar.edges= edges;
            end
        end
    end
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function p = remove_redundant(p);

if isempty(p.F_struc)
    return
end

b = p.F_struc(1+p.K.f:p.K.l+p.K.f,1);
A = -p.F_struc(1+p.K.f:p.K.l+p.K.f,2:end);

redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));

if length(redundant)>1
    p.InequalityConstraintState(redundant) = inf;
end

if p.options.bmibnb.lpreduce
    b = p.lpcuts(:,1);
    A = -p.lpcuts(:,2:end);
    redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));
    if length(redundant)>1
        p.lpcuts(redundant,:) = [];
        p.cutState(redundant) = [];
    end
end

if p.K.f > 0
    b = p.F_struc(1:p.K.f,1);
    A = -p.F_struc(1:p.K.f,2:end);
    s1 = ((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <1e-6);
    s2 = ((-A>0).*(-A)*(p.ub-p.lb) - ((-b)-(-A)*p.lb) <1e-6);
    redundant = find(s1 & s2);
    if length(redundant)>1
        p.EqualityConstraintState(redundant) = inf;
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

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

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = sum(X)./size(X,1);
murep = repmat(mu, [size(X,1) 1]);
sigma2 = sum((X-murep).^2)/size(X,1);

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

end

repo: machine-learning-coursera


text:
function vu = burgers_solution ( nu, vxn, vx, vtn, vt )

%*****************************************************************************80
%
%% BURGERS_SOLUTION evaluates a solution to the Burgers equation.
%
%  Discussion:
%
%    The form of the Burgers equation considered here is
%
%      du       du        d^2 u
%      -- + u * -- = nu * -----
%      dt       dx        dx^2
%
%    for -1.0 < x < +1.0, and 0 < t.
%
%    Initial conditions are u(x,0) = - sin(pi*x).  Boundary conditions
%    are u(-1,t) = u(+1,t) = 0.  The viscosity parameter nu is taken
%    to be 0.01 / pi, although this is not essential.
%
%    The authors note an integral representation for the solution u(x,t),
%    and present a better version of the formula that is amenable to
%    approximation using Hermite quadrature.
%
%    This program library does little more than evaluate the exact solution
%    at a user-specified set of points, using the quadrature rule.
%    Internally, the order of this quadrature rule is set to 8, but the
%    user can easily modify this value if greater accuracy is desired.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2011
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Claude Basdevant, Michel Deville, Pierre Haldenwang, J Lacroix,
%    J Ouazzani, Roger Peyret, Paolo Orlandi, Anthony Patera,
%    Spectral and finite difference solutions of the Burgers equation,
%    Computers and Fluids,
%    Volume 14, Number 1, 1986, pages 23-41.
%
%  Parameters:
%
%    Input, real NU, the viscoscity.
%
%    Input, integer VXN, the number of spatial grid points.
%
%    Input, real VX(VXN), the spatial grid points.
%
%    Input, integer VTN, the number of time grid points.
%
%    Input, real VT(VTN), the time grid points.
%
%    Output, real VU(VXN,VTN), the solution of the Burgers
%    equation at each space and time grid point.
%
  qn = 8;
%
%  Compute the rule.
%
  [ qx, qw ] = hermite_ek_compute ( qn );
%
%  Evaluate U(X,T) for later times.
%
  vu = zeros ( vxn, vtn );

  for vti = 1 : vtn

    if ( vt(vti) == 0.0 )

      vu(1:vxn,vti) = - sin ( pi * vx(1:vxn) );

    else

      for vxi = 1 : vxn

        top = 0.0;
        bot = 0.0;

        for qi = 1 : qn

          c = 2.0 * sqrt ( nu * vt(vti) );

          top = top - qw(qi) * c * sin ( pi * ( vx(vxi) - c * qx(qi) ) ) ...
            * exp ( - cos ( pi * ( vx(vxi) - c * qx(qi)  ) ) ...
            / ( 2.0 * pi * nu ) );

          bot = bot + qw(qi) * c ...
            * exp ( - cos ( pi * ( vx(vxi) - c * qx(qi)  ) ) ...
            / ( 2.0 * pi * nu ) );

          vu(vxi,vti) = top / bot;

        end

      end

    end

  end

  return
end


text:
%MIT IAP Radar Course 2011
%Resource: Build a Small Radar System Capable of Sensing Range, Doppler,
%and Synthetic Aperture Radar Imaging
%
%Gregory L. Charvat

%SAR algorithm from:
%Range Migration Algorithm from ch 10 of Spotlight Synthetic Aperture Radar
%Signal Processing Algorithms, Carrara, Goodman, and Majewski

%NOTE: set up-ramp sweep from 2-3.2V to stay within ISM band
%change fstart and fstop bellow when in ISM band

%-------------------------------------------%
%Process raw data here
clear all;
close all;

%read the raw data .wave file here
[Y,FS,NBITS] = wavread('towardswarehouse.wav');

%constants
c = 3E8; %(m/s) speed of light

%radar parameters
Tp = 20E-3; %(s) pulse time
Trp = 0.25; %(s) min range profile time duration
N = Tp*FS; %# of samples per pulse
fstart = 2260E6; %(Hz) LFM start frequency
fstop = 2590E6; %(Hz) LFM stop frequency
%fstart = 2402E6; %(Hz) LFM start frequency for ISM band
%fstop = 2495E6; %(Hz) LFM stop frequency for ISM band
BW = fstop-fstart; %(Hz) transmti bandwidth
f = linspace(fstart, fstop, N/2); %instantaneous transmit frequency

%the input appears to be inverted
trig = -1*Y(:,1);
s = -1*Y(:,2);
clear Y;

%parse data here by position (silence between recorded data)
rpstart = abs(trig)>mean(abs(trig));
count = 0;
Nrp = Trp*FS; %min # samples between range profiles

for ii = Nrp+1:size(rpstart,1)-Nrp
    if rpstart(ii) == 1 & sum(rpstart(ii-Nrp:ii-1)) == 0
        count = count + 1;
        RP(count,:) = s(ii:ii+Nrp-1);
        RPtrig(count,:) = trig(ii:ii+Nrp-1);
    end
end

%parse data by pulse
count = 0;
thresh = 0.08;
clear ii;
for jj = 1:size(RP,1)
    %clear SIF;
    SIF = zeros(N,1);
    start = (RPtrig(jj,:)> thresh);
    count = 0;
    jj
    for ii = 12:(size(start,2)-2*N)
        [Y I] =  max(RPtrig(jj,ii:ii+2*N));
        if mean(start(ii-10:ii-2)) == 0 & I == 1
            count = count + 1;
            SIF = RP(jj,ii:ii+N-1)' + SIF;
        end
    end
    %hilbert transform
    q = ifft(SIF/count);
    sif(jj,:) = fft(q(size(q,1)/2+1:size(q,1)));
end
sif(find(isnan(sif))) = 1E-30; %set all Nan values to 0

%SAR data should be ready here
clear s;
s = sif;
save routsidewarehouse2 s; %for image data

%-------------------------------------------%
%load additional varaibles and setup constants for radar here
clear all;
c = 3E8; %(m/s) speed of light

%load IQ converted data here
load routsidewarehouse2 s; %load variable sif %for image data

for ii = 1:size(s,1)
    s(ii,:) = s(ii,:) - mean(s,1);
end

%sif = s-sif_sub; %perform coherent background subtraction
%sif = sif_sub; %image just the background
sif = s; %image without background subtraction
clear s;
clear sif_sub;

%***********************************************************************
%radar parameters
fc = (2590E6 - 2260E6)/2 + 2260E6; %(Hz) center radar frequency
B = (2590E6 - 2260E6); %(hz) bandwidth
cr = B/20E-3; %(Hz/sec) chirp rate
Tp = 20E-3; %(sec) pulse width
%VERY IMPORTANT, change Rs to distance to cal target
%Rs = (12+9/12)*.3048; %(m) y coordinate to scene center (down range), make this value equal to distance to cal target
Rs = 0;
Xa = 0; %(m) beginning of new aperture length
delta_x = 2*(1/12)*0.3048; %(m) 2 inch antenna spacing
L = delta_x*(size(sif,1)); %(m) aperture length
Xa = linspace(-L/2, L/2, (L/delta_x)); %(m) cross range position of radar on aperture L
Za = 0;
Ya = Rs; %THIS IS VERY IMPORTANT, SEE GEOMETRY FIGURE 10.6
t = linspace(0, Tp, size(sif,2)); %(s) fast time, CHECK SAMPLE RATE
Kr = linspace(((4*pi/c)*(fc - B/2)), ((4*pi/c)*(fc + B/2)), (size(t,2)));

%Save background subtracted and callibrated data
save sif sif delta_x Rs Kr Xa;
%clear all;

%run IFP
SBAND_RMA_IFP;

author: lukeweston

repo: SimpleFMCWRadar


text:
%kEllipse 'Ellipse Object Menuform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Ellipse.pane file
%
% Parameters:
%
% Example: kEllipse( {})
%
% Khoros helpfile follows below:
% .begin tagged
%
% .item "The Ellipse Visual Object"
% An ellipse visual object supports the display of an ellipse annotation.
% The (x, y) location of the ellipse, its width and height may be specified in
% device coordinates or world coordinates by the application.  The menuform,
% however, displays the world coordinates of the ellipse, regardless of whether
% the application used world or device coordinates to specify that location.
% Attributes of the circle object that can be interactively specified in addition
% to its world coordinates include its border width, border type, border color,
% whether or not it is filled, and (if it is filled), its fill color.
%
% .item "Filled"
% When set to true, this attribute indicates that the ellipse should be filled
% with the fill color (background color).
%
% .item "Border Width"
% This list selection allows you to set the line width used by
% the ellipse object. Line widths range from "Extra Fine" to "Extra Wide".
%
% .item "Border Type"
% This list selection lets you set the line type used on
% the ellipse object. There are seven line types, including "Solid",
% "Dotted", "Dot Dash", "Short Dash", "Long Dash", "Odd Dash", and
% "Grid Dotted".
%
% .item "Border Color"
% This stringlist selection lets you set the border color (foreground color)
% of the ellipse object, using the color name.
%
% .item "Fill Color"
% This stringlist selection lets you set the fill color (background color)
% of the ellipse object, using the color name.
%
% .item "X Circle, Y Circle"
% This pair of double values specifies the (x, y) location of the ellipse object
% on its parent in world coordinates.  The world coordinate range of the ellipse
% is specified by the application, and thus varies depending on context.
%
% .item "Width"
% This double value specifies the width of the ellipse in world coordinates.
%
% .item "Height"
% This double value specifies the height of the ellipse in world coordinates.
%
% .end tagged

function varargout = kEllipse(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kEllipse(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'Ellipse" '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% prtUtilSumExp(x)
%   returns log(sum(exp(x),1)) while avoiding underflow issues
%
% Notes: This only sums down, thus the sum( * , 1)
%        This only accepts real doubles
%        If there is a large spread between the min and max values
%           some underflow is still possible. Although unlikely to
%           matter in the end due to the large spread in your values.
%
%        Error checking is minimal and seg faults are possible

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function plotIPDF(SO3F,r,varargin)
% plot inverse pole figures
%
% Input
%  odf - @SO3FunHarmonic
%  r   - @vector3d specimen directions
%
% Options
%  RESOLUTION - resolution of the plots
%
% Flags
%  antipodal - include <VectorsAxes.html antipodal symmetry>
%  complete  - plot entire (hemi)--sphere
%
% See also
% S2Grid/plot savefigure Plotting Annotations_demo ColorCoding_demo PlotTypes_demo
% SphericalProjection_demo

if numel(SO3F)>1
  warning(['You try to plot an multivariate function. Plot the desired components ' ...
    'manually. In the following the first component is plotted.'])
end

plotIPDF@SO3Fun(SO3F.subSet(1),r,varargin{:});

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
classdef StressNormPduringIterAsPostprocess < handle

    properties (Access = public)

    end

    properties (Access = private)
        mesh
        dataRes
        iterC
        cost
        fileNameIter
        finalIter
    end

    properties (Access = private)
        microCase
        path
        fileResName
        fileNameMesh
        pNorm
        valueT
        caseName
    end

    methods (Access = public)

        function obj = StressNormPduringIterAsPostprocess()
            obj.init();
            obj.obtainCost();
            value = zeros(length(obj.pNorm),obj.finalIter);
            for iter = 1:obj.finalIter
                fNameIter = [obj.fileResName,num2str(iter)];
                obj.wrapResAndMesh(fNameIter);
                value(:,iter) = obj.computeStressNorm();
                if mod(iter,1)== 0
                obj.plotValue(value,iter)
                end
            end
           obj.plotFinal(value);
        end

    end

    methods (Access = private)

        function init(obj)
            obj.valueT = 0.18324;
            %obj.valueT  = 34.6295;
            %obj.valueT = 0.5646;
            %obj.microCase = 'Rectangle';
            obj.microCase = 'SuperEllipse';
            meshType = 'Small';
            %meshType  = 'Medium';
            obj.finalIter = 73;%607;%20;%46;%607;%42;18;405;42;296;57;46;250;485;250;
%             fCase = [obj.microCase,meshType];
%             obj.caseName = ['ExperimentingPlot'];
%             obj.path = ['/media/alex/My Passport/LatticeResults/StressNorm',fCase,'/'];
%             obj.fileResName  = ['ExperimentingPlot'];

            fCase = [obj.microCase,meshType];
            obj.caseName = ['ExperimentingPlot7',fCase];
            obj.path = ['/media/alex/My Passport/LatticeResults/',obj.caseName,'/'];
            obj.fileResName  = ['ExperimentingPlot',obj.microCase];
            obj.fileNameMesh = ['CantileverSquare',meshType];
            obj.pNorm = 2:2:16;
        end

        function obtainCost(obj)
            fNameMon = fullfile(obj.path,'Monitoring.fig');
            h = openfig(fNameMon,'invisible');
            handles = findobj(h,'Type','line');
            iter = get(handles(5),'Xdata');
            c = get(handles(5),'Ydata');
            close all
            obj.iterC = iter;
            obj.cost  = c;
        end

        function plotValue(obj,value,iter)
            figure(1)
            clf
            plot(obj.iterC,obj.cost*obj.valueT)
            hold on
            plot(1:iter,value(:,1:iter))
            hold off
            drawnow
        end

        function plotFinal(obj,value)
            f = figure(1);
            clf
            p{1} = plot(obj.iterC,obj.cost*obj.valueT);
            hold on
            for i = 1:size(value,1)
              p{i+1} = plot(value(i,:));
            end
           pP = plotPrinter(f,p);
           fName = ['/home/alex/Dropbox/GregoireMeeting7Decembre/',obj.caseName];
           pP.print(fName)
        end

        function wrapResAndMesh(obj,fileNameIter)
            s.folderPath = obj.path;
            s.fileName   = fileNameIter;
            w = WrapperMshResFiles(s);
            w.compute();
            obj.mesh    = w.mesh;
            obj.dataRes = w.dataRes;
        end

         function value = computeStressNorm(obj)
            s.m1            = obj.dataRes.DesignVar1;
            s.m2            = obj.dataRes.DesignVar2;
            s.alpha         = obj.dataRes.AlphaGauss';
            switch obj.microCase
                case 'Rectangle'
                    s.vademecumName = 'SuperEllipseQMax';
                case 'SuperEllipse'
                    s.vademecumName = 'SuperEllipseQOptAnalytic';
            end
            s.mesh          = obj.mesh;
            s.fileName      = obj.fileNameMesh;
            s.pNorm         = obj.pNorm;
            sC = StressNormFromVademecumComputer(s);
            value = sC.compute();
         end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% RES = corrDn(IM, FILT, EDGES, STEP, START, STOP)
%
% Compute correlation of matrices IM with FILT, followed by
% downsampling.  These arguments should be 1D or 2D matrices, and IM
% must be larger (in both dimensions) than FILT.  The origin of filt
% is assumed to be floor(size(filt)/2)+1.
%
% EDGES is a string determining boundary handling:
%    'circular' - Circular convolution
%    'reflect1' - Reflect about the edge pixels
%    'reflect2' - Reflect, doubling the edge pixels
%    'repeat'   - Repeat the edge pixels
%    'zero'     - Assume values of zero outside image boundary
%    'extend'   - Reflect and invert (continuous values and derivs)
%    'dont-compute' - Zero output when filter overhangs input boundaries
%
% Downsampling factors are determined by STEP (optional, default=[1 1]),
% which should be a 2-vector [y,x].
%
% The window over which the convolution occurs is specfied by START
% (optional, default=[1,1], and STOP (optional, default=size(IM)).
%
% NOTE: this operation corresponds to multiplication of a signal
% vector by a matrix whose rows contain copies of the FILT shifted by
% multiples of STEP.  See upConv.m for the operation corresponding to
% the transpose of this matrix.

% Eero Simoncelli, 6/96, revised 2/97.

function res = corrDn(im, filt, edges, step, start, stop)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

% fprintf(1,'WARNING: You should compile the MEX version of "corrDn.c",\n         found in the MEX subdirectory of matlabPyrTools, and put it in your matlab path.  It is MUCH faster, and provides more boundary-handling options.\n');

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('edges') == 1)
  if (strcmp(edges,'reflect1') ~= 1)
    warning('Using REFLECT1 edge-handling (use MEX code for other options).');
  end
end

if (exist('step') ~= 1)
	step = [1,1];
end

if (exist('start') ~= 1)
	start = [1,1];
end

if (exist('stop') ~= 1)
	stop = size(im);
end

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

% Reverse order of taps in filt, to do correlation instead of convolution
filt = filt(size(filt,1):-1:1,size(filt,2):-1:1);

tmp = rconv2(im,filt);
res = tmp(start(1):step(1):stop(1),start(2):step(2):stop(2));

author: tyshiwo


text:
function f = flipud(f)
%FLIPUD   Flip/reverse a TRIGTECH object.
%   G = FLIPUD(F) returns G such that G(x) = F(-x) for all x in [-1,1].

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

% Flip the values:
f.values = [ f.values(1,:); flipud(f.values(2:end,:)) ];

% Flip the coefficients taking into account where f is odd or even
if mod(size(f.coeffs,1),2)
    % Odd length is easy, just flip the coefficients
    f.coeffs = flipud(f.coeffs);
else
    % Even length requires keeping the first coefficient in place and
    % flipping the remaining ones.  This follows since we interpret the
    % first coefficient to correspond to the 1/2*cos(-N/2 x) mode.
    f.coeffs(1) = conj(f.coeffs(1));
    f.coeffs(2:end,:) = flipud(f.coeffs(2:end,:));
end

end


text:
function [tgraph,tedge,tloc]= createtimegraph(Graph,edge,loc,t)
%% Written by Muhammet Balcilar, France
% all rights reverved

nn=size(Graph,1);
tgraph=zeros(t*nn);
tedge=[];
tloc=loc;
eyeI=eye(nn);

for i=1:t-1
    tloc=[tloc;loc];
    tgraph((i-1)*nn+1:i*nn,i*nn+1:(i+1)*nn)=Graph+eyeI;
    tmp1=edge(:,1)+(i-1)*nn;
    tmp2=edge(:,2)+(i)*nn;
    tedge=[tedge; tmp1 tmp2];
    tmp=[1:nn]';
    tmp3=tmp+(i-1)*nn;
    tmp4=tmp+(i)*nn;
    tedge=[tedge; tmp3 tmp4];
end

author: balcilar

repo: Multi-Robot-Path-Planning-on-Graphs


function dx = l2lossBackward(x,r,p )

dx = 2*p*(x-r);

dx = dx / (size(x,1)*size(x,2));

end


text:
function value = i4_characteristic ( q )

%*****************************************************************************80
%
%% I4_CHARACTERISTIC gives the characteristic for an integer.
%
%  Discussion:
%
%    For any positive integer Q, the characteristic is:
%
%    Q, if Q is a prime;
%    P, if Q = P^N for some prime P and some integer N;
%    0, otherwise, that is, if Q is negative, 0, 1, or the product
%       of more than one distinct prime.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    FORTRAN77 original version by Bratley, Fox, Niederreiter.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Harald Niederreiter,
%    Algorithm 738:
%    Programs to Generate Niederreiter's Low-Discrepancy Sequences,
%    ACM Transactions on Mathematical Software,
%    Volume 20, Number 4, pages 494-495, 1994.
%
%  Parameters:
%
%    Input, integer Q, the value to be tested.
%
%    Output, integer VALUE, the characteristic of Q.
%
  if ( q <= 1 )
    value = 0;
    return
  end
%
%  If Q is not prime, then there is at least one prime factor
%  of Q no greater than SQRT(Q)+1.
%
%  A faster code would only consider prime values of I,
%  but that entails storing a table of primes and limiting the
%  size of Q.  Simplicity and flexibility for now!
%
  i_max = floor ( sqrt ( q ) ) + 1;

  for i = 2 : i_max

    if ( mod ( q, i ) == 0 )

      while ( mod ( q, i ) == 0 )
        q = q / i;
      end

      if ( q == 1 )
        value = i;
      else
        value = 0;
      end

      return

    end

  end
%
%  If no factor was found, then Q is prime.
%
  value = q;

  return
end

repo: jburkardt-m


text:
function diff = getDiffuseness_TV(i_vecs)
%GETDIFFUSENESS_TV Compute diffuseness from temporal variation of intensity
%vectors
%
%   This routine estimates diffuseness using the temporal variation of
%   intensity vectors, as proposed in
%
%       Ahonen, J. and Pulkki, V., 2009.
%       Diffuseness estimation using temporal variation of intensity vectors.
%       In 2009 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA).
%
%   also termed coefficient of variation in other works. It is based on the
%   relation sqrt( 1 - ||<Ia>||/<||Ia||> ), where Ia is the active intensity
%   vector. It is similar to the spherical variance estimator, but it takes
%   into account the magnitude of the DoA vectors, hence giving more weight
%   to strong sound estimates instead of low noisy ones. It is also flexible
%   in the sense that alternative DoA vectors with magnitude that relates
%   to the energy of the sound-field can be used instead of the intensity
%   vector.
%
%   Inputs:
%       i_vecs:  Kx3 matrix of DoA vectors, for K temporal observations
%
%   Outputs:
%       diff:       diffuseness value
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% GETDIFFUSENESS_TV.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% mean of active intensity vectors
mean_i_vecs = mean(i_vecs, 1);
% norm of mean vector
norm_mean_i_vecs = sqrt(sum(mean_i_vecs.^2, 2));

% norm of intensity vectors
norm_i_vecs = sqrt(sum(i_vecs.^2, 2));
% mean of norms
mean_norm_i_vecs = mean(norm_i_vecs);

% diffuseness
diff = sqrt( 1 - norm_mean_i_vecs/mean_norm_i_vecs );

end

author: polarch

repo: Spherical-Array-Processing


text:
function evaluate_poses_comparison

opt = globals();
% your results
resdir = '/var/Projects/PoseRNN/output/lov/lov_keyframe/vgg16_fcn_color_single_frame_2d_lov_iter_80000/mat';

% read class names
fid = fopen('classes.txt', 'r');
C = textscan(fid, '%s');
object_names = C{1};
fclose(fid);

% load model points
num_objects = numel(object_names);
models = cell(num_objects, 1);
for i = 1:num_objects
    filename = fullfile(opt.root, 'models', object_names{i}, 'points.xyz');
    disp(filename);
    models{i} = load(filename);
end

% load the keyframe indexes
fid = fopen('keyframe.txt', 'r');
C = textscan(fid, '%s');
keyframes = C{1};
fclose(fid);

% save results
distances_sys = zeros(100000, 2);
distances_non = zeros(100000, 2);
errors_rotation = zeros(100000, 2);
errors_translation = zeros(100000, 2);
results_seq_id = zeros(100000, 1);
results_frame_id = zeros(100000, 1);
results_object_id = zeros(100000, 1);
results_cls_id = zeros(100000, 1);

% for each image
count = 0;
for i = 1:numel(keyframes)

    % parse keyframe name
    name = keyframes{i};
    pos = strfind(name, '/');
    seq_id = str2double(name(1:pos-1));
    frame_id = str2double(name(pos+1:end));

    % load PoseCNN result
    filename = sprintf('results_PoseCNN/%04d.mat', i - 1);
    result = load(filename);

    % load 3D coordinate regression result
    filename = sprintf('%s/%04d.mat', resdir, i - 1);
    result_new = load(filename);

    % load gt poses
    filename = fullfile(opt.root, 'data', sprintf('%04d/%06d-meta.mat', seq_id, frame_id));
    disp(filename);
    gt = load(filename);

    % for each gt poses
    for j = 1:numel(gt.cls_indexes)
        count = count + 1;
        cls_index = gt.cls_indexes(j);
        RT_gt = gt.poses(:, :, j);

        results_seq_id(count) = seq_id;
        results_frame_id(count) = frame_id;
        results_object_id(count) = j;
        results_cls_id(count) = cls_index;

        % network result
        roi_index = find(result.rois(:, 2) == cls_index);
        if isempty(roi_index) == 0
            RT = zeros(3, 4);

            % pose from network
            RT(1:3, 1:3) = quat2rotm(result.poses(roi_index, 1:4));
            RT(:, 4) = result.poses(roi_index, 5:7);
            distances_sys(count, 1) = adi(RT, RT_gt, models{cls_index}');
            distances_non(count, 1) = add(RT, RT_gt, models{cls_index}');
            errors_rotation(count, 1) = re(RT(1:3, 1:3), RT_gt(1:3, 1:3));
            errors_translation(count, 1) = te(RT(:, 4), RT_gt(:, 4));
        else
            distances_sys(count, 1) = inf;
            distances_non(count, 1) = inf;
            errors_rotation(count, 1) = inf;
            errors_translation(count, 1) = inf;
        end

        % your result
        roi_index = find(result_new.rois(:, 2) == cls_index);
        if isempty(roi_index) == 0
            RT = zeros(3, 4);
            RT(1:3, 1:3) = quat2rotm(result_new.poses(roi_index, 1:4));
            RT(:, 4) = result_new.poses(roi_index, 5:7);
            distances_sys(count, 2) = adi(RT, RT_gt, models{cls_index}');
            distances_non(count, 2) = add(RT, RT_gt, models{cls_index}');
            errors_rotation(count, 2) = re(RT(1:3, 1:3), RT_gt(1:3, 1:3));
            errors_translation(count, 2) = te(RT(:, 4), RT_gt(:, 4));
        else
            distances_sys(count, 2) = inf;
            distances_non(count, 2) = inf;
            errors_rotation(count, 2) = inf;
            errors_translation(count, 2) = inf;
        end
    end
end
distances_sys = distances_sys(1:count, :);
distances_non = distances_non(1:count, :);
errors_rotation = errors_rotation(1:count, :);
errors_translation = errors_translation(1:count, :);
results_seq_id = results_seq_id(1:count);
results_frame_id = results_frame_id(1:count);
results_object_id = results_object_id(1:count, :);
results_cls_id = results_cls_id(1:count, :);
save('results_comparison.mat', 'distances_sys', 'distances_non', 'errors_rotation', 'errors_translation',...
    'results_seq_id', 'results_frame_id', 'results_object_id', 'results_cls_id');

function pts_new = transform_pts_Rt(pts, RT)
%     """
%     Applies a rigid transformation to 3D points.
%
%     :param pts: nx3 ndarray with 3D points.
%     :param R: 3x3 rotation matrix.
%     :param t: 3x1 translation vector.
%     :return: nx3 ndarray with transformed 3D points.
%     """
n = size(pts, 2);
pts_new = RT * [pts; ones(1, n)];

function error = add(RT_est, RT_gt, pts)
%     """
%     Average Distance of Model Points for objects with no indistinguishable views
%     - by Hinterstoisser et al. (ACCV 2012).
%
%     :param R_est, t_est: Estimated pose (3x3 rot. matrix and 3x1 trans. vector).
%     :param R_gt, t_gt: GT pose (3x3 rot. matrix and 3x1 trans. vector).
%     :param model: Object model given by a dictionary where item 'pts'
%     is nx3 ndarray with 3D model points.
%     :return: Error of pose_est w.r.t. pose_gt.
%     """
pts_est = transform_pts_Rt(pts, RT_est);
pts_gt = transform_pts_Rt(pts, RT_gt);
diff = pts_est - pts_gt;
error = mean(sqrt(sum(diff.^2, 1)));

function error = adi(RT_est, RT_gt, pts)
%     """
%     Average Distance of Model Points for objects with indistinguishable views
%     - by Hinterstoisser et al. (ACCV 2012).
%
%     :param R_est, t_est: Estimated pose (3x3 rot. matrix and 3x1 trans. vector).
%     :param R_gt, t_gt: GT pose (3x3 rot. matrix and 3x1 trans. vector).
%     :param model: Object model given by a dictionary where item 'pts'
%     is nx3 ndarray with 3D model points.
%     :return: Error of pose_est w.r.t. pose_gt.
%     """
pts_est = transform_pts_Rt(pts, RT_est);
pts_gt = transform_pts_Rt(pts, RT_gt);

% Calculate distances to the nearest neighbors from pts_gt to pts_est
MdlKDT = KDTreeSearcher(pts_est');
[~, D] = knnsearch(MdlKDT, pts_gt');
error = mean(D);

function error = re(R_est, R_gt)
%     """
%     Rotational Error.
%
%     :param R_est: Rotational element of the estimated pose (3x1 vector).
%     :param R_gt: Rotational element of the ground truth pose (3x1 vector).
%     :return: Error of t_est w.r.t. t_gt.
%     """

error_cos = 0.5 * (trace(R_est * inv(R_gt)) - 1.0);
error_cos = min(1.0, max(-1.0, error_cos));
error = acos(error_cos);
error = 180.0 * error / pi;

function error = te(t_est, t_gt)
% """
% Translational Error.
%
% :param t_est: Translation element of the estimated pose (3x1 vector).
% :param t_gt: Translation element of the ground truth pose (3x1 vector).
% :return: Error of t_est w.r.t. t_gt.
% """
error = norm(t_gt - t_est);

author: yuxng


text:
function dicomrt_plotcontourc(data)
% dicomrt_plotcontourc(data)
%
% Plot contours created by contourc.
%
% See also: contourc
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% scan and plot data
start=1;
highestpair=0;
while (length(data)-highestpair)~=0
    npair=data(2,start);
    lowestpair=start+1;
    highestpair=start+npair;
    plot(data(1,lowestpair:highestpair),data(2,lowestpair:highestpair),'w','LineWidth',2);
    start=highestpair+1;
end

author: cerr

repo: CERR


text:
function exportZtdLittleR(ztd_map, lat_map, lon_map, h_map, time, out_dir, file_name, qc_ztd)
% Export a ZTD map into little_R format (input of 3Dvar)
%
% EXAMPLE
%   exportZtdLittleR(ztd_map, lat_map, lon_map, h_map, time, 'reserved/tmp', 'sar_map', 0.015)
%
% SYNTAX
%   exportZtdLittleR(ztd_map, lat_map, lon_map, h_map, time, out_dir, file_name, qc_ztd)
%
%
% info on LITTLE R export
% reference: http://www2.mmm.ucar.edu/wrf/users/wrfda/OnlineTutorial/Help/littler.html#FM

% qc_ztd = 0.015;
% ztd_map = ones(1,1);
n_row = size(ztd_map, 1);
n_col = size(ztd_map, 2);
% lat_map = randn(n_row,n_col);
% lon_map = 10 + randn(n_row,n_col);
% h_map = 100 + randn(n_row,n_col);
% time = GPS_Time(now);
% out_dir = 'reserved/tmp';
% file_name = 'sar_map';
qc_ztd = uint32(qc_ztd * 1e4); % expressed in tenths of millimeters

log = Core.getLogger();
save_on_disk = true;
no_data = -888888.00000;
qc_no_data = -888888;
end_data = -777777.00000;
end_qc = 0;
empty_data_line = sprintf('%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d\n', ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data, ...
    no_data, qc_no_data);
end_data_line = sprintf('%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d\n', ...
    end_data, end_qc, ...
    end_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc, ...
    no_data, end_qc);
tail_integers = sprintf('%7d%7d%7d\n',1,0,0);

% LITTLE R export
txt = '';
if save_on_disk
    if ~exist(out_dir, 'dir')
        log.newLine();
        log.addWarning(sprintf('The out folder does not exists\n Creating %s', out_dir));
        mkdir(out_dir);
    end
end

if save_on_disk
    file_name = sprintf('%s%c%s', out_dir, filesep, file_name);
    fid = fopen(file_name,'w');
end
i = 0;
fprintf('\n');
wb = Go_Wait_Bar.getInstance(n_col * n_row);
wb.createNewBar(sprintf('Exporting ZTD map into LITTLE_R format (%s)', time.toString('yyyy - mm - dd')));
for c = 1 : n_col
    % get the closest epoch to the sar passage
    for r = 1 : n_row
        i = i+1;
        % reference: http://www2.mmm.ucar.edu/wrf/users/wrfda/OnlineTutorial/Help/littler.html#FM
        lat = lat_map(r,c);
        lon = lon_map(r,c);
        id_dtr = n_row * (c-1) + (r-1);
        name = sprintf('CALSAR_%05d', id_dtr);
        platform = 'FM-114 GPSZD';
        source = 'goGPS software';
        h_ortho = h_map(r, c);
        obs_time = time.toString('yyyymmddHHMMSS');
        ztd = ztd_map(r,c) * 100; % expressed in cm

        if ~isnan(ztd)
            %Latitude	Longitude	ID	Name	Platform (FM�Code)note	Source	Elevation	Valid fields	Num. errors	Num. warnings	Sequence number	Num. duplicates	Is sounding?	Is bogus?	Discard?	Unix time	Julian day	Date	SLP, QC	Ref Pressure, QC	Ground Temp, QC	SST,   QC	SFC Pressure, QC	Precip, QC	Daily Max T, QC	Daily Min T, QC	Night Min T, QC	3hr Pres Change, QC	24hr Pres Change, QC	Cloud cover, QC	Ceiling, QC	Precipitable water, QC (see note)
            %F20.5	    F20.5	    A40	A40	    A40	                    A40	    F20.5	    I10            	I10	        I10	            I10	            I10	            L	            L	        L	        I10	        I10	        A20	    F13.5,  I7	F13.5,    I7	F13.5,       I7	F13.5, I7	F13.5,        I7	F13.5,  I7	F13.5,       I7	F13.5,       I7	F13.5,       I7	F13.5,           I7	F13.5,            I7	F13.5,       I7	F13.5,   I7	F13.5,              I7

            txt = sprintf('%s%20.5f%20.5f%-40s%-40s%-40s%-40s%20.5f%10d%10d%10d%10d%10d%10s%10s%10s%10d%10d%20s%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d%13.5f%7d\n%s%s%s', ...
                txt, ...
                lat, ...
                lon, ...
                num2str(id_dtr), ...
                name, ...
                platform, ...
                source, ...
                h_ortho, ...
                1, ...
                0, ...
                0, ...
                i, ...
                0, ...
                'F', ...
                'T', ...
                'F', ...
                uint32(time.getUnixTime), ...
                uint32(time.getJD), ...
                obs_time, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                no_data, qc_no_data, ...
                ztd, qc_ztd, ...
                empty_data_line, ...
                end_data_line, ...
                tail_integers);
        end
        if save_on_disk && (length(txt) > (4 * 1024))
            wb.go(r + (c-1) * n_row);
            fprintf(fid,'%s', txt);
            txt = [];
        end
    end
    wb.go(r + (c-1) * n_row);
    if save_on_disk
        fprintf(fid,'%s', txt);
        txt = [];
    end
end
wb.close();
if save_on_disk
    fclose(fid);
end

author: goGPS-Project


text:
function [H, xx, xw]= plot_barWithAntenna(y, y_sem, varargin)
%PLOT_BARWITHANTENNA - Bar plot with antennas, e.g., to show the SEM.
%
%Synopsis:
% [H, X]= plot_barWithAntenna(Y, Y_SEM, <OPT>)
%
%Input:
% Y:     Heights of the bars, DOUBLE [MxN]
% Y_SEM: Length of the antenna, DOUBLE [MxN]
% OPT:   Structure or property/value list of optional properties:
%  .Clf          - Clear figure before plotting, BOOL (default 1)
%  .ColOrder     - Colormap of N colors, default blue to white, DOUBLE [Nx3]
%                  (generated by cmap_bluewhite(N, 'min_val',0.15)
%  .AntennaType  - Specifies how anntennas are drawn, options are 'o','oo',
%                  'T','I', CHAR (default 'o')
%  .AntennaColor - Color of antennas, with 'auto' meaning same colors as
%                  bars, DOUBLE [Nx3] or 'auto' (default)
%  .AntennaSpec  - Line specifications for antennas
%
%Output:
% H: handles to graphic objects
% X: horizonal position of the middle of the bars
%
%Description:
% This function makes a bar plot in which a (symmetric) interval around
% each value is visualized by antennas. The typically application is
% that the values itself are the means of some observed variable and
% the intervals around each mean are given as the standard deviation
% (STD) or the standard error of the mean (SEM).
%
%Examples:
% Y= diag([5 6 2])*rand(3,6)+1;
% Y_SEM= rand(3,6)+0.2;
% H= bar_with_antenna(Y, Y_SEM);
% set(gca, 'XTickLabel','Group 1|Group 2|Group 3');
% legend(H.bar, str_cprintf('bar %d',1:6));
%
% [H, xx]= bar_with_antenna(Y(1,:), Y_SEM(1,:), 'ColOrder',cmap_rainbow(6), 'AntennaColor',[0 0 0]);
% set(gca, 'XTick',xx, 'XTickLabel',1:6);
% H= bar_with_antenna(Y, Y_SEM, 'AntennaType','oo');
%
% H= bar_with_antenna(Y(1,:), Y_SEM(1,:), 'AntennaType','I', 'AntennaColor','auto', 'DrawMarker',1);
% delete(H.bar);

% Benjamin Blankertz
% Modified Laura Acqualagna 04.15: Adjusted for handling new graphcs of Matlab versions later than 8.4

[nGroups, nCol]= size(y);
cmap_default= cmap_bluewhite(nCol, 'MinVal',0.15);

props = {'Clf',             1,              '!BOOL';
         'ColOrder',        cmap_default,   '!DOUBLE[- 3]';
         'AntennaType',     'o',            '!CHAR(o oo T TT I none)';
         'AntennaColor',    'auto',         '!CHAR(auto)|DOUBLE[- 3]';
         'AntennaSpec',     {},             'CELL';
         'AntennaWidth',    .5,             'DOUBLE[1]';
         'DrawMarker',      0,              '!BOOL';
         'MarkerSpec',      {},             'CELL';
         };

if nargin==0,
  H= props; return
end

misc_checkType(y,'!DOUBLE[- -]');
misc_checkType(y_sem,'!DOUBLE[- -]');
opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if nargin<2 || isempty(y_sem),
  opt.AntennaType= 'none';
end

if ismember(opt.AntennaType, {'oo','TT','I'},'legacy'),
  [opt, isdefault]= opt_overrideIfDefault(opt, isdefault, ...
                                          'AntennaColor',0.5*[1 1 1]);
end
if isequal(opt.AntennaColor, 'auto'),
  opt.AntennaColor= opt.ColOrder;
end

if opt.Clf,
  clf;
end

if nGroups==1,  %% work-around for case nGroups=1
  H.bar= bar([y; y]);
else
  H.bar= bar(y);
end

hold on;
for c= 1:nCol,
 if verLessThan('matlab', '8.4') % old Matlab versions before 2014b
 xdata= get(get(H.bar(c),'Children'),'XData');
 xw= xdata(3,1)-xdata(1,1);
 xtmp= mean(xdata);
 else
 set(H.bar(c), 'FaceColor',opt.ColOrder(c,:));
 xtmp =  H.bar(c).XData + [H.bar(c).XOffset];
 xw=H.bar(c).BarWidth/(length(H.bar)*2-1); % BarWidth refers to the overall width of the group.
 %To get width of single bar, divide BarWidth for number of bars + number of spaces between bars.
 end
 if nGroups==1,  %% work-around for case nGroups=1
   xtmp= xtmp(1);
 end
 xx(c,:)= xtmp;
 xtw= xw*opt.AntennaWidth/2;
 switch(opt.AntennaType),
  case 'o',
   h1= line([xx(c,:); xx(c,:)], ...
            [y(:,c), y(:,c)+y_sem(:,c)]');
   h2= line([xx(c,:); xx(c,:)], ...
            [y(:,c)+y_sem(:,c), y(:,c)+y_sem(:,c)]');
   H.antenna(:,c)= [h1; h2];
  case 'T',
   h1= line([xx(c,:); xx(c,:)], ...
            [y(:,c), y(:,c)+y_sem(:,c)]');
   h2= line([xx(c,:)-xtw; xx(c,:)+xtw], ...
            [y(:,c)+y_sem(:,c), y(:,c)+y_sem(:,c)]');
   H.antenna(:,c)= [h1; h2];
  case 'oo',
   H.antenna(:,c)= line([xx(c,:); xx(c,:)], ...
                        [y(:,c)-y_sem(:,c), y(:,c)+y_sem(:,c)]');
  case {'I','TT'},
   h1= line([xx(c,:); xx(c,:)], ...
            [y(:,c)-y_sem(:,c), y(:,c)+y_sem(:,c)]');
   h2= line([xx(c,:)-xtw; xx(c,:)+xtw], ...
            [y(:,c)+y_sem(:,c), y(:,c)+y_sem(:,c)]');
   h3= line([xx(c,:)-xtw; xx(c,:)+xtw], ...
            [y(:,c)-y_sem(:,c), y(:,c)-y_sem(:,c)]');
   H.antenna(:,c)= [h1; h2; h3];
  case 'none',
   H.antenna(:,c)= NaN;
 end
 if verLessThan('matlab', '8.4')
     if ~any(isnan(H.antenna(:,c))),
       ci= 1+mod(c-1, size(opt.AntennaColor,1));
       set(H.antenna(:,c), 'Color',opt.AntennaColor(ci,:));
     end
 else %in the new version the previous check cannot be done in one line, if condition has been removed
       ci= 1+mod(c-1, size(opt.AntennaColor,1));
       set(H.antenna(:,c), 'Color',opt.AntennaColor(ci,:));
 end
end

switch(opt.AntennaType),
 case 'o',
  ii= nGroups+[1:nGroups];
  set(H.antenna(ii,:), 'Marker','o', 'MarkerSize',8, 'MarkerFaceColor','w');
 case 'oo',
  set(H.antenna, 'Marker','o', 'MarkerSize',8, 'MarkerFaceColor','w');
end
 if verLessThan('matlab', '8.4')
     if ~any(isnan(H.antenna(:,c))),
      set(H.antenna, 'LineWidth',4, opt.AntennaSpec{:});
     end
 else
     set(H.antenna, 'LineWidth',4, opt.AntennaSpec{:});
 end
set(gca, 'YGrid','on', 'TickLength',[0 0]);

if opt.DrawMarker,
  for c= 1:nCol,
    H.marker(:,c)= plot(xx(c,:)', y(:,c), 'd');
    ci= 1+mod(c-1, size(opt.AntennaColor,1));
    set(H.marker(:,c), 'Color',opt.AntennaColor(ci,:));
  end
  set(H.marker, 'MarkerSize',12, 'LineWidth',4, 'MarkerFaceColor','w', ...
                opt.AntennaSpec{:}, opt.MarkerSpec{:});
end

% Adjust XLim such that there is one bar-width space at the borders,
set(gca, 'XLim',reshape(xx([1 end]),[1 2])+[-1 1]*xw*1.5);
hold off;

author: bbci


text:
clc;
pzstat = input('enter a puzzle in 3 x 3 form (use 0 for space) ..........:   ');
pzrstat = [1 2 3;4 5 6;7 8 0];
intv000 = size(pzstat);
if not(all(intv000 == [3 3]))
    msgbox('Input is not in correct format.............. The format must be 3 x 3 matrix. The blank space should be replaced by 0.','8-Puzzle suresh Kumar Gadi');
else
    pzinl = [];
    pzhis = [];
    pzrhis = [];
    intv00g = 0;
    [intv0h1,intv0h2] = gskcostastar(pzstat,pzrstat);
    intv00h = intv0h1 + intv0h2;
    intv00f = intv00h + intv00g;
    pzinl = gskpzjoin(pzinl,pzstat,intv00f,intv00g,intv00h);
    pzhis = gskpzjoin(pzhis,pzstat,intv00f,intv00g,intv00h);
    pzhis(:,13)=0;
    pzhis(1,14)=0;
    pzhis(1,15)=1;
    intv001 = 1;
    intv0ct = 0;
    intv0nd = 1;
    intv1nd = 1;
    while intv001 == 1;
        intv0ct = intv0ct + 1;
        [pzout, intv00f, intv00g, intv00h] =  gskpzget(pzinl,1);
        [intv002, intv003] = find(pzout==0);
        pzinl = [];
        if intv00h ~= 0
% % % % % % left
            if intv003 > 1
                pzout1 = pzout;
                pzout1(intv002,intv003) = pzout(intv002,intv003 - 1);
                pzout1(intv002,intv003 - 1) = 0;
                [intv0h1,intv0h2] = gskcostastar(pzout1,pzrstat);
                intv00h = intv0h1 + intv0h2;
                intv0g1 = intv00g + 1;
                intv00f = intv0g1 + intv00h;
                [pzinl] =  gskpzjoin(pzinl,pzout1,intv00f,intv0g1,intv00h);
            end
% % % % % % top
            if intv002 > 1
                pzout1 = pzout;
                pzout1(intv002,intv003) = pzout(intv002-1,intv003);
                pzout1(intv002-1,intv003) = 0;
                [intv0h1,intv0h2] = gskcostastar(pzout1,pzrstat);
                intv00h = intv0h1 + intv0h2;
                intv0g1 = intv00g + 1;
                intv00f = intv0g1 + intv00h;
                [pzinl] =  gskpzjoin(pzinl,pzout1,intv00f,intv0g1,intv00h);
            end
% % % % % % right
            if intv003 < 3
                pzout1 = pzout;
                pzout1(intv002,intv003) = pzout(intv002,intv003 + 1);
                pzout1(intv002,intv003 + 1) = 0;
                [intv0h1,intv0h2] = gskcostastar(pzout1,pzrstat);
                intv00h = intv0h1 + intv0h2;
                intv0g1 = intv00g + 1;
                intv00f = intv0g1 + intv00h;
                [pzinl] =  gskpzjoin(pzinl,pzout1,intv00f,intv0g1,intv00h);
            end
% % % % % % down
            if intv002 < 3
                pzout1 = pzout;
                pzout1(intv002,intv003) = pzout(intv002 + 1,intv003);
                pzout1(intv002 + 1,intv003) = 0;
                [intv0h1,intv0h2] = gskcostastar(pzout1,pzrstat);
                intv00h = intv0h1 + intv0h2;
                intv0g1 = intv00g + 1;
                intv00f = intv0g1 + intv00h;
                [pzinl] =  gskpzjoin(pzinl,pzout1,intv00f,intv0g1,intv00h);
            end
% % % % % % compare with history
            [intv003, intv004] = size(pzinl);
            for intv004 = 1 : intv003
                intv005 = gsksearch(pzhis,pzinl(intv004,:));
                if ~isempty(intv005)
%                     pzhis(intv005,13) = 1;
                else
                    [intv007 , intv008] = size(pzhis);
                    intv008 = pzinl(intv004,1:12);
                    intv008(13) = 0;
                    intv008(14) = intv0nd;
                    intv1nd = intv1nd + 1;
                    intv008(15) = intv1nd;
                    pzhis(intv007+1,:)=intv008;
                end
            end
% % % % % % selecting new node
            intv005 = find(pzhis(:,13) == 1);
            intv006 = length(intv005);
            pzrhis = pzhis;
            for intv007 = 1:intv006
                pzrhis(intv005(intv007)-(intv007-1),:)=[];
            end
            intv003 = min(pzrhis,[],1);
            [pzout, intv00f, intv00g, intv00h] =  gskpzget(intv003,1);
            [intv002, intv003] = find((pzrhis(:,10)==intv00f) & (pzrhis(:,13) == 0));
%             intv002 = intv002(1,1);
            [pzout, intv00f, intv00g, intv00h] =  gskpzget(pzrhis,intv002(1,1));
            intv0nd = pzrhis(intv002(1,1),15);
            pzinl = [];
            [pzinl] =  gskpzjoin(pzinl,pzout,intv00f,intv00g,intv00h);
            intv005 = gsksearch(pzhis,pzinl);
            pzhis(intv005,13) = 1;
            intv001 = 1;
        else
            intv001 = 2;
        end
    end
% % Display of steps
    [intv000, intv001] = size(pzhis);
    intv002 = 1;
    pzshow = [];
    intvcnt = 0;
    while intv002 ==1
        intvcnt = intvcnt + 1;
        pzshow(intvcnt) = intv000;
        intv000 = pzhis(intv000,14);
        if intv000 == 0
            intv002 = 2;
        else
            intv002 = 1;
        end
    end
    for intv000 = 1 : intvcnt
        [pzout, intv00f, intv00g, intv00h] =  gskpzget(pzhis,pzshow(intvcnt-intv000+1));
        intv000
        pzout
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function breadth = polyhedronMeanBreadth(vertices, edges, faces)
%POLYHEDRONMEANBREADTH Mean breadth of a convex polyhedron
%
%   BREADTH = polyhedronMeanBreadth(V, E, F)
%   Return the mean breadth (average of polyhedron caliper diameter over
%   all direction) of a convex polyhedron.
%
%   The mean breadth is computed using the sum, over the edges of the
%   polyhedron, of the edge dihedral angles multiplied by the edge length,
%   the final sum being divided by (4*PI).
%
%   Note: the function assumes that the faces are correctly oriented. The
%   face vertices should be indexed counter-clockwise when considering the
%   supporting plane of the plane, with the outer normal oriented outwards
%   of the polyhedron.
%
%   Typical values for classical polyhedra are:
%     cube side a               breadth = (3/2)*a
%     cuboid sides a, b, c      breadth = (a+b+c)/2
%     tetrahedron side a        breadth = 0.9123*a
%     octaedron side a          beradth = 1.175*a
%     dodecahedron, side a      breadth = 15*arctan(2)*a/(2*pi)
%     icosaehdron, side a       breadth = 15*arcsin(2/3)*a/(2*pi)
%
%   Example
%   [v e f] = createCube;
%   polyhedronMeanBreadth(v, e, f)
%   ans =
%       1.5
%
%   See also
%   meshes3d, meshEdgeFaces, meshDihedralAngles, checkMeshAdjacentFaces
%   trimeshMeanBreadth
%
%   References
%   Stoyan D., Kendall W.S., Mecke J. (1995) "Stochastic Geometry and its
%       Applications", John Wiley and Sons, p. 26
%   Ohser, J., Muescklich, F. (2000) "Statistical Analysis of
%       Microstructures in Materials Sciences", John Wiley and Sons, p.352

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% compute dihedral angle of each edge
alpha = meshDihedralAngles(vertices, edges, faces);

% compute length of each edge
lengths = meshEdgeLength(vertices, edges);

% compute product of length by angles
breadth = sum(alpha.*lengths)/(4*pi);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% LORENZ system
% Train ensemble of models (e.g. for varying training data length, noise
% level, etc.)
% System identification: SINDYc

clear all, close all, clc

SystemModel = 'LORENZ';
DERIV_NOISE = 0;
TRACK_MODEL_BEST = 0;
%% Set Paths
figpath = '../FIGURES/EX_LORENZ_Dependencies/SINDYc/'; mkdir(figpath)
datapath = '../DATA/EX_LORENZ_Dependencies/SINDYc/'; mkdir(datapath)

% Overwrite if necessary
if exist('DERIV_NOISE')
    if DERIV_NOISE == 1
        figpath = '../FIGURES/EX_LORENZ_Dependencies/SINDYc/TVRegDiff/'; mkdir(figpath)
        datapath = '../DATA/EX_LORENZ_Dependencies/SINDYc/TVRegDiff/'; mkdir(datapath)
    end
end
addpath('../utils');

%% Parameters
ModelName = 'SINDYc';
polyorder = 2;
usesine = 0;
lambda0 = 0.1;     % lambda is our sparsification knob.
dep_trainlength = 1;
dep_noise = 0;

MOD_VAL = 10;

%% Select case 1 or 2
% 1) Dependency on training length
Ntrain_vec = [5:15,20:5:95,100:100:1000];%,1250,1500,2000,3000];%,1500:500:3000];
eta_vec = 0;%0.05;
Nr = 1;%50;

% 2) Noise dependency
% Ntrain_vec = 3000;
% eta_vec = [0.01 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5];
% Nr = 50;

N_LENGTHS = length(Ntrain_vec);
N_ETA = length(eta_vec);

%% Parameters
SAVE_MODEL = 0;
SHOW_RESULTS = 0;
SHOW_STATS = 1;
ONLY_TRAINING_LENGTH = 1; % if 0 : 3000 time unit, otherwise 1000

if ONLY_TRAINING_LENGTH == 1
    PostName = ['_TrainLength'];
else
    PostName = [];
end
%% Generate Data %{'sine2', 'chirp','prbs', 'sphs'}
InputSignalType = 'sphs'; %prbs; chirp; noise; sine2; sphs; mixed
getTrainingData

DataTrain.x = x;
DataTrain.t = t;
DataTrain.tspan = tspan;
DataTrain.u = u;

xstd = std(xv(:,3));

close all
rng(0,'twister')

errBest = inf*ones(N_LENGTHS,N_ETA);
% BestModels(1:N_LENGTHS,1:N_ETA) = struct('name', [], 'polyorder', [], 'usesine', [], 'Xi', [], ...
%     'dt', [], 'N', [], 'Ttraining', [], 'Err', [], 'ErrM', []);
BestModelsList = zeros(N_LENGTHS,N_ETA);
Lambda = zeros(N_LENGTHS,N_ETA);
Nvar = 3;
%% SINDYc
if ONLY_TRAINING_LENGTH == 0
    for iN = 1:N_LENGTHS
        for iNoise = 1:N_ETA
            disp(['Running for noise case ', num2str(iNoise), ' of ', num2str(N_ETA)])

            Results = struct('err', zeros(Nr,1), 'errM', zeros(Nr,1), 'xA', zeros(length(tv),Nvar), 'xB', zeros(length(tv),Nvar,Nr),'Ttraining',zeros(Nr,1));

            for iR = 1:Nr

                % Setup data
                x = DataTrain.x(1:Ntrain_vec(iN),:);
                u = DataTrain.u(1:Ntrain_vec(iN));
                t = DataTrain.t(1:Ntrain_vec(iN));
                tspan = DataTrain.tspan(1:Ntrain_vec(iN));

                % Add noise
                eps = eta_vec(iNoise)*xstd;
                x = x + eps*randn(size(x));

                % Train model
                lambda = lambda0;
                tic
                trainSINDYc
                telapsed = toc

                % Prediction over training phase
                p.ahat = Xi(:,1:3);
                p.polyorder = polyorder;
                p.usesine = usesine;
                p.dt = dt;
                [N,Ns] = size(DataTrain.x);
                xSINDYc = zeros(Ns,N); xSINDYc(:,1) = x0';
                for ct=1:N-1
                    xSINDYc(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xSINDYc(:,ct),DataTrain.u(ct),dt,1,[],p);
                end
                xSINDYc = xSINDYc';

                % Show validation
                SHOW_PREDICTION_FOR_TRAINING_PHASE

                %% Prediction
                % Reference
                xA      = [xv];
                tA      = tv;

                % Model
                [N,Ns] = size(xA);
                xB = zeros(Ns,N); xB(:,1) = DataTrain.x(end,:)';
                for ct=1:N
                    xB(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xB(:,ct),uv(ct),dt,1,[],p);
                end
                xB = xB(:,1:N+1)';
                tB = [tv(1)-dt;tv];

                % Show training and prediction
                SHOW_PREDICTION_FOR_VALIDATION_PHASE

                %% Error
                err = mean(sum((xB(2:end,:)-xA).^2,2));
                errM = mean(sum((xB(2:250+1,:)-xA(1:250,:)).^2,2));

                %% Save Data
                Model.name = 'SINDYc';
                Model.polyorder = polyorder;
                Model.usesine = usesine;
                Model.Xi = Xi;
                Model.dt = dt;
                Model.N  = Ntrain_vec(iN);
                Model.Ttraining = telapsed;
                Model.Err = err;
                Model.ErrM = errM;

                if SAVE_MODEL == 1
                    if mod(iR,MOD_VAL) == 0 || iR == 1

                        save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_iR',num2str(iR),'.mat']),'Model')
                    end
                end

                Results.err(iR) = err;
                Results.errM(iR) = errM;
                Results.xA = xA;
                Results.xB(:,1:3,iR) = xB(2:end,:);
                Results.Ttraining(iR) = telapsed;

                %% Track best model
                % errBest = 10^12*ones(N_LENGTHS,N_ETA)
                if TRACK_MODEL_BEST == 1
                    if Results.err(iR)<errBest(iN,iNoise) || iR == 1
                        errBest(iN,iNoise) = Results.err(iR);

                        BestModel = Model;
                        save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_BEST_MODEL.mat']),'Model')

                        %BestModels(iN,iNoise) = Model;
                        BestModelsList(iN,iNoise) = iR;
                        Lambda(iN,iNoise) = lambda;

                        clear ph
                        f1 = figure('visible','off');box on, hold on,
                        ccolors = get(gca,'colororder');
                        plot([tB(1),tB(1)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
                        plot([t(end),t(end)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
                        ylim([-25 65])
                        text(5,55,'Training', 'FontSize',12)
                        %                     text(10+tA(1),230,'Validation', 'FontSize',12)
                        text(3+tA(1),55,'Validation', 'FontSize',12)

                        if eps~=0
                            ph(4) = plot(tspan,x(:,1),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(1,:)+[0.15 0.3 0.25]
                            plot(tspan,x(:,2),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(2,:)+[0.15 0.3 0.25]
                            plot(tspan,x(:,3),'-','Color',0.7*ones(1,3),'LineWidth',1);
                            ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',0.5);
                            ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',0.5);
                            ph(3) = plot([DataTrain.t;tA],[DataTrain.x(:,3);xA(:,3)],'-','Color',ccolors(3,:),'LineWidth',0.5);
                        else
                            ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',1);
                            ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',1);
                            ph(3) = plot(t,x(:,1),'--','Color',[0 1 0],'LineWidth',1); % Training data
                            plot(t,x(:,2),'--','Color',[0 1 0],'LineWidth',1);
                        end

                        ph(5) = plot(tB,xB(:,1),'-.','Color',ccolors(1,:)-[0 0.2 0.2],'LineWidth',2);
                        ph(6) = plot(tB,xB(:,2),'-.','Color',ccolors(2,:)-[0.1 0.2 0.09],'LineWidth',2);
                        ph(7) = plot(tB,xB(:,3),'-.','Color',ccolors(3,:)-[0.1 0.2 0.09],'LineWidth',2);
                        grid off
                        xlim([0 tv(end)])
                        xlabel('Time')
                        ylabel('xi')
                        set(gca,'LineWidth',1, 'FontSize',14)
                        set(gcf,'Position',[100 100 300 200])
                        set(gcf,'PaperPositionMode','auto')
                        print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_BEST_MODEL_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_noleg.eps']);

                        if eps~=0
                            lh = legend(ph([1,4,5]),'Truth','Training',ModelName,'Location','NorthWest');
                            %                         lh.Position = [lh.Position(1)+0.13,lh.Position(2)-0.2,lh.Position(3:4)];
                            lh.Position = [lh.Position(1)+0.02,lh.Position(2)-0.06,lh.Position(3:4)];
                        else
                            lh = legend(ph([1,4,5]),'Truth','Training',ModelName,'Location','NorthWest');
                            lh.Position = [lh.Position(1)+0.13,lh.Position(2)-0.2,lh.Position(3:4)];
                        end
                        print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_BEST_MODEL_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'.eps']);
                        close(f1);
                    end
                end
                lambda = lambda0;
            end
            if SHOW_STATS == 1
                %%
                xBmin = min(Results.xB(:,1:3,:),[],3);
                xBmax = max(Results.xB(:,1:3,:),[],3);
                clear ph
                figure('visible','off'),box on, hold on,
                ccolors = get(gca,'colororder');
                plot([tB(1),tB(1)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
                plot([t(end),t(end)],[-25 65],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
                ylim([-25 65])
                t1 = text(5,55,'Training', 'FontSize',12);
                t2 = text(5+tA(1),55,'Validation', 'FontSize',12);

                X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
                Y=[xBmin(:,1)',flipud(xBmax(:,1))'];              %#create y values for out and then back
                fillh1 = fill(X,Y,ccolors(1,:)-[0 0.2 0.2]);                  %#plot filled area
                fillh1.EdgeColor = ccolors(1,:)-[0 0.2 0.2]; fillh1.FaceAlpha = 0.5;

                X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
                Y=[xBmin(:,2)',flipud(xBmax(:,2))'];              %#create y values for out and then back
                fillh = fill(X,Y,ccolors(2,:)-[0.1 0.2 0.09]);                  %#plot filled area
                fillh.EdgeColor = ccolors(2,:)-[0.1 0.2 0.09]; fillh.FaceAlpha = 0.5;

                X=[tB(2:end)',fliplr(tB(2:end)')];                %#create continuous x value array for plotting
                Y=[xBmin(:,3)',flipud(xBmax(:,3))'];              %#create y values for out and then back
                fillh = fill(X,Y,ccolors(3,:)-[0.1 0.2 0.09]);                  %#plot filled area
                fillh.EdgeColor = ccolors(3,:)-[0.1 0.2 0.09]; fillh.FaceAlpha = 0.5;

                if eps~=0
                    ph(4) = plot(tspan,x(:,1),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(1,:)+[0.15 0.3 0.25]
                    plot(tspan,x(:,2),'-','Color',0.7*ones(1,3),'LineWidth',1); %ccolors(2,:)+[0.15 0.3 0.25]
                    plot(tspan,x(:,3),'-','Color',0.7*ones(1,3),'LineWidth',1);
                    ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',1);
                    ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',1);
                    ph(3) = plot([DataTrain.t;tA],[DataTrain.x(:,3);xA(:,3)],'-','Color',ccolors(3,:),'LineWidth',1);
                else
                    ph(1) = plot([DataTrain.t;tA],[DataTrain.x(:,1);xA(:,1)],'-','Color',ccolors(1,:),'LineWidth',1);
                    ph(2) = plot([DataTrain.t;tA],[DataTrain.x(:,2);xA(:,2)],'-','Color',ccolors(2,:),'LineWidth',1);
                    ph(3) = plot([DataTrain.t;tA],[DataTrain.x(:,3);xA(:,3)],'-','Color',ccolors(3,:),'LineWidth',1);
                    ph(4) = plot(t,x(:,1),'--','Color',[0 1 0],'LineWidth',1); % Training data
                    plot(t,x(:,2),'--','Color',[0 1 0],'LineWidth',1);
                    plot(t,x(:,3),'--','Color',[0 1 0],'LineWidth',1);
                end

                grid off
                xlim([0 tv(end)])
                xlabel('Time')
                ylabel('xi')
                set(gca,'LineWidth',1, 'FontSize',14)
                set(gcf,'Position',[100 100 300 200])
                set(gcf,'PaperPositionMode','auto');

                if eps~=0
                    lh = legend([ph([1,4]),fillh1],'Truth','Training',ModelName,'Location','NorthWest');
                    lh.Position = [lh.Position(1)+0.02,lh.Position(2)-0.06,lh.Position(3:4)];
                else
                    lh = legend(ph([1,4]),'Truth','Training',ModelName,'Location','NorthWest');
                    lh.Position = [lh.Position(1)+0.02,lh.Position(2)-0.06,lh.Position(3:4)];
                end

                print('-depsc2', '-painters','-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.eps']);

                delete(lh), delete(t1), delete(t2)
                print('-depsc2', '-painters','-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS_noleg.eps']);
            end
            %%
            save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.mat']),'Results')
        end
    end
    % % save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_BEST_MODELS.mat']),'BestModels')
    save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_BEST_MODELS_LIST',PostName,'.mat']),'BestModelsList')
    save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_BEST_MODELS_Lambda',PostName,'.mat']),'Lambda')
end

%%
if ONLY_TRAINING_LENGTH == 1
    length_vec = 5:2:Ntrain-1;
    Results = struct('err', zeros(length(length_vec),2),'PredHor_Ball', zeros(length(length_vec),2), 'PredLength', zeros(length(length_vec),2),'RelErr', zeros(length(length_vec),2),'RelErrMax', zeros(length(length_vec),2), 'Ttraining', zeros(length(length_vec),1));
    Results.Ntrain_vec = length_vec;
    for iL = 1:length(length_vec)
        iN = length_vec(iL);
        if mod(iN,50)
            disp(['Running for length ', num2str(iN), 'of ', num2str(Ntrain-1)])
        end
        % Setup data
        x = DataTrain.x(1:iN,:);
        u = DataTrain.u(1:iN);
        t = DataTrain.t(1:iN);

        % Train model
        lambda = lambda0;
        tic
        trainSINDYc
        telapsed = toc

        % Prediction over training phase
        p.ahat = Xi(:,1:3);
        p.polyorder = polyorder;
        p.usesine = usesine;
        p.dt = dt;
        [N,Ns] = size(DataTrain.x);
        xSINDYc = zeros(Ns,N); xSINDYc(:,1) = x0';
        for ct=1:N-1
            xSINDYc(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xSINDYc(:,ct),DataTrain.u(ct),dt,1,[],p);
        end
        xSINDYc = xSINDYc';

        %% Error over training phase
        Results.err(iL,1) = mean(sum((xSINDYc(1:end,:)-DataTrain.x).^2,2));
        Results.RelErr(iL,1) = mean( sum( abs( (DataTrain.x - xSINDYc)./DataTrain.x ),2) );
        Results.RelErrMax(iL,1) = mean( max( abs( (DataTrain.x - xSINDYc)./DataTrain.x ),[],2) );

        tmp = abs( (DataTrain.x - xSINDYc)./DataTrain.x );
        TF = tmp>0.1;
        idx = [];
        for iVar = 1:Nvar
            idx = [idx ; find(TF(:,iVar)==1,1,'first')];
        end
        Results.PredLength(iL,1) = max(idx);

        tmp2 = sqrt(sum(abs(DataTrain.x - xSINDYc).^2,2));
        idx = find(tmp2>3,1,'first');
        Results.PredHor_Ball(iL,1) = idx-1;
        %% Prediction
        % Reference
        xA      = [xv];
        tA      = tv;

        % Model
        [N,Ns] = size(xA);
        xB = zeros(Ns,N); xB(:,1) = DataTrain.x(end,:)';
        for ct=1:N
            xB(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xB(:,ct),uv(ct),dt,1,[],p);
        end
        xB = xB(:,1:N+1)';
        tB = [tv(1)-dt;tv];

        %% Error over prediction phase
        xB = xB(2:end,:);
        Results.err(iL,2) = mean(sum((xB-xA).^2,2));
        Results.RelErr(iL,2) = mean( sum( abs( (xA - xB)./xA ),2) );
        Results.RelErrMax(iL,2) = mean( max( abs( (xA - xB)./xA ),[],2) );
        Results.Ttraining(iL) = telapsed;

        tmp = abs( (xA - xB)./xA );
        TF = tmp>0.1;
        idx = [];
        for iVar = 1:Nvar
            idx = [idx ; find(TF(:,iVar)==1,1,'first')];
        end
        Results.PredLength(iL,2) = max(idx);

        tmp2 = sqrt(sum(abs(xA - xB).^2,2));
        idx = find(tmp2>3,1,'first');
        Results.PredHor_Ball(iL,2) = idx-1;
    end

    Results.DataTrain = DataTrain;
    Results.DataValid.x = xv;
    Results.DataValid.u = uv;
    Results.DataValid.t = tv;
    Results.DataValid.tspan = tspanv;
    save(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_Nevol','_STATS.mat']),'Results')
end

author: eurika-kaiser

repo: SINDY-MPC