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

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

% Last Modified by GUIDE v2.5 08-Jan-2019 21:37:20

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

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

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

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

    centerUIFigure(hObject);

    elem = varargin{1};
    setappdata(hObject,'elem',elem);

    allowFuPctEdit = varargin{2};
    setappdata(hObject,'allowFuPctEdit',allowFuPctEdit);

    curve = varargin{3};
    setappdata(hObject,'curve',curve);

    % Update handles structure
    guidata(hObject, handles);

    populateUI(elem, allowFuPctEdit, curve, handles);

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

function populateUI(elem, allowFuPctEdit, curve, handles)
    handles.titleLabel.String = sprintf('Edit %s Element', curve.getCurveName());

    handles.indepValLabel.String = elem.getIndepVarName();
    handles.fuRemainPctText.String = fullAccNum2Str(elem.indepVar);
    handles.indepUnitLbl.String = elem.getIndepVarUnit();

    handles.depValLabel.String = elem.getDepVarName();
    handles.throtModifierText.String = fullAccNum2Str(elem.depVar);
    handles.depUnitLbl.String = elem.getDepVarUnit();

    if(allowFuPctEdit == false)
        handles.fuRemainPctText.Enable = 'off';
    end

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

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        elem = getappdata(hObject,'elem');
        elem.indepVar = str2double(handles.fuRemainPctText.String);
        elem.depVar = str2double(handles.throtModifierText.String);

        varargout{1} = true;
        close(handles.lvd_editThrottleModifierProfileElementGUI);
    end

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

    if(isempty(errMsg))
        uiresume(handles.lvd_editThrottleModifierProfileElementGUI);
    else
        msgbox(errMsg,'Invalid Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    elem = getappdata(handles.lvd_editThrottleModifierProfileElementGUI,'elem');

    val = str2double(get(handles.fuRemainPctText,'String'));
    enteredStr = get(handles.fuRemainPctText,'String');
    numberName = elem.getIndepVarName();
    lb = elem.minIndepValue;
    ub = elem.maxIndepValue;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

    if(isempty(errMsg))
        if((val == elem.minIndepValue || val == elem.maxIndepValue) && elem.indepVar ~= val)
            errMsg = {sprintf('Cannot add points at %.1f %s or %.1f %s.  Edit these points from the existing list if necessary.', elem.minIndepValue, elem.maxIndepValue, elem.getIndepVarUnit(), elem.getDepVarUnit())};
        end
    end

    val = str2double(get(handles.throtModifierText,'String'));
    enteredStr = get(handles.throtModifierText,'String');
    numberName = elem.getDepVarName();
    lb = elem.minDepValue;
    ub = elem.maxDepValue;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

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

function fuRemainPctText_Callback(hObject, eventdata, handles)
% hObject    handle to fuRemainPctText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of fuRemainPctText as text
%        str2double(get(hObject,'String')) returns contents of fuRemainPctText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

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

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

function throtModifierText_Callback(hObject, eventdata, handles)
% hObject    handle to throtModifierText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of throtModifierText as text
%        str2double(get(hObject,'String')) returns contents of throtModifierText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

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

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

% --- Executes on key press with focus on lvd_editThrottleModifierProfileElementGUI or any of its controls.
function lvd_editThrottleModifierProfileElementGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_editThrottleModifierProfileElementGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_editThrottleModifierProfileElementGUI);
    end

author: Arrowstar

repo: ksptot


text:
function APPdownload(username, destImages, destAnnotations)
% Downloads the images and annotations collected with the LabelMe app.
%
% For more details visit:
% http://labelme2.csail.mit.edu/Release3.0/browserTools/php/iPhoneHelp.php
%
% Usage:
%
%   APPdownload(username)
%
%   or, you can specify the destination folders where the images and
%   annotations will be saved:
%
%   APPdownload(username, destImages, destAnnotations)
%
% Example:
%   APPdownload(username)
%   D=LMdatabase('Annotations');
%   LMdbshowscenes(D, 'Images');

if nargin == 1
    destImages = 'Images';
    destAnnotations = 'Annotations';
end

HOMEANNOTATIONS = ['http://labelme.csail.mit.edu/Annotations/users/' username];
HOMEIMAGES = ['http://labelme.csail.mit.edu/Images/users/' username];

% folder index
folder = {'iPhoneCollection'};

% create list of folders to copy, copy and rename folders in XML:
LMinstall(folder, destImages, destAnnotations, HOMEIMAGES, HOMEANNOTATIONS)

% make sure folder names are coherent
ADMINrenamefolder(destAnnotations);

author: CSAILVision

repo: LabelMeToolbox


text:
function [GAmp,GTime]=GxUser(p)

global VVar;

load(p.GxFile); % load Gx file, Gx file needs to contain GAmp, GTime
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

GAmp=GAmp(min(end,VVar.TRCount),:);
GTime=GTime(min(end,VVar.TRCount),:);

GAmp(1)=0;
GAmp(end)=0;

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function write_tsurf_to_obj_and_png(fileprefix, tsh)

  % Check that plot is set up for per-vertex data
  assert(strcmp(tsh.FaceColor,'interp'));
  % Even though you may pass it  with 'CData' seems it gets moved to
  % 'FaceVertexCData'
  assert(size(tsh.FaceVertexCData,2) == 1);
  assert(size(tsh.FaceVertexCData,1) == size(tsh.Vertices,1));

  % get color map
  CM = tsh.Parent.Colormap;
  im = permute(flipud(CM),[1 3 2]);
  % make it a square
  im = repmat(im,[1 size(im,1), 1]);
  % use at least 1024
  im = imresize(im,ceil(1024/size(im,1)),'nearest');
  imwrite(im,[fileprefix '.png']);

  V = tsh.Vertices;
  F = tsh.Faces;
  CA = caxis(tsh.Parent);
  D = min( max( (tsh.FaceVertexCData-CA(1))/(CA(2)-CA(1)) , 0+1/size(im,1)),1-1/size(im,1));

  save([fileprefix '.mat'],'V','F','CA','D','CM');

  writeOBJ([fileprefix '.obj'],V,F,[D D],F);
  text = fileread([fileprefix '.obj']);
  fp = fopen([fileprefix '.obj'],'w');
  fprintf(fp,'mtllib %s.mtl\n',fileprefix);
  fprintf(fp,'usemtl Textured\n');
  fprintf(fp,'%s',text);
  fclose(fp);

  fp = fopen([fileprefix '.mtl'],'w');
  fprintf(fp,'newmtl Textured\n');
  fprintf(fp,'Ka 1.000 1.000 1.000\n');
  fprintf(fp,'Kd 1.000 1.000 1.000\n');
  fprintf(fp,'Ks 0.000 0.000 0.000\n');
  fprintf(fp,'d 1.0\n');
  fprintf(fp,'illum 2\n');
  fprintf(fp,'# the diffuse texture map\n');
  fprintf(fp,'map_Kd %s.png\n',fileprefix);
  fclose(fp);

end

author: alecjacobson

repo: gptoolbox


text:
function yp = p07_fun ( neqn, t, y )

%*****************************************************************************80
%
%% P07_FUN evaluates the function for problem P07.
%
%  Discussion:
%
%    y1' = -y1 +   y2
%    y2' =  y1 - 2 y2 + y3
%    y3' =         y2 - y3
%    y1(0) = 2
%    y2(0) = 0
%    y3(0) = 1
%
%    3 equations.
%    Enright and Pryce nonstiff problem #B2.
%    Autonomous.
%
%    Note that the quantity (y1+y2+y3) is conserved by the exact solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wayne Enright, John Pryce,
%    Algorithm 648,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 1, pages 28-34.
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the derivative
%    function.
%
%    Output, real YP(NEQN), the value of the derivative function.
%
  yp = zeros ( neqn, 1 );

  yp(1) = - y(1) +       y(2);
  yp(2) =   y(1) - 2.0 * y(2) + y(3);
  yp(3) =                y(2) - y(3);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hrf,p] = spm_hrf(RT,P,T)
% Haemodynamic response function
% FORMAT [hrf,p] = spm_hrf(RT,p,T)
% RT   - scan repeat time
% p    - parameters of the response function (two Gamma functions)
%
%                                                           defaults
%                                                          {seconds}
%        p(1) - delay of response (relative to onset)          6
%        p(2) - delay of undershoot (relative to onset)       16
%        p(3) - dispersion of response                         1
%        p(4) - dispersion of undershoot                       1
%        p(5) - ratio of response to undershoot                6
%        p(6) - onset {seconds}                                0
%        p(7) - length of kernel {seconds}                    32
%
% T    - microtime resolution [Default: 16]
%
% hrf  - haemodynamic response function
% p    - parameters of the response function
%__________________________________________________________________________
%
% The parameters p(1:4) correspond to the shape and scale parameters of two
% probability density functions of the Gamma distribution (see spm_Gpdf.m),
% one corresponding to the main response and the other one to the
% undershoot.
% Note that the mean of the Gamma distribution is shape*scale and its mode
% is (shape-1)*scale.  This means that with the default values of the
% parameters the peak of the heamodynamic response function will be around
% 5 seconds.
%__________________________________________________________________________
% Copyright (C) 1996-2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_hrf.m 7721 2019-11-27 13:03:32Z guillaume $

%-Parameters of the response function
%--------------------------------------------------------------------------
try
    p = spm_get_defaults('stats.fmri.hrf');
catch
    p = [6 16 1 1 6 0 32];
end
if nargin > 1
    p(1:length(P)) = P;
end

%-Microtime resolution
%--------------------------------------------------------------------------
if nargin > 2
    fMRI_T = T;
else
    fMRI_T = spm_get_defaults('stats.fmri.t');
end

%-Modelled haemodynamic response function - {mixture of Gammas}
%--------------------------------------------------------------------------
dt  = RT/fMRI_T;
u   = [0:ceil(p(7)/dt)] - p(6)/dt;
hrf = spm_Gpdf(u,p(1)/p(3),dt/p(3)) - spm_Gpdf(u,p(2)/p(4),dt/p(4))/p(5);
hrf = hrf([0:floor(p(7)/RT)]*fMRI_T + 1);
hrf = hrf'/sum(hrf);

author: spm


text:
% given 3D point set of face, draw structure of face in 3D space, or
% given 2D reprojection point set of face, draw it in corresponding image
function drawFace(KeyPointSet, showKpsIdx)
%   KeyPointSet - 66x2 or 66x3

    dim = size(KeyPointSet, 2);

    % 9 sets of keypoints(outlines) to be draw using line function
    sets = {0:16; 17:21; 22:26; 27:30; 31:35;... % face, eyebrows(r,l), nose(t,b)
        [36:41 36]; [42:47 42]; [48:59 48]; [48, 60:62, 54, 63:65, 48]}; % eyes(r,l), mouth(o,i)
    % notice: 0 indexed here

    % draw connections
    for i = 1:size(sets,1)
        indices = sets{i} + 1; % convert to 1-indexed
        if dim == 3
            line(KeyPointSet(indices, 1), KeyPointSet(indices, 2), KeyPointSet(indices, 3)...
                , 'Color', [0 150 230]/255, 'LineWidth', 3);
        elseif dim == 2
            line(KeyPointSet(indices, 1), KeyPointSet(indices, 2), 'Color', [0 150 230]/255, 'LineWidth', 2);
        end
    end

    hold on

    % draw kps
    if dim == 3
        scatter3(KeyPointSet(:,1), KeyPointSet(:,2), KeyPointSet(:,3), 500, 'r.');
        if showKpsIdx
            text(KeyPointSet(:,1), KeyPointSet(:,2), KeyPointSet(:,3), num2str((0 : 65)'));
        end
    elseif dim == 2
        scatter(KeyPointSet(:,1), KeyPointSet(:,2), 200, 'r.');
        if showKpsIdx
            text(KeyPointSet(:,1), KeyPointSet(:,2), num2str((0 : 65)'));
        end
    end
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
function x = utl_check_fingerprint(x,opts,ctx,exp)
% Check whether the given argument is an inconsistent impure expression.
% Data = utl_check_fingerprint(Data,Options,Context,Expressions)
%
% The remaining arguments are used only when the function is used as an argstep in exp_beginfun.
% There, it serves as an argstep for 'filter'/'editing' functions, to check, fix up and warn about the value of
% inconsistent impure expressions (expressions referring to signals, in particular) and data set values
% (without any notion of expressions).
%
% See also:
%   hlp_fingerprint, exp_beginfun, exp_endfun
%
%                                        Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                        2010-04-15

if nargin < 2
    opts.fingerprint_check = true;
    exp = [];
    ctx = [];
end

if isfield(x,'tracking') && all(isfield(x.tracking,{'expression','fingerprint'}))
    % get the fingerprint checking expression (either from the options or from the dynamic context)
    if isempty(opts.fingerprint_check)
        opts.fingerprint_check = hlp_resolve('fingerprint_check',true,ctx); end
    % this check can be selectively overridden with a block-scoped symbol 'fingerprint_check', which may be an expression that involves @expression
    % (the expression for which we are about to execute the check)
    if hlp_microcache('fprint_lookup',@is_enabled,opts.fingerprint_check,exp)
        if ~isequal(hlp_fingerprint(rmfield(x,'tracking')),x.tracking.fingerprint)
            error('The fingerprint of the given data structure does not match that of its expression.'); end
    end
end

% find out whether a given reference expression yields true if @expression is substituted with some substitution expression
function res = is_enabled(ref_exp,subs_exp)
if exp_eval(utl_releasehold(utl_replacerepeated(ref_exp,{exp_rule(@expression,subs_exp)})),inf)
    res = true;
else
    res = false;
end

repo: REST


text:
function CrowdDis = ConvergenceScore(front, p)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Annibale Panichella

[m,~] = size(front);
CrowdDis = zeros(1,m);

for i=1:m
    CrowdDis(i) = -norm(front(i,:),p);
end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% Select Grain Boundaries
%
%%
% In this section we explain how to extract specific grain boundaries.
% Therefore we start by importing some EBSD data and reconstructing the
% grain structure.

close all; plotx2east

% import the data
mtexdata forsterite silent

% restrict it to a subregion of interest.
ebsd = ebsd(inpolygon(ebsd,[5 2 10 5]*10^3));

[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));

% remove very small grains
ebsd(grains(grains.grainSize <= 5)) = [];

% and recompute grains
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));

% smooth the grains a bit
grains = smooth(grains,4);

% visualize as a phase map
plot(ebsd)
hold on
plot(grains.boundary,'linewidth',2)
hold off

%%
% The output of

grains.boundary

%%
% tells us the number of boundary segments between the different phsaes.
% Those segments with notIndexed phase include also those boundary segments
% where the grains are cutted by the scanning boundary. To restrict the
% grain boundaries to a specific phase transistion you shall do

hold on
plot(grains.boundary('Fo','Fo'),'lineColor','blue','micronbar','off','lineWidth',2)
hold off

%%
% Similarly we may select all Forsterite to enstatite boundary segements.

hold on
plot(grains.boundary('Fo','En'),'lineColor','darkgreen','micronbar','off','lineWidth',2)
hold off

%%
% Note, that the order of the phase names matter when considering the
% corresponding misorintations

grains.boundary('Fo','En').misorientation(1)
grains.boundary('En','Fo').misorientation(1)

%%
% In the fist case the misorientation returned is from Forsterite to
% Enstatite and in the second case its exactly the inverse

%%
% The selection of grain boundaries according to specific misorientationsm
% according to twist / tild character or twinning is explained in linked
% sections.
%

author: mtex-toolbox

repo: mtex

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


text:
%FEATSELO Branch and bound feature selection
%
%   [W,R] = FEATSELO(A,CRIT,K,T)
%   [W,R] = A*FEATSELO([],CRIT,K,T)
%   [W,R] = A*FEATSELO(CRIT,K,T)
%   [W,R] = FEATSELO(A,CRIT,K,N)
%   [W,R] = A*FEATSELO([],CRIT,K,N)
%   [W,R] = A*FEATSELO(CRIT,K,N)
%
% INPUT
%   A     Input dataset
%   CRIT  String name of the criterion or untrained mapping
%           (optional, def= 'maha-s')
%   K     Numner of features to select (optional, def: K=2)
%   T     Validation set (optional)
%   N     Number of cross-validations (optional)
%
% OUTPUT
%   W     Output feature selection mapping
%   R     Matrix with step-by-step results
%
% DESCRIPTION
% Backward selection of K features by baktracking using the branch
% and bound procedure on the data set A. CRIT sets the criterion
% used by the feature evaluation routine FEATEVAL. If the data set T
% is given, it is used as test set for FEATEVAL. Alternatively a number
% of cross-validations N may be supplied. The resulting W can be used for
% the selecting features of a dataset B by B*W.
% The selected features are stored in W.DATA and can be found by +W.
%
% This procedure finds the optimum feature set if a monotoneous
% criterion is used. The use of a testset does not guarantee that.
%
% REFERENCE
% P. M. Narendra and K. Fukunaga
% A Branch and Bound Algorithm for Feature Subset Selection,
% IEEE Trans. Computer, 26(9), pp. 917-922, September 1977
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELF, FEATSELB, FEATSELI,
% FEATSEL, FEATSELP, FEATSELM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: featselo.m,v 1.7 2009/07/01 09:33:23 duin Exp $

function [W,R] = featselo(varargin)

  varargin = shiftargin(varargin,{'char','prmapping'});
  argin = setdefaults(varargin,[],'maha-s',2,[],[]);
  if mapping_task(argin,'definition')
    W = define_mapping(argin,'untrained','B&B FeatSel');
    return
  end

  [A,crit,kmin,T,fid] = deal(argin{:});

	isvaldfile(A,1,2); % at least 1 object per class, 2 classes
	A = testdatasize(A);
	if isdataset(T), iscomdset(A,T); end

	[m,k,c] = getsize(A);
	featlist = getfeatlab(A);
  A = setprior(A,getprior(A));

	if ((kmin < 1) | (kmin >= k))
		error('The desired feature size should be > 0 and < dataset feature size')
	end

	% space for criteria values
	feat = zeros(1,k);

   % Get performance of the individual features:
	if isempty(T)
		for j=1:k
			feat(j) = feateval(A(:,j),crit);
		end
	elseif is_scalar(T)
		for j=1:k
			feat(j) = feateval(A(:,j),crit,T);
		end
	else
		for j=1:k
			feat(j) = feateval(A(:,j),crit,T(:,j));
		end
	end

   % Get the kmin worst(?) individual features according to their
   % individual performance:
	[F,S] = sort(feat);

	%sometimes the above line is bad compared to the following two
	%w = featselb(A,crit,[]);
	%S = fliplr(+w);

	Iopt = [k-kmin+1:k];

	I = [1:k];
	J = [zeros(1,kmin),1:(k-kmin-1),k-kmin+1,k+1];
	level = k;

   % Get the performance of Iopt
	if isdataset(T)
		bound = feateval(A(:,S(Iopt)),crit,T(:,S(Iopt)));
	elseif is_scalar(T)
		bound = feateval(A(:,S(Iopt)),crit,T);
	else
		bound = feateval(A(:,S(Iopt)),crit);
	end

	C = inf;
	prwaitbar(100,'Branch & Bound Feature Selection')
	iter = 0;
	while numel(I) > 0 && J(k+1) == k+1;
		iter = iter+1;
		prwaitbar(100,100-100*exp(-iter/25),['Branch & Bound Feature Selection: ' num2str(iter)]);
		if J(level) == J(level+1) | level <= kmin | C <= bound
			J(level) = level - kmin;
			level = level + 1;
			I = sort([I,J(level)]);
			J(level) = J(level) + 1;
			C = inf;
		else
			I(J(level)) = [];
			level = level - 1;
			if J(level+1) < 3 & level == kmin+1 & 0 % never happens ??
				;
			else
				if isdataset(T)
					C = feateval(A(:,S(I)),crit,T(:,S(I)));
				elseif is_scalar(T)
					C = feateval(A(:,S(I)),crit,T);
				else
					C = feateval(A(:,S(I)),crit);
        end
				if level == kmin & C > bound
					bound = C;
					Iopt = I;
          disp([bound,iter,numel(I)])
        end
			end
		end
  end
	prwaitbar(0);

   % Store the optimal features in the mapping:
	W = featsel(k,S(Iopt));
  W = setmapping_type(W,'trained');
  W = setsize(W,[k length(S(Iopt))]);
	if ~isempty(featlist)
		W = setlabels(W,featlist(S(Iopt),:));
	end
	W = setname(W,'B&B FeatSel');

	R = [];  %DXD I'm still not sure what to return

return

author: marianux

repo: ecg-kit

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


text:
function jed = ymd_to_jed_gregorian ( y, m, d )

%*****************************************************************************80
%
%% YMD_TO_JED_GREGORIAN converts a Gregorian YMD date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, M, D, the YMD date.
%
%    Output, real JED, the corresponding JED.
%

%
%  Check the date.
%
  [ y, m, d, ierror ] = ymd_check_gregorian ( y, m, d );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Account for the missing year 0 by moving negative years up one.
%
  y2 = y_common_to_astronomical ( y );
%
%  Convert the calendar date to a computational date.
%
  y_prime = y2 + 4716 - floor ( ( 14 - m ) / 12 );
  m_prime = mod ( m + 9, 12 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 1461 * y_prime ) / 4 );

  j2 = floor ( ( 153 * m_prime + 2 ) / 5 );

  g = floor ( ( y_prime + 184 ) / 100 );

  g = floor ( 3 * g / 4 ) - 38;

  jed = j1 + j2 + d_prime - 1401 - g - 0.5;

  return
end


text:
function b = ne(q1,q2)
% q1 ~= q2 ?

b = ~ (q1 == q2);


text:
%Navigation frame=Local geodetic Frame.
%Position is mechanized in Cen (essentially Ceg) and h.
%Attitude is mechanized in quat

%Compare this with strapdown_Cen_dcm and strandown_wander_quat. Both of
%these strapdowns implement wander frame mechanization. However, only the
%2nd one is non-singular. Furthermore, also note that the only difference
%between this script and the strandown_wander_quat is the 1st argument of
%geoparam function.

function [qbn_new, Vn_new, Cen_new, h_new]=strapdown_Cen_quat(qbn, Vn, Cen, h, velinc, anginc, dt)
%Compute the transport rate of geodetic frame (also the gravity)
[Fc, wen_n, wie_n, g]=geoparam_v001(2, Cen(:,3), h, Vn); %Compute the curvature matrix for local geodetic frame

%% Update the velocity
vel_inc1=quatrot_v000(qbn,velinc,0); %Note: 'a' is assumed to be output of a sculling module.
vel_inc2=(cross(Vn,2*wie_n+wen_n)+[0;0;g])*dt;
Vn_new=Vn+vel_inc1+vel_inc2;

%%update the attitude
%Body frame updates
qb=rvec2quat_v000(anginc);    %Note: 'w' is assumed to be the output of coning module
qbn_new=quatmult_v000(qbn,qb);

%Navigation frame updates
qn=rvec2quat_v000(-(wen_n+wie_n)*dt);
qbn_new=quatmult_v000(qn,qbn_new);

%%Position update
Vmid=(Vn+Vn_new)/2;
wen_n=Fc*[Vmid(2);-Vmid(1);0];
Cn=rot2dcm_v000(-wen_n*dt);
Cen_new=Cn*Cen;

h_new=h-Vmid(3)*dt;

repo: ekfmonoslam


text:
% add the path of RBM code
addpath('..');
addpath('~/work/Algorithms/liblinear-1.7/matlab');

use_whitening = 1;

% load natural image patches
load 'bsds500bw_patches_8.mat';
X = (Xbw / 255);

% shuffle the training data
perm_idx = randperm (size(X,1));

n_all = size(X, 1);
n_train = ceil(n_all * 3 / 4);
n_valid = floor(n_all /4);

X_valid = X(perm_idx(n_train+1:end), :);
X = X(perm_idx(1:n_train), :);

if use_whitening
    %% ZCA
    %[Z, Wsep, Wmix, mX] = zca(X, 0.1);
    load patch8_whiten.mat;
    X = zca_whiten(X, Wsep, Wmix, mX);
    X_valid = zca_whiten(X_valid, Wsep, Wmix, mX);
end

% construct RBM and use default configurations
D = default_dae (size(X, 2), 320);

if use_whitening
    D.do_normalize = 0;
    D.do_normalize_std = 0;
end

D.data.binary = 0;
D.hidden.binary = 0;

D.learning.lrate = 1e-1;
D.learning.lrate0 = 10000;
%D.learning.momentum = 0.5;
D.learning.weight_decay = 0;

D.adagrad.use = 1;

D.rica.cost = 0.01;

D.noise.drop = 0;
D.noise.level = 0;

D.sparsity.cost = 0;
D.sparsity.target = 0;

% max. 100 epochs
D.iteration.n_epochs = 500;

% set the stopping criterion
D.stop.criterion = 0;
D.stop.recon_error.tolerate_count = 1000;

% save the intermediate data after every epoch
D.hook.per_epoch = {@save_intermediate, {'rica_patch8.mat'}};

% print learining process
D.verbose = 0;

% display the progress
D.debug.do_display = 0;

fprintf(1, 'Training rICA\n');
tic;
D = dae (D, X, X_valid, 0.1);
fprintf(1, 'Training is done after %f seconds\n', toc);

author: kyunghyuncho

repo: deepmat


text:
function [nrm] = lognrm(tt)
%Frobenius norm of the TT-tensor
%   [NRM]=LOGNRM(TT) Computes log10 of the Frobenius norm of a tensor
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

d=tt.d;
n=tt.n;
r=tt.r;
pos=tt.ps;
cr=tt.core;
pos1=1;
nrm=zeros(d,1);
core0=cr(1:r(1)*n(1)*r(2));
%Orthogonalization from left-to-tight
for i=1:d-1
   core0=reshape(core0,[r(i)*n(i),r(i+1)]);
   [core0,ru]=qr(core0,0); nrm(i)=norm(ru,'fro');
   nrm(i)=max(nrm(i),1e-308);
   ru=ru./nrm(i);
   core1=cr(pos(i+1):pos(i+2)-1);
   core1=reshape(core1,[r(i+1),n(i+1)*r(i+2)]);
   core1=ru*core1;
   r(i+1)=size(core0,2);
   cr(pos1:pos1-1+r(i)*n(i)*r(i+1))=core0(:);
   cr(pos1+r(i)*n(i)*r(i+1):pos1+r(i)*n(i)*r(i+1)+r(i+1)*n(i+1)*r(i+2)-1)=core1(:);
   core0=core1;
   pos1=pos1+r(i)*n(i)*r(i+1);
end
pos1=pos1+r(d)*n(d)*r(d+1)-1;
nrm(d)=norm(core0(:));
nrm=sum(log10(nrm));
return
end

author: oseledets

repo: TT-Toolbox


text:
function sphere_triangle_quad_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests SPHERE01_TRIANGLE_QUAD_00.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';
  n_mc1 = 1000;
  n_mc2 = 10000;
  n_mc3 = 100000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Approximate an integral on a random spherical triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  QUAD_MC1 uses a Monte Carlo method with %d points.\n', n_mc1 );
  fprintf ( 1, '  QUAD_MC2 uses a Monte Carlo method with %d points.\n', n_mc2 );
  fprintf ( 1, '  QUAD_MC3 uses a Monte Carlo method with %d points.\n', n_mc3 );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  QUAD_MC1      QUAD_MC2      QUAD_MC3\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

    [ result_mc1, seed ] = sphere01_triangle_quad_00 ( n_mc1, v1, v2, v3, ...
      @polyterm_value_3d, seed );
    [ result_mc2, seed ] = sphere01_triangle_quad_00 ( n_mc2, v1, v2, v3, ...
      @polyterm_value_3d, seed );
    [ result_mc3, seed ] = sphere01_triangle_quad_00 ( n_mc3, v1, v2, v3, ...
      @polyterm_value_3d, seed );

    fprintf ( 1, '  %14.6g  %14.6g  %14.6g\n', ...
      result_mc1, result_mc2, result_mc3 );

  end

  return
end


text:
function res = isequal(this, that)
% Method to check if 2 MEEG objects are the same
% FORMAT res = isequal(this, that)
% _______________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id: isequal.m 5025 2012-10-31 14:44:13Z vladimir $

res = isequal(struct(this), struct(that));


text:
function imgOut = normalizeFromAnything(img)
%
%       imgOut = normalizeFromAnything(img)
%
%
%       input:
%         - img: an image
%
%       output:
%         - imgOut: an image with value in [0,1]
%
%     Copyright (C) 2015  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Video Matching"
% 	  by Peter Sand and Seth Teller, MIT Computer Science and Aritficial
% 	  Intelligence Laboratory, Technical Report, November 11 2004.
%

imgOut = [];

if(isa(img, 'uint8'))
    imgOut = single(img) / 255.0;
end

if(isa(img, 'uint16'))
    imgOut = single(img) / 65535.0;
end

if(isa(img, 'double') | isa(img, 'single'))
    max_val = max(img(:));
    if(max_val > 1.0)
        imgOut = img / max_val;
    end
end

if(isempty(imgOut))
    imgOut = img;
end

end

author: banterle


text:
%% Define the objective function.
addpath(genpath(pwd()));
cd AC_blackbox_eval/blackbox_eval_source/
config_file = 'config_file.txt';
numRun = 1;
scenario_file = 'scenario-lpsolve-CORLAT-inst.txt';
[func, obj] = define_func_handle(config_file, scenario_file, numRun);
cd ../..

%% Optimize.
yaml = 'lpsolve.yaml';
model = highdim_opt(yaml, obj);

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function result = torus_square_14c ( func, r1, r2 )

%*****************************************************************************80
%
%% TORUS_SQUARE_14C approximates an integral in a "square" torus in 3D.
%
%  Discussion:
%
%    A 14-th degree 960 point formula is used.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      R1 - R2 <= SQRT ( X**2 + Y**2 ) <= R1 + R2,
%       -R2 <= Z <= R2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function of three variables which is to be integrated, of the form:
%      function value = func ( x, y, z )
%
%    Input, real R1, R2, the radii that define the torus.
%
%    Output, real RESULT, the approximate integral of the function.
%
  norder = 8;

  [ rtab, weight ] = legendre_set ( norder );

  w = 1.0E+00 / ( 60.0E+00 * r1 );
  quad = 0.0E+00;

  for n = 1 : 15

    angle = 2.0E+00 * pi * n / 15.0E+00;
    cth = cos ( angle );
    sth = sin ( angle );

    for i = 1 : norder

      u = r1 + rtab(i) * r2;
      x = u * cth;
      y = u * sth;

      for j = 1 : norder
        z = rtab(j) * r2;
        quad = quad + u * w * weight(i) * weight(j) * feval ( func, x, y, z );
      end

    end

  end

  volume = torus_square_volume_3d ( r1, r2 );
  result = quad * volume;

  return
end


text:
function d = stoi(x, y, fs_signal)
% The Short-Time Objective Intelligibility measure
%   d = stoi(x, y, fs_signal) returns the output of the short-time
%   objective intelligibility (STOI) measure described in [1, 2], where x
%   and y denote the clean and processed speech, respectively, with sample
%   rate fs_signal in Hz. The output d is expected to have a monotonic
%   relation with the subjective speech-intelligibility, where a higher d
%   denotes better intelligible speech. See [1, 2] for more details.
%
%   References:
%      [1] C.H.Taal, R.C.Hendriks, R.Heusdens, J.Jensen 'A Short-Time
%      Objective Intelligibility Measure for Time-Frequency Weighted Noisy
%      Speech', ICASSP 2010, Texas, Dallas.
%
%      [2] C.H.Taal, R.C.Hendriks, R.Heusdens, J.Jensen 'An Algorithm for
%      Intelligibility Prediction of Time-Frequency Weighted Noisy Speech',
%      IEEE Transactions on Audio, Speech, and Language Processing, 2011.
%
%
% Copyright 2009: Delft University of Technology, Signal & Information
% Processing Lab. The software is free for non-commercial use. This program
% comes WITHOUT ANY WARRANTY.
%
%
%
% Updates:
% 2011-04-26 Using the more efficient 'taa_corr' instead of 'corr'

if length(x)~=length(y)
    error('x and y should have the same length');
end

% initialization
x           = x(:);                             % clean speech column vector
y           = y(:);                             % processed speech column vector

fs          = 10000;                            % sample rate of proposed intelligibility measure 10000
N_frame    	= 256;                              % window support 256
K           = 512;                              % FFT size 512
J           = 15;                               % Number of 1/3 octave bands
mn          = 150;                              % Center frequency of first 1/3 octave band in Hz.
H           = thirdoct(fs, K, J, mn);           % Get 1/3 octave band matrix
N           = 30;                               % Number of frames for intermediate intelligibility measure (Length analysis window)
Beta        = -15;                           	% lower SDR-bound -15
dyn_range   = 40;                               % speech dynamic range 40

% resample signals if other samplerate is used than fs
if fs_signal ~= fs
    x	= resample(x, fs, fs_signal);
    y 	= resample(y, fs, fs_signal);
end

% remove silent frames
[x y] = removeSilentFrames(x, y, dyn_range, N_frame, N_frame/2);

% apply 1/3 octave band TF-decomposition
x_hat     	= stdft(x, N_frame, N_frame/2, K); 	% apply short-time DFT to clean speech
y_hat     	= stdft(y, N_frame, N_frame/2, K); 	% apply short-time DFT to processed speech

x_hat       = x_hat(:, 1:(K/2+1)).';         	% take clean single-sided spectrum
y_hat       = y_hat(:, 1:(K/2+1)).';        	% take processed single-sided spectrum

X           = zeros(J, size(x_hat, 2));         % init memory for clean speech 1/3 octave band TF-representation
Y           = zeros(J, size(y_hat, 2));         % init memory for processed speech 1/3 octave band TF-representation

for i = 1:size(x_hat, 2)
    X(:, i)	= sqrt(H*abs(x_hat(:, i)).^2);      % apply 1/3 octave bands as described in Eq.(1) [1]
    Y(:, i)	= sqrt(H*abs(y_hat(:, i)).^2);
end

% loop al segments of length N and obtain intermediate intelligibility measure for all TF-regions
d_interm  	= zeros(J, length(N:size(X, 2)));                               % init memory for intermediate intelligibility measure
c           = 10^(-Beta/20);                                                % constant for clipping procedure

for m = N:size(X, 2)
    X_seg  	= X(:, (m-N+1):m);                                              % region with length N of clean TF-units for all j
    Y_seg  	= Y(:, (m-N+1):m);                                              % region with length N of processed TF-units for all j
    alpha   = sqrt(sum(X_seg.^2, 2)./sum(Y_seg.^2, 2));                     % obtain scale factor for normalizing processed TF-region for all j
    aY_seg 	= Y_seg.*repmat(alpha, [1 N]);                               	% obtain \alpha*Y_j(n) from Eq.(2) [1]
    for j = 1:J
      	Y_prime             = min(aY_seg(j, :), X_seg(j, :)+X_seg(j, :)*c); % apply clipping from Eq.(3)
        d_interm(j, m-N+1)  = taa_corr(X_seg(j, :).', Y_prime(:));          % obtain correlation coeffecient from Eq.(4) [1]
    end
end

d = mean(d_interm(:));                                                      % combine all intermediate intelligibility measures as in Eq.(4) [1]

%%
function  [A cf] = thirdoct(fs, N_fft, numBands, mn)
%   [A CF] = THIRDOCT(FS, N_FFT, NUMBANDS, MN) returns 1/3 octave band matrix
%   inputs:
%       FS:         samplerate
%       N_FFT:      FFT size
%       NUMBANDS:   number of bands
%       MN:         center frequency of first 1/3 octave band
%   outputs:
%       A:          octave band matrix
%       CF:         center frequencies

f               = linspace(0, fs, N_fft+1);
f               = f(1:(N_fft/2+1));
k               = 0:(numBands-1);
cf              = 2.^(k/3)*mn;
fl              = sqrt((2.^(k/3)*mn).*2.^((k-1)/3)*mn);
fr              = sqrt((2.^(k/3)*mn).*2.^((k+1)/3)*mn);
A               = zeros(numBands, length(f));

for i = 1:(length(cf))
    [a b]                   = min((f-fl(i)).^2);
    fl(i)                   = f(b);
    fl_ii                   = b;

	[a b]                   = min((f-fr(i)).^2);
    fr(i)                   = f(b);
    fr_ii                   = b;
    A(i,fl_ii:(fr_ii-1))	= 1;
end

rnk         = sum(A, 2);
numBands  	= find((rnk(2:end)>=rnk(1:(end-1))) & (rnk(2:end)~=0)~=0, 1, 'last' )+1;
A           = A(1:numBands, :);
cf          = cf(1:numBands);

%%
function x_stdft = stdft(x, N, K, N_fft)
%   X_STDFT = X_STDFT(X, N, K, N_FFT) returns the short-time
%	hanning-windowed dft of X with frame-size N, overlap K and DFT size
%   N_FFT. The columns and rows of X_STDFT denote the frame-index and
%   dft-bin index, respectively.

frames      = 1:K:(length(x)-N);
x_stdft     = zeros(length(frames), N_fft);

w           = hanning(N);
x           = x(:);

for i = 1:length(frames)
    ii              = frames(i):(frames(i)+N-1);
	x_stdft(i, :) 	= fft(x(ii).*w, N_fft);
end

%%
function [x_sil y_sil] = removeSilentFrames(x, y, range, N, K)
%   [X_SIL Y_SIL] = REMOVESILENTFRAMES(X, Y, RANGE, N, K) X and Y
%   are segmented with frame-length N and overlap K, where the maximum energy
%   of all frames of X is determined, say X_MAX. X_SIL and Y_SIL are the
%   reconstructed signals, excluding the frames, where the energy of a frame
%   of X is smaller than X_MAX-RANGE

x       = x(:);
y       = y(:);

frames  = 1:K:(length(x)-N);
w       = hanning(N);
msk     = zeros(size(frames));

for j = 1:length(frames)
    jj      = frames(j):(frames(j)+N-1);
    msk(j) 	= 20*log10(norm(x(jj).*w)./sqrt(N));
end

msk     = (msk-max(msk)+range)>0;
count   = 1;

x_sil   = zeros(size(x));
y_sil   = zeros(size(y));

for j = 1:length(frames)
    if msk(j)
        jj_i            = frames(j):(frames(j)+N-1);
        jj_o            = frames(count):(frames(count)+N-1);
        x_sil(jj_o)     = x_sil(jj_o) + x(jj_i).*w;
        y_sil(jj_o)  	= y_sil(jj_o) + y(jj_i).*w;
        count           = count+1;
    end
end

x_sil = x_sil(1:jj_o(end));
y_sil = y_sil(1:jj_o(end));

%%
function rho = taa_corr(x, y)
%   RHO = TAA_CORR(X, Y) Returns correlation coeffecient between column
%   vectors x and y. Gives same results as 'corr' from statistics toolbox.
xn    	= x-mean(x);
xn  	= xn/sqrt(sum(xn.^2));
yn   	= y-mean(y);
yn    	= yn/sqrt(sum(yn.^2));
rho   	= sum(xn.*yn);

author: jtkim-kaist

repo: Speech-enhancement


text:
function plotSuperVsRectangleOptim
mCase = 'Middle';
folderR = ['/home/alex/Desktop/ExperimentingPlotRectangle',mCase,'/'];
[iter,cost,maxStress] = obtainCostShapes(folderR);
f = figure();
h{1} = plot(iter,[cost;maxStress]','b');
hline = findobj(gcf, 'type', 'line');
set(hline(1),'LineStyle','--')

folderS = ['/home/alex/Desktop/ExperimentingPlotSuperEllipse',mCase,'/'];
[iter,costS,maxStressS] = obtainCostShapes(folderS);
hold on
h{2} = plot(iter,[costS;maxStressS]','r');
hline = findobj(gcf, 'type', 'line');
set(hline(1),'LineStyle','--')

leg = {'Rectangle Amplified Stress p', 'Amplified Stress max',...
       'Superellipse Amplified Stress p', 'Superellipse Stress max'};
legend(leg)
printer = plotPrinter(f,h);
path = '/home/alex/Dropbox/GregoireMeetings/GregoireMeeting25Janvier';
output = fullfile(path,['Convergence',mCase]);
printer.print(output)

end

function [iter,stress,maxStress] = obtainCostShapes(folder)
fNameMon = fullfile(folder,'Monitoring.fig');
h = openfig(fNameMon);
handles = findobj(h,'Type','line');
iter = get(handles(5),'Xdata');
stress = get(handles(11),'Ydata');
maxStress = get(handles(4),'Ydata');
close(h)
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ a_lu, pivot, rcond, z ] = r8ge_co ( n, a )

%*****************************************************************************80
%
%% R8GE_CO factors a R8GE matrix and estimates its condition number.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%    For the system A * X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    EPSILON/RCOND.
%
%    If RCOND is so small that the logical expression
%      1.0E+00 + rcond == 1.0E+00
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate
%    underflows.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979
%
%  Parameters:
%
%    Input, integer N, the order of the matrix A.
%
%    Input, real A(N,N), a matrix to be factored.
%
%    Output, real A_LU(N,N), the LU factorization of the matrix.
%
%    Output, integer PIVOT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal condition number of A.
%
%    Output, real Z(N), a work vector whose contents are usually unimportant.
%    If A is close to a singular matrix, then Z is an approximate null vector
%    in the sense that
%      norm ( A * Z ) = RCOND * norm ( A ) * norm ( Z ).
%

%
%  Compute the L1 norm of A.
%
  anorm = 0.0E+00;
  for j = 1 : n
    anorm = max ( anorm, sum ( abs ( a(1:n,j) ) ) );
  end
%
%  Compute the LU factorization.
%
  [ a, pivot, info ] = r8ge_fa ( n, a );
%
%  RCOND = 1 / ( norm(A) * (estimate of norm(inverse(A))) )
%
%  estimate of norm(inverse(A)) = norm(Z) / norm(Y)
%
%  where
%    A * Z = Y
%  and
%    A' * Y = E
%
%  The components of E are chosen to cause maximum local growth in the
%  elements of W, where U'*W = E.  The vectors are frequently rescaled
%  to avoid overflow.
%
%  Solve U' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( z(k) ~= 0.0 )
      ek = - r8_sign ( z(k) ) * abs ( ek );
    end

    if ( abs ( a(k,k) ) < abs ( ek - z(k) ) )
      s = abs ( a(k,k) ) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );

    if ( a(k,k) ~= 0.0 )
      wk = wk / a(k,k);
      wkm = wkm / a(k,k);
    else
      wk = 1.0;
      wkm = 1.0;
    end

    if ( k + 1 <= n )

      for j = k + 1 : n
        sm = sm + abs ( z(j) + wkm * a(k,j) );
        z(j) = z(j) + wk * a(k,j);
        s = s + abs ( z(j) );
      end

      if ( s < sm )
        t = wkm - wk;
        wk = wkm;
        z(k+1:n) = z(k+1:n) + t * a(k,k+1:n);
      end

    end

    z(k) = wk;

  end

  t = sum ( abs ( z(1:n) ) );
  z(1:n) = z(1:n) / t;
%
%  Solve L' * Y = W
%
  for k = n : -1 : 1

    z(k) = z(k) + a(k+1:n,k)' * z(k+1:n)';

    t = abs ( z(k) );

    if ( 1.0E+00 < t )
      z(1:n) = z(1:n) / t;
    end

    l = pivot(k);

    t = z(l);
    z(l) = z(k);
    z(k) = t;

  end

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) );

  ynorm = 1.0E+00;
%
%  Solve L * V = Y.
%
  for k = 1 : n

    l = pivot(k);

    t = z(l);
    z(l) = z(k);
    z(k) = t;

    z(k+1:n) = z(k+1:n) + t * a(k+1:n,k)';

    if ( 1.0E+00 < abs ( z(k) ) )
      ynorm = ynorm / abs ( z(k) );
      z(1:n) = z(1:n) / abs ( z(k) );
    end

  end

  s = sum ( abs ( z(1:n) ) );
  z(1:n) = z(1:n) / s;
  ynorm = ynorm / s;
%
%  Solve U * Z = V.
%
  for k = n : -1 : 1

    if ( abs ( a(k,k) ) < abs ( z(k) ) )
      s = abs ( a(k,k) ) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    if ( a(k,k) ~= 0.0E+00 )
      z(k) = z(k) / a(k,k);
    else
      z(k) = 1.0E+00;
    end

    z(1:k-1) = z(1:k-1) - z(k) * a(1:k-1,k)';

  end
%
%  Normalize Z in the L1 norm.
%
  s = 1.0E+00 / sum ( abs ( z(1:n) ) );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0E+00 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0E+00;
  end

  return
end


text:
function c = acorr2(x,maxlag)
%ACORR Estimate autocorrelation function of time series
%
%   C = ACORR(X,MAXLAG) returns normalized autocorrelation
%   sequences for each column of X computing correlation via FFT
%

% Copyright (C) 2000-2008 Aki Vehtari
%
% This software is distributed under the GNU General Public
% Licence (version 3 or later); please refer to the file
% Licence.txt, included with the software, for details.

if nargin < 1
  error('Not enough input arguments.');
end
if nargin < 2
  maxlag=length(x)-1;
end
[m,n]=size(x);
c=zeros(maxlag,n);
for i1=1:n
  xn = x(:,i1)-mean(x(:,i1));
  xf = fft(xn,2^nextpow2(2*m-1));
  cf = ifft(abs(xf).^2);
  c(:,i1) = cf(2:(maxlag+1))./cf(1);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% This file will generate the guess of the SINDy-PI left hand side.
% Last Updated: 2019/04/21
% Coded By: K
function [P_Data,P_sym]=GuessLib(X,dX,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order)
%% First get the size of the X matrix, determin the data length and the number of variables we have.
[Data_Length,Variable_Number]=size(X);
[~,Variable_Number_dX]=size(dX);
[~,Variable_Number_u]=size(u);
% Also create the symbolic variable
Symbol=sym('z',[Variable_Number,1]);
Symbol_dX=sym('dz',[Variable_Number,1]);
Symbol_u=sym('u',[Variable_Number_u,1]);

%% Now according the Highest Polynomial Order entered, we will calculate the data matrix.
Data=[];
Index=1;

%% First calculate the polynomial term

%Order Zero:
Data(:,Index)=ones(Data_Length,1);
Sym_Struct{1,Index}=sym(1);

%Order One:
if Highest_Poly_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=X(:,i);
        Sym_Struct{1,Index}=Symbol(i,1);
    end
end

%Order Two:
if Highest_Poly_Order>=2
    for i=1:Variable_Number
        for j=1:Variable_Number
            Index=Index+1;
            Data(:,Index)=X(:,i).*X(:,j);
            Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1);
        end
    end
end

%Order Three:
if Highest_Poly_Order>=3
    for i=1:Variable_Number
        for j=1:Variable_Number
            for k=1:Variable_Number
                Index=Index+1;
                Data(:,Index)=X(:,i).*X(:,j).*X(:,k);
                Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1);
            end
        end
    end
end

%Order Four:
if Highest_Poly_Order>=4
    for i=1:Variable_Number
        for j=1:Variable_Number
            for k=1:Variable_Number
                for p=1:Variable_Number
                    Index=Index+1;
                    Data(:,Index)=X(:,i).*X(:,j).*X(:,k).*X(:,p);
                    Sym_Struct{1,Index}=Symbol(i,1)*Symbol(j,1)*Symbol(k,1)*Symbol(p,1);
                end
            end
        end
    end
end

%% Then add the Trigonometric Function in the output data:

%Order One:
if Highest_Trig_Order>=1
    for i=1:Variable_Number
        Index=Index+1;
        Data(:,Index)=sin(X(:,i));
        Sym_Struct{1,Index}=sin(Symbol(i,1));
    end
end

pin=Index;
j=0;

%% Frome here, we add the right hand side to our data
%If the Highest_U_Order is zero, we assume that our system does not have
%control input, vice versa.
if Highest_U_Order==0
    for k=1:pin
        j=j+1;
        P_Data(:,j)=dX(:,1).*Data(:,k);
        P_sym{1,j}=Symbol_dX(1,1)*(Sym_Struct{1,k});
    end
else
    for i=1:Variable_Number_u
        for k=1:pin
            j=j+1;
            P_Data(:,j)=u(:,1).*Data(:,k);
            P_sym{1,j}=Symbol_u(1,1)*(Sym_Struct{1,k});
        end
    end
end

author: dynamicslab

repo: SINDy-PI


text:
function [a,varargout] = Insets(location,s,varargin)

%Insets - Create insets in current axes.
%
%  Create insets (= small figures) in current axes. Insets are laid out side
%  by side from left to right. The number of insets is determined by the
%  number of output parameters.
%
%  USAGE
%
%    [a,b,...] = Insets(location,size,<options>)
%
%    location       'topleft', 'top', 'topright', 'left', 'center', 'right',
%                   'bottomleft', 'bottom', or 'bottomright'
%    size           total [h v] size of insets, as proportions of current axes
%                   (after subtracting margins)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'margin'      margin around insets, as a proportion of figure size
%                   (default = 0.1)
%    =========================================================================
%
%  EXAMPLE
%
%    % Create a figure and add the main plot
%    figure;plot(x,y);
%    % Create 3 insets in the top right corner, leaving 10% of the available
%    % width and height for margins (default). In total, the insets will occupy
%    % half of the remaining width, and one quarter of the remaining height.
%    [a,b,c] = Insets('topright',[0.5 0.25]);
%    % Plot data in each of the insets
%    axes(a);plot(x,z1);
%    axes(b);plot(x,z2);
%    axes(c);plot(x,z3);

% Copyright (C) 2009-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Default values
margin = 0.1;

% Check number of parameters
if nargin < 2,
  error('Incorrect number of parameters (type ''help Insets'' for details).');
end

% Check parameters
if any(s>1|s<0),
  error('Sizes must be expressed as proportion [0..1] of total figure size (type ''help Insets'' for details).');
end
location = lower(location);
if ~isstring_FMAT(location,'topleft','top','topright','left','center','right','bottomleft','bottom','bottomright'),
  error('Incorrect location (type ''help Insets'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help Insets'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'margin',
			margin = lower(varargin{i+1});
			if ~isdscalar(margin,'>=0','<=1'),
				error('Incorrect value for property ''margin'' (type ''help Insets'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help Insets'' for details).']);
	end
end

% Position and size of current axes
p = get(gca,'position');
axesBottom = p(2);
axesTop = p(2)+p(4);
axesLeft = p(1);
axesRight = p(1)+p(3);
axesWidth = p(3);
axesHeight = p(4);

% Inset margins
hMargin = margin*axesWidth;
vMargin = margin*axesHeight;

% Size of each mini figure
n = nargout;
dx = s(1)/n*(axesWidth-2*hMargin);
dy = s(2)*(axesHeight-2*vMargin);

% Start point
switch(lower(location)),
	case 'topleft',
		x0 = axesLeft+hMargin;
		y0 = axesTop-dy-vMargin;
	case 'top',
		x0 = axesLeft+hMargin+dx*n/2;
		y0 = axesTop-dy-vMargin;
	case 'topright',
		x0 = axesRight-hMargin-dx*n;
		y0 = axesTop-dy-vMargin;
	case 'left',
		x0 = axesLeft+hMargin;
		y0 = axesBottom+vMargin+dy/2;
	case 'center',
		x0 = axesLeft+hMargin+dx*n/2;
		y0 = axesBottom+vMargin+dy/2;
	case 'right',
		x0 = axesRight-hMargin-dx*n;
		y0 = axesBottom+vMargin+dy/2;
	case 'bottomleft',
		x0 = axesLeft+hMargin;
		y0 = axesBottom+vMargin;
	case 'bottom',
		x0 = axesLeft+hMargin+dx*n/2;
		y0 = axesBottom+vMargin;
	case 'bottomright',
		x0 = axesRight-hMargin-dx*n;
		y0 = axesBottom+vMargin;
end

% Create all axes
for i = 1:n,
	varargout{i} = axes('position',[x0+(i-1)*dx y0 dx dy]);
end

% Matlab cannot return only varargout, it needs at least one explicit parameter...
a = varargout{1};
varargout = {varargout{2:end}};

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Flush out the MATLAB.
close all;
clc;
clear all;

% Read the desired image file.
ImageData=imread('BC.jpg');

% Calculate the number of pixels.
NoOfPixel=size(ImageData,1)*size(ImageData,2);

% Display the original image.
figure,imshow(ImageData);
title(' Original Image without Histogram: ');

% Histogram image implementation.
HistogramImage=uint8(zeros(size(ImageData,1),size(ImageData,2)));

Frequency=zeros(256,1);

ProbabilityFrequency=zeros(256,1);

ProbabilityC=zeros(256,1);

Cumulative=zeros(256,1);

OutputZeros=zeros(256,1);

% Cumulative histogram for each and every pixel og image.
for a=1:size(ImageData,1)

    for b=1:size(ImageData,2)
        FinalValue=ImageData(a,b);

        Frequency(FinalValue+1)=Frequency(FinalValue+1)+1;

        ProbabilityFrequency(FinalValue+1)=Frequency(FinalValue+1)/NoOfPixel;
    end

end

FinalSum=0;

NumberBin=255;

% Distributation probability for each and every pixel.
for a=1:size(ProbabilityFrequency)
   FinalSum=FinalSum+Frequency(a);

   Cumulative(a)=FinalSum;

   ProbabilityC(a)=Cumulative(a)/NoOfPixel;

   OutputZeros(a)=round(ProbabilityC(a)*NumberBin);
end

% Final calculation.
for a=1:size(ImageData,1)

    for b=1:size(ImageData,2)
            HistogramImage(a,b)=OutputZeros(ImageData(a,b)+1);
    end

end

% Display the output image.
figure,imshow(HistogramImage);
title(' Final Image with Histogram: ');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function spatial_dist_extrema = spatial_distance_extrema_light(virtual_topography, virtual_chanlocs)
% spatial_dist_extrema calculates the log of the 2-norm of the distance
% between the two electrodes with the highest and lowest activations
% as described in
% Automatic Classification of Artifactual ICA-Components for Artifact
% Removal in EEG Signals by Irene Winkler, Stefan Haufe and Michael
% Tangermann http://www.behavioralandbrainfunctions.com/content/7/1/30
%
% Input:
% eeg: EEGLab data structure with additional fields
% virtual_topography and virtual_chanlocs. See the readme file for the
% toolbox for details on these fields.
%
% Output:
% spatial_dist_extrema: A vector containing the log of the 2-norm of the
% distance between the electrodes with highest and lowest activation in the
% scalp map of each IC in the EEGLab data structure given as input.

% Copyright (C) 2013  Laura Froelich (laura.frolich@gmail.com)
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

nics = size(virtual_topography, 1);
spatial_dist_extrema = NaN(nics,1);
for ic =1:nics
    channel_activation = virtual_topography(ic,:);
    [dummy, maxindex] = max(channel_activation);
    [dummy, minindex] = min(channel_activation);
    maxlocation = [virtual_chanlocs(maxindex).X virtual_chanlocs(maxindex).Y virtual_chanlocs(maxindex).Z];
    minlocation = [virtual_chanlocs(minindex).X virtual_chanlocs(minindex).Y virtual_chanlocs(minindex).Z];

    spatial_dist_extrema(ic) = log(norm(maxlocation-minlocation,2));
end

end


text:
% SYNTAX:
% mlActAuto = hmrR_PruneChannels(data, probe, mlActMan, tInc, dRange, SNRthresh, SDrange)
%
% UI NAME:
% Prune_Channels
%
% DESCRIPTION:
% Prune channels from the measurement list if their signal is too weak, too
% strong, or their standard deviation is too great. This function
% updates MeasListAct based on whether data 'd' meets these conditions
% as specified by'dRange' and 'SNRthresh'.
%
% INPUTS:
% d - SNIRF object containing time course data (nTpts x nChannels )
% probe - SNIRF object describing the probe - optode positions and wavelengths.
% mlActMan -
% dRange - if mean(d) < dRange(1) or > dRange(2) then it is excluded as an
%      active channel
% SNRthresh - if mean(d)/std(d) < SNRthresh then it is excluded as an
%      active channel
% SDrange - will prune channels with a source-detector separation <
%           SDrange(1) or > SDrange(2)
%
% OUTPUTS:
% mlAct - cell array of all data blocks - each data block is an array
%         of true/false for all channels in the contanining data block
%         specifying active/inactive status. (# of data blocks x # of Channels)
%
% USAGE OPTIONS:
% Prune_Channels: mlActAuto = hmrR_PruneChannels(data, probe, mlActMan, tIncMan, dRange, SNRthresh, SDrange)
%
% PARAMETERS:
% dRange: [1e4, 1e7]
% SNRthresh: 2
% SDrange: [0.0, 45.0]
%
% TO DO:
% consider Conc as well as wavelength data
%
function mlActAuto = hmrR_PruneChannels(data, probe, mlActMan, tIncMan, dRange, SNRthresh, SDrange)

% Init output
mlActAuto = cell(length(data),1);

% Check input args
if nargin<7
    disp( 'USAGE: hmrR_PruneChannels(data, probe, mlActMan, tIncMan, dRange, SNRthresh, SDrange)' )
    return
end
if isempty(tIncMan)
    tIncMan = cell(length(data),1);
end
if isempty(mlActMan)
    mlActMan = cell(length(data),1);
end

for iBlk = 1:length(data)

    d        = data(iBlk).GetDataTimeSeries();
    t        = data(iBlk).GetTime();
    MeasList = data(iBlk).GetMeasList();
    Lambda   = probe.GetWls();
    SrcPos   = probe.GetSrcPos();
    DetPos   = probe.GetDetPos();

    mlActMan{iBlk} = mlAct_Initialize(mlActMan{iBlk}, MeasList);

    if isempty(tIncMan{iBlk})
        tIncMan{iBlk} = ones(length(t),1);
    end
    tInc = tIncMan{iBlk};

    lstInc = find(tInc==1);
    d = d(lstInc,:);

    % check for dRange and SNRthresh
    dmean = mean(d,1);
    dstd = std(d,[],1);

    idxs1 = [];
    idxs2 = [];
    idxs3 = [];
    nLambda = length(Lambda);
    lst1 = find(MeasList(:,4)==1);

    % Start by including all channels
    chanList = ones(size(MeasList,1),1);

    for ii = 1:nLambda
        lst = [];
        rhoSD = [];
        for jj = 1:length(lst1)
            lst(jj) = find(MeasList(:,1)==MeasList(lst1(jj),1) & ...
                           MeasList(:,2)==MeasList(lst1(jj),2) & ...
                           MeasList(:,4)==ii);
            rhoSD(jj) = norm( SrcPos(MeasList(lst1(jj),1),:) - DetPos(MeasList(lst1(jj),2),:) );
        end

        % dRange exclusion criteria
        idxs1 = [idxs1, find(dmean(lst)<=dRange(1) | dmean(lst)>=dRange(2))];

        % SNRthresh exclusion criteria
        idxs2 = [idxs2, find((dmean(lst)./dstd(lst)) <= SNRthresh)];

        % SDrange exclusion criteria
        idxs3 = [idxs3, find(rhoSD<SDrange(1) | rhoSD>SDrange(2))];

        idxsExcl = unique([idxs1(:)', idxs2(:)', idxs3(:)']);

        chanList(lst(idxsExcl)) = 0;
    end

    idxs4 = find(mlActMan{iBlk}(:,3) == 0);

    if ~isempty(idxs1)
        fprintf('hmrR_PruneChannels:  excluded channels [ %s ] based on dRange=[ %s ] at wavelength %d\n', num2str(idxs1(:)'), num2str(dRange), ii);
    end
    if ~isempty(idxs2)
        fprintf('hmrR_PruneChannels:  excluded channels [ %s ] based on SNRthresh=%0.1f at wavelength %d\n', num2str(idxs2(:)'), SNRthresh, ii);
    end
    if ~isempty(idxs3)
        fprintf('hmrR_PruneChannels:  excluded channels [ %s ] based on SDrange=[ %s ] at wavelength %d\n', num2str(idxs3(:)'), num2str(SDrange), ii);
    end
    if ~isempty(idxs4)
        fprintf('hmrR_PruneChannels:  manually excluded channels [ %s ]\n', num2str(idxs4(:)'));
    end

    chanList = chanList(:) & mlActMan{iBlk}(:,3);

    % update MeasListAct
    mlActAuto{iBlk} = mlAct_Initialize(chanList, MeasList);
end

author: BUNPC


text:
function sys = assignschur(AConstraint,thecompiler,varargin)

sys = sdpvar(1);
sys.typeflag = 30;
Schurfun_defined = 0;
Var_defined = 0;

for i = 1:nargin
    if isa(varargin{i},'sdpvar')
        sys.extra.variables{i} = varargin{i};
        Var_defined = 1;
    elseif isa(varargin{i},'double')
        sys.extra.data{i} = varargin{i};
    elseif isa(varargin{i},'char')
        sys.extra.Schurfun = varargin{i};
        Schurfun_defined = 1;
    end
end
sys.extra.negated = 0;

if ~Schurfun_defined
    error('A Schur-function must be defined.')
end

sys = constraint(sys,'>=',0);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function h = addAsterisks(Y,E,sel,varargin);
% h = addAsterisks(Y,E,sel,[plot options]);
%
% Creates an errorbar plot, with asterisks over selected
% points. Y is a matrix of column vectors containing the
% data to be plotted, E is a matrix the same size as Y containing
% error lengths for each point. (In this respect, it is the same
% as calling ERRORBAR(Y,E).) sel can either be:
%      1) a matrix the same size as Y and E, with a 1 at each
%         location which will get an asterisk and a 0 otherwise, or
%
%      2) a vector equal to the # of rows or columns, selecting entire
%         rows or columns for asterisks.
%
% Any fourth or later input arguments are passed on to the ERRORBAR command
% (e.g., formatting arguments -- see HELP PLOT for more details.), except
% for the following:
%
% 'X',[vals]: provide the values to plot on the X-axis as the next arg.
%
% 'leg',{'labels'}: add a legend using the labels in the next argument
% (as a cell-of-strings)

% 03/03 ras
% 01/04 ras: fixed bug in which asterisks were never centered
if ~exist('sel','var')  | isempty(sel)
    sel = zeros(size(Y));
end

% figure out if entire rows/colums are selected for asterisks
if length(sel)==size(Y,1) & isequal(unique(sel),[0 1])
    sel = repmat(sel,size(Y,2),1)';
elseif length(sel)==size(Y,2) & isequal(unique(sel),[0 1])
    sel = repmat(sel,size(Y,1),1);
end

%%%%% defaults
X = repmat(1:size(Y,1),size(Y,2),1)';
plotOptions = {};

%%%%% parse the option flags
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch lower(varargin{i})
            case {'leg','legend'}
                leg = varargin{i+1};
            case 'x',
                X = varargin{i+1};
                if size(X,1)==1
                    X = repmat(X,size(Y,2),1)';
                end
            otherwise,
                plotOptions{end+1} = varargin{i};
        end
    end
end

if isempty(plotOptions)
    h = errorbar(X,Y,E);
else
    cmd = ['h = errorbar(X,Y,E'];
    for i = 1:length(plotOptions)
        cmd = [cmd ',''' plotOptions{i} ''''];
    end
    cmd = [cmd ');'];
    eval(cmd);
end
AX = axis;
xSz = AX(2) - AX(1);
ySz = AX(4) - AX(3);

whichPoints = find(sel);

for i = 1:length(whichPoints)
    pt = whichPoints(i);
    xLoc = X(pt);
    yLoc = Y(pt) + E(pt) + 0.05*ySz;
    h2 = text(xLoc,yLoc,'*','FontSize',24,'HorizontalAlignment','center');
    %set(h2,'HorizontalAlignment','center');
    %set(h2,'Fontsize',20);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function xhpf = HPF(u)
%
%
persistent prevX
persistent prevU
persistent dt tau
persistent firstRun

if isempty(firstRun)
  prevX = 0;
  prevU = 0;
  dt    = 0.01;
  tau   = 0.0233;

  firstRun = 1;
end

alpha = tau / (tau + dt);
xhpf = alpha*prevX + alpha*(u - prevU);

prevX = xhpf;
prevU = u;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function fx = jdfunc (jdin)

% objective function for tdb2utc

% input

%  jdin = current value for UTC julian date

% output

%  fx = delta julian date

% Orbital Mechanics with MATLAB

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

global jdsaved

tai_utc = findleap(jdin);

fx = utc2tdb (jdin, tai_utc) - jdsaved;

end

author: Sable

repo: mcbench-benchmarks


text:
function [min] = ns2min(ns)
% Convert time from nanoseconds to minutes.
% Chad Greene 2012
min = ns/60*1e-9 ;

% Dear Postmaster General,
%
% The following is a list of individuals who do not live in my house:
%
%     Current Resident
%     Smart Shopper
%     Valued Customer
%     Our Neighbor
%     Deal Seeker
%     Postal Patron
%     Thrifty Consumer
%     Box Holder
%     Home Owner
%     Tax Payer
%     Texas Voter
%     Pet Lover
%     Bargain Buyer
%     Sports Fan
%     Safe Driver
%     Future Client
%     Home Chef
%     Steak Lover
%     Darryl Jenkins
%
% Please refrain from delivering their mail to my address.


text:
function x = ncc_compute_points ( n )

%*****************************************************************************80
%
%% NCC_COMPUTE_POINTS: points of a Newton-Cotes Closed quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
  x_min = -1.0;
  x_max =  1.0;

  if ( n == 1 )

    x(1) = ( x_max + x_min ) / 2.0;

  else

    for i = 1 : n
      x(i) = ( ( n - i     ) * x_min   ...
             + (     i - 1 ) * x_max ) ...
             / ( n     - 1 );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c = diag(a,k)
%DIAG         Implements  diag(a,k)  for hessians
%
%   c = diag(a,k)
%
% functionality as Matlab function diag for matrices
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_HESSIAN_NUMVAR = getappdata(0,'INTLAB_HESSIAN_NUMVAR');

  if nargin==1
    k = 0;
  end

  [m n] = size(a.x);
  c.x = diag(a.x,k);
  if m==1 | n==1
    index = diag(ones(1,m*n),k);
    if issparse(a.hx)
      c.dx = sparse([],[],[],INTLAB_HESSIAN_NUMVAR,prod(size(index)),0);
      c.hx = sparse([],[],[],INTLAB_HESSIAN_NUMVAR^2,prod(size(index)),0);
    else
      c.dx = zeros(INTLAB_HESSIAN_NUMVAR,prod(size(index)));
      c.hx = zeros(INTLAB_HESSIAN_NUMVAR^2,prod(size(index)));
    end
    if isa(a.x,'intval')
      c.dx = intval(c.dx);
      c.hx = intval(c.hx);
    end
    c.dx(:,index~=0) = a.dx;
    c.hx(:,index~=0) = a.hx;
  else
    index = diag( reshape( 1:(m*n) , m , n ) , k );
    c.dx = a.dx(:,index(:));
    c.hx = a.hx(:,index(:));
  end

  c = class(c,'hessian');

repo: OpenMAS


text:
%RECOVERPOSE  Recover relative camera rotation and translation from an estimated essential matrix and the corresponding points in two images, using cheirality check
%
%     [R, t, good] = cv.recoverPose(E, points1, points2)
%     [R, t, good, mask, triangulatedPoints] = cv.recoverPose(...)
%     [...] = cv.recoverPose(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __E__ The input essential matrix, 3x3.
% * __points1__ Cell array of N 2D points from the first image, or numeric
%   array Nx2/Nx1x2/1xNx2. The point coordinates should be floating-point
%   (single or double precision).
% * __points2__ Cell array or numeric array of the second image points of the
%   same size and format as `points1`.
%
% ## Output
% * __R__ Recovered relative rotation, 3x3 matrix.
% * __t__ Recovered relative translation, 3x1 vector.
% * __good__ the number of inliers which pass the cheirality check.
% * __mask__ Output mask for inliers in `points1` and `points2`. In the output
%   mask only inliers which pass the cheirality check. Vector of length N, see
%   the `Mask` input option.
% * __triangulatedPoints__ 3D points which were reconstructed by triangulation,
%   see cv.triangulatePoints
%
% ## Options
% * __CameraMatrix__ Camera matrix `K = [fx 0 cx; 0 fy cy; 0 0 1]`. Note that
%   this function assumes that `points1` and `points2` are feature points from
%   cameras with the same camera matrix. default `eye(3)`.
% * __DistanceThreshold__ threshold distance which is used to filter out far
%   away points (i.e. infinite points). default 50.0
% * __Mask__ Input mask of length N for inliers in `points1` and `points2`
%   (0 for outliers and to 1 for the other points (inliers). If it is not
%   empty, then it marks inliers in `points1` and `points2` for then given
%   essential matrix `E`. Only these inliers will be used to recover pose.
%   Not set by default.
%
% This function decomposes an essential matrix using cv.decomposeEssentialMat
% and then verifies possible pose hypotheses by doing cheirality check. The
% cheirality check basically means that the triangulated 3D points should have
% positive depth. Some details can be found in [Nister03].
%
% This function can be used to process output `E` and `mask` from
% cv.findEssentialMat. In this scenario, `points1` and `points2` are the same
% input for cv.findEssentialMat.
%
% ## Example
%
%     % Estimation of fundamental matrix using the RANSAC algorithm
%     point_count = 100;
%     points1 = cell(1, point_count);
%     points2 = cell(1, point_count);
%     % initialize the points here ...
%     for i=1:point_count
%         points1{i} = ...;  % [x,y]
%         points2{i} = ...;  % [x,y]
%     end
%
%     % cametra matrix with both focal lengths = 1, and principal point = [0 0]
%     cameraMatrix = eye(3,3);
%
%     [E, mask] = cv.findEssentialMat(points1, points2, ...
%         'CameraMatrix',cameraMatrix, 'Method','Ransac');
%     [R, t, ~, mask] = cv.recoverPose(E, points1, points2, ...
%         'CameraMatrix',cameraMatrix, 'Mask',mask);
%
% ## References
% [Nister03]:
% > David Nister. "An efficient solution to the five-point relative pose
% > problem". Pattern Analysis and Machine Intelligence, IEEE Transactions on,
% > 26(6):756-770, 2004.
%
% See also: cv.findEssentialMat, cv.decomposeEssentialMat,
%  cv.triangulatePoints, relativeCameraPose
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [a,C] = biconnected_components(A,varargin)
% BICONNECTED_COMPONENTS Compute the biconnected components and
% articulation points for a symmetric graph A.
%
% [a C] = biconnected_components(A) returns a list of articulation points
% a and the component graph C where each non-zero indicates the connected
% component of the edge.  That is, C is a matrix with the same non-zero
% structure as A, but with the values replaced with the index of the
% biconnected component of that edge.  The vector a is a list of
% articulation points in the graph.  Articulation points are vertices that
% belong to more than one biconnected component.  Removing an articulation
% point disconnects the graph.
%
% If C is not requested, it is not built.
%
% This method works on undirected graphs.
% The runtime is O(V+E), the algorithm is just depth first search.
%
% ... = biconnected_components(A,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   There are no additional options for this function.
%
% Note: the input to this function must be symmetric, so this function
% ignores the 'notrans' default option and never transposes the input.
%
% Note: this function does not depend upon the non-zero values of A, but
% only uses the non-zero structure of A.
%
% Example:
%    load graphs/tarjan-biconn.mat
%    biconnected_components(A)
%
% See also COMPONENTS

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

%% History
%  2006-04-19: Initial version
%  2006-05-31: Added full2sparse check
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

if trans, end

if check
    % make sure the matrix is symmetric
    check_matlab_bgl(A,struct('sym',1));
end;

% the graph has to be symmetric, so trans doesn't matter.

if (nargout > 1)
    [a ci] = biconnected_components_mex(A);

    % convert the indices from the graph back into a new matrix.
    [i j] = find(A);
    C = sparse(i,j,ci,size(A,1),size(A,1));

    C = max(C,C');
else
    a = biconnected_components_mex(A);
end;

% 0 a indicates it isn't an articulation point.
a = a(a > 0);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [net, info] = cnn_imagenet(varargin)
%CNN_IMAGENET   Demonstrates training a CNN on ImageNet
%  This demo demonstrates training the AlexNet, VGG-F, VGG-S, VGG-M,
%  VGG-VD-16, and VGG-VD-19 architectures on ImageNet data.

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

% opts.dataDir = fullfile(vl_rootnn, 'data','ILSVRC2012') ;
opts.dataDir = fullfile('/data','jiangqy','ImageNet');
opts.modelType = 'alexnet' ;
opts.networkType = 'simplenn' ;
opts.batchNormalization = false ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.modelType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
sfx = [sfx '-' opts.networkType] ;
% opts.expDir = fullfile(vl_rootnn, 'data', ['imagenet12-' sfx]) ;
opts.expDir = fullfile('/data','jiangqy','ImageNet');
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.numFetchThreads = 12 ;
opts.lite = false ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = [12]; end;

% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------

net = cnn_imagenet_init('model', opts.modelType, ...
                        'batchNormalization', opts.batchNormalization, ...
                        'weightInitMethod', opts.weightInitMethod, ...
                        'networkType', opts.networkType) ;

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
else
  imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% Set the class names in the network
net.meta.classes.name = imdb.classes.name ;
net.meta.classes.description = imdb.classes.description ;

% Compute image statistics (mean, RGB covariances, etc.)
imageStatsPath = fullfile(opts.expDir, 'imageStats.mat') ;
if exist(imageStatsPath)
  load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else
  [averageImage, rgbMean, rgbCovariance] = getImageStats(opts, net.meta, imdb) ;
  save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end

% Set the image average (use either an image or a color)
%net.meta.normalization.averageImage = averageImage ;
net.meta.normalization.averageImage = rgbMean ;

% Set data augmentation statistics
[v,d] = eig(rgbCovariance) ;
net.meta.augmentation.rgbVariance = 0.1*sqrt(d)*v' ;
clear v d ;

% -------------------------------------------------------------------------
%                                                                     Learn
% -------------------------------------------------------------------------

% switch opts.networkType
%   case 'simplenn', trainFn = @cnn_train ;
%   case 'dagnn', trainFn = @cnn_train_dag ;
% end

[net, info] = cnn_train(net, imdb, getBatchFn(opts, net.meta), ...
                      'expDir', opts.expDir, ...
                      net.meta.trainOpts, ...
                      opts.train) ;

% -------------------------------------------------------------------------
%                                                                    Deploy
% -------------------------------------------------------------------------

net = cnn_imagenet_deploy(net) ;
modelPath = fullfile(opts.expDir, 'net-deployed.mat')

switch opts.networkType
  case 'simplenn'
    save(modelPath, '-struct', 'net') ;
  case 'dagnn'
    net_ = net.saveobj() ;
    save(modelPath, '-struct', 'net_') ;
    clear net_ ;
end

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------
useGpu = numel(opts.train.gpus) > 0 ;

bopts.numThreads = opts.numFetchThreads ;
bopts.imageSize = meta.normalization.imageSize ;
bopts.border = meta.normalization.border ;
bopts.averageImage = meta.normalization.averageImage ;
bopts.rgbVariance = meta.augmentation.rgbVariance ;
bopts.transformation = meta.augmentation.transformation ;

switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(bopts,x,y) ;
  case 'dagnn'
    fn = @(x,y) getDagNNBatch(bopts,useGpu,x,y) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getSimpleNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal
  % training
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0) ;
else
  % validation: disable data augmentation
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0, ...
                              'transformation', 'none') ;
end

if nargout > 0
  labels = imdb.images.label(batch) ;
end

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, useGpu, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
isVal = ~isempty(batch) && imdb.images.set(batch(1)) ~= 1 ;

if ~isVal
  % training
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0) ;
else
  % validation: disable data augmentation
  im = cnn_imagenet_get_batch(images, opts, ...
                              'prefetch', nargout == 0, ...
                              'transformation', 'none') ;
end

if nargout > 0
  if useGpu
    im = gpuArray(im) ;
  end
  labels = imdb.images.label(batch) ;
  inputs = {'input', im, 'label', labels} ;
end

% -------------------------------------------------------------------------
function [averageImage, rgbMean, rgbCovariance] = getImageStats(opts, meta, imdb)
% -------------------------------------------------------------------------
train = find(imdb.images.set == 1) ;
train = train(1: 101: end);
bs = 256 ;
opts.networkType = 'simplenn' ;
fn = getBatchFn(opts, meta) ;
avg = {}; rgbm1 = {}; rgbm2 = {};

for t=1:bs:numel(train)
  batch_time = tic ;
  batch = train(t:min(t+bs-1, numel(train))) ;
  fprintf('collecting image stats: batch starting with image %d ...', batch(1)) ;
  temp = fn(imdb, batch) ;
  z = reshape(permute(temp,[3 1 2 4]),3,[]) ;
  n = size(z,2) ;
  avg{end+1} = mean(temp, 4) ;
  rgbm1{end+1} = sum(z,2)/n ;
  rgbm2{end+1} = z*z'/n ;
  batch_time = toc(batch_time) ;
  fprintf(' %.2f s (%.1f images/s)\n', batch_time, numel(batch)/ batch_time) ;
end
averageImage = mean(cat(4,avg{:}),4) ;
rgbm1 = mean(cat(2,rgbm1{:}),2) ;
rgbm2 = mean(cat(3,rgbm2{:}),3) ;
rgbMean = rgbm1 ;
rgbCovariance = rgbm2 - rgbm1*rgbm1' ;

author: jiangqy


text:
function Y=repmat(varargin)
%REPMAT (overloaded)

Y = varargin{1};
dim = Y.dim;
resize = [varargin{2:end}];
index = (1:size(Y.basis,1))';
index = reshape(index,Y.dim);
index = repmat(index,resize);
Y.basis = Y.basis(index(:),:);
Y.dim = size(index);
% Reset info about conic terms
Y.conicinfo = [0 0];


text:
function [mFMDC, mFMD] = calc_FMDMag(mCatalog)
% function [mFMDC, mFMD] = calc_FMD(mCatalog)
% -------------------------------------------
% Calculates the cumulative and non-cumulative frequency magnitude distribution
%   for given magnitudes of an earthquake catalog
%
% Input parameter:
%   vMagnitudes    Either EQ catalog in ZMAP format or vector of magnitudes of an earthquake catalog
%
% Output parameters:
%   mFMDC       cumulative frequency magnitude distribution
%               mFMDC(1,:) = magnitudes (x-axis)
%               mFMDC(2,:) = number of events (y-axis)
%   mFMD        non-cumulative frequency magnitude distribution
%
% J. Woessner; woessner@seismo.ifg.ethz.ch
% last update: 12.07.02

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

[nXSize, nYSize] = size(mCatalog);
if nXSize > 1
  % Use magnitude column from ZMAP data catalog format
  vMagnitudes = mCatalog(:,6);
else
  % Use one column magnitude vector
  vMagnitudes = mCatalog;
end

% Determine the magnitude range
fMaxMagnitude = ceil(10 * max(vMagnitudes)) / 10;
fMinMagnitude = floor(min(vMagnitudes));
if fMinMagnitude > 0
  fMinMagnitude = 0;
end

% Naming convention:
%   xxxxR : Reverse order
%   xxxxC : Cumulative number

% Do the calculation
[vNumberEvents] = hist(vMagnitudes, (fMinMagnitude:0.1:fMaxMagnitude));
vNumberEventsR  = vNumberEvents(length(vNumberEvents):-1:1);
vNumberEventsCR = cumsum(vNumberEvents(length(vNumberEvents):-1:1));

% Create the x-axis values
vXAxis = (fMaxMagnitude:-0.1:fMinMagnitude);

% Merge the x-axis values with the FMDs and return them
mFMD  = [vXAxis; vNumberEventsR];
mFMDC = [vXAxis; vNumberEventsCR];

author: CelsoReyes


text:
function varargout=potsband(varargin)
%
% For calling details please see v_potsband.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_potsband.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: potsband.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_potsband(varargin{:});
else
    v_potsband(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
  function hh = fig_text(varargin)
%|function hh = fig_text(...)
%| add text to figure, using coordinates relative to entire window
%| options:
%|	'-date'		prepend date
%|	'-tex'		tex parse string
%|	x,y		coordinates relative to [0 0 1 1]
%|	{args}		style arguments for text() command
%| jeff fessler

comment = '';
thedate = '';
interpret = 'none';
textarg = {'fontsize', 11};
x = 0.01;
y = 0.01;
while length(varargin)
	arg = varargin{1};
	if ischar(arg)
		if streq(arg, '-date')
			thedate = [mydate ' '];
		elseif streq(arg, '-tex')
			interpret = 'tex';
		else
			comment = arg;
		end
		varargin = {varargin{2:end}};

	elseif isnumeric(arg)
		if length(varargin) < 2, error 'need x, y', end
		x = varargin{1};
		y = varargin{2};
		varargin = {varargin{3:end}};

	elseif iscell(arg)
		textarg = {textarg{:}, arg{:}};
		varargin = {varargin{2:end}};

	else
		help(mfilename), error 'bad args'
	end
end

hh = gca;

hold on
h = axes('position', [0 0 1 1]); axis off
h = text(x, y, [thedate comment], 'interpreter', interpret, textarg{:});
hold off

if nargout
	hh = h;
else
	axes(hh) % reset axes to the one before the text added
	clear hh
end

function s = mydate
%s = datestr(now, 'yyyy-mm-dd HH:MM:SS');
%s = datestr(now, 'yy-mm-dd HH:MM:SS');
t1 = datestr(now, 'yy');
t2 = datestr(now, 'mm');
t3 = datestr(now, 'dd');
t4 = datestr(now, 'HH:MM:SS');
s = sprintf('%s-%s-%s %s', t1, t2, t3, t4);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef WrapperResFile < FileReader

    properties (Access = private)
       dataBase
       stress
       exponent
       alpha
       lineNumber
       linesEndValue
    end

    properties (Access = private)
       nElem
       nNodes
       dimension
    end

    methods (Access = public)

        function obj = WrapperResFile(cParams)
            obj.init(cParams);
        end

        function read(obj)
            obj.openFile();
            obj.readStress();
            obj.readStrain();
            obj.readDisplacement();
            obj.readStress();
            obj.readStrain();
            obj.readDisplacement();
            obj.readAmplifiedStressNormP();
            %obj.readCompliance();
            obj.readAngle();
            obj.readAngle();
            obj.readAmplifiedStressNormP();
            obj.readDensity();
            obj.readScalarInNodes();
            obj.readScalarInNodes();
            obj.readScalarInGauss();
            obj.readScalarInNodes();
            obj.readScalarInNodes();
            obj.closeFile();
        end

        function d = getDataBase(obj)
            d = obj.dataBase;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.filePath  = cParams.filePath;
            obj.nElem     = cParams.nElem;
            obj.nNodes    = cParams.nNodes;
            obj.dimension = cParams.dimension;
            obj.lineNumber = 14;
            obj.linesEndValue = 2;
        end

        function readStress(obj)
            obj.readTensorFieldInGauss()
        end

        function readStrain(obj)
            obj.readTensorFieldInGauss()
        end

        function readDisplacement(obj)
            obj.readVectorFieldInNodes();
        end

        function readCompliance(obj)
            obj.readScalarInGauss();
        end

        function readAmplifiedStressNormP(obj)
            obj.readScalarInGauss();
        end

        function readAngle(obj)
            obj.readVectorFieldInGauss();
        end

        function readDensity(obj)
            obj.readScalarInGauss();
        end

        function readLevelSet(obj)
            obj.readScalarInNodes();
        end

        function readVectorFieldInGauss(obj)
            obj.readVectorInGauss(obj.dimension)
        end

        function readVectorFieldInNodes(obj)
            obj.readVectorInNodes(obj.dimension)
        end

        function readTensorFieldInGauss(obj)
            switch obj.dimension
                case 2
                    nComp = 3;
                case 3
                    nComp = 6;
            end
            obj.readVectorInGauss(nComp);
        end

        function readScalarInGauss(obj)
            nLines = obj.nElem;
            linesJump = 2;
            nComp = 1;
            obj.readField(obj.lineNumber,nLines,linesJump,nComp);
            obj.updateLineNumber(nLines + linesJump + obj.linesEndValue);
        end

        function readVectorInGauss(obj,nComp)
            nLines = obj.nElem;
            linesJump = 3;
            obj.readField(obj.lineNumber,nLines,linesJump,nComp);
            obj.updateLineNumber(nLines + linesJump + obj.linesEndValue);
        end

        function readVectorInNodes(obj,nComp)
            nLines = obj.nNodes;
            linesJump = 3;
            obj.readField(obj.lineNumber,nLines,linesJump,nComp);
            obj.updateLineNumber(nLines + linesJump + obj.linesEndValue);
        end

        function readScalarInNodes(obj)
            nLines = obj.nNodes;
            linesJump = 2;
            nComp = 1;
            obj.readField(obj.lineNumber,nLines,linesJump,nComp);
            obj.updateLineNumber(nLines + linesJump + obj.linesEndValue);
        end

        function [fName,fValue] = readField(obj,lineNumber,nLines,linesJump,nComp)
            fName = obj.readFieldName(lineNumber);
            lineNumber = lineNumber + linesJump;
            fValue = obj.readlFieldValues(lineNumber,nLines,nComp);
            obj.dataBase.(fName) = fValue;
        end

        function name = readFieldName(obj,lineNumber)
            lineText = obj.readLines(lineNumber,1);
            lSplit = split(lineText{1});
            name = lSplit(2);
            name = obj.eraseQuationMarks(name{1});
        end

        function fValue = readlFieldValues(obj,lineNumber,nLines,nComp)
            lines = obj.readLines(lineNumber,nLines);
            namesSplit = split(lines{:});
            lines = str2double(namesSplit);
            index = (1:nComp) + 1;
            fValue = lines(:,index);
        end

        function tLine = readLines(obj,lineNumber,nLines)
            fseek(obj.fid,0,'bof');
            tLine = textscan(obj.fid,'%s',nLines,'delimiter','\n', 'headerlines',lineNumber-1);
        end

        function updateLineNumber(obj,lines)
            obj.lineNumber = obj.lineNumber + lines;
        end

    end

    methods (Access = private, Static)

        function str = eraseQuationMarks(str)
            str = erase(str,"""");
        end

    end

end

author: SwanLab