text:
function varargout = contourfcmap(x,y,z,clev,cmap,lo,hi,cbarloc)
%CONTOURFCMAP Filled contour plot with specified colors
%
% h = contourfcmap(x,y,z,clev,cmap,lo,hi,cbarloc)
% h = contourfcmap(x,y,z,clev,cmap)
% h = contourfcmap(x,y,z,clev,cmap,lo)
% h = contourfcmap(x,y,z,clev,cmap,lo,hi)
%
% This function creates a shaded contour map, similar to that created by
% the contourf function.  However, the relationship between a contourf plot
% and it's colormap (i.e. exactly which color corresponds to each contour
% interval), can often be confusing and inconsistent, in my opinion.  This
% function instead allows the user to specify exactly which colors to use
% in each interval, and also to choose colors for regions that exceed the
% contour line limits.
%
% Input variables
%
%   x:          x-coordinates of grid, following size restrictions of surf
%               command
%
%   y:          y-coordinates of grid, following size restrictions of surf
%               command
%
%   z:          array of data to be contoured, following size restritions
%               of surf command
%
%   clev:       vector of length n, contour levels, must be monotonically
%               increasing
%
%   cmap:       n-1 x 3 colormap array, specifying colors used to shade
%               each contour interval
%
%   lo:         1 x 3 colormap array, specifying color used for all data
%               falling below the first contour level.  If not included or
%               empty, default will be to white.
%
%   hi:         1 x 3 colormap array, specifying color used for all data
%               falling above the last contour level.  If not included or
%               empty, default will be to white.
%
%   cbarloc:    string specifying colorbar location (see colorbar),  or a
%               1 x 4 position vector for colorbar.  If not included, no
%               colorbar will be created.  Note that the colorbar created
%               is not a colorbar object, but just an axis with plotted
%               patches; it is for labeling only and is not linked to the
%               "peer axis" in any way.
%
% Output variables:
%
%   h:          1 x 1 structure
%
%               c:      contour matrix for filled contour plot
%
%               h:      handle to contourgroup
%
%               cbax:   handle to axis of colorbar
%
% Example:
%
% [x,y] = meshgrid(linspace(0,1,100));
% z = peaks(100);
% contourfcmap(x,y,z,[-5 -3 -2:.5:2 3 5],jet(12), ...
%              [.8 .8 .8], [.2 .2 .2], 'eastoutside')
%

% Copyright 2010 Kelly Kearney

%------------------------
% Parse input
%------------------------

error(nargchk(5,8,nargin));

% Check clev

if ~isvector(clev)
    error('clev must be a vector');
end
nlev = length(clev);

% Check cmap

if size(cmap,2) ~=3 || size(cmap,1) ~= (nlev-1) || any(cmap(:)<0|cmap(:)>1)
    error('cmap must be nlev-1 x 3 colormap array');
end

% Colors for too-low and too-high patches

if nargin < 6 || isempty(lo)
    lo = [1 1 1];
end

if nargin < 7 || isempty(hi)
    hi = [1 1 1];
end

if ~isequal(size(lo),[1 3])
    error('lo must be 1 x 3 colormap array');
end

if ~isequal(size(hi),[1 3])
    error('hi must be 1 x 3 colormap array');
end

% Colorbar location

if nargin == 8 && ~isempty(cbarloc)
    pos = {'north', 'south', 'east', 'west', 'northoutside', 'southoutside', 'eastoutside', 'westoutside'};
    if ischar(cbarloc)
        if ~any(strcmp(lower(cbarloc), pos))
            error('Unrecognizd colorbar position');
        end
    elseif ~isequal(size(cbarloc), [1 4]) || any(cbarloc > 1 | cbarloc < 0)
        error('cbarloc must be position string or  1 x 4 normalized position');
    end
    showcb = true;
else
    showcb = false;
end

% Axis

ax = gca;

%------------------------
% Create filled contour
% plot
%------------------------

[c,h] = contourf(x,y,z,clev);

hpatch = get(h, 'children');

cdata = cell2mat(get(hpatch, 'cdata'));

% Mark too-high contours

isabove = cdata == max(clev);

% Distinguish between too-lo contours and NaN contours

if ~any(isnan(z(:)))
    isbelow = isnan(cdata);
else
    warning('NaN in data; haven''t figured NaN vs lo yet');
    isbelow = isnan(cdata);
end

level = get(h, 'LevelList');
if ~isequal(level, clev)
    error('oops, something new with levels, check this');
end

[tf, idx] = ismember(cdata, clev);

isweirdextra = idx == 0 & ~isbelow & ~isabove; % Why do these show up?!

for ip = 1:length(hpatch)
    if isbelow(ip)
            set(hpatch(ip), 'facecolor', lo);
    elseif isabove(ip)
            set(hpatch(ip), 'facecolor', hi);
    elseif isweirdextra(ip)
        dist = cdata(ip) - clev;
        dist(dist<0) = Inf;
        [blah, imin] = min(dist);
        set(hpatch(ip), 'facecolor', cmap(imin,:));
    else
        set(hpatch(ip), 'facecolor', cmap(idx(ip),:));
    end
end

%------------------------
% Create colorbar
%------------------------

if showcb
    clevcbar = reshape(clev,1,[]);
    height = 0.1 * (max(clevcbar) - min(clevcbar));
    y1 = [clev(1)-height; clev(1); clev(1); clev(1)-height; clev(1)-height];
    y2 = [clev(end); clev(end)+height; clev(end)+height; clev(end); clev(end)];

    yp = [y1 [clev(1:end-1); clev(2:end); clev(2:end); clev(1:end-1); clev(1:end-1)] y2];
    xp = [0;0;1;1;0] * ones(1,nlev+1);
    cp = [lo; cmap; hi];
    cp = permute(cp, [3 1 2]);

    if ~ischar(cbarloc)
        cbarcoord = cbarloc;
        hascbcoord = true;
        if cbarcoord(3)-cbarcoord(1) < cbarcoord(4)-cbarcoord(2)
            isvert = true;
            cbarloc = 'east';
        else
            isvert = false;
            cbarloc = 'south';
        end
    else
        hascbcoord = false;
        if any(strcmp(pos([3 4 7 8]), lower(cbarloc)))
            isvert = true;
        else
            isvert = false;
        end
    end

    if ~isvert
        tmp = xp;
        xp = yp;
        yp = tmp;
    end

    cbax = colorbar(cbarloc);
    axpos = get(ax, 'position');
    cbpos = get(cbax, 'position');

    delete(cbax);
    set(ax, 'position', axpos);

    cbax = axes('position', cbpos);
    set(cbax, 'userdata', 'contourfcolorbar');

    patch(xp,yp,cp);
    if isvert
        set(cbax, 'ytick', clev, 'ylim', minmax(yp), 'xlim', [0 1], 'xtick', []);
    else
        set(cbax, 'xtick', clev, 'xlim', minmax(xp), 'ylim', [0 1], 'ytick', []);
    end
end

if showcb && hascbcoord
    set(cbax, 'position', cbarcoord);
end

%------------------------
% Output
%------------------------

hndl.c = c;
hndl.h = h;
if showcb
    hndl.cbax = cbax;
end

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

% Return focus to axis (needed if colorbar created)

axes(ax);

function a = minmax(b)
a = [min(b(:)) max(b(:))];

author: Sable

repo: mcbench-benchmarks


text:
function [ozf] = dyn2ozf(dyn)
% Convert force from dyne to ounces (force).
% Chad A. Greene 2012
ozf = dyn*0.000035969431019;


text:
function [LOSS, CSP, LDA] = MI_calibration_2(file,band,fs,interval,varargin)

opt=opt_cellToStruct(varargin{:});
cls=opt.nClass;
switch cls
    case 2
        marker={'1','right';'2','left'};
    case 3
        marker={'1','right';'2','left';'3','foot'};
    case 4
        marker={'1','right';'2','left';'3','foot';'4','rest'};
end
field={'x','t','fs','y_dec','y_logic','y_class','class', 'chan'};
[EEG.data, EEG.marker, EEG.info]=Load_EEG(file,{'device','brainVision';'marker',marker;'fs',fs});
CNT_=opt_eegStruct({EEG.data, EEG.marker, EEG.info}, field);

if ~isfield(opt,'channel'),opt.channel=1:EEG.data.nCh;end
if ~isfield(opt,'erd'),opt.erd=0;end

% cross validation parameters
CV.var.band=band;
CV.var.interval=interval;
CV.var.fv='logvar';        % temp
CV.var.classfier='LDA';    % temp
CV.var.evaluation='KFold'; % temp
CV.var.k=5;                % temp
CV.prep={
    'CNT=prep_filter(CNT, {"frequency", band})'
    'SMT=prep_segmentation(CNT, {"interval", interval})'
    };
CV.train={
    '[SMT, CSP_W, CSP_D]=func_csp(SMT,{"nPatterns", [3]})'
    'FT=func_featureExtraction(SMT, {"feature",fv})'
    '[CF_PARAM]=func_train(FT,{"classifier",classfier})'
    };
CV.test={
    'SMT=func_projection(SMT, CSP_W)'
    'FT=func_featureExtraction(SMT, {"feature",fv})'
    '[cf_out]=func_predict(FT, CF_PARAM)'
    };
CV.option={
    'evaluation' , 'k'
    };

n=size(CNT_.class,1);
pairIndex=[];
for i=1:(n-1)
    pairIndex=[pairIndex,[1:i;(1:i)+(n-i)]];
end
[~,I]=sort(sum(pairIndex));
pairIndex=pairIndex(:,I);
pairClass=cell(size(pairIndex));
for i=1:length(pairClass(:))
    pairClass{i}=CNT_.class(pairIndex(i),2);
end
pairNum=size(pairIndex,2);
LOSS=cell(pairNum,2);
CSP=cell(pairNum,2);
LDA=cell(pairNum,2);

for pair=1:pairNum

    pairstr=sprintf('%s vs %s',char(pairClass{1,pair}),char(pairClass{2,pair}));
    disp(pairstr);

    CNT=prep_selectClass(CNT_,{'class',{pairClass{1,pair}, pairClass{2,pair}}});
    CNT=prep_selectChannels(CNT,{'Index',opt.channel});

    %% Pre-processing module
    CNT=prep_filter(CNT, {'frequency', band});
    SMT=prep_segmentation(CNT, {'interval', interval});

    %% Spatial-frequency optimization module
    [SMT, CSP_W, CSP_D]=func_csp(SMT,{'nPatterns', [3]});
    FT=func_featureExtraction(SMT, {'feature',CV.var.fv});

    %% Classifier module
    [CF_PARAM]=func_train(FT,{'classifier',CV.var.classfier});

    %% Evaluation
    [loss]=eval_crossValidation_(CNT, CV);

    CSP{pair,1}=CSP_W; CSP{pair,2}=pairstr;
    LDA{pair,1}=CF_PARAM; LDA{pair,2}=pairstr;
    LOSS{pair,1}=loss;LOSS{pair,2}=pairstr;

end

author: PatternRecognition

repo: OpenBMI


text:
function [volume, flipvecout] = volumeflip(volume, flipvecin)

% VOLUMEFLIP
%
% See also VOLUMEPERMUTE, ALIGN_IJK2XYZ, ALIGN_XYZ2IJK

if nargin<2
  flipvecin = 'auto';
end

% do a low-level check on the input data
if ~isfield(volume, 'transform'), ft_error('the input volume needs a transformation matrix'); end
if ~isfield(volume, 'dim'),       ft_error('the input volume needs a dim field');             end

isrighthanded = det(volume.transform(1:3,1:3))>0;

% check the requested behavior
if ischar(flipvecin)
  switch flipvecin
    case {'left' 'lefthanded'}
      if isrighthanded
        ft_warning('left-handed axes system is requested, performing flip');
        flipvecin = [1 0 0];
      else
        flipvecin = [0 0 0];
      end
    case {'right' 'righthanded'}
      if ~isrighthanded
        ft_warning('right-handed axes system is requested, performing flip');
        flipvecin = [1 0 0];
      else
        flipvecin = [0 0 0];
      end
    case 'auto'
    otherwise
      ft_error('unsupported option specified');
  end
end

% get the largest values on the main diagonal of the transformation matrix
[volume, P] = volumepermute(volume);

% define some variable locally
T   = volume.transform;
dim = volume.dim;

% determine which fields can be flipped
fn = parameterselection('all', volume);

% pre-allocate
flipvecout = false(1,3);

% do the flipping on the numeric data and transformation matrix
for m = 1:3
  if ischar(flipvecin)
    flipvecout(m) = T(m,m)<0;
  else
    flipvecout(m) = flipvecin(m);
  end

  if flipvecout(m)
    % get the reflection matrix
    flipmat = eye(4); flipmat(m,m) = -1; flipmat(m,4) = dim(m)+1;
    for k = 1:numel(fn)
      tmp = volume.(fn{k});
      volume.(fn{k}) = flipdim(tmp, m);
    end
    T = T*flipmat;
  else
    % do nothing
  end
end

% update the transformation matrix
volume.transform = T;

% permute back
volume = volumepermute(volume, P);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%SF_LAGMP Constant lagrange multiplier.
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_LAGMP( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates constant lagrange multiplier.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1-3            Number of space dimensions
%       n_vert      scalar: 2-8            Number of vertices per cell
%       i_dof       scalar: 1-n_ldof       Local basis function to evaluate
%       xi          [n_sdim(+1)]           Local coordinates of evaluation point
%       aInvJac     [n,n_sdim(+1)*n_sdim]  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [n_sdim,n_ldof]        Local coordinates of local dofs
%       sfun        string                 Function name of called shape function

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function h = prtUtilCellPlot(xCell,yCell)
%h = prtUtilCellPlot(xCell,yCell)
% Plot data in cell arrays xCell and yCell

xCell = xCell(:);
yCell = yCell(:);
fullCell = cat(1,xCell',yCell');
fullCell = fullCell(:)';
h = plot(fullCell{:});

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function I = percentile2i(h,P)
%PERCENTILE2I Computes an intensity value given a percentile.
%   I = PERCENTILE2I(H,P) Given a percentile, P, and a histogram, H,
%   this function computes an intensity, I, representing the Pth
%   percentile and returns the value in I. P must be in the range [0, 1]
%   and I is returned as a value in the range [0, 1] also.
%
%   Example:
%
%   Suppose that h is a uniform histogram of an 8-bit image.  Typing
%
%       I = percentile2i(h,0.5)
%
%   would output I = 0.5. To convert to the (integer) 8-bit range
%   [0, 255], we let I = floor(255*I).
%
%   See also I2PERCENTILE.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Check value of P.
if P < 0 || P > 1
   error('The percentile must be in the range [0, 1].')
end

% Normalized the histogram to unit area. If it is already normalized
% the following computation has no effect.
h = h/sum(h);

% Cumulative distribution.
C = cumsum(h);

% Calculations.
idx = find(C >= P,1,'first');
% Subtract 1 from idx because indexing starts at 1, but intensities
% start at 0. Also, normalize to the range [0, 1].
I = (idx - 1)/(numel(h) - 1);

author: dipum

repo: dipum-toolbox

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


text:
filename = 'CantileverSquareSym';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
% m1 = 0.0101;
% m2 = 0.0101;
m1 = 0.8;
m2 = 0.8;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
%vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 50;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 500;

line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 2;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function test_bug2375

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_prepare_headmodel ft_headmodel_localspheres

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2375/localspheres_bug.mat'));

vol = ft_prepare_headmodel(cfg, headshape);

% they should all have the same sphere, since the cfg is faulty
assert(all(vol.r==vol.r(1)));

% fix the cfg and try again
cfg = rmfield(cfg, 'radius');
cfg = rmfield(cfg, 'maxradius');
vol = ft_prepare_headmodel(cfg, headshape);

% they should not all have the same sphere any more
assert(~all(vol.r==vol.r(1)));


text:
function [r,varargout]=circfit(x,y)
%CIRCFIT  Least squares fit of X-Y data to a circle.
%   R = CIRCFIT(X,Y) returns scalar radius R of a fitted circle. X and Y are 1-D
%   arrays of position data in a rectilinear coordinate system. X and Y must be
%   the same length and must contain at least three non-colinear points in order
%   for a valid solution to be found.
%
%   [R,ERR] = CIRCFIT(X,Y) additionally returns the scalar root mean squared
%   error of the fitted circle radius and center relative to the position data.
%
%   [R,XC,YC] = CIRCFIT(X,Y) additionally returns the scalar positions, XC and
%   YC, of center of the fitted circle.
%
%   [R,XC,YC,ERR] = CIRCFIT(X,Y) returns both the center positions of the circle
%   as well as the root mean squared error.
%
%   Examples:
%       % Fit of just five noisy points
%       x1=[1 0 -1 0 1]+0.05*randn(1,5); y1=[0 1 0 -1 0]+0.05*randn(1,5);
%       r1=circfit(x1,y1)
%
%       % CIRCFIT can sometimes perfom poorly with less than 180-degree arc
%       t=0:0.1:pi; lt=length(t);
%       x2=cos(t)+0.04*randn(1,lt); y2=sin(t)+0.04*randn(1,lt);
%       r2_90deg=circfit(x2(1:floor(lt/2)),y2(1:floor(lt/2)))
%       r2_180deg=circfit(x2,y2)

%   Andrew D. Horchler, adh9 @ case . edu, Created 5-12-7
%   Revision: 1.3, 4-6-16

% Check inputs
if nargout > 4
    error('circfit:circfit:TooManyOutputs','Too many output arguments.');
end

if ~isvector(x) || ~isfloat(x) || ~isreal(x) || ~all(isfinite(x))
    error('circfit:circfit:NonFiniteRealVectorX',...
          'X must be a finite real vector of floating point numbers.');
end
if ~isvector(y) || ~isfloat(y) || ~isreal(y) || ~all(isfinite(y))
    error('circfit:circfit:NonFiniteRealVectory',...
          'Y must be a finite real vector of floating point numbers.');
end

lx=length(x);
if lx ~= length(y)
    error('circfit:circfit:LengthMismatch',...
          'The vectors X and Y must have the same length.');
end
if lx < 3
    error('circfit:circfit:Min3Points',...
          'The vectors X and Y must contain at least three points.');
end
x=x(:);
y=y(:);

% Check collinearity, assume with sufficient points, some will be non-collinear
if rank(diff([x([1:min(50,lx) 1]) y([1:min(50,lx) 1])])) == 1
    if lx <= 50 || rank(diff([x y;x(1) y(1)])) == 1
        error('circfit:circfit:Colinearity',...
             ['The points in vectors X and Y must not all be collinear, or '...
              'nearly collinear, with each other.']);
    end
end

xx=x.*x;
yy=y.*y;
xy=x.*y;
xxyy=xx+yy;
sx=sum(x);
sy=sum(y);
sxx=sum(xx);
syy=sum(yy);
sxy=sum(xy);

% Solve linear system without inverting
% a=[sx sy lx;sxy syy sy;sxx sxy sx]\[sxx+syy;sum(xxyy.*y);sum(xxyy.*x)];
[L,U]=lu([sx sy lx;sxy syy sy;sxx sxy sx]);
a=U\(L\[sxx+syy;sum(xxyy.*y);sum(xxyy.*x)]);

xc=0.5*a(1);          	% X-position of center of fitted circle
yc=0.5*a(2);          	% Y-position of center of fitted circle
r=sqrt(xc^2+yc^2+a(3));	% Radius of fitted circle

% Set variable outputs
if nargout > 2
    varargout{1}=xc;
    varargout{2}=yc;
end
if nargout == 2 || nargout == 4
    varargout{nargout-1}=sqrt(mean((sqrt((x-xc).^2+(y-yc).^2)-r).^2));	% RMSE
end

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function [w] = projectSimplex(v)
% Computest the minimum L2-distance projection of vector v onto the probability simplex
nVars = length(v);
mu = sort(v,'descend');
sm = 0;
for j = 1:nVars
    sm = sm+mu(j);
   if mu(j) - (1/j)*(sm-1) > 0
       row = j;
       sm_row = sm;
   end
end
theta = (1/row)*(sm_row-1);
w = max(v-theta,0);

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
clear all; close all;

dt = 0.02;
sim_t = 20;
x0 = [0; 20; 100];

%% Parameters are from
% Aaron Ames et al. Control Barrier Function based Quadratic Programs
% with Application to Adaptive Cruise Control, CDC 2014, Table 1.

params.v0 = 14;
params.vd = 24;
params.m  = 1650;
params.g = 9.81;
params.f0 = 0.1;
params.f1 = 5;
params.f2 = 0.25;
params.ca = 0.3;
params.cd = 0.3;
params.T = 1.8;

params.u_max = params.ca * params.m * params.g;
params.u_min  = -params.cd * params.m * params.g;

params.clf.rate = 5;
params.cbf.rate = 5;

params.weight.input = 2/params.m^2;
params.weight.slack = 2e-2;

%%
accSys = ACC(params);

odeFun = @accSys.dynamics;
controller = @accSys.ctrlCbfClfQp;
odeSolver = @ode45;

total_k = ceil(sim_t / dt);
x = x0;
t = 0;
% initialize traces.
xs = zeros(total_k, 3);
ts = zeros(total_k, 1);
us = zeros(total_k-1, 1);
slacks = zeros(total_k-1, 1);
hs = zeros(total_k-1, 1);
Vs = zeros(total_k-1, 1);
xs(1, :) = x0';
ts(1) = t;
for k = 1:total_k-1
    t
    Fr = accSys.getFr(x);
    % Determine control input.
    [u, slack, h, V] = controller(x, Fr);
    us(k, :) = u';
    slacks(k, :) = slack;
    hs(k) = h;
    Vs(k) = V;

    % Run one time step propagation.
    [ts_temp, xs_temp] = odeSolver(@(t, s) odeFun(t, s, u), [t t+dt], x);
    x = xs_temp(end, :)';

    xs(k+1, :) = x';
    ts(k+1) = ts_temp(end);
    t = t + dt;
end

plot_results(ts, xs, us, slacks, hs, Vs, params)

function plot_results(ts, xs, us, slacks, hs, Vs, params)
    fig_sz = [10 15];
    plot_pos = [0 0 10 15];
    yellow = [0.998, 0.875, 0.529];
    blue = [0.106, 0.588, 0.953];
    navy = [0.063, 0.075, 0.227];
    magenta = [0.937, 0.004, 0.584];
    orange = [0.965, 0.529, 0.255];
    grey = 0.01 *[19.6, 18.8, 19.2];

    figure(1);
    subplot(6,1,1);
    p = plot(ts, xs(:, 2));
    p.Color = blue;
    p.LineWidth = 1.5;
    hold on;
    plot(ts, params.vd*ones(size(ts, 1), 1), 'k--');
    ylabel("v (m/s)");
    title("State - Velocity");
    set(gca,'FontSize',14);
    grid on;

    subplot(6,1,2);
    p = plot(ts, xs(:, 3));
    p.Color = magenta;
    p.LineWidth = 1.5;
    ylabel("z (m)");
    title("State - Distance to lead vehicle");
    set(gca, 'FontSize', 14);
    grid on;

    subplot(6,1,3);
    p = plot(ts(1:end-1), us); hold on;
    p.Color = orange;
    p.LineWidth = 1.5;
    plot(ts(1:end-1), params.u_max*ones(size(ts, 1)-1, 1), 'k--');
    plot(ts(1:end-1), params.u_min*ones(size(ts, 1)-1, 1), 'k--');
    ylabel("u(N)");
    title("Control Input - Wheel Force");
    set(gca, 'FontSize', 14);
    grid on;

    subplot(6,1,4);
    p = plot(ts(1:end-1), slacks); hold on;
    p.Color = magenta;
    p.LineWidth = 1.5;
    ylabel("slack");
    title("Slack variable");
    set(gca, 'FontSize', 14);
    grid on;

    subplot(6,1,5);
    p = plot(ts(1:end-1), hs);
    p.Color = navy;
    p.LineWidth = 1.5;
    ylabel("CBF (h(x))");
    title("CBF");
    set(gca, 'FontSize', 14);
        grid on;

    subplot(6,1,6);
    set(gca, 'FontSize', 14);
    p = plot(ts(1:end-1), Vs);
    p.Color = navy;
    p.LineWidth = 1.5;
    xlabel("t(s)");
    ylabel("CLF (V(x))");
    title("CLF");
    set(gca, 'FontSize', 14);
    grid on;

    set(gcf, 'PaperSize', fig_sz);
    set(gcf, 'PaperPosition', plot_pos);
end

author: HybridRobotics

repo: CBF-CLF-Helper


text:
function scale = p30_scale ( neqn )

%*****************************************************************************80
%
%% P30_SCALE returns scale factors for problem p30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
  function xrs = xray_read_spectra(stype, varargin)
%|function xrs = xray_read_spectra(stype, [options])
%|
%| Read X-ray spectra data and initialize a structure that describes "M"
%| piecewise constant polyenergetic X-ray spectra, where M=2 for dual-kVp case.
%|
%| in
%|	stype		char	which spectrum model:
%|				'mono,60,100' dual mono-energetic
%|				'poly1,kvp1[,kvp2][,...]' polyenergetic
%|		or:	cell	{en [ne 1], sp [ne M]} sampled spectra
%|
%| option
%|	'en'		[ne 1]	specify energy sampling
%|	'filters'	cell{M} optional additional filtration of each spectrum
%|		{{mtype11, thick11, mtype12, thick12, ...}, ...
%|		 {mtypeM1, thickM1, mtypeM2, thickM2, ...}}
%|	'units'		''	'cm' (default) or 'mm' for filter thicknesses
%|	'show'		1|0	plot?
%|	'name'		{char}	names for each spectrum
%|	'kvp'		[M 1]
%|	'dir_spectra'	char	directory with spectra.  default:
%|				[path_find_dir('ct') filesep 'xray-spectra'];
%| out
%|	xrs		strum
%| data:
%|	xrs.en		[ne 1]	energy list
%|	xrs.sp		[ne M]	M spectra, for M kVp settings: I_m(E)
%|	xrs.Ide		[ne M]	differential array for "integrating": I_m(E) dE
%|	xrs.I		[M 1]	total intensity: sum(Ide)
%|	xrs.kvp		[M 1]	tube potentials (voltage in kVp)
%|	xrs.eff		[M 1]	mean energies (incident)
%| methods:
%|		.plot	plot spectra
%|
%| Note: xrs.eff perhaps should be named xrs.mean because it is a mean energy
%| of the incident spectrum.  It is not the "effective energy" as defined
%| e.g., on p56 of 4th edition of W. R. Hendee's "Medical Imaging Physics"
%| because that definition depends on the attenuation of a particular material.
%|
%| Copyright 2001-04-27, Jeff Fessler, University of Michigan

if ~nargin, ir_usage, end
if streq(stype, 'test'), xray_read_spectra_test, return, end

% defaults
arg.en = [];
arg.filters = {};
arg.show = false;
arg.name = {};
arg.kvp = [];
arg.units = 'cm';
arg.dir_spectra = [path_find_dir('ct') filesep 'xray-spectra'];
arg = vararg_pair(arg, varargin);

xrs.kvp = arg.kvp;

if ischar(stype)
	xrs = xray_read_spectra_char(stype, arg.en, arg.filters, arg.dir_spectra);
	xrs.type = stype;

elseif iscell(stype) && numel(stype) == 2
	xrs.en = stype{1};
	xrs.sp = stype{2};
	xrs.type = '';
	xrs.filters = arg.filters;
	if ~isequal(size(xrs.en,1), size(xrs.sp,1))
		error 'bad stype as cell usage'
	end
else
	error 'bad stype'
end

MM = size(xrs.sp, 2);

%
% apply filtration, if any
%
if ~isempty(xrs.filters)
	if numel(xrs.filters) ~= 1 && numel(xrs.filters) ~= MM
		error 'should be 1 or M sets of filters'
	end
	for mm=1:MM
		xrs.sp(:,mm) = xray_apply_filters(xrs.sp(:,mm), xrs.en, ...
			xrs.filters{min(mm,end)}, arg.units);
	end
end

%
% precompute some aspects of the spectra
%
xrs.Ide = xrs.sp .* repmat(difff(xrs.en), 1, MM); % [N M] I(E) dE
xrs.I = sum(xrs.Ide); % [1 M] spectrum integral
xrs.eff = xrs.en' * xrs.Ide ./ xrs.I;
for mm=1:MM
	xrs.sp_at_eff(mm) = interp1(xrs.en, xrs.sp(:,mm), xrs.eff(mm));
end

if isempty(arg.name)
	if ischar(stype)
		xrs.name = cat(2, {stype}, cell(1,MM-1));
	else
		xrs.name = '';
	end
else
	xrs.name = arg.name;
end
xrs = strum(xrs, {'plot', @xray_plot_spectra});

if arg.show
	xrs.plot;
end

%
% xray_read_spectra_char()
%
function xrs = xray_read_spectra_char(stype, en, filters, dir_spectra)
xrs.en = en;
xrs.filters = filters;

%
% simplest option is mono-energetic (for testing)
% usage: mono,kvp1,kvp2,...
%
if streq(stype, 'mono', 4)
	xrs.kvp = str2num(strrep(stype(6:end), ',', ' '));
	for mm=1:numel(xrs.kvp)
		xrs.en = [20:140]';
		xrs.sp(:,mm) = xrs.en == xrs.kvp(mm);
	end

%
% polyenergetic spectra
% usage: poly1,kvp1,kvp2,...
%
elseif streq(stype, 'poly1', 5)
	if ~exist(dir_spectra, 'dir')
		fail('spectra dir "%s" not in path', dir_spectra)
	end
	xrs.kvp = str2num(strrep(stype(7:end), ',', ' ')); % [M]

	% read raw data
	MM = numel(xrs.kvp);
	for mm=1:MM
		kvp = xrs.kvp(mm);
		raw = sprintf('spectra.%d', kvp);
		raw = [dir_spectra filesep raw];
		tmp = load_ascii_skip_header(raw);
		sr.enc{mm} = tmp(:,1);
		sr.spc{mm} = tmp(:,2);

		% The Wilderman/Sukovic spectra must be scaled by energy!
		sr.spc{mm} = sr.spc{mm} .* sr.enc{mm};
	end

	% interpolate onto same energy sampling
	[xrs.en xrs.sp] = xray_sp_interp(xrs.en, sr.enc, sr.spc);

%
% 1st spectra from predrag sukovic
%
elseif streq(stype, 'ps1')
	if ~exist(dir_spectra, 'dir')
		fail('spectra dir "%s" not in path', dir_spectra)
	end

	xrs.kvp = [80 140];

	for mm=1:numel(xrs.kvp)
		file = sprintf('xray%03d.mat', xrs.kvp(mm));
		file = [dir_spectra filesep file];
		if ~exist(file, 'file')
			error(sprintf('file "%s" not found', file))
		end
		raw = load(file);
		ie = raw.energy >= 20 & raw.energy <= 140;
		sr.enc{mm} = raw.energy(ie);
		sr.spc{mm} = raw.spe(ie) .* raw.energy(ie);
	end

	[xrs.en xrs.sp] = xray_sp_interp(xrs.en, sr.enc, sr.spc);

%
% spectra used for 2002 SPIE talk (fix: or did i just use 'ps1' ???)
%
elseif streq(stype, 'spie02')
	xrs = xray_read_spectra('poly1,80,140');
%	xrs.filters = { {{'aluminum', 0.25}, {'copper', 0.05}} };
	xrs.filters = {{'aluminum', 0.25, 'copper', 0.05}};
%	xrs.filters = {{'aluminum', 0.10}, {'copper', 0.05}};
%	xrs.filters = {{'aluminum', 0.5, 'copper', 0.04, 'csi', 0.1}};
%	xrs.filters = {{'aluminum', 0.5, 'copper', 0.04, 'gadolinium', 0.1}};

else
	fail('bad stype "%s"', stype)
end

%
% interpolate onto same energy sampling
%
function [en, sp] = xray_sp_interp(en, enc, spc)
MM = numel(spc);
if isempty(en)
	tmp = zeros(MM,1);
	for mm=1:MM
		tmp(mm) = max(enc{mm});
	end
	mm = imax(tmp); % find col with largest max
	en = enc{mm};
end
sp = zeros(numel(en), MM);
for mm=1:MM
	sp(:,mm) = interp1(enc{mm}, spc{mm}, ...
		en, 'linear', 0); % extrapolate with zeros
end

%
% xray_apply_filters()
%
function sp = xray_apply_filters(sp, en, filters, units);
% {mtype1, thick1, mtype2, thick2, ...}, ...
if ~iscell(filters)
	error 'filtration arguments must be cells'
end
for ii=1:2:numel(filters)
	sp = sp .* xray_filters(filters{ii}, filters{ii+1}, en, 'units', units);
end

%
% xray_plot_spectra()
% plot routine
%
function dummy = xray_plot_spectra(xrs, varargin)
dummy = [];

arg.title = xrs.type;
arg.kmin = min(xrs.en);
arg.kmax = max(xrs.en);
arg.spike = false;
arg.threshold = 1.2;
arg.hold = false;
arg.linetype = 'y.-';
arg = vararg_pair(arg, varargin);

MM = size(xrs.sp,2);
if ~im, return, end

if ~arg.hold, clf, end
pl=@(m) subplot(MM*100+10+m);
for mm=1:MM
	pl(mm)
	if arg.hold, hold on, end
	plot(xrs.en, xrs.sp(:,mm), arg.linetype, ...
		xrs.eff(mm) * [1 1], [0 xrs.sp_at_eff(mm)], 'm--')
	if arg.hold
		hold off
	else
		axis([arg.kmin arg.kmax [-0.00 1.05]*max(xrs.sp(:,mm))])
	end
	xt = [arg.kmin round(xrs.eff(mm)) arg.kmax];
	if arg.spike
		where = xrs_find_spike(xrs.sp(:,mm), arg.threshold);
		xt = [xt col(xrs.en(where))'];
	end
	xtick(sort(xt)), ytick([0])
	ylabelf('I%d(E)', mm)
	if mm == 1
		title(arg.title)
	end
	if ~isempty(xrs.name) && ~isempty(xrs.name{mm})
		title(xrs.name{mm})
	end

	if ~isempty(xrs.kvp)
		t = sprintf('%g kVp Spectrum', xrs.kvp(mm));
		texts(0.5, 0.8, t, 'color', 'green')
	end
end
xlabel 'Energy [keV]'

%
% xrs_find_spike()
% find spikes in spectrum
%
function where = xrs_find_spike(sp, threshold)
local = (sp(1:end-2) + sp(3:end)) / 2;
where = (sp(2:end-1) - local) > threshold * local;
where = 1 + find(where);

%
% xray_read_spectra_test()
% test routine
%
function xray_read_spectra_test

stype = 'mono,60,100';
stype = 'poly1,80,100';
stype = 'spie02';
stype = 'ps1';
stype = 'poly1,160,100,60';
filts = {{}, {}, {}};
xrs1 = xray_read_spectra(stype, 'filters', filts);
filts = {{'aluminum', 0.1}, {'copper', 0.2}, {}};
xrs2 = xray_read_spectra(stype, 'filters', filts);
xrs1.plot;
xrs2.plot('hold', true, 'linetype', 'c-');
%clf, plot(xrs1.en, xrs1.sp, '-', xrs2.en, xrs2.sp, '--')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function c = tapas_ehgf_binary_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the enhanced Hierarchical Gaussian Filter (eHGF)
% for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_ehgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_ehgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_ehgf_binary_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'ehgf_binary';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -3,   2];
c.omsa = [NaN,   4,   4];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 3*c.n_levels+2*(c.n_levels-1)+1;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_ehgf_binary;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_ehgf_binary_transp;

end

author: translationalneuromodeling

repo: tapas


text:
function positions = getFigurePositions(tile, screenBorder)
    numRows   = tile(1);
    numCols   = tile(2);
    figBorder = [20, 60];

    screenSize = get(0,'screensize'); % Determine terminal size in pixels

    % Adjust to actual visible screen size
    screenSize(3:4) = screenSize(3:4) - (nargin == 2) * 2 * screenBorder;

    % Set figure width and height with border
    figWidth  = fix(screenSize(3)/numCols);
    figHeight = fix(screenSize(4)/numRows);

    % Calculate horizontal and vertical position for each figure
    horPos = screenSize(3) - (numCols:-1:1)*figWidth  + figBorder(1);
    verPos = screenSize(4) - (1:numRows)*figHeight    + figBorder(2);
    horPosAll = repmat(horPos', [numRows, 1]);
    verPosAll = repmat(verPos,  [numCols, 1]);

    % Return position
    positions = zeros(prod(tile), 4);
    positions(:,1) = horPosAll;
    positions(:,2) = verPosAll(:);
    positions(:,3) = figWidth  - 2 * figBorder(1);
    positions(:,4) = figHeight - 2 * figBorder(2);
end

author: analogdevicesinc


text:
%%%%  Demodulator

function [data_out]=demod_d(data_in,rate_id)
% define parameter
switch (rate_id)
    case 0
       M=2;  %Size of signal constellation for BFSK
    case {1,2}
       M=4;  %Size of signal constellation QPSK
    case {3,4}
       M=16;  %Size of signal constellation 16QAM
    case {5,6}
       M=64;  %Size of signal constellation 64QAM
    otherwise
       display('error in constellation modulator give proper rate_id')
end

k=log2(M); % no. of bits per symbol
%% Received Signal
yrx = data_in;
% %%
% scatterplot(yrx)

% Demodulate signal % result in column vector containing the value 0 to M-1
switch (rate_id)
    case {0,1,2}
      zsym = pskdemod(yrx,M);
    case {3,4,5,6}
       zsym = qamdemod(yrx,M);
    otherwise
       display('error in demodulation give proper rate_id')
end

% Symbol-to-Bit Mapping
z = de2bi(zsym,'left-msb');
data_out = reshape(z.',prod(size(z)),1);


text:
% DEMO  --  Chebyshev Polynomial Interpolation Order vs Accuracy
% UPDATED  --  October 28, 2013
% Written by Matthew Kelly, Cornell University
%
%
% This script demonstrates how the accuracy of a chebyshev approximation is
% dependant on the order of the underlying polynomial.
%
% RESULTS:
%   The accuracy is very bad if the order is too small, but then rapidly
%   improves until it reaches machine precision. The accuracy of the
%   derivative will get slowly worse for very large order apprimations due
%   to the compounding of rounding errors in the matrix multiplication.
%

%% General Settings
clear; clc;

%What order should the approximation be?
orderRange = [2,1000];
orderNum = 20;
orderList = unique(round(logspace(...
    log10(orderRange(1)),...
    log10(orderRange(2)),...
    orderNum)));

%How many points should be used for error calculations and plotting?
nTime = 1000;

%What domain should we be looking at?
d = [-2,1];

%Time for use in plots
time = linspace(d(1),d(2),nTime);

%Set up the input/output struct:
IO.domain = d;
IO.userFunc = @testFunction;

%Set up the data logging:
N = length(orderList);
cpuTime = zeros(N,1);
maxError = zeros(N,2);
meanError = zeros(N,2);

%Run the calculations:
for i=1:N
    disp(['Order: ' num2str(orderList(i))])
    tic;
    f = chebyshevFit(IO, orderList(i));
    [y, Dy] = chebyshevInterpolate(f,time,d);
    cpuTime(i) = toc;
    [g,Dg] = testFunction(time);
    e = abs(g-y);
    De = abs(Dg-Dy);
    maxError(i,:) = [max(e),max(De)];
    meanError(i,:) = [mean(e),mean(De)];
end

%Make nice plots:
    figure(403); clf;
    subplot(3,2,1)
        loglog(orderList,cpuTime)
        title('CPU Time')
        xlabel('Order of Chebyshev Polynomial')
        ylabel('Time (s)')
    subplot(3,2,2)
        plot(time,testFunction(time))
        title('Function')
        xlabel('Input')
        ylabel('Output')
    subplot(3,2,3)
        loglog(orderList,maxError(:,1))
        title('Max Error in Function')
        xlabel('Order of Chebyshev Polynomial')
    subplot(3,2,5)
        loglog(orderList,meanError(:,1))
        title('Mean Error in Function')
        xlabel('Order of Chebyshev Polynomial')
    subplot(3,2,4)
        loglog(orderList,maxError(:,2))
        title('Max Error in Derivative')
        xlabel('Order of Chebyshev Polynomial')
    subplot(3,2,6)
        loglog(orderList,meanError(:,2))
        title('Mean Error in Derivative')
        xlabel('Order of Chebyshev Polynomial')

author: MatthewPeterKelly

repo: dscTutorials


text:
function eclipVec=ecliptic2ICRS(vec,TT1,TT2,method)
%%ECLIPTIC2ICRS Convert a location vector from the ecliptic coordinate
%               system to the International Celestial Reference System
%               (ICRS) either using the IAU 2006 precession model or the
%               Vondrak 400 millennia precession model.
%
%INPUTS: x The NXnumVec collection of vectors in the ecliptic coordinate
%          system to convert. N can be 2, or 3. If the vectors are 2D,
%          then they are assumed to be azimuth and elevation in radians.
%          3D vectors are assumed to be Cartesian position.
% Jul1, Jul2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of
%          both terms. The date is broken into two parts to provide more
%          bits of precision. It does not matter how the date is split.
%   method An optional parameter specifying which algorithm is to be used.
%          Possible values are
%          0 (The default if omitted or an empty matrix is passed) Use
%            the IAU 2006 precession model.
%          1 Use the long-term (Vondrak) precession model.
%
%OUTPUTS: eclipVec The vectors rotated into the ICRS. If the input was 2D
%                  azimuth and elevation, the output will be the same. If
%                  the input is Cartesian, then the output will be
%                  Cartesian.
%
%This function is a Matlab interface for the relevant functions in the
%International Astronomical Union's (IAU) Standard's of Fundamental
%Astronomy library.
%
%The ecliptic is defined in the IERS Conventions [1] to be the "the
%plane perpendicular to the mean heliocentric orbital angular momentum
%vector of the Earth-Moon barycentre in the BCRS".
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%eclipVec=ecliptic2ICRS(vec,TT1,TT2,method);
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%
%July 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function G = Grid2GlinearROI(Grid,Grid_dec,npxdtc,ctrIdx,block,interp,nblk,vmsk,nnzmsk)
% Grid2Glinear
% input Grid/s generated in fbkp_geometry
% -Grid is a matrix of integers that are coordinate of the detector pixel
% -Grid_dec is [] is the 'nearest neighbour' interp option is used
% -block (integer & scalar) sets the projection undergoing backprojection
% -nblk number of blocks (i.e. projections) used to make up the (rows) sparse matrix
% i.e. the the large block goes from block (i.e. projection) to block + nblk
% -vmsk is the mask as a vector (logical)
% -npxdtc number of detector pixels ( nb in Glinear)
% output one sparse matrix G as by Glinear (see Fessler code)
%
% Gianni Schena March  2003

G=[];
block
nargin

Nx=size(Grid(:,:,1),1); Ny=size(Grid(:,:,1),2); % get the size of the problem
nm=Nx*Ny
%

if nargin < 8
	vmsk=logical(ones(1,nm)); nnzmsk=nm;
	% define default mask - if the mask is not an input parameter
end

if nargin < 7, nblk=1, end

%SPALLOC(M,N,NZMAX) creates an M-by-N all zero sparse matrix with room to eventually hold NZMAX nonzeros.

if strcmp(interp, 'nearest neighbor')
	G = spalloc(nblk*npxdtc,nnzmsk, nblk*nnzmsk);
else
	G = spalloc(nblk*npxdtc,nnzmsk, 2*nblk*nnzmsk); % pre-allocation
end

jc=[1:1:nnzmsk]; % index of columns in the matrix G

for ib=1:nblk % for nblk projections and starting from block

	Gv=Grid(:,:,block+ib-1) ; Gv=Gv(:) ; % starts from block
	ir = double(Gv) + ctrIdx  ; % the coordinate (saved as 'round') + centre
	ir=[ir+npxdtc*(ib-1)]; ir=ir(vmsk)';

	% S = SPARSE(i,j,s,m,n,nzmax) uses the rows of [i,j,s] to generate an
	% m-by-n sparse matrix with space allocated for nzmax nonzeros.

	if strcmp(interp, 'nearest neighbor')
		v=1.;
		T=sparse(ir,jc,double(v),npxdtc*nblk,nnzmsk);

	else % i.e. interp == linear
		dec=double(Grid_dec(:,:,block+ib-1))/100; % use also the decimal saved as integer
		vc = dec(:); % weight for (floor +1) == weight of 'ceil' term
		vc=vc(vmsk);
		T= [sparse(ir,jc,(1-vc),npxdtc*(nblk),nnzmsk) + sparse(ir+1,jc,vc,npxdtc*(nblk),nnzmsk)];
		%T= [sparse(ir,jc',(1-vc)) + sparse([ir+1],jc',vc)];
	end

	G=[G+T];
end % end ib


text:
%--------------------------------------------------------------------------
%----------------------  Validation of the result -------------------------
%--------------------------------------------------------------------------

% With this function you can re-estimated the CCR for the optimal Feature
% Subset found by a Feature Selection method

 function [ResultMat, UpLimitCOD,  LowLimitCOD] = ...
                    CCRForOptSet(ResultMat, DatasetToUse, ErrorEstMethod)
 format short g
 axis([0 length(ResultMat(:,1)) 0 1])
 axis manual
%----------- Feed from Feature Selection method----------------------------
%Cross-val
%PoolSelectFeat = [18;19;45;14;28;40;6;90;11;68;84;22;69;2;44;70;25;7;65;1;5;4;8;57;9;21;13;10;29;3;53;20;27;15;47;23;16;33;26;30;35;89;24;82;31;59;62;79;32;50;34;36;51;55;60;37;72;43;12;38;17;58;39;41;86;67;42;46;48;49;52;54;56;64;88;61;66;63;75;71;85;73;74;76;80;77;78;81;87;83;];
%Resub
PoolSelectFeat =ResultMat(:,3);
%--------------------------------------------------------------------------

 if strcmp(DatasetToUse, 'finalvecDESsubset')
     [Patterns, Targets] = Preprocessing('finalvecDES');
     Patterns = Patterns(1:360,:);
     Targets  = Targets(1:360);
 elseif strcmp(DatasetToUse, 'finalvecSUSASsubset')
     [Patterns, Targets] = Preprocessing('finalvecSUSAS');
     rand('twister',260); % To have always the same set
     IndexSubset = randperm(length(Targets));
     IndexSubset = IndexSubset(1:end);
     Patterns = Patterns(IndexSubset,:);
     Targets  = Targets(IndexSubset);
     size(Patterns)
 else
     [Patterns, Targets] = Preprocessing(DatasetToUse);
 end

[NPatterns, KFeatures] = size(Patterns);
CClasses = max(Targets);
%  ------------------------- sffs settings -------------
PercTest             = 10;             % Percentage of data to use for testing. Options 5,...,50.
GammaParam = 0;                        % confidence interval to control number of repetitions.
                                       % Options:0<GammaParam <1. The lowest, the better.

if strcmp(ErrorEstMethod,'ProposedAB')
    NRepThres = [];                   % Only for programming. Do not change.
elseif strcmp(ErrorEstMethod,'ProposedA') || strcmp(ErrorEstMethod,'Standard')
   NRepThres = 1;                  % Number of repetitions for the standard SFFS. Options: >10,
                                               % The greater, the better.
   GammaParam = [];                          % Only for programming. Do not change.
elseif strcmp(ErrorEstMethod,'Resubstitution')
    NRepThres  =1;
end

ResultMat           = [];          		        	     % Result matrix with selection history.

for FeatIndx = 1:length(PoolSelectFeat)
    ResultMat(end+1,1)  = FeatIndx;
    ResultMat(end,2) = PoolSelectFeat(FeatIndx);

    [Critval(FeatIndx)] = EnsembCriterEval(...
                               Patterns(:,PoolSelectFeat(1:FeatIndx)), ...
                                        Targets, PercTest, ErrorEstMethod);

    ResultMat(end,4) = 1;
    ResultMat(end,3) = Critval(FeatIndx);

    NDc = NPatterns/CClasses;

    [InfoLossResub]   = CalcInfoLoss(FeatIndx, floor(NDc), 1);
    [InfoLossCross]   = CalcInfoLoss(FeatIndx, floor(NDc), 0);
    AverInfoLoss      = (InfoLossCross + InfoLossResub)/2;

    CCR                   = ResultMat(end,3);
    UpLimitCOD(FeatIndx)  = CCR  + AverInfoLoss*(1-CCR);
    LowLimitCOD(FeatIndx) = CCR  - AverInfoLoss*(CCR-1/CClasses);

    hold on
    plot(ResultMat(1:end,3))
    plot(   UpLimitCOD, 'g.-')
    plot(   LowLimitCOD, 'r.-')

    drawnow
end

return

%--------------------------------------------------------------------------------------------------------------
%--------------------------------------------------------------------------------------------------------------
%--------------------------------------------------------------------------------------------------------------
% -------------------------------  Criterion evaluation method ---------------------------------------
%--------------------------------------------------------------------------------------------------------------
%--------------------------------------------------------------------------------------------------------------

function [AverCorr] = EnsembCriterEval(Patterns, Targets, PercTest,...
               ErrorEstMethod)

NPatterns = size(Patterns,1);
Kfeatures = size(Patterns,2);
CClasses  = max(Targets);

GaussConst = (2*pi)^(-Kfeatures/2);

CCRTable = zeros(1, 2000);
Prior         = zeros(1,CClasses);
NRep = 500;

rep=0;
while rep  <  NRep,
        rep=rep+1;
        if strcmp(ErrorEstMethod,'ProposedAB') || strcmp(ErrorEstMethod,'ProposedA') ...
                || strcmp(ErrorEstMethod,'Standard'),
           IndexAll  = randperm(NPatterns);
           IndexTrain = IndexAll((floor(NPatterns*PercTest/100)+1):end);
           TrainSetPatterns = Patterns(IndexTrain,:);
           TrainSetTargets = Targets(IndexTrain);
           IndexTest = IndexAll(1:floor(NPatterns*PercTest/100));
           TestSetPatterns   = Patterns(IndexTest,:);
           TestSetTargets = Targets(IndexTest);
        elseif strcmp(ErrorEstMethod,'Resubstitution')
            TrainSetPatterns = Patterns;
            TrainSetTargets = Targets;
            TestSetPatterns = Patterns;
            TestSetTargets = Targets;
        end

           NTestPatterns = size(TestSetPatterns,1);
           NTrainPatterns = size(TrainSetPatterns,1);
           TargetsProbs = zeros(NTestPatterns, CClasses);

            for IndexClass = 1:CClasses
                    TrainSetIndexPerClass = find(TrainSetTargets == IndexClass);
                    TrainSetPatternPerClass = TrainSetPatterns(TrainSetIndexPerClass,: );
                    ResCl.mean(IndexClass,:)         = mean(TrainSetPatternPerClass);
                    ResCl.cov(:,:,IndexClass)         = cov(TrainSetPatternPerClass);
                    ResCl.covDet(IndexClass)        = det(ResCl.cov(:,:,IndexClass));
                    if ResCl.covDet(IndexClass) ~= 0
                       ResCl.covInv(:,:,IndexClass)     = inv(ResCl.cov(:,:,IndexClass));
%                    [ResCl.covInv(:,:,IndexClass), ResCl.covDet(IndexClass)] =...
   %                                                                               InvDet(ResCl.cov(:,:,IndexClass));
%                        NTestPatternsPerClass  =  sum(TestSetTargets == IndexClass);
%                       Prior(IndexClass) = NTestPatternsPerClass /length(TestSetTargets);
                          Prior(IndexClass)  = length(TrainSetIndexPerClass) / NTrainPatterns;
%                              Prior(IndexClass)  = 1/CClasses; % all equiprobable
                          ResCl.prior(IndexClass) = Prior(IndexClass);
                          ClassConst = GaussConst * ResCl.prior(IndexClass)*ResCl.covDet(IndexClass)^(-0.5);
                          DistMeanVectorTest = TestSetPatterns - ...
                                                                   repmat( ResCl.mean(IndexClass,:), NTestPatterns,1);
                          A = DistMeanVectorTest * ResCl.covInv(:,:,IndexClass);
                          G = sum(A.*DistMeanVectorTest,2);
                          TargetsProbs(1:NTestPatterns, IndexClass) = ClassConst * exp(-0.5*G);
                    else
                       TargetsProbs(1:NTestPatterns, IndexClass) = eps*rand(NTestPatterns,1);
                    end
            end

           [TargetMaxProb, TargetPrediction] = max(TargetsProbs');
           TargetPrediction = TargetPrediction';
           CorrClassified = length(find(TargetPrediction - TestSetTargets == 0));
           CCRTable(rep) = CorrClassified/ NTestPatterns;

 end                 % repetetion
AverCorr = sum(CCRTable(1:NRep))/NRep;
return


text:
function out = getMinuteDataFromDB(tableName)

% Copyright 2011, The MathWorks, Inc.
% All rights reserved.

% Set preferences with setdbprefs.
s.DataReturnFormat = 'numeric';
s.ErrorHandling = 'store';
s.NullNumberRead = 'NaN';
s.NullNumberWrite = 'NaN';
s.NullStringRead = 'null';
s.NullStringWrite = 'null';
s.JDBCDataSourceFile = '';
s.UseRegistryForSources = 'yes';
s.TempDirForRegistryOutput = 'C:\Temp';
s.DefaultRowPreFetch = '10000';
setdbprefs(s)

% Make connection to database.  Note that the password has been omitted.
% Using ODBC driver.
conn = database('Minute Bars','','password');

% Read data from database.
e = exec(conn,['SELECT ALL Dates,High,Low,Close FROM ' tableName]);
e = fetch(e);
close(e)

% Assign data to output variable.
out = e.Data;

% Close database connection.
close(conn)


text:
function [nota1, nota2]=check_my_kinematics(robot)

   %q=zeros(1, robot.DOF) + 0.05*rand(1, robot.DOF);
   q=[0.1 0.1 0.1 0.1 0.1 0.1] + 0.05*rand(1, robot.DOF);

   %Now compute direct kinematics for this position q
   T = directkinematic(robot, q);

   %Set to zero if you want to see the robot transparent
   robot.graphical.draw_transparent=1;
   robot.graphical.draw_axes=1;

   qinv = inversekinematic(robot, T);

resultado_q = [];
resultado_T = [];

for i=1:size(qinv,2),
    Ti = directkinematic(robot, qinv(:,i)); %Ti is constant for the different solutions

    resultado_q= [resultado_q sum(qinv(:,i)-q')];
    resultado_T= [resultado_T sum(sum(T-Ti))];

    %now draw the robot to see the solution
    drawrobot3d(robot, qinv(:,i))

    pause(1)
end

indice1 = find(abs(resultado_q) < 0.05);
indice2 = sum(abs(resultado_T));

if size(indice1,2)==0
   nota1=0;
else
   nota1=10;
end

if indice2>0.01
    nota2=0;
else
    nota2=10;
end

repo: ARTE


text:
%% generate eye inf.
% author: Shuang Wang
% email: shw070@ucsd.edu
% Division of Biomedical Informatics, University of California, San Diego.
function a = eyeInf(len)
    a = zeros(len);
    a(1:length(a)+1:numel(a)) = Inf;
end


text:
function [outMat,procInfo,dcmInfo]=ReadDicomDir(inDir,debugLevel)
% [outMat,procInfo,dcmInfo]=ReadDicomDir(inDir [,debugLevel])
%   Read DICOM files from directory 'inDir'.
%   If mag and phase images in same directory: read and sort these images;
%   phase images will be at 2nd half of outMat
%
% inDir: directory where the DICOM files are stored
% debugLevel: if >0, print debug message; debugLevel=0, 1, 2
%
% outMat: a 2D matrix for a single slice, 3D matrix for multiple slices
% procInfo: abbreviated DICOM information
% dcmInfo: DICOM information from 1st DICOM file
%
% Example:
% d1 = 'C:\Data\TAG7MM-92_11'
% [a,pname,pp]=ReadDicomDir(d1);
% figure; imagesc(abs(a(:,:,20))
%
% check current directory
% [a]=ReadDicomDir('.');
%
% use GUI to choose directory
% [a]=ReadDicomDir;
%
% Tested using Siemens dicom images and Philips mag, phase dicom images.
% Tested using Philips mag and phase images
% Tested using GE dicom images: mag, phase, mag, phase, mag, phase ....
%
% For Siemens dicom in mosaic format: default to be magnitude.
% Works for Siemens PET dicom images.
% For Philips real, imag images: check ImageType and revise flags
%
% Wen-Tung Wang, 201106
% Revised to read and sort mag and phase images (currently Siemens images only). 201107
% Revised to read Siemens dicom in mosaic format (all magnitude images).201205
% Revised to read GE dicom images. 201205

% msg = nargchk(1,2,nargin);
% if ~isempty(msg)
%    fprintf(1,'\tusage: [outMat,procInfo,dcmInfo]=ReadDicomDir(inDir [, debugLevel])\n');
%    fprintf(1,'\tRead multiple single-slice DICOM data (*.dcm, not *.ima)\n');
%    fprintf(1,'\tcheck "help ReadDicomDir" for more details\n');
%    return;
% end
% input from ui
    if (nargin<1) || (isempty(inDir)),
        inDir=uigetdir('','Select DICOM Directory');
        if ~ischar(inDir)
            disp('no valid Directory selected.');
            outMat=[]; procInfo=[]; dcmInfo=[];
            return;
        end
    end
% debugLevel
    if nargin<2, debugLevel=0; end
% check if directory exist
    if (debugLevel), tic; end
    if (exist(inDir,'dir'))~=7,
	    fprintf(2,'%s \nnot exist or not a folder\n',inDir);
	    outMat=[]; procInfo=[]; dcmInfo=[];
	    return
    end
	aaDir = strfind(inDir,filesep);
	lenDir = length(inDir);
    if (~isempty(aaDir)),
        if (aaDir(end)==lenDir), % remove last filesep: either / or \
            inDir = inDir(1:end-1);
        end
    end
% read all files
    ffext = sprintf('*');
    fdir = fullfile(inDir,ffext);
    if (debugLevel), fprintf(1,'%s\n',fdir); end
	ff = dir(fdir);
    flen0 = numel(ff);
	flen = length(ff);
% exclude '.' and '..', and other files if needed
    ff(cell2mat({ff(:).isdir})>0)=[];
	flen = length(ff);
    ff0 = ff; %filenames without directory
% add path to filenames
    if (debugLevel),
        fprintf(1,'%d files. 1st file = %s\n',flen,ff(1).name);
    end
	for k=1:flen
		ff(k).name = [inDir,filesep,ff(k).name];
	end
% GE or Siemens or Philips DICOM
% Use magmatch, phamatch, realmatch, imagmatch to distinguish magnitude,
% phase, real, and imaginary dicom images
	display(ff(1).name);
	dinfo=dicominfo(ff(1).name);
    GEflag = 0; Siemensflag = 0; Philipsflag = 0;
    dcmInfo = dinfo;
    if (strfind(dinfo.Manufacturer,'SIEMENS')),
        Siemensflag = 1;
        vendor = 'Siemens';
        magmatch = 'ORIGINAL\PRIMARY\M\ND\NORM';
        magmatch = 'ORIGINAL\PRIMARY\M\';
        magmatch2= 'ORIGINAL\PRIMARY\M\ND\MOSAIC'; %mosaic format: magnitude
        phamatch = 'ORIGINAL\PRIMARY\P\';
    elseif (strfind(dinfo.Manufacturer,'GE')),
        % GE mag and phase images alternates.
        % e.g. i1110030.MRDC.1 is magnitude; i1110031.MRDC.2 is phase
        GEflag = 1;
        vendor = 'GE';
        magmatch = 'ORIGINAL\PRIMARY\OTHER'; %
        phamatch = 'ORIGINAL\PRIMARY\M'; % need to be revised
        if (mod(flen,2)~=0),
            fprintf(1,'GE data: #files=%d. Need even number of files.\n',flen);
            outMat=[]; procInfo=[]; dcmInfo=[]; return;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %fprintf(1,'CHECK IMAGETYPE for GE Dicom images (%s), stop\n');
	    %outMat=[]; procInfo=[]; dcmInfo=[]; return;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif (strfind(dinfo.Manufacturer,'Philips')),
        Philipsflag = 1;
        vendor = 'Philips';
        magmatch = '\M\';
        phamatch = '\PHASE MAP\';
        realmatch = '\Real\'; % may need to be revised
        imagmatch = '\Imag\';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %fprintf(1,'CHECK IMAGETYPE for Philips Dicom images, stop\n');
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else
        fprintf(1,'not recognized Manufacturer\n');
	    outMat=[]; procInfo=[]; dcmInfo=[]; return;
    end
    if (debugLevel>1), fprintf(1,'vendor: %s\n',vendor); end
% name strings: dicom tags to read. Add if need more.
    namestring = {'FileModDate','SeriesDescription','EchoTime',...
        'SeriesNumber','RepetitionTime','EchoTrainLength',...
        'PixelBandwidth','SequenceName','PercentPhaseFieldOfView',...
        'StudyDate','Manufacturer','InstanceNumber',...
        'SliceLocation','SliceThickness','FlipAngle','EchoNumber',...
        'Rows','Columns','Width','Height','ProtocolName',...
        'ImageType'};
    lenstr = length(namestring);
% assign parameters
    %matching variables
    for kk=1:lenstr
        str = cell2mat(namestring(kk));
        ttxt = sprintf('if(isfield(dinfo,''%s'')), %s = dinfo.%s; end',str,str,str);
        %disp(ttxt); %for checking
        eval([ttxt]);
    end
    if (debugLevel>1), fprintf(1,'variable matching done\n'); end
    % trim
    if (exist('SequenceName','var')),
        ind1 = strfind(SequenceName,'*');
        SequenceName = SequenceName(ind1+1:end);
    else
        ind1 = [];
    end
% display procInfo
    if (exist('ProtocolName')),
        line1 = sprintf('%s Protocol: %s',vendor,ProtocolName);
    else
        line1 = sprintf('%s Protocol: %s',vendor,SequenceName);
    end
% for multiple images: (1) read in all images. (2) sort mag and phase
% images, if exist. Phase images will be at 2nd half of outMat
	slc1=dicomread(dinfo);
    slcIndex = InstanceNumber;
    if (exist('EchoNumber','var')),
        EchoNumberMax = EchoNumber;
    else
        EchoNumberMax = 0;
    end
    if (debugLevel>1), fprintf(1,'start looping through files\n'); end
	if flen==1,
        outMat = double(slc1);
        % for procInfo
        magFlag = 0;
        phaFlag = 0;
		line2 = sprintf('  TE/TR/FA/SlcThk=%.2fms/%.0fms/%.1fdeg/%.1fmm',EchoTime,RepetitionTime,FlipAngle,SliceThickness);
		line3 = sprintf('  matrix=%dx%d, FoV=%dx%d',Rows,Columns,Height,Width);
        if (strfind(ImageType,magmatch)),
            magFlag = 1;
        elseif (strfind(ImageType,phamatch2)),
            magFlag = 1;
        elseif (strfind(ImageType,phamatch)),
            phaFlag = 1;
            if (magFlag==1),
                slcIndex = slcIndex + skipNumber;
            end
        else
            fprintf(1,'neither mag nor phase. ImageType=%s\n',ImageType);
            procInfo=[];
        end
		if (magFlag==1)&&(phaFlag==0),
			line3 = [line3,', Magnitude images'];
		elseif (magFlag==0)&&(phaFlag==1),
			line3 = [line3,', Phase images'];
		elseif (magFlag==1)&&(phaFlag==1),
			line3 = [line3,', Magnitude+Phase images'];
		else
        end
		if (EchoNumberMax>1),
            ETline=sprintf('  %d TEs=',EchoNumberMax);
            for kk=1:EchoNumberMax-1, ETline=[ETline,sprintf('%.2f/',ETsort(kk))]; end
            ETline=[ETline,sprintf('%.1f ms',ETsort(EchoNumberMax))];
            procInfo = sprintf('%s\n%s\n%s\n%s',line1,line2,line3,ETline);
        else
            procInfo = sprintf('%s\n%s\n%s',line1,line2,line3);
        end
		disp(procInfo);
    else
        kk = 1;
        magFlag = 0;
        phaFlag = 0;
        magSlices = [];
        phaSlices = [];
        skipNumber = flen/2;
        if (strfind(ImageType,magmatch)),
            magFlag = 1;
            magSlices = [magSlices; kk InstanceNumber];
        elseif (strfind(ImageType,magmatch2)),
            magFlag = 1;
            magSlices = [magSlices; kk InstanceNumber];
        elseif (strfind(ImageType,phamatch)),
            phaFlag = 1;
            phaSlices = [phaSlices; kk InstanceNumber+skipNumber];
        else
            fprintf(1,'neither mag nor phase. Default: mag\n');
            magFlag = 1;
            magSlices = [magSlices; kk InstanceNumber];
            %procInfo=[];
        end
        % philips may output mag, phase, real, imag
        % still use magSlices and phaSlice to store real and imag data
        if (Philipsflag),
            if (strfind(ImageType,realmatch)), %using strfind maybe better than strcmp
                magFlag = 1;
                magSlices = [magSlices; kk InstanceNumber];
            elseif (strfind(ImageType,imagmatch)),
                phaFlag = 1;
                phaSlices = [phaSlices; kk InstanceNumber+skipNumber];
            else
                fprintf(1,'neither real nor imag\n');
                procInfo=[];
            end
        end
        if (debugLevel>2), fprintf(1,'%3d ',kk); end
        outMat=zeros([Rows,Columns,flen]);
        outMat(:,:,kk) = double(slc1);
        InstanceNumbers=zeros([flen,1]);
        if (exist('EchoTime','var')),
            EchoTimes=zeros([flen,1]);
            EchoTimes(kk)=EchoTime;
        end
        InstanceNumbers(kk)=InstanceNumber;
		for kk=2:flen,
			dinfo = dicominfo(ff(kk).name);
            %EchoTrainLength = dinfo.EchoTrainLength;
            %SliceLocation = dinfo.SliceLocation;
            InstanceNumber = dinfo.InstanceNumber;
            slcIndex = InstanceNumber;
            if (exist('EchoTime','var')),
                EchoTime = dinfo.EchoTime;
            end
            if (exist('EchoNumber','var')),
                EchoNumber = dinfo.EchoNumber;
                if (EchoNumber>EchoNumberMax), EchoNumberMax = EchoNumber; end
            end
            ImageType = dinfo.ImageType;
            if (debugLevel>4),
                fprintf(1,'%s ',ff0(kk).name);
                if (debugLevel>5),
                    fprintf(1,' type=%s ',ImageType);
                end
            end
			if (strfind(ImageType,magmatch)),
                magFlag = 1;
                magSlices = [magSlices; kk slcIndex];
            elseif (strfind(ImageType,magmatch2)),
                magFlag = 1;
                magSlices = [magSlices; kk slcIndex];
			elseif (strfind(ImageType,phamatch)),
                phaFlag = 1;
                phaSlices = [phaSlices; kk slcIndex+skipNumber];
				if (magFlag==1),
                    slcIndex = slcIndex + skipNumber;
                end
            else %default: mag
                magFlag = 1;
                magSlices = [magSlices; kk slcIndex];
            end
            InstanceNumbers(kk)=slcIndex;
            if (exist('EchoTime','var')),
                EchoTimes(kk)=EchoTime;
            end
            % output
			outMat(:,:,kk) = double(dicomread(dinfo));
            if (debugLevel>2), fprintf(1,'%3d ',kk); end
            if (debugLevel>2)&&((mod(kk,20)==0)||(mod(kk,flen)==0)), fprintf(1,'\n'); end
        end
        % swap order according to InstanceNumber
        % for mag+phase images: phase images at 2nd half of outMat
        tmp = zeros(size(outMat));
        lenMag = length(magSlices);
        lenPha = length(phaSlices);
        if (debugLevel>1), fprintf(1,'#(magintude, phase)images = (%d, %d)\n',lenMag,lenPha); end
        if ((lenMag<1)&&(lenPha<1)), procInfo=[]; fprintf(1,'check magSlices (%s) and phasSlices (%s)\n',magmatch,phamatch); return; end
        % Philips data are in order of mag+phase already, they don't mixed
        if (~Philipsflag),
            if (lenMag*lenPha)>0, % mag + phase images
                tmp(:,:,magSlices(:,2))=outMat(:,:,magSlices(:,1));
                tmp(:,:,phaSlices(:,2))=outMat(:,:,phaSlices(:,1));
            else
                if (lenMag>0), % magnitude images only
                    tmp(:,:,magSlices(:,2))=outMat(:,:,magSlices(:,1));
                else % phase images only: subtract skipNumber
                    phaSlices(:,2) = phaSlices(:,2) - skipNumber;
                    tmp(:,:,phaSlices(:,2))=outMat(:,:,phaSlices(:,1));
                end
            end
            outMat=tmp;
        end

        % find echo times
        if (exist('EchoTime','var')),
            uEchoTimes = unique(EchoTimes);
            [ETsort] = sort(uEchoTimes,'ascend');
        end
        lenETsort = length(ETsort);

	   	%display procInfo
        %if (debugLevel>1), fprintf(1,'\n'); end
        if (exist('EchoTime','var')),
            line2 = sprintf('  TE/TR/FA/SlcThk=%.2fms/%.0fms/%.1fdeg/%.1fmm',EchoTime,RepetitionTime,FlipAngle,SliceThickness);
        else
            line2 = sprintf('  SlcThk=%.1fmm',SliceThickness);
        end
		line3 = sprintf('  matrix=%dx%d, FoV=%dx%d',Rows,Columns,Height,Width);
		if (magFlag==1)&&(phaFlag==0),
			line3 = [line3,', Magnitude images'];
		elseif (magFlag==0)&&(phaFlag==1),
			line3 = [line3,', Phase images'];
		elseif (magFlag==1)&&(phaFlag==1),
			line3 = [line3,', Magnitude+Phase images'];
		else
		end
		if ((EchoNumberMax>1)&(lenETsort>1)),
            ETline=sprintf('  %d TEs=',EchoNumberMax);
            for kk=1:EchoNumberMax-1, ETline=[ETline,sprintf('%.2f/',ETsort(kk))]; end
            ETline=[ETline,sprintf('%.1f ms',ETsort(EchoNumberMax))];
            procInfo = sprintf('%s\n%s\n%s\n%s',line1,line2,line3,ETline);
        else
            procInfo = sprintf('%s\n%s\n%s',line1,line2,line3);
		end
		disp(procInfo);
		if (flen>1),
			fprintf(1,'total %d files/images\n',flen);
		end
    end
    if ((GEflag)&(flen>1)),
        reply = input('Place odd slices as 2nd half of 3D matrix? Y/N [Y]:','s');
        if isempty(reply)
            reply = 'Y';
        end
        if ((reply=='Y')|(reply=='y')),
            tmpOdd=outMat(:,:,1:2:end-1);
            tmpEven=outMat(:,:,2:2:end);
            outMat = cat(3, tmpOdd, tmpEven);
        end
    end
% total elapsed time
    if (debugLevel),
        fprintf(1,'elapsed time = %f\n',toc);
    end
return;


text:
function CNet(Net)

format compact
format long e
theta = linspace(0,2*pi,length(Net)+1);
xy = zeros(length(Net)+1,2);
x = cos(theta);
y = sin(theta);
xy(1:length(Net)+1,1) = x(1:length(Net)+1);
xy(1:length(Net)+1,2) = y(1:length(Net)+1);
figure, gplot(Net,xy,'.-');
set(gcf, 'Color', [1 1 1]);
axis('equal');
xlim([-1.1 1.1]);
ylim([-1.1 1.1]);
axis off;


text:
function [ fplanes,cplanes,info ] = tfjigsawsep( f, varargin )
%TFJIGSAWSEP Time frequency jigsaw puzzle tonal-transient separation
%   Usage:  fplanes = tfjigsawsep(f);
%           fplanes = tfjigsawsep(f,r1,r2);
%           fplanes = tfjigsawsep(f,r1,r2,p);
%           [fplanes, cplanes, info] = tfjigsawsep(...);
%
%   Input parameters:
%            f        : Input signal
%            r1       : Significance level of the tonal layer refered to
%                       a white noise reference
%            r2       : Same for the transient layer
%            p        : Proportionfactor of the supertilesizes relative
%                       to the time-, and frequency stepsize
%
%   Output parameters:
%           fplanes   : signallength-by-3 array containing the 3 produced
%                       layers, tonal in fplanes(:,1), transient in
%                       fplanes(:,2) and the noisy residual in fplanes(:,3).
%           cplanes   : 3-by-1 cellarray containing the Gabor coefficients
%                       for the individual layers
%
%   `tfjigsawsep(f)` applies the separation algorithm on the input signal *f*
%   and returns the tonal, the transient and the residual parts.
%   The following default values are used, *r1=r2=0.95*, *p=4* and for the
%   3 Gabor systems:
%
%       "Tonal" system:     g1 = {'hann',4096}; a1 = 512; M1 = 4096;
%
%       "Transient" system: g2 = {'hann',256};  a2 = 32;  M2 = 256;
%
%       "Residual" system:  g3 = {'hann',2048}; a3 = 512; M3 = 2048;
%
%   `tfjigsawsep(f,r1,r2)` works as before, but allows changing threshold
%   parameters *r1* and *r2*. Good values are in range [0.85,0.95].
%   *t2* sometimes has to be chosen larger (~ 1.05), eg. for
%   percussions in music signals.
%
%   `tfjigsawsep(f,r1,r2,p)` (recommended) additionally allows changing the
%   size of the supertiles to a1*p in timesamples and b2*p in
%   frequencysamples. The choice of this particular proportion is
%   reasonable since it provides equal numbers of coefficients of the two
%   Gabor systems in each supertile. Good values are in the range of [1,10],
%   but it depends very much on the type of signal.
%   E.g. for speech signals, higher values yield better results.
%
%   `[fplanes, cplanes, info] = tfjigsawsep(...)` additionally returns
%   a 3 element cell array *cplanes* containing |dgtreal| coefficients
%   of the respective separated layers and a structure *info*, with the
%   parameters used in the algorithm. The relationship between *fplanes*
%   and *cplanes* is the following:

%   Additional parameters:
%   ----------------------
%
%   The function accepts the following flags:
%
%       'ver2' Uses the second version of the algorithm.
%
%       'plot' Plots the separated waveforms and the coefficients.
%
%       'verbose' Information about the boundary condition.
%
%
%   and the following key-value pairs:
%
%       'wintype' Requested windowtype
%
%       'winsize1' and 'winsize2' Windowlengths
%
%       'a1' and 'a2' Time stepsizes
%
%       'M1' and 'M2' Number of frequency channels
%
%       'T' and 'F' Supertile sizes in samples manually - alternative to p!
%
%       'maxit' Maximum number of iterations. The default value is 15.
%
%   Algorithm:
%   ----------
%
%   The algorithm is based on [1]. It transforms a signal into a two-windowed
%   Gabor expansion such that one wide window shall lead to a high frequency
%   resolution (tonal layer is represented well) and a narrow one to a high
%   time resolution (transient layer is repr. well). The resulting Gabor
%   coefficients are respectively considered within rectangular 'supertiles'
%   in the time-frequency plane. An entropy criterion chooses those tiles
%   respectively, where tonal and transient parts of the signal are
%   represented better and are below a estimated threshold. The rest is set
%   to zero. The leftover Gabor coefficients are transformed back and
%   subtracted from the original signal. By applying this procedure
%   iteratively on the residual, tonal and transient layers emerge.
%
%   Examples:
%   ---------
%
%   The following example shows the basic usage of the function:::
%
%     % Load the glockenspiel test signal and add some noise
%     [f,fs] = gspi; f = f + 0.001*randn(size(f));
%     % Setup the parameters
%     p = 2; r1 = 0.92; r2 = 0.93;
%     [fplanes,cplanes,info] = tfjigsawsep(f,r1,r2,p,'plot','ver2','fs',fs);
%
%   See also: dgtreal plottfjigsawsep
%
%   References: jato07

%AUTHOR: Daniel Haider, 2017

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Remarks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   The residual condition is computed very heuristically, designing it   %
%   more flexible would be a nice improvement of the algorithm.           %
%   It would also be useful to provide parameter settings for specific    %
%   types of signals (ie. speech, music,...) .                            %
%                                                                         %
%   Version 1 of the algorithm (default) works particularly well for      %
%   speech signals, but also for depicting the transient layer            %
%   (eg. percussive elements) nicely in musical signals.                  %
%   Version 2 works particularly well for depicting a nice tonal layer.   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,upper(mfilename),0);

if W>1
    error('%s: Multichannel inputs are not supported.',upper(mfilename));
end

definput.keyvals.r1 = [];
definput.keyvals.r2 = [];
definput.keyvals.p = [];
definput.keyvals.T = [];
definput.keyvals.F = [];
definput.keyvals.wintype = 'hann';
definput.keyvals.winsize1 = 4096;
definput.keyvals.a1 = 512;
definput.keyvals.M1 = 4096;
definput.keyvals.winsize2 = 512;
definput.keyvals.a2 = 64;
definput.keyvals.M2 = 512;
definput.keyvals.maxit = 15;
definput.keyvals.fs = [];
definput.flags.algver = {'ver1','ver2'};
definput.flags.plot = {'noplot','plot'};
definput.flags.verbose = {'noverbose','verbose'};
[flags,kv] = ltfatarghelper({'r1','r2','p','T','F','wintype','winsize1','a1','M1','winsize2','a2','M2'},definput,varargin);

% significance level
r1 = kv.r1;
r2 = kv.r2;
% Gabor system settings
wintype = kv.wintype;
winsize1 = kv.winsize1;
a1 = kv.a1;
M1 = kv.M1;
winsize2 = kv.winsize2;
a2 = kv.a2;
M2 = kv.M2;
% supertile sizes
if isempty(kv.T) && isempty(kv.F)
    if isempty(kv.p)
        p = 4;
    else
        p = kv.p;
    end
    T = a1*p;
    F = dgtlength(Ls,a2,M2)/M2*p;
elseif isempty(kv.p)
    T = kv.T;
    F = kv.F;
else
    error('%s: Use EITHER the proportional setting OR set T and F manually.',upper(mfilename));
end

if winsize1 < winsize2
    error('%s: The tonal system uses a shorter window than the transient system!',upper(mfilename));
end

% why this and not setting the kv on top?
if xor(isempty(r1),isempty(r2))
    error('%s: Both r1 and r2 must be defined.',upper(mfilename));
else
    if isempty(r1), r1 = 0.95; end
    if isempty(r2), r2 = 0.95; end
end

% windows
g1 = {wintype,winsize1};
g2 = {wintype,winsize2};

L1 = dgtlength(Ls,a1,M1);
L2 = dgtlength(Ls,a2,M2);
b1 = L1/M1;
b2 = L2/M2;

complainif_notposint(T,'T',mfilename);
if ~( T < min([L1,L2]) )
    error('%s: Supertile length must be smaller than the signal length [%i]',...
          upper(mfilename),min([L1,L2]));
elseif ~( F < min([L1/2+1,L2/2+1]) )
    error('%s: Supertile height must be smaller than the frequency range [%i]',...
      upper(mfilename),min([L1/2+1,L2/2+1]));
elseif ~( T > max([a1,a1]) )
    error('%s: Supertile length must be larger than the time stepsizes [%i]',...
          upper(mfilename),max([a1,a1]));
elseif ~( F > max([b1,b2]) )
    error('%s: Supertile height must be larger than the frequency stepsizes [%i]',...
          upper(mfilename),max([b1,b2]));
end

% entropy reference from noise signal
% thresholds tau1,tau2 for tonal resp. transient layer are chosen to have
% a certain significance with respect to the estimated reference
[ref1,ref2] = noisest(Ls,g1,g2,a1,a2,M1,M2,T,F);
tau1 = ref1*r1;
tau2 = ref2*r2;

% initialization of residual, layers, min and max number of iterations
% and epsilon, the upper limit for the residual condition, which is
% computed at the kmin-th iteration
% R = postpad(f,L);
R = f;
l1 = zeros(Ls,1);
l2 = zeros(Ls,1);
k = 1;
kmin = 5;
kmax = kv.maxit;
epsilon = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% main loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% runs until all values in R are below epsilon to exclude single peaks
while all(abs(R) < epsilon) ~= 1

    switch flags.algver
       case 'ver1'
        [x1,x2]=jigsaw1(R,g1,g2,a1,a2,M1,M2,T,F,tau1,tau2);
       case 'ver2'
        [x1,x2]=jigsaw2(R,g1,g2,a1,a2,M1,M2,T,F,tau1,tau2);
    end

    % residual parts of the signal
    R = R-x1-x2;
    l1 = l1+x1;
    l2 = l2+x2;

    if k == kmax
        break
    end

    if k == kmin
        % epsilon is computed as upper limit, corresponding to
        % an empirical p-quantile
        epsilon = sort(abs(R),'ascend');
        epsilon = 5/2*epsilon(round(0.998*numel(epsilon)));
        if flags.do_verbose
            disp(['The current maximum in the residual is ',num2str(max(abs(R)))])
            disp(['An upper bound is estimated to be ',num2str(epsilon)])
        end
    end

    k = k+1;
end

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

if flags.do_verbose
    disp('max. number of iteration reached: condition for residual is not fulfilled!')
    disp('Try to increase t2 slightly or change the size of the supertiles.')
end

% fplanes as Ls-by-3 array containing the layers
fplanes = [l1,l2,R];
fplanes = postpad(fplanes,Ls);

% info structure
info.g1 = {wintype, winsize1};
info.g2 = {wintype, winsize2};
info.g3 = {wintype, 2048};
info.M1 = M1;
info.M2 = M2;
info.M3 = 2048;
info.a1 = a1;
info.a2 = a2;
info.a3 = 512;
info.supertilesizes = [F,T];
info.noiseentropy_tonal = ref1;
info.threshold_tonal = tau1;
info.noiseentropy_transient = ref2;
info.threshold_transient = tau2;

% cplanes as 3-by-1 cell array containing
% the gabor coefficients corresp. to the layers
if nargout > 1 || flags.do_plot
    cplanes = cell(3,1);
    cplanes{1} = dgtreal(l1,g1,a1,M1);
    cplanes{2} = dgtreal(l2,g2,a2,M2);
    cplanes{3} = dgtreal(R,info.g3,info.a3,info.M3);
end

% option for plots
if flags.do_plot
    plottfjigsawsep(fplanes,cplanes,info,'fs',kv.fs,'showbuttons','equalyrange');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% compiling functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ r ] = renyi( t,alpha )

% computes the Renyi entropy for an array
% yields high values for peaky data and
% low values for almost constant

if isempty(t)
    r = inf;
elseif norm(t) == 0
    r = inf;
else
    switch nargin
        case 1
            alpha = 2.4;
            r = (1/(1-alpha))*log2(sum(sum(abs(t).^(2*alpha)))*(sum(sum(abs(t).^2)).^(-alpha)));
        case 2
            if alpha <= 0 || alpha == 1
                error('alpha must be chosen positive and unequal to 1.');
            else
                r = (1/(1-alpha))*log2(sum(sum(abs(t).^(2*alpha))).*norm(t(:),2)^(-2*alpha));
            end
        otherwise
            error('This function takes at least one input argument.');
    end
end

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

function [ref1,ref2] = noisest (Ls,g1,g2,a1,a2,M1,M2,T,F)

% computes the entropy for a white noise signal within one supertile
% as estimation reference for the tresholds tau1,tau2

n = noise(Ls,'white');
n1 = dgtreal(n,g1,a1,M1);
n2 = dgtreal(n,g2,a2,M2);
r1 = n1(1:floor(F*M1/Ls),1:floor(T/a1));
r2 = n2(1:floor(F*M2/Ls),1:floor(T/a2));
ref1 = renyi(r1);
ref2 = renyi(r2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Version 1 of the jigsaw puzzle algorithm %%%%%%%%%%%%%%%%%%

function [x1,x2] = jigsaw1(R,g1,g2,a1,a2,M1,M2,T,F,tau1,tau2)

% error check?

% signal lengths
Ls = length(R);
L1 = dgtlength(Ls,a1,M1);
L2 = dgtlength(Ls,a2,M2);

% gabor transformations
c1 = dgtreal(R,g1,a1,M1); % M1/2+1-by-L/a1
c2 = dgtreal(R,g2,a2,M2); % M2/2+1-by-L/a2

% frequency hop sizes
b1 = L1/M1;
b2 = L2/M2;

% indices on the TF plane (L/2+1 x L), where the coefficients belong to
aa1 = 1:a1:L1;
aa2 = 1:a2:L2;
bb1 = 1:b1:L1/2+1;
bb2 = 1:b2:L2/2+1;

L = max([L1,L2]);

% find the indices for the coefficients in every single TF-supertile
for m=0:floor((L/2+1)/F)-1
    for n=0:floor(L/T)-1
        f1 = find(bb1<=(m+1)*F & bb1>m*F);
        f2 = find(bb2<=(m+1)*F & bb2>m*F);
        t1 = find(aa1<=(n+1)*T & aa1>n*T);
        t2 = find(aa2<=(n+1)*T & aa2>n*T);
        % look for parts of c1,c2 where tonal/transient parts are repr well
        [c1,c2] = decision1(c1,c2,f1,f2,t1,t2,tau1,tau2);
    end
end

% last column of remaining supertiles
for m=0:floor((L/2+1)/F-1)
    f1 = find(bb1<=(m+1)*F & bb1>m*F);
    f2 = find(bb2<=(m+1)*F & bb2>m*F);
    t1 = find(aa1<=L & aa1>floor(L/T)*T);
    t2 = find(aa2<=L & aa2>floor(L/T)*T);

    [c1,c2] = decision1(c1,c2,f1,f2,t1,t2,tau1,tau2);
end

% upper row of remaining supertiles
for n=0:floor(L/T-1)
    f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
    f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
    t1 = find(aa1<=(n+1)*T & aa1>n*T);
    t2 = find(aa2<=(n+1)*T & aa2>n*T);

    [c1,c2] = decision1(c1,c2,f1,f2,t1,t2,tau1,tau2);
end

% right upper remaining supertiles
f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
t1 = find(aa1<=L & aa1>floor(L/T)*T);
t2 = find(aa2<=L & aa2>floor(L/T)*T);

[c1,c2] = decision1(c1,c2,f1,f2,t1,t2,tau1,tau2);

% synthesis with the canonical dual windows
x1 = idgtreal(c1,{'dual',g1},a1,M1,Ls);
x2 = idgtreal(c2,{'dual',g2},a2,M2,Ls);

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

function [c1,c2] = decision1(c1,c2,f1,f2,t1,t2,tau1,tau2)

% decision procedure for version 1
%   case 1: both entropies are above their thresholds -> set them to zero
%   case 2: g1 has a better repr. and the entropy is below its
%           threshold -> set the tile corr. to g2 to zero
%   case 3: vice versa
%   case 4: g1 has a better repr. and its entropy is above its
%           threshold tau1 but despite of that, the entropy corr. to
%           g2 is below its threshold tau2 -> set the tile corr. to
%           g1 to zero
%   case 5: vice versa

E1 = renyi(c1(f1,t1));
E2 = renyi(c2(f2,t2));

if E1 > tau1 && E2 > tau2
    c1(f1,t1) = 0;
    c2(f2,t2) = 0;
else
    if (min([E1,E2]) == E1 && E1 < tau1) || (min([E1,E2]) == E2 && E2 > tau2 && E1 < tau1)
        c2(f2,t2) = 0;
    elseif (min([E1,E2]) == E2 && E2 < tau2) || (min([E1,E2]) == E1 && E1 > tau1 && E2 < tau2)
        c1(f1,t1) = 0;
    else
        warning('something went wrong with the decision criteria!');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Version 2 of the jigsaw puzzle algorithm %%%%%%%%%%%%%%%%%%

function [x1,x2] = jigsaw2(R,g1,g2,a1,a2,M1,M2,T,F,tau1,tau2)

% signal lengths
Ls = length(R);
L1 = dgtlength(Ls,a1,M1);
L2 = dgtlength(Ls,a2,M2);

% frequency hop sizes
b1 = L1/M1;
b2 = L2/M2;

% gabor transformation
c1 = dgtreal(R,g1,a1,M1); % M1/2+1-by-L/a1
c2 = dgtreal(R,g2,a2,M2); % M2/2+1-by-L/a2

% indices on the TF plane (L/2+1 x L), where the coefficients belong to
aa1 = 1:a1:L1;
aa2 = 1:a2:L2;
bb1 = 1:b1:L1/2+1;
bb2 = 1:b2:L2/2+1;

L = max([L1,L2]);

for m=0:floor((L/2+1)/F)-1
    for n=0:floor(L/T)-1
        % find the indices for the coefficients in each TF-supertile
        f1 = find(bb1<=(m+1)*F & bb1>m*F);
        f2 = find(bb2<=(m+1)*F & bb2>m*F);
        t1 = find(aa1<=(n+1)*T & aa1>n*T);
        t2 = find(aa2<=(n+1)*T & aa2>n*T);
        % keep the tonals
        [c1,c2] = decision2ton(c1,c2,f1,f2,t1,t2,tau1);
    end
end

% last column of remaining supertiles
for m=0:floor((L/2+1)/F-1)
    f1 = find(bb1<=(m+1)*F & bb1>m*F);
    f2 = find(bb2<=(m+1)*F & bb2>m*F);
    t1 = find(aa1<=L & aa1>floor(L/T)*T);
    t2 = find(aa2<=L & aa2>floor(L/T)*T);
    % keep the tonals
    [c1,c2] = decision2ton(c1,c2,f1,f2,t1,t2,tau1);
end

% upper row of remaining supertiles
for n=0:floor(L/T-1)
    f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
    f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
    t1 = find(aa1<=(n+1)*T & aa1>n*T);
    t2 = find(aa2<=(n+1)*T & aa2>n*T);
    % keep the tonals
    [c1,c2] = decision2ton(c1,c2,f1,f2,t1,t2,tau1);
end

% right upper remaining supertile
f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
t1 = find(aa1<=L & aa1>floor(L/T)*T);
t2 = find(aa2<=L & aa2>floor(L/T)*T);
% keep the tonals
[c1,~] = decision2ton(c1,c2,f1,f2,t1,t2,tau1);

% synthesis of the tonal parts with the canonical dual window of g1
x1 = idgtreal(c1,{'dual',g1},a1,M1,Ls);
RR = R-x1;

% the same procedure is now applied with respect to the narrow window g2

% gabor transformation
c1 = dgtreal(RR,g1,a1,M1); % M1/2+1-by-L/a1
c2 = dgtreal(RR,g2,a2,M2); % M2/2+1-by-L/a2

for m=0:floor((L/2+1)/F)-1
    for n=0:floor(L/T)-1
        % find the indices for the coefficients in each TF-supertile
        f1 = find(bb1<=(m+1)*F & bb1>m*F);
        f2 = find(bb2<=(m+1)*F & bb2>m*F);
        t1 = find(aa1<=(n+1)*T & aa1>n*T);
        t2 = find(aa2<=(n+1)*T & aa2>n*T);
        % keep the transients
        [~,c2] = decision2trans(c1,c2,f1,f2,t1,t2,tau2);
    end
end

% last column of remaining supertiles
for m=0:floor((L/2+1)/F-1)
    f1 = find(bb1<=(m+1)*F & bb1>m*F);
    f2 = find(bb2<=(m+1)*F & bb2>m*F);
    t1 = find(aa1<=L & aa1>floor(L/T)*T);
    t2 = find(aa2<=L & aa2>floor(L/T)*T);
    % keep the transients
    [~,c2] = decision2trans(c1,c2,f1,f2,t1,t2,tau2);
end

% upper row of remaining supertiles
for n=0:floor(L/T-1)
    f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
    f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
    t1 = find(aa1<=(n+1)*T & aa1>n*T);
    t2 = find(aa2<=(n+1)*T & aa2>n*T);
    % keep the transients
    [~,c2] = decision2trans(c1,c2,f1,f2,t1,t2,tau2);
end

% right upper remaining supertile
f1 = find(bb1<=L/2+1 & bb1>floor((L/2+1)/F)*F);
f2 = find(bb2<=L/2+1 & bb2>floor((L/2+1)/F)*F);
t1 = find(aa1<=L & aa1>floor(L/T)*T);
t2 = find(aa2<=L & aa2>floor(L/T)*T);
% keep the transients
[~,c2] = decision2trans(c1,c2,f1,f2,t1,t2,tau2);

% synthesis of the transient parts with the canonical dual window of g2
x2 = idgtreal(c2,{'dual',g2},a2,M2,Ls);

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

function [c1,c2] = decision2ton(c1,c2,f1,f2,t1,t2,tau1)
% decision procedure for tonals

E1 = renyi(c1(f1,t1));
E2 = renyi(c2(f2,t2));

if min([E1,E2]) == E2 || E1 > tau1
    c1(f1,t1) = 0;
end

function [c1,c2] = decision2trans(c1,c2,f1,f2,t1,t2,tau2)
% decision procedure for transients

E1 = renyi(c1(f1,t1));
E2 = renyi(c2(f2,t2));

if min([E1,E2]) == E1 || E2 > tau2
    c2(f2,t2) = 0;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%%  SVM FOR FAULT DIAGNOSIS
% Script containing the SVM models of the various faults, and the measures data you wish to test,
% and get results of the fault diagnosis method.By Nassim laouti, Date 13.11.2010.
% Latest modification 16.02.2012.
%%

clear all
close all
clc

global st1 st2 st3 st5 st6 st7 st8 st9 st10 st11 st12 st13 st14 st15 st16 st17 st18 theta

%% Parameters
% Constants using in the preprocessing block
c3=10;c4=1000;wr0=1;wg0=100;omega_n=11.11;xi=0.6;P_r=4.8e6;
T1=0.6;T2=0.02;T3=0.06;T4=0.08;T5=0.04;

% Simple time and time
Ts=1/100
Time=Ts:Ts:4400;

%% Load the data_test that contains all measures and the controller output
% data_test=[v_hub_m;omega_r_m1;omega_r_m2;omega_g_m1;omega_g_m2;Tau_g_m;P_g_m;Beta_1_m1;Beta_1_m2;Beta_2_m1;Beta_2_m2;Beta_3_m1;Beta_3_m2;Tau_g_r;Beta_r]

load data_test_3;
measures=measures';
wr0=measures(1,3);wg0=measures(1,5);

%% Load SVM models of the various faults
% Load the svm model of the sensor fault beta1m1 type fixed value
load model_cap1_beta1m1

% Load the svm model of the sensor fault beta1m2 type fixed value
load model_cap1_beta1m2

% Load the svm model of the sensor fault beta2 type gain factor
load model_cap_beta2_t2

% Load the svm model of the sensor fault omega_r_m1 type fixed value
load model_cap1_wr1

% Load the svm model of the sensor fault omega_r_m2 type fixed value
load model_cap1_wr2

% Load the svm model of the sensor fault omega_r_m2 type gain factor
load model_cap_wr2_t2

% Load the svm model of the sensor fault omega_r_m1 type gain factor
load model_cap_wr1_t2

% Load the svm model of the sensor fault omega_g_m1 type fixed value
load model_cap1_wg1

% Load the svm model of the sensor fault omega_g_m2 type fixed value
load model_cap1_wg2

% Load the svm model of the sensor fault omega_g_m2 type gain factor
load model_cap_wg2_t2

% Load the svm model of the sensor fault omega_g_m1 type gain factor
load model_cap_wg1_t2

% Load the svm model of the actuator fault Tau_g type offset
load model_act_couple

% Load the svm model of the system fault
load model_system

% Load the svm model of the sensor fault omega_g type gain factor
load model_cap_wg_t2

% Load the svm model of the sensor fault omega_r type gain factor
load model_cap_wr_t2

% Load the svm model of the sensor fault beta2m2 type gain factor
load model_cap_beta2m2_t2

% Load the svm model of the sensor fault beta2m2 type gain factor
load model_cap_beta2m1_t2

%% Simulation

 open_system('FDI_SVM_WT.mdl');

 tic
 sim('FDI_SVM_WT.mdl',Time(end));
 toc

author: Sable

repo: mcbench-benchmarks