text:
function [taskReport, essentialRxns, taskStructure, essentialFluxes]=checkTasks(model,inputFile,printOutput,printOnlyFailed,getEssential,taskStructure)
% checkTasks
%   Performs a set of simulations as defined in a task file.
%
%   Input:
%   model           a model structure
%   inputFile       a task list in Excel format. See the function
%                   parseTaskList for details (opt if taskStructure is
%                   supplied)
%   printOutput     true if the results of the test should be displayed
%                   (opt, default true)
%   printOnlyFailed true if only tasks that failed should be displayed
%                   (opt, default false)
%   getEssential    true if the essential reactions should be calculated for
%                   all the tasks. This option is used with runINIT (opt,
%                   default false)
%   taskStructure   structure with the tasks, as from parseTaskList. If
%                   this is supplied then inputFile is ignored (opt)
%
%
%   Output:
%   taskReport          structure with the results
%       id              cell array with the id of the task
%       description     cell array with the description of the task
%       ok              boolean array with true if the task was successful
%   essentialRxns       MxN matrix with the essential reactions (M) for each
%                       task (N). An element is true if the corresponding
%                       reaction is essential in the corresponding task.
%                       Failed tasks and SHOULD FAIL tasks are ignored.
%                       This is used by the INIT algorithm (if tasks
%                       are supplied). If getEssential=false then
%                       essentialRxns=false(nRxns,nTasks)
%   taskStructure       structure with the tasks, as from parseTaskList
%   essentialFluxes     The fluxes of the essential rxns - same structure as essentialRxns
%
%
%   This function is used for defining a set of tasks for a model to
%   perform. The tasks are defined by defining constraints on the model,
%   and if the problem is feasible, then the task is considered successful.
%   In general, each row can contain one constraint on uptakes, one
%   constraint on outputs, one new equation, and one change of reaction
%   bounds. If more bounds are needed to define the task, then several rows
%   can be used for each task.
%
%   Usage: [taskReport, essentialReactions, taskStructure]=checkTasks(model,inputFile,...
%           printOutput,printOnlyFailed,getEssential,taskStructure)

if nargin<3 || isempty(printOutput)
    printOutput=true;
end
if nargin<4 || isempty(printOnlyFailed)
    printOnlyFailed=false;
end
if nargin<5 || isempty(getEssential)
    getEssential=false;
end

%Prepare the input model a little
model.b=zeros(numel(model.mets),2);

modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
if ~isfield(model,'unconstrained')
    EM='Exchange metabolites should normally not be removed from the model when using checkTasks. Inputs and outputs are defined in the task file instead. Use importModel(file,false) to import a model with exchange metabolites remaining';
    dispEM(EM,false);
end

%Parse the task file
if nargin<6
    taskStructure=parseTaskList(inputFile);
end

essentialRxns=false(numel(model.rxns),numel(taskStructure));
essentialFluxes = NaN(numel(model.rxns),numel(taskStructure));

tModel=model;
taskReport=[];
for i=1:numel(taskStructure)
    taskReport.id{i,1}=taskStructure(i).id;
    taskReport.description{i,1}=taskStructure(i).description;
    %Set the inputs
    if ~isempty(taskStructure(i).inputs)
        [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
        J=J(I); %Only keep the ones with matches
        K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
        L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
        %Check that all metabolites are either real metabolites or
        %ALLMETS/ALLMETSIN
        if ~all(I|K|L)
            fprintf(['ERROR: Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
            taskReport.ok(i,1)=false;
            tModel=model;
            continue;
        end
        if numel(J)~=numel(unique(J))
            EM=['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
            dispEM(EM);
        end
        %If all metabolites should be added
        if any(K)
            %Check if ALLMETS is the first metabolite. Otherwise print a
            %warning since it will write over any other constraints that
            %are set
            if K(1)==0
                EM=['ALLMETS is used as an input in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
                dispEM(EM,false);
            end
            %Use the first match of ALLMETS. There should only be one, but
            %still..
            tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
        end
        %If metabolites in a specific compartment should be used
        if any(L)
            L=find(L);
            for j=1:numel(L)
                %The compartment defined
                compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
                %Check if it exists in the model
                C=find(ismember(upper(model.comps),compartment));
                if any(C)
                    %Match to metabolites
                    tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
                else
                    EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
                    dispEM(EM);
                end
            end
        end
        %Then add the normal constraints
        if any(J)
            tModel.b(J,1)=taskStructure(i).UBin(I)*-1;
            tModel.b(J,2)=taskStructure(i).LBin(I)*-1;
        end
    end
    %Set the outputs
    if ~isempty(taskStructure(i).outputs)
        [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
        J=J(I); %Only keep the ones with matches
        K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
        L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
        %Check that all metabolites are either real metabolites or
        %ALLMETS/ALLMETSIN
        if ~all(I|K|L)
            fprintf(['ERROR: Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
            taskReport.ok(i,1)=false;
            tModel=model;
            continue;
        end
        if numel(J)~=numel(unique(J))
            EM=['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
            dispEM(EM);
        end
        %If all metabolites should be added
        if any(K)
            %Check if ALLMETS is the first metabolite. Otherwise print a
            %warning since it will write over any other constraints that
            %are set
            if K(1)==0
                EM=['ALLMETS is used as an output in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
                dispEM(EM,false);
            end
            %Use the first match of ALLMETS. There should only be one, but
            %still..
            tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
        end
        %If metabolites in a specific compartment should be used
        if any(L)
            L=find(L);
            for j=1:numel(L)
                %The compartment defined
                compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
                %Check if it exists in the model
                C=find(ismember(upper(model.comps),compartment));
                if any(C)
                    %Match to metabolites
                    tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
                else
                    EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
                    dispEM(EM);
                end
            end
        end
        %Then add the normal constraints
        if any(J)
            %Verify that IN and OUT bounds are consistent. Cannot require
            %that a metabolite is simultaneously input AND output at some
            %nonzero flux.
            I = find(I);  % otherwise indexing becomes confusing
            nonzero_LBin = tModel.b(J,2) < 0;
            nonzero_LBout = taskStructure(i).LBout(I) > 0;
            if any(nonzero_LBin & nonzero_LBout)
                EM=['The IN LB and OUT LB in "[' taskStructure(i).id '] ' taskStructure(i).description '" cannot be nonzero for the same metabolite'];
                dispEM(EM);
            end
            tModel.b(J(nonzero_LBout),1)=taskStructure(i).LBout(I(nonzero_LBout));
            tModel.b(J,2)=taskStructure(i).UBout(I);
        end
    end
    %Add new rxns
    if ~isempty(taskStructure(i).equations)
        rxn.equations=taskStructure(i).equations;
        rxn.lb=taskStructure(i).LBequ;
        rxn.ub=taskStructure(i).UBequ;
        rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
        %Allow for new metabolites to be added. This is because it should
        %be possible to add, say, a whole new pathway
        tModel=addRxns(tModel,rxn,3,[],true);
    end
    %Add changed bounds
    if ~isempty(taskStructure(i).changed)
        tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
        tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
    end

    %Solve and print
    sol=solveLP(tModel);
    if ~isempty(sol.x)
        %assign the fluxes
        essentialFluxes(:,i) = sol.x(1:numel(model.rxns));

        if ~taskStructure(i).shouldFail
            taskReport.ok(i,1)=true;
            if printOnlyFailed==false && printOutput==true
                fprintf(['PASS: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
            %Calculate the essential reactions
            if getEssential==true
                [~, taskEssential]=getEssentialRxns(tModel);
                %This is because there could be more reactions in tModel
                %than in model
                essentialRxns(taskEssential(taskEssential<=numel(model.rxns)),i)=true;
            end
        else
            taskReport.ok(i,1)=false;
            if printOutput==true
                fprintf(['PASS (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        end
    else
        if ~taskStructure(i).shouldFail
            taskReport.ok(i,1)=false;
            if printOutput==true
                fprintf(['FAIL: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        else
            taskReport.ok(i,1)=true;
            if printOnlyFailed==false && printOutput==true
                fprintf(['FAIL (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        end
    end
    if taskStructure(i).printFluxes && ~isempty(sol.x)
        sol=solveLP(tModel,1);
        if ~isempty(sol.x)
            printFluxes(tModel,sol.x,false,10^-6,[],'%rxnID (%eqn):%flux\n');
            fprintf('\n');
        end
    end
    tModel=model;
end
end

author: SysBioChalmers

repo: RAVEN


text:
function run_experiment(imdb_video)
%RUN_EXPERIMENT
%	contains the parameters that are specific to the current experiment.
% 	These are the parameters that change frequently and should not be committed to
% 	the repository but should be saved with the results of the experiment.

	% Parameters that should have no effect on the result.
	opts.prefetch = false;
	opts.gpus = 1;

	% Parameters that should be recorded.
	% opts.foo = 'bla';

	if nargin < 1
	    imdb_video = [];
    end
    opts.expDir = 'data_gray';
    opts.augment.grayscale = 0.25;
	experiment(imdb_video, opts);

end

author: shenjianbing

repo: TripletTracking


text:
function out = pedfremap(x)
% PEDFREMAP Piecewise extended density format remap
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

    out = amplitudetodensity(x);
    out(out > 128) = .5 * (out(out>128) + 128);
    out = uint8(out);
end
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function varargout = lvd_EditConstraintsGUI(varargin)
% LVD_EDITCONSTRAINTSGUI MATLAB code for lvd_EditConstraintsGUI.fig
%      LVD_EDITCONSTRAINTSGUI, by itself, creates a new LVD_EDITCONSTRAINTSGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITCONSTRAINTSGUI returns the handle to a new LVD_EDITCONSTRAINTSGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITCONSTRAINTSGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITCONSTRAINTSGUI.M with the given input arguments.
%
%      LVD_EDITCONSTRAINTSGUI('Property','Value',...) creates a new LVD_EDITCONSTRAINTSGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditConstraintsGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditConstraintsGUI_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_EditConstraintsGUI

% Last Modified by GUIDE v2.5 03-Dec-2018 17:19:14

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditConstraintsGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditConstraintsGUI_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_EditConstraintsGUI is made visible.
function lvd_EditConstraintsGUI_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_EditConstraintsGUI (see VARARGIN)

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

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

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, lvdData)
    set(handles.constraintsListBox,'String',lvdData.optimizer.constraints.getListboxStr());

    setTooltipStrings(handles, lvdData);

function setTooltipStrings(handles, lvdData)
    warning('off','MATLAB:ui:javaframe:PropertyToBeRemoved');

    hListbox = handles.constraintsListBox;
    jScrollPane = findjobj(hListbox);
    jListbox = jScrollPane.getViewport.getComponent(0);
    jListbox = handle(jListbox, 'CallbackProperties');

    tooltipStrs = lvdData.optimizer.constraints.getToolboxStrs();
    set(jListbox, 'MouseMovedCallback', @(src, evt) mouseMovedCallback(src, evt, hListbox, tooltipStrs));

    warning('on','MATLAB:ui:javaframe:PropertyToBeRemoved');

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditConstraintsGUI_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
        varargout{1} = true;
        close(handles.lvd_EditConstraintsGUI);
    end

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

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

% Hints: contents = cellstr(get(hObject,'String')) returns constraintsListBox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from constraintsListBox
    if(strcmpi(get(handles.lvd_EditConstraintsGUI,'SelectionType'),'open'))
        lvdData = getappdata(handles.lvd_EditConstraintsGUI,'lvdData');

        constraintSet = lvdData.optimizer.constraints;
        selConstraint = get(handles.constraintsListBox,'Value');
        constraint = constraintSet.getConstraintForInd(selConstraint);

        constraint.openEditConstraintUI(lvdData);

        listBoxStr = constraintSet.getListboxStr();
        set(handles.constraintsListBox,'String',listBoxStr);

        setTooltipStrings(handles, lvdData);
    end

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

% Hint: listbox 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 button press in addConstraintButton.
function addConstraintButton_Callback(hObject, eventdata, handles)
% hObject    handle to addConstraintButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditConstraintsGUI,'lvdData');

    listBoxStr = ConstraintEnum.getListBoxStr();

    [Selection,ok] = listdlg('PromptString','Select the constraint type:',...
                    'SelectionMode','single',...
                    'Name','Constraint Type',...
                    'ListString',listBoxStr, ...
                    'ListSize',[300 300]);

    if(ok == 1)
        [~, enum] = ConstraintEnum.getIndForName(listBoxStr{Selection});
        constClass = enum.class;
        input1 = enum.constructorInput1;

        newConstraint = eval(sprintf('%s.getDefaultConstraint(%s, %s)', constClass, 'input1', 'lvdData'));

        addConstraintTf = newConstraint.openEditConstraintUI(lvdData);

        if(addConstraintTf)
            constraintSet = lvdData.optimizer.constraints;
            constraintSet.addConstraint(newConstraint);

            listBoxStr = constraintSet.getListboxStr();
            set(handles.constraintsListBox,'String',listBoxStr);

            if(handles.constraintsListBox.Value <= 0)
                handles.constraintsListBox.Value = 1;
            elseif(handles.constraintsListBox.Value > length(listBoxStr))
                handles.constraintsListBox.Value = length(listBoxStr);
            end
        end
    end

    setTooltipStrings(handles, lvdData);

% --- Executes on button press in removeConstraintButton.
function removeConstraintButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeConstraintButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditConstraintsGUI,'lvdData');
    constraintSet = lvdData.optimizer.constraints;

    selConstraint = get(handles.constraintsListBox,'Value');
    constraint = constraintSet.getConstraintForInd(selConstraint);

    constraintSet.removeConstraint(constraint);

	listBoxStr = constraintSet.getListboxStr();
	set(handles.constraintsListBox,'String',listBoxStr);

	numConstraints = length(listBoxStr);
    if(selConstraint > numConstraints)
        set(handles.constraintsListBox,'Value',numConstraints);
    end

    setTooltipStrings(handles, lvdData);

% --- Executes on key press with focus on lvd_EditConstraintsGUI or any of its controls.
function lvd_EditConstraintsGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditConstraintsGUI (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'
            uiresume(handles.lvd_EditConstraintsGUI);
        case 'enter'
            uiresume(handles.lvd_EditConstraintsGUI);
        case 'escape'
            uiresume(handles.lvd_EditConstraintsGUI);
    end

% Mouse-movement callback
function mouseMovedCallback(jListbox, jEventData, hListbox, tooltipStrs)
   % Get the currently-hovered list-item
   mousePos = java.awt.Point(jEventData.getX, jEventData.getY);
   hoverIndex = jListbox.locationToIndex(mousePos) + 1;
%    listValues = get(hListbox,'string');
%    hoverValue = listValues{hoverIndex};

   % Modify the tooltip based on the hovered item
   if(hoverIndex > 0 && hoverIndex <= length(tooltipStrs))
       msgStr = tooltipStrs{hoverIndex};
       set(hListbox, 'Tooltip',msgStr);
   else
       set(hListbox, 'Tooltip','');
   end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function a = stirling_inverse ( n )

%*****************************************************************************80
%
%% STIRLING_INVERSE returns the inverse of the STIRLING matrix.
%
%  Comments:
%
%    The inverse of S1, the matrix of Stirling numbers of the first kind,
%    is S2, the matrix of Stirling numbers of the second kind.
%
%    S2(I,J) represents the number of distinct partitions of I elements
%    into J nonempty sets.  For any I, the sum over J of the Stirling
%    numbers S2(I,J) is represented by B(I), called "Bell's number",
%    and represents the number of distinct partitions of I elements.
%
%    For example, with 4 objects, there are:
%
%    1 partition into 1 set:
%
%      (A,B,C,D)
%
%    7 partitions into 2 sets:
%
%      (A,B,C) (D)
%      (A,B,D) (C)
%      (A,C,D) (B)
%      (A) (B,C,D)
%      (A,B) (C,D)
%      (A,C) (B,D)
%      (A,D) (B,C)
%
%    6 partitions into 3 sets:
%
%      (A,B) (C) (D)
%      (A) (B,C) (D)
%      (A) (B) (C,D)
%      (A,C) (B) (D)
%      (A,D) (B) (C)
%      (A) (B,D) (C)
%
%    1 partition into 4 sets:
%
%      (A) (B) (C) (D)
%
%    So S2(4,1) = 1, S2(4,2) = 7, S2(4,3) = 6, S2(4,4) = 1, and B(4) = 15.
%
%
%  First terms:
%
%    I/J: 1    2    3    4    5    6    7    8
%
%    1    1    0    0    0    0    0    0    0
%    2    1    1    0    0    0    0    0    0
%    3    1    3    1    0    0    0    0    0
%    4    1    7    6    1    0    0    0    0
%    5    1   15   25   10    1    0    0    0
%    6    1   31   90   65   15    1    0    0
%    7    1   63  301  350  140   21    1    0
%    8    1  127  966 1701 1050  266   28    1
%
%  Recursion:
%
%    S2(I,1) = 1 for all I.
%    S2(I,I) = 1 for all I.
%    S2(I,J) = 0 if I < J.
%
%    S2(I,J) = J * S2(I-1,J) + S2(I-1,J-1)
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is lower triangular.
%
%    A is nonnegative.
%
%    det ( A ) = 1.
%
%    A is unimodular.
%
%    LAMBDA(1:N) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the  matrix.
%
  a = zeros ( n, n );

  a(1,1) = 1.0;
  a(1,2:n) = 0.0;

  for i = 2 : n

    a(i,1) = 1.0;

    for j = 2 : n
      a(i,j) = j * a(i-1,j) + a(i-1,j-1);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
%   FILE NAME:    getRxnsFromMNX.m
%
%   PURPOSE: Generate the data structure for MetaNetX reactions
%

% Move to the target MNX folder

% Load the MNX reactions
T=readtable('reac_prop.xlsx','ReadVariableNames',1);
MNXRxns=table2struct(T,'ToScalar',true);
% Refine the equations by removing compartment suffix
MNXRxns.equations=regexprep(MNXRxns.Equation,'\@MNXD\d+','');
MNXRxns.equations=regexprep(MNXRxns.equations,'=','<=>');
MNXRxns.Description=regexprep(MNXRxns.Description,'1 `','');
MNXRxns.Description=regexprep(MNXRxns.Description,'`','');

% Construct the S matrix and list of metabolites
[S, mets, badRxns, reversible]=constructS(MNXRxns.equations);  %time-consuming
MNXRxns.S=S;
MNXRxns.mets=mets;
MNXRxns.badRxns=badRxns;
% Count metabolite number for reactions
num=numel(MNXRxns.EC);
MNXRxns.metNum=zeros(num,1);
for i=1:num
	MNXRxns.metNum(i,1)=numel(find(MNXRxns.S(:,i)));
end

% Construct the full S matrix and list of metabolites
MNXRxns.fullequations=regexprep(MNXRxns.Equation,'=','<=>');
[fullS, fullMets, badRxns, reversible]=constructS(MNXRxns.fullequations);  %time-consuming
MNXRxns.fullS=fullS;
MNXRxns.fullMets=fullMets;
MNXRxns.fullSbadRxns=badRxns;
% Count metabolite number for reactions
MNXRxns.fullMetNum=zeros(num,1);
for i=1:num
	MNXRxns.fullMetNum(i,1)=numel(find(MNXRxns.fullS(:,i)));
end

save('MNXRxns.mat','MNXRxns');

repo: Human-GEM


text:
function [eddies, tracks] = process_eddies_and_tracks_tolerance(type, dates, eddies_path, tracks, eddies_save_dir, tracks_save_dir)
% This function is for preparing eddy and track data for being shown in the
% eddy viewer. Simply run this function with the required parameters, then
% replace the data in the folders inside tracks_viewer with your own data.

%% function to load eddies.mat files, add flagged track eddies to the .mat
%%file, and then give that eddy all the required fields/properties with
%%values that are the average of the two eddies on either side of it.
if ~strcmp(eddies_save_dir(end), '/')
    eddies_save_dir = strcat(eddies_save_dir, '/');
end
if ~strcmp(tracks_save_dir(end), '/')
    tracks_save_dir = strcat(tracks_save_dir, '/');
end
eddies_names = get_eddies_names(eddies_path, type);
eddy_length = length(eddies_names);
fake_eddy_count = ones(eddy_length, 1);
orig_eddy_count = zeros(eddy_length, 1);
eddies = cell(eddy_length, 1);
iterations_since_needed = zeros(eddy_length, 1);
has_been_loaded = false(eddy_length, 1);
max_iterations = round(length(tracks) / 10);

first_time_slice = num2str(dates(1));

if ~exist(eddies_save_dir, 'dir')
    mkdir(eddies_save_dir);
end

if isempty(tracks)
    error('Please provide a non-empty input for tracks');
end

%% now go through and add the flagged eddies to eddies_cells and alter eddy indices in tracks, and resave them

%disp('altering eddies and tracks')
for i = 1:length(tracks)
    iterations_since_needed = iterations_since_needed + 1;
    if size(tracks{i},2) == 5
        for j = 1:size(tracks{i},1)
            if tracks{i}(j,5) ~= 0
                eddy_before_index = tracks{i}(j-1, 4);
                date_index = tracks{i}(j,3);
                iterations_since_needed(date_index - 1:date_index) = 0;
                has_been_loaded(date_index - 1:date_index) = true;
                serial_date = date_index + datenum(first_time_slice, 'yyyymmdd') - 1;
                date = datestr(serial_date, 'yyyymmdd');
                [eddies, orig_eddy_count] = fetch_eddies(eddies, eddies_names, date_index, orig_eddy_count, eddies_save_dir);
                [eddies, orig_eddy_count] = fetch_eddies(eddies, eddies_names, date_index - 1, orig_eddy_count, eddies_save_dir);
                fake_eddy_index = orig_eddy_count(date_index) + fake_eddy_count(date_index);
                fake_eddy_count(date_index) = fake_eddy_count(date_index) + 1;
                eddies{date_index}(fake_eddy_index) = eddies{date_index - 1}(eddy_before_index);
                eddies{date_index}(fake_eddy_index).Date = str2num(date);
                tracks{i}(j,4) = length(eddies{date_index});
            end
        end
    end
    for k = 1:length(iterations_since_needed)
        if iterations_since_needed(k) > max_iterations && ~isempty(eddies{k})
            %disp(['Saving eddies on day ', num2str(k), ' back to file.']);
            path_cell = strsplit(eddies_names{k}, '/');
            file_name = path_cell{end};
            if ~exist(eddies_save_dir, 'dir')
                mkdir(eddies_save_dir);
            end
            personal_save([eddies_save_dir, file_name], eddies{k});
            fake_eddy_count(k) = 1;
            eddies{k} = [];
        end
    end
end

for i = 1:length(eddies)
    if ~isempty(eddies{i})
        %disp(['Saving eddies on day ', num2str(i), ' back to file.']);
        path_cell = strsplit(eddies_names{i}, '/');
        file_name = path_cell{end};
        if ~exist(eddies_save_dir, 'dir')
            mkdir(eddies_save_dir);
        end
        personal_save([eddies_save_dir, file_name], eddies{i});
        eddies{i} = [];
    end
    if ~has_been_loaded(i)
        %disp(['Saving non-loaded eddies on day ', num2str(i), ' back to file.']);
        vars = load(eddies_names{i});
        s = fieldnames(vars);
        path_cell = strsplit(eddies_names{i}, '/');
        file_name = path_cell{end};
        if ~exist(eddies_save_dir, 'dir')
            mkdir(eddies_save_dir);
        end
        personal_save([eddies_save_dir, file_name], vars.(s{1}));
    end
end

%% now save the altered tracks:

if strcmp(type, 'cyclonic') == 1
    cyclonic_tracks = tracks;%#ok
    save([tracks_save_dir 'cyclonic_tracks_processed.mat'], 'cyclonic_tracks');
else
    anticyc_tracks = tracks;%#ok
    save([tracks_save_dir 'anticyc_tracks_processed.mat'], 'anticyc_tracks');
end

end
%% wrapper function to save eddies in parfor:

function personal_save(eddies_path, eddies)%#ok
save(eddies_path, 'eddies');
end

function [eddies_names] = get_eddies_names(path, type)
% path is the path to the eddies directory
% type is anticyclonic or cyclonic
if ~strcmp(path(end), '/')
    path = strcat(path, '/');
end
files = dir(path);
x = 0;
for i = 1:length(files)
    if ~isempty(strfind(files(i).name, [type, '_']))
        x = x + 1;
    end
end
eddies_names = cell(x, 1);
x = 1;
for i = 1:length(files)
    if files(i).isdir && ~isequal(files(i).name, '.') && ~isequal(files(i).name, '..')
        rec_names = get_eddies_names([path, files(i).name, '/'], type);
        for j = 1:length(rec_names)
            eddies_names{x} = rec_names{j};
            x = x + 1;
        end
        continue;
    end
    file = files(i).name;
    [~, name, ext] = fileparts([path, file]);
    if ~isempty(strfind(name, [type, '_'])) && strcmp(ext, '.mat')
        eddies_names{x} = [path, file];
        x = x + 1;
    end
end
end

function [eddies, orig_eddy_count] = fetch_eddies(eddies, eddies_names, index, orig_eddy_count, eddies_save_dir)
if index < 1 || index > length(eddies)
    disp(['Date index is ', num2str(index)]);
    error('Date index outside of eddy bounds.');
end
if isempty(eddies{index})
    %disp(['Loading eddies on day ', num2str(index)]);
    path_cell = strsplit(eddies_names{index}, '/');
    file_name = path_cell{end};
    if exist([eddies_save_dir, file_name], 'file')
        vars = load([eddies_save_dir, file_name]);
        s = fieldnames(vars);
        eddies{index} = vars.(s{1});
    else
        vars = load(eddies_names{index});
        s = fieldnames(vars);
        eddies{index} = vars.(s{1});
    end
    orig_eddy_count(index) = length(eddies{index});
end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function colin = points_colin_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% POINTS_COLIN_2D estimates the colinearity of 3 points in 2D.
%
%  Discussion:
%
%    The estimate of colinearity that is returned is the ratio
%    of the area of the triangle spanned by the points to the area
%    of the equilateral triangle with the same perimeter.
%
%    This estimate is 1 if the points are maximally noncolinear, 0 if the
%    points are exactly colinear, and otherwise is closer to 1 or 0 depending
%    on whether the points are far or close to colinearity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2), P2(2), P3(2), the points.
%
%    Output, real COLIN, the colinearity estimate.
%
  dim_num = 2;

  t(1:dim_num,1:3) = [ p1(1:dim_num); p2(1:dim_num); p3(1:dim_num) ]';

  area_triangle = triangle_area_2d ( t );

  if ( area_triangle == 0.0 )

    colin = 0.0;

  else

    perim = sqrt ( sum ( ( p2(1:dim_num) - p1(1:dim_num) ).^2 ) ) ...
          + sqrt ( sum ( ( p3(1:dim_num) - p2(1:dim_num) ).^2 ) ) ...
          + sqrt ( sum ( ( p1(1:dim_num) - p3(1:dim_num) ).^2 ) );

    side = perim / 3.0;

    area2 = 0.25 * sqrt ( 3.0 ) * side * side;

    colin = abs ( area_triangle ) / area2;

  end

  return
end


text:
% Example of fixed lag smoothing

rand('state', 1);
S = 2;
O = 2;
T = 7;
data = sample_discrete([0.5 0.5], 1, T);
transmat = mk_stochastic(rand(S,S));
obsmat = mk_stochastic(rand(S,O));
obslik = multinomial_prob(data, obsmat);
prior = [0.5 0.5]';

[alpha0, beta0, gamma0, ll0, xi0] = fwdback(prior, transmat, obslik);

w = 3;
alpha1 = zeros(S, T);
gamma1 = zeros(S, T);
xi1 = zeros(S, S, T-1);
t = 1;
b = obsmat(:, data(t));
olik_win = b; % window of conditional observation likelihoods
alpha_win = normalise(prior .* b);
alpha1(:,t) = alpha_win;
for t=2:T
  [alpha_win, olik_win, gamma_win, xi_win] = ...
      fixed_lag_smoother(w, alpha_win, olik_win, obsmat(:, data(t)), transmat);
  alpha1(:,max(1,t-w+1):t) = alpha_win;
  gamma1(:,max(1,t-w+1):t) = gamma_win;
  xi1(:,:,max(1,t-w+1):t-1) = xi_win;
end

e = 1e-1;
%assert(approxeq(alpha0, alpha1, e));
assert(approxeq(gamma0(:, T-w+1:end), gamma1(:, T-w+1:end), e));
assert(approxeq(xi0(:,:,T-w+1:end), xi1(:,:,T-w+1:end), e));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function h = r8mat_house_pre ( n, a, row, col )

%*****************************************************************************80
%
%% R8MAT_HOUSE_PRE computes a Householder pre-multiplier matrix.
%
%  Discussion:
%
%    H(ROW,COL) has the property that the COL-th column of
%    H(ROW,COL)*A is zero from entry ROW+1 to the end.
%
%    In the most common case, where a QR factorization is being computed,
%    ROW = COL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrices.
%
%    Input, real A(N,N), the matrix whose Householder matrix
%    is to be computed.
%
%    Input, integer ROW, COL, specify the location of the
%    entry of the matrix A which is to be preserved.  The entries in
%    the same column, but higher rows, will be zeroed out if A is
%    premultiplied by H.
%
%    Output, real H(N,N), the Householder matrix.
%

%
%  Set up the vector V.
%
  a_col(1:row-1,1) = 0.0;
  a_col(row:n,1) = a(row:n,col);

  v = r8vec_house_column ( n, a_col, row );
%
%  Form the matrix H(V).
%
  h = r8mat_house_form ( n, v );

  return
end


text:
function [data,times] = fiff_read_raw_segment_times(raw,from,to,sel)
%
% [data,times] = fiff_read_raw_segment_times(raw,from,to)
%
% Read a specific raw data segment
%
% raw    - structure returned by fiff_setup_read_raw
% from   - starting time of the segment in seconds
% to     - end time of the segment in seconds
% sel    - optional channel selection vector
%
% data   - returns the data matrix (channels x samples)
% times  - returns the time values corresponding to the samples (optional)
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.3  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.2  2006/04/21 14:23:16  msh
%   Further improvements in raw data reading
%
%

me='MNE:fiff_read_raw_segment_times';

if nargin == 3
   sel = [];
elseif nargin ~= 4
   error(me,'Incorrect number of arguments');
end
%
%   Convert to samples
%
from = floor(from*raw.info.sfreq);
to   = ceil(to*raw.info.sfreq);
%
%   Read it
%
[ data, times ] = fiff_read_raw_segment(raw,from,to,sel);

return;

end


text:
function tests = test_lang
  tests = functiontests(localfunctions);
end

function test_1(testCase)
    verifyTrue(testCase, spx.lang.is_class('spx.lang'));
    verifyFalse(testCase, spx.lang.is_class('spx.lang.noop'));
    spx.lang.noop;
end

author: indigits

repo: sparse-plex


text:
function [R, fv] = tnsre_applyClassifier(epo, Out, csp_w)

%     epo_proc = proc_filtButter(epo, 3, bandpass_filter);
%     epo_proc = proc_selectChannels(epo_proc, sub_channel); % Channel Selection
%     epo_proc = proc_commonAverageReference(epo_proc); % Applying Common Average Reference (CAR)

    fv = proc_linearDerivation(epo, csp_w);
    fv = proc_variance(fv);
    fv = proc_logarithm(fv);

%     fv.x = fv.x';

    R = applyClassifier(fv, 'wr_multiClass', Out.C);
    R = out2label(R);

    i = 0;
    for i=1:size(R, 2)
        switch R(i)
            case 1
                R(i) = 0;
            case 2
                R(i) = 1;

        end
    end
end

author: PatternRecognition

repo: OpenBMI


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function varargout = overlayImage2D(T,R,omega,m,varargin)
%
% A 2D image viewer for overlaying images R and T
%
% Input:
%   T           discretized image
%   R           discretized image
%   omega       describing the domain
%    m          number of discretization points
%   varargin    optional parameters like {'colorT','r'}
%
% Output:
%  i1,i2        image handles for T and R
%==============================================================================

function varargout = overlayImage2D(T,R,omega,m,varargin)

if nargin==0
    runMinimalExample;
	return;
end

colorT = [1,0,0];
colorR = [0,1,1];
scale  = 1;
alphadata = 0.5;

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

h  = (omega(2:2:end)-omega(1:2:end))./m;
xi = @(i) (omega(2*i-1)+h(i)/2:h(i):omega(2*i)-h(i)/2)';

if scale == 1,
  T = 255/max(T(:))*T;
  R = 255/max(R(:))*R;
end;

if length(size(T)) == 2,
  T = permute(reshape(uint8(T*colorT),[m,3]),[2,1,3]);
  R = permute(reshape(uint8(R*colorR),[m,3]),[2,1,3]);
else
  11
end;

cla;
i1  = image(xi(1),xi(2),R); axis xy image; hold on;
i2  = image(xi(1),xi(2),T); axis xy image; hold off;
set(i2,'alphaData',alphadata);

if nargout == 1, varargout = {[i1,i2]};  end;

%------------------------------------------------------------------------------
function runMinimalExample

help(mfilename)
setup2DhandData;
overlayImage2D(dataT(:),dataR(:),omega,m);
%==============================================================================


text:
function imdb = alot_get_database(alotDir, varargin)
opts.seed = 0 ;
opts.version = 'grey2';
opts = vl_argparse(opts, varargin) ;

rng(opts.seed, 'twister');

imdb.imageDir = fullfile(alotDir, opts.version) ;
imdb.maskDir = fullfile(alotDir, 'mask'); % DO NOT USE

cats = dir(imdb.imageDir);
cats = cats([cats.isdir] & ~ismember({cats.name}, {'.','..'})) ;
imdb.classes.name = {cats.name} ;
imdb.images.id = [] ;

for c=1:numel(cats)
  ims = dir(fullfile(imdb.imageDir, imdb.classes.name{c}, '*.png'));
  %imdb.images.name{c} = fullfile(imdb.classes.name{c}, {ims.name}) ;
  imdb.images.name{c} = cellfun(@(S) fullfile(imdb.classes.name{c}, S), ...
    {ims.name}, 'Uniform', 0);
  imdb.images.label{c} = c * ones(1,numel(ims)) ;
  if numel(ims) ~= 100, error('ops') ; end
  % http://cmp.felk.cvut.cz/~sulcmila/papers/Sulc-TR-2014-12.pdf
  % uses 20 for training, 80 for testing, per class
  sets = [1 * ones(1,20), 3 * ones(1,80)];
  imdb.images.set{c} = sets(randperm(100)) ;
end
imdb.images.name = horzcat(imdb.images.name{:}) ;
imdb.images.label = horzcat(imdb.images.label{:}) ;
imdb.images.set = horzcat(imdb.images.set{:}) ;
imdb.images.id = 1:numel(imdb.images.name) ;

imdb.segments = imdb.images ;
imdb.segments.imageId = imdb.images.id ;
imdb.segments.mask = strrep(imdb.images.name, 'image', 'mask') ;

% make this compatible with the OS imdb
imdb.meta.classes = imdb.classes.name ;
imdb.meta.inUse = true(1,numel(imdb.meta.classes)) ;
imdb.segments.difficult = false(1, numel(imdb.segments.id)) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function [e, edata, eprior] = rbferr(net, x, t)
%RBFERR	Evaluate error function for RBF network.
%
%	Description
%	E = RBFERR(NET, X, T) takes a network data structure NET together
%	with a matrix X of input vectors and a matrix T of target vectors,
%	and evaluates the appropriate error function E depending on
%	NET.OUTFN.  Each row of X corresponds to one input vector and each
%	row of T contains the corresponding target vector.
%
%	[E, EDATA, EPRIOR] = RBFERR(NET, X, T) additionally returns the data
%	and prior components of the error, assuming a zero mean Gaussian
%	prior on the weights with inverse variance parameters ALPHA and BETA
%	taken from the network data structure NET.
%
%	See also
%	RBF, RBFFWD, RBFGRAD, RBFPAK, RBFTRAIN, RBFUNPAK
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check arguments for consistency
switch net.outfn
case 'linear'
   errstring = consist(net, 'rbf', x, t);
case 'neuroscale'
   errstring = consist(net, 'rbf', x);
otherwise
   error(['Unknown output function ', net.outfn]);
end
if ~isempty(errstring);
  error(errstring);
end

switch net.outfn
case 'linear'
   y = rbffwd(net, x);
   edata = 0.5*sum(sum((y - t).^2));
case 'neuroscale'
   y = rbffwd(net, x);
   y_dist = sqrt(dist2(y, y));
   % Take t as target distance matrix
   edata = 0.5.*(sum(sum((t-y_dist).^2)));
otherwise
   error(['Unknown output function ', net.outfn]);
end

% Compute Bayesian regularised error
[e, edata, eprior] = errbayes(net, edata);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Parameters

F_update = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1];

Npop_particles = 4000;

Xstd_rgb = 50;
Xstd_pos = 25;
Xstd_vec = 5;

Xrgb_trgt = [255; 0; 0];

%% Loading Movie

vr = VideoReader('Person.wmv');

Npix_resolution = [vr.Width vr.Height];
Nfrm_movie = floor(vr.Duration * vr.FrameRate);

%% Object Tracking by Particle Filter

X = create_particles(Npix_resolution, Npop_particles);

for k = 1:Nfrm_movie

    % Getting Image
    Y_k = read(vr, k);

    % Forecasting
    X = update_particles(F_update, Xstd_pos, Xstd_vec, X);

    % Calculating Log Likelihood
    L = calc_log_likelihood(Xstd_rgb, Xrgb_trgt, X(1:2, :), Y_k);

    % Resampling
    X = resample_particles(X, L);

    % Showing Image
    show_particles(X, Y_k);
%    show_state_estimated(X, Y_k);

end

author: Sable

repo: mcbench-benchmarks


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function [v, a, b] = subsref(r, s)

	if s(1).type  ~= '.'
		%error('only .field supported')
	end

	% NOTE WELL:  the following code can't use getfield() since
	% getfield()  uses this, and Matlab will crash!!

	el = char(s(1).subs);
	switch el,
%-------------------------------------------------------------------------------------
% External Functions
	case 'fkine',
		q = s(2).subs;
		v = fkine(r,q{:});

	case 'ikine6s',
		q = s(2).subs;
		v = ikine6s(r,q{:});

	case 'ikine',
		q = s(2).subs;
		v = ikine(r,q{:});

	case 'maniplty',
		q = s(2).subs;
		[v, a] = maniplty(r,q{:});

	case 'rne'
		q = s(2).subs;
		[v, a] = rne(r,q{:});

	case 'isspherical',
		v = isspherical(r);

	case 'gravload'
		q = s(2).subs;
		v = gravload(r,q{:});

	case 'inertia'
		q = s(2).subs;
		v = inertia(r,q{:});

	case 'cinertia',
		q = s(2).subs;
		v = cinertia(r,q{:});

	case 'coriolis'
		q = s(2).subs;
		v = coriolis(r,q{:});

	case 'friction'
		q = s(2).subs;
		v = friction(r,q{:});

	case 'nofriction'
		if numel(s(2).subs) < 1
		v = nofriction(r);
		else
		q = s(2).subs;
		v = nofriction(r,q{:});
		end

	case 'islimit',
		q = s(2).subs;
		v = islimit(r,q{:});

	case 'payload',
		q = s(2).subs;
		l = r.links;
		l(r.n).m = q{1};
		l(r.n).r = q{2};
		r.links = l;
		v = SerialLink(r);

	case 'jacob0',
		q = s(2).subs;
		v = jacob0(r,q{:});

	case 'accel',
		q = s(2).subs;
		v = accel(r,q{:});

	case 'plot',
		q = s(2).subs;
		plot(r,q{:})

	case 'fdyn',
		q = s(2).subs;
		[v, a, b] = fdyn(r,q{:});

	case 'jacobn',
		q = s(2).subs;
		v = jacobn(r,q{:});

	case 'jtraj',
		q = s(2).subs;
		v = jtraj(r,q{:});

	case 'perturb',
		q = s(2).subs;
		v = perturb(r,q{:});

	case 'itorque',
		q = s(2).subs;
		v = itorque(r,q{:});

	case 'jacob_dot',
		q = s(2).subs;
		v = jacob_dot(r,q{:});

	case 'dyn',
		q = s(2).subs;
		v = dyn(r,q{:});

	case 'showlink',
		q = s(2).subs;
		v = showlink(r,q{:});

%-------------------------------------------------------------------------------------
	case 'links',

		if length(s) == 1,
			v = r.links;

		elseif length(s) == 2,
			if s(2).type == '()'
				j = s(2).subs;
				j = j{:};
				if (j < 1) | (j > r.n)
					error('link index out of bounds')
				end
			end
			v = r.links(j);

		elseif length(s) >= 3,
			if s(2).type == '()'
				j = s(2).subs;
				j = j{:};
				if (j < 1) | (j > r.n)
					error('link index out of bounds')
				end
			end
			if s(3).subs = 'dyn',
				link = r.links(j);
				v = link.dyn();
			end
		end
	case 'offset',
		L = r.links;
		v = [];
		for i=1:r.n,
			v = [v; L(i).offset];
		end
	case 'qlim',
		L = r.links;
		v = [];
		for i=1:r.n,
			v = [v; L(i).qlim];
		end

	case 'n',
		v = r.n;
	case 'name',
		v = r.name;
	case 'dh',
		v = rdh(r);
	case 'dyn'
		v = rdyn(r);
	case 'gravity'
		v = r.gravity;
	case 'tool'
		v = r.tool;
	case 'base'
		v = r.base;
	case 'mdh',
		v = r.mdh;
	case 'q',
		v = r.q;
	case 'plotopt',
		v = r.plotopt;
	case 'lineopt'
		v = r.lineopt;
	case 'shadowopt'
		v = r.shadowopt;
	case {'show', 'handle'}
		v = r.handle';
	otherwise, error('Unknown method in subsref')
	end

author: petercorke

repo: robotics-toolbox-matlab


text:
%RESIZE  Resizes an image
%
%     dst = cv.resize(src, dsize)
%     dst = cv.resize(src, fx, fy)
%     dst = cv.resize(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ input image.
% * __dsize__ output image size `[w,h]`.
% * __fx__ Scale factors along the horizontal axis.
% * __fy__ Scale factors along the vertical axis.
%
% When scale factors are specified, the output image size is computed as:
% `dsize = round([fx*size(src,2), fy*size(src,1)])`.
%
% ## Output
% * __dst__ output image. It has the size `dsize` or the size computed from
%   `size(src)` and `fx` and `fy`. The type of `dst` is the same as of `src`.
%
% ## Options
% * __Interpolation__ interpolation method, default 'Linear'. One of:
%   * __Nearest__ a nearest-neighbor interpolation
%   * __Linear__ a bilinear interpolation (used by default)
%   * __Cubic__ a bicubic interpolation over 4x4 pixel neighborhood
%   * __Area__ resampling using pixel area relation. It may be a preferred
%     method for image decimation, as it gives moire-free results. But when
%     the image is zoomed, it is similar to the 'Nearest' method.
%   * __Lanczos4__ a Lanczos interpolation over 8x8 pixel neighborhood
%   * __LinearExact__ a bit exact bilinear interpolation
%
% The function cv.resize resizes the image `src` down to or up to the
% specified size. The size and type of `dst` are derived from `src`, `dsize`,
% `fx` and `fy`. If you want to resize `src` so that it fits a specified size,
% you may call the function as follows:
%
%     dst = cv.resize(src, [w,h], 'Interpolation',interp)
%
% If you want to decimate the image by factor of 2 in each direction, you can
% call the function this way:
%
%     dst = cv.resize(src, 0.5, 0.5, 'Interpolation',interp)
%
% To shrink an image, it will generally look best with 'Area' interpolation,
% whereas to enlarge an image, it will generally look best with 'Cubic' (slow)
% or 'Linear' (faster but still looks OK).
%
% See also: cv.warpAffine, cv.warpPerspective, cv.remap, imresize
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sigmaVal = fwhm2sigma(fwhmVal)
%FWHM2SIGMA
%   Conversion factor proved here:
%   https://brainder.org/2011/08/20/gaussian-kernels-convert-fwhm-to-sigma/

    conversionFactor = sqrt(8*log(2));

    sigmaVal = fwhmVal./conversionFactor;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [ edge_pointer, edge_data, xy ] = grf_example ( node_num, edge_num )

%*****************************************************************************80
%
%% GRF_EXAMPLE sets up a GRF example.
%
%  Discussion:
%
%    The example is known as the Coxeter graph.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Skiena,
%    Implementing Discrete Mathematics,
%    Combinatorics and Graph Theory with Mathematica,
%    Addison-Wesley, 1990.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer EDGE_NUM, the number of edges.
%
%    Output, integer EDGE_POINTER(NODE_NUM+1), pointers to
%    the beginning of edge data for each node.
%
%    Output, integer EDGE_DATA(EDGE_NUM), the edge data.
%
%    Output, real XY(2,NODE_NUM), the node coordinates.
%
  edge_pointer = [ ...
    1,  4,  7, 10, 13, 16, 19, 22, 25, 28, ...
   31, 34, 37, 40, 43, 46, 49, 52, 55, 58, ...
   61, 64, 67, 70, 73, 76, 79, 82, 85 ]';

  edge_data = [ ...
     8,   2,   3, ...
    14,   1,   5, ...
     9,   4,   1, ...
    10,   7,   3, ...
    13,   2,   6, ...
    12,   5,   7, ...
    11,   6,   4, ...
    25,  20,   1, ...
    24,  21,   3, ...
    23,  15,   4, ...
    22,  16,   7, ...
    28,  17,   6, ...
    27,  18,   5, ...
    26,  19,   2, ...
    10,  18,  19, ...
    11,  19,  20, ...
    12,  21,  20, ...
    13,  15,  21, ...
    14,  16,  15, ...
     8,  17,  16, ...
     9,  18,  17, ...
    11,  27,  24, ...
    10,  28,  25, ...
     9,  26,  22, ...
     8,  23,  27, ...
    14,  24,  28, ...
    13,  25,  22, ...
    12,  26,  23 ]';

  xy = [ ...
    0.412,   0.984; ...
    0.494,   0.984; ...
    0.366,   0.926; ...
    0.388,   0.862; ...
    0.546,   0.926; ...
    0.518,   0.860; ...
    0.458,   0.818; ...
    0.152,   0.684; ...
    0.264,   0.682; ...
    0.354,   0.680; ...
    0.458,   0.670; ...
    0.554,   0.672; ...
    0.658,   0.668; ...
    0.774,   0.692; ...
    0.164,   0.450; ...
    0.228,   0.448; ...
    0.274,   0.390; ...
    0.242,   0.330; ...
    0.194,   0.278; ...
    0.146,   0.328; ...
    0.102,   0.390; ...
    0.668,   0.472; ...
    0.638,   0.416; ...
    0.656,   0.334; ...
    0.714,   0.270; ...
    0.798,   0.326; ...
    0.830,   0.408; ...
    0.754,   0.466 ]';

  return
end


text:
sim('pend');

a = logsout.angle;
f = logsout.force;

figure(3)
plot(a.time, a.data*180/pi)
xlabel('time')
ylabel('angle [degree]')
grid on;

figure(4)
plot(f.time, f.data)
xlabel('time')
ylabel('force [N]')
grid on;


text:
% Test file for ShortTimeFourierTransformFrameComputer computer
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef test_compute < matlab.unittest.TestCase
    properties
        computer
    end

    properties (TestParameter)
        bank = {GaborFilterBank(MelScaling())}
        frame_length_ms = {25}
        frame_shift_ms = {10}
        frame_style = {'causal'} %{'causal', 'centered'}
        include_energy = {true}  %{true, false}
        pad_to_nearest_power_of_two = {true} % {true, false}
        use_log = {true} %{true, false}
        use_power = {true} %{true, false}
        buff_len = {0, 1, 2 ^ 8, 2^10} %{ 0, 1, 2 ^ 8, 2 ^ 10} % empty buffer', 'length 1 buffer', 'medium buffer', 'large buffer'
    end

    methods (Test)
        % Test validility of the function
        function test_validate(testCase,...
                bank, ...
                frame_length_ms, ...
                frame_shift_ms, ...
                frame_style, ...
                include_energy, ...
                pad_to_nearest_power_of_two, ...
                use_log, ...
                use_power)

            testCase.computer = ShortTimeFourierTransformFrameComputer(...
                                        bank, ...
                                        frame_length_ms, ...
                                        frame_shift_ms, ...
                                        frame_style, ...
                                        include_energy, ...
                                        pad_to_nearest_power_of_two, ...
                                        use_log, ...
                                        use_power);
            testCase.verifyTrue(testCase.computer.isvalid);
        end

        function test_framewise_matches_full(testCase,...
                bank, ...
                frame_length_ms, ...
                frame_shift_ms, ...
                frame_style, ...
                include_energy, ...
                pad_to_nearest_power_of_two, ...
                use_log, ...
                use_power,...
                buff_len)
            % This test might fail due to minor numerical error

            % Create buffer
            buff = rand(buff_len,1);

            % Create Computer
            testCase.computer = ShortTimeFourierTransformFrameComputer(...
                                        bank, ...
                                        frame_length_ms, ...
                                        frame_shift_ms, ...
                                        frame_style, ...
                                        include_energy, ...
                                        pad_to_nearest_power_of_two, ...
                                        use_log, ...
                                        use_power);

            feats_full = testCase.computer.compute_full(buff);
            feats_framewise = Util.frame_by_frame_calculation(testCase.computer, buff);
            testCase.verifyEqual(feats_full, feats_framewise, 'RelTol', 1);
        end

        function test_zero_samples_generate_zero_features(testCase,...
                bank, ...
                frame_length_ms, ...
                frame_shift_ms, ...
                frame_style, ...
                include_energy, ...
                pad_to_nearest_power_of_two, ...
                use_log, ...
                use_power)
            % Create Computer
            testCase.computer = ShortTimeFourierTransformFrameComputer(...
                                        bank, ...
                                        frame_length_ms, ...
                                        frame_shift_ms, ...
                                        frame_style, ...
                                        include_energy, ...
                                        pad_to_nearest_power_of_two, ...
                                        use_log, ...
                                        use_power);

            testCase.verifyEqual(size(testCase.computer.compute_full(double.empty(0))), [0 testCase.computer.num_coeffs]);
            testCase.verifyEqual(size(testCase.computer.compute_chunk(double.empty(0))), [0  testCase.computer.num_coeffs]);
            testCase.verifyEqual(size(testCase.computer.finalize()), [0 testCase.computer.num_coeffs]);
        end

        function test_finalize_twice_generates_no_coefficients(testCase,...
                bank, ...
                frame_length_ms, ...
                frame_shift_ms, ...
                frame_style, ...
                include_energy, ...
                pad_to_nearest_power_of_two, ...
                use_log, ...
                use_power)
            % Create Computer
            testCase.computer = ShortTimeFourierTransformFrameComputer(...
                                        bank, ...
                                        frame_length_ms, ...
                                        frame_shift_ms, ...
                                        frame_style, ...
                                        include_energy, ...
                                        pad_to_nearest_power_of_two, ...
                                        use_log, ...
                                        use_power);
            buff = rand(testCase.computer.frame_length * 2, 1);
            coeffs = [...
                testCase.computer.compute_chunk(buff); ...
                testCase.computer.finalize()...
                ];

            shape_coeffs = size(coeffs);
            testCase.verifyGreaterThanOrEqual(shape_coeffs(1), 1);
            shape_finalize = size(testCase.computer.finalize());
            testCase.verifyEqual(shape_finalize, [0 testCase.computer.num_coeffs]);
        end

        function test_repeated_calls_generate_same_results(testCase,...
                bank, ...
                frame_length_ms, ...
                frame_shift_ms, ...
                frame_style, ...
                include_energy, ...
                pad_to_nearest_power_of_two, ...
                use_log, ...
                use_power, ...
                buff_len)

            % Create buffer
            buff = rand(buff_len,1);

            % Create Computer
            testCase.computer = ShortTimeFourierTransformFrameComputer(...
                                        bank, ...
                                        frame_length_ms, ...
                                        frame_shift_ms, ...
                                        frame_style, ...
                                        include_energy, ...
                                        pad_to_nearest_power_of_two, ...
                                        use_log, ...
                                        use_power);
            testCase.verifyEqual(testCase.computer.compute_full(buff), ...
                testCase.computer.compute_full(buff));
            testCase.verifyEqual(...
                Util.frame_by_frame_calculation(testCase.computer, buff), ...
                Util.frame_by_frame_calculation(testCase.computer, buff));
        end

    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function vol = read_asa_vol(fn)

% READ_ASA_VOL reads an ASA volume conductor file
%
% all data is converted to the following units
%   vertices        mm
%   conductivities  S/m

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

Nbnd  = read_ini(fn, 'NumberBoundaries=', '%d');
UnitC = read_ini(fn, 'UnitConduct', '%s');
UnitP = read_ini(fn, 'UnitPosition', '%s');
cond  = read_ini(fn, 'Conductivities', '%f');
radii = read_ini(fn, 'Radii', '%f');
pos   = read_ini(fn, 'Positions', '%f');
bnd1  = read_ini(fn, 'Boundary1', '%s');
bnd2  = read_ini(fn, 'Boundary2', '%s');
bnd3  = read_ini(fn, 'Boundary3', '%s');
bnd4  = read_ini(fn, 'Boundary4', '%s');

if ~isempty(radii) || ~isempty(pos)
  % this appears to be a spherical volume conductor
  if strcmpi(UnitP,'mm')
    radii = 1*radii;
    pos   = 1*pos;
  elseif strcmpi(UnitP,'cm')
    radii = 100*radii;
    pos   = 100*pos;
  elseif strcmpi(UnitP,'m')
    radii = 1000*radii;
    pos   = 1000*pos;
  else
    ft_error(sprintf('Unknown unit of distance for volume (%s)', UnitP));
  end
end

if strcmpi(UnitC,'s/m')
  cond = cond/1;
elseif strcmpi(UnitC,'s/cm')
  cond = cond/100;
elseif strcmpi(UnitC,'s/mm')
  cond = cond/1000;
else
  ft_error(sprintf('Unknown unit of conductivity for volume (%s)', UnitC));
end

if ~isempty(radii)
  % this appears to be a spherical volume conductor
  vol.radius = radii;
  vol.cond   = cond;
  vol.center = pos;
else
  % this appears to be a realistical volume conductor
  [path, name, ext] = fileparts(fn);
  if Nbnd>=1
    vol.bnd(1) = read_asa_bnd(fullfile(path, bnd1));
  end
  if Nbnd>=2
    vol.bnd(2) = read_asa_bnd(fullfile(path, bnd2));
  end
  if Nbnd>=3
    vol.bnd(3) = read_asa_bnd(fullfile(path, bnd3));
  end
  if Nbnd>=4
    vol.bnd(4) = read_asa_bnd(fullfile(path, bnd4));
  end
  if Nbnd>=5
    ft_error('cannot read more than 4 boundaries');
  end

  % if there is a precomputed matrix, read it from an external file
  mat_file = read_ini(fn, 'Matrix', '%s');
  if ~isempty(mat_file)
    nr = read_ini(fullfile(path, mat_file), 'NumberRows=', '%d');
    nc = read_ini(fullfile(path, mat_file), 'NumberColumns=', '%d');
    mab_file = read_ini(fullfile(path, mat_file), 'Matrix', '%s');
    fid = fopen_or_error(fullfile(path, mab_file), 'rb', 'ieee-le');
    vol.mat = fread(fid, [nr nc], 'float32');
    fclose(fid);
    % remove the factor 2 that ASA assumes in the system matrix
    vol.mat = vol.mat/2;
    % scale the system matrix corresponding to vertex coordinates in mm
    vol.mat = vol.mat*100;
  end

  vol.cond   = cond;
end

% remove all empty fields
field=fieldnames(vol);
for i=1:length(field)
  if isempty(getfield(vol, field{i}))
    vol = rmfield(vol, field{i});
  end
end


text:
function Hbar = computeHbar(Cbar, Qbar, Rbar)
    Hbar = Cbar'*Qbar*Cbar + Rbar;
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
% EEG_EVENTTYPES  - return a list of event or urevent types in a dataset and
%                     the respective number of events of each type. Output event
%                     types are sorted in reverse order of their number. If no
%                     outputs, print this list on the commandline instead.
%
% Usage:
%        >> [types,numbers] = eeg_eventtypes(EEG);
% Inputs:
%        EEG        - EEGLAB dataset structure
% Outputs:
%        types      - cell array of event type strings
%        numbers    - vector giving the numbers of each event type in the data
%
% Example:
%           >> eeg_eventtypes(EEG);       % print number of each event types
%
% Author: Scott Makeig, SCCN/INC/UCSD, April 28, 2004-

%        >> [types,numbers] = eeg_eventtypes(EEG,types);
%        >> [types,numbers] = eeg_eventtypes(EEG,'urevents',types);
% Inputs:
%        EEG        - EEGLAB dataset structure
%        'urevents' - return event information for the EEG.urevent structure
%        types      - {cell array} of event types to return or print.
% Outputs:
%        types      - cell array of event type strings
%        numbers    - vector giving the numbers of each event type in the data
%
% Note:  Numeric (ur)event types are converted to strings, so, for example,
%        types {13} and {'13'} are not distinguished.
%
% Example:
%           >> eeg_eventtypes(EEG);       % print number of each event types
%
% Currently disabled:
%           >> eeg_eventtypes(EEG,'urevent');  % print hist. of urevent types
%           >> eeg_eventtypes(EEG,{'rt'});% print number of 'rt' events
%           >> eeg_eventtypes(EEG,'urevent',{'rt','break'});
%                                         % print numbers of 'rt' and 'break'
%                                         % type urevents

% Copyright (C) 2004 Scott Makeig, SCCN/INC/UCSD, smakeig@ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
%
%
% event types can be numbers, Stefan Debener, 05/12/2006,
% added 2nd and 3rd args, sorted outputs by number, Scott Makeig, 09/12/06

function [types,numbers] = eeg_eventtypes(EEG,arg2,arg3)

if nargin< 1
   help eeg_eventtypes
   return
end
if ~isstruct(EEG)
   error('EEG argument must be a dataset structure')
end

if ~isfield(EEG,'event')
   error('EEG.event field not found');
end
if nargin > 1
    error('Multiple input arguments are currently disabled');
end

UREVENTS = 0; % flag returning info for urevents instead of events
typelist = [];
if nargin>1
   if ischar(arg2)
       if strcmp(arg2,'urevent') || strcmp(arg2,'urevents')
            UREVENTS = 1; % change flag
       else
            error('second argument string not understood')
       end
       if nargin>2
            if iscell(arg3)
                  typelist = arg3;
            end
       end
   elseif iscell(arg2)
       typelist = arg2;
   end
end
if ~isempty(typelist)     % cast to cell array of strings
   for k=1:length(typelist)
       if isnumeric(typelist{k})
          typelist{k} = num2str(typelist{k});
       end
   end
end

if ~UREVENTS
   nevents = length(EEG.event);
   alltypes = cell(nevents,1);
   for k=1:nevents
       if isnumeric(EEG.event(k).type)
          alltypes{k} = num2str(EEG.event(k).type);
       else
          alltypes{k} = EEG.event(k).type;
       end
   end
else
   nevents = length(EEG.urevent);
   alltypes = cell(nevents,1);
   for k=1:nevents
       if isnumeric(EEG.urevent(k).type)
          alltypes{k} = num2str(EEG.urevent(k).type);
       else
          alltypes{k} = EEG.urevent(k).type;
       end
   end
end

[types i j] = unique_bc(alltypes);

istypes = 1:length(types);
notistypes = [];
if ~isempty(typelist)
  notistypes = ismember_bc(typelist,types);
  istypes = ismember_bc(types,typelist(find(notistypes==1))); % types in typelist?
  notistypes = typelist(find(notistypes==0));
end

types(~istypes) = []; % restrict types to typelist

ntypes = length(types);
numbers = zeros(ntypes + length(notistypes),1);
for k=1:ntypes
  numbers(k) = length(find(j==k));
end
types = [types(:); notistypes(:)]; % concatenate the types not found
ntypes = length(types);
for j = 1:length(notistypes)
   numbers(k+j) = 0;
end

% sort types in reverse order of event numbers
[numbers nsort] = sort(numbers);
numbers = numbers(end:-1:1);
types = types(nsort(end:-1:1));

%
% print output on commandline %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargout < 1
  fprintf('\n');
  if UREVENTS
    fprintf('EEG urevent types:\n\n')
  else
    fprintf('EEG event types:\n\n')
  end

  maxx = 0;
  for k=1:ntypes
    x = length(types{k});
    if x > maxx
        maxx = x; % find max type name length
    end
  end
  for k=1:ntypes
    fprintf('  %s',types{k});
    for j=length(types{k})+1:maxx+4
       fprintf(' ');
    end
    fprintf('%d\n',numbers(k));
 end
 fprintf('\n');
 clear types % no return variables
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%% Tire comparison
% Comparison between tire models.
%
% <html>
% <script src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
% </html>
%
%% Description
% The typical relation between the lateral force and the slip angle can be observed in the figure below (Adapted from [1]). Besides, its possible to verify the definition of slip angle.
%
% <<CurvaCaracteristica.svg>>
%
%% Equivalence
% Given a reference Pacejka tire model it is possible to obtain an equivalent linear and polynomial model. The cornering stiffness of all models must be equal and the maximal lateral force of the Pacejka and Polynomial models must be the same.
%
% The Pacejka model depends on the parameters \(a_0\), \(a_1\), \(a_2\), \(a_3\), \(a_4\), \(a_5\), \(a_6\) e \(a_7\) that defines the constants \(B\), \(C\), \(D\) and \(E\) wich can be used to define the constants of the equivalent models.
%
% The equivalent linear tire model has cornering stiffness \(K\) given by
%
% \[ K = B C D \]
%
% The equivalent polynomial tire model has coeficients \(k_1\) and \(k_2\) given by
%
% \[ k_1 = B C D \]
%
% \[ k_2 = (4 k_1^3)/(27 F_{y, Max}^2) \]
%
% where \(F_{y, Max}\) is the maximal lateral force of the reference characteristic curve.
%
%% Model comparison
%

% Code start

clear ; close all ; clc

deriva = (0:0.1:15)*pi/180;         % ngulo de deriva [rad]

a0  = 1.3;
a1  = 2.014156;
a2  = 710.5013;
a3  = 5226.341;
a4  = 78.87699;
a5  = 0.01078379;
a6  = -0.004759443;
a7  = -1.8572;
a8  = 0;
a9  = 0;
a10 = 0;
a11 = 0;
a12 = 0;
a13 = 0;

TirePac = VehicleDynamicsLateral.TirePacejka();

Fz = 4e+03;
camber = 0;
TirePac.a0 = a0;
TirePac.a1 = a1;
TirePac.a2 = a2;
TirePac.a3 = a3;
TirePac.a4 = a4;
TirePac.a5 = a5;
TirePac.a6 = a6;
TirePac.a7 = a7;
TirePac.a8 = a8;
TirePac.a9 = a9;
TirePac.a10 = a10;
TirePac.a11 = a11;
TirePac.a12= a12;
TirePac.a13 = a13;

muy0 = TirePac.a1 * Fz/1000 + TirePac.a2;
D = muy0 * Fz/1000;
BCD = TirePac.a3 * sin(2 * atan(Fz/1000/TirePac.a4))*(1-TirePac.a5 * abs(camber));

% Linear tire model

K = BCD * 180/pi;

TireLin = VehicleDynamicsLateral.TireLinear();
TireLin.k = K;

% Polynomial tire MODEL

k1 = BCD * 180/pi;
k2 = (4 * k1^3)/(27 * D^2);

TirePol = VehicleDynamicsLateral.TirePolynomial();
TirePol.k1 = k1;
TirePol.k2 = k2;

% Lateral force
FyPac = TirePac.Characteristic(deriva, Fz, muy0/1000);
FyLin = TireLin.Characteristic(deriva);
FyPol = TirePol.Characteristic(deriva);

% Graphics
g = VehicleDynamicsLateral.Graphics(TirePac);

figure(1)
ax = gca;
set(ax, 'NextPlot', 'add', 'Box', 'on', 'XGrid', 'on', 'YGrid', 'on')
p = plot(deriva * 180/pi,-FyLin, 'Color', 'g', 'Marker', 's', 'MarkerFaceColor', 'g', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
p = plot(deriva * 180/pi,-FyPol, 'Color', 'b', 'Marker', '^', 'MarkerFaceColor', 'b', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
p = plot(deriva * 180/pi,-FyPac, 'Color', 'r', 'Marker', 'o', 'MarkerFaceColor', 'r', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
xlabel('\(\alpha\) [grau]', 'Interpreter', 'Latex')
ylabel('\(F_y\) [N]', 'Interpreter', 'Latex')
l = legend('Linear', 'Polynomial', 'Pacejka');
set(l, 'Interpreter', 'Latex', 'Location', 'NorthWest')

%%
% In the above figure the 3 characteristic curves are plotted. For small slip angles all models are equivalents.
%
%% Comparison treatment slip angle
%

deriva180 = (0:0.1:180)*pi/180;     % ngulo de deriva de 0  180 graus [rad]

% Sem tratamento
ALPHA = deriva180 * 180/pi;
C = a0;
muy = muy0;
E = a6 * Fz/1000 + a7;
B = BCD/(C * D);
Sh = a8 * camber + a9 * Fz/1000 + a10;
Sv = a11 * Fz/1000 * camber + a12 * Fz/1000 + a13;
ALPHAeq = muy0/muy*(ALPHA + Sh);
% Reference characteristics
fy = D * sin(C * atan(B * ALPHAeq - E*(B * ALPHAeq - atan(B * ALPHAeq))));
% Lateral force
FyPacSem180 = -muy/muy0*(fy + Sv);

% Com tratamento
FyPacCom180 = TirePac.Characteristic(deriva180, Fz, muy0/1000);

figure(2)
ax = gca;
set(ax, 'NextPlot', 'add', 'Box', 'on', 'XGrid', 'on', 'YGrid', 'on')
p = plot(deriva180 * 180/pi,-FyPacSem180, 'Color', 'm', 'Marker', 'd', 'MarkerFaceColor', 'm', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
p = plot(deriva180 * 180/pi,-FyPacCom180, 'Color', 'r', 'Marker', 'o', 'MarkerFaceColor', 'r', 'MarkeredgeColor', 'k', 'MarkerSize', 7);
g.changeMarker(p, 10);
plot([90 90],[0 3000],'--k')    % Linha vertical de simetria
xlabel('\(\alpha\) [grau]', 'Interpreter', 'Latex')
ylabel('\(F_y\) [N]', 'Interpreter', 'Latex')
l = legend('Pacejka without treatment', 'Pacejka with treatment ');
set(l, 'Interpreter', 'Latex', 'Location', 'SouthEast')

%%
% In the above figure we can see that the curve is symmetric to a vertical line positioned at \(\alpha = 90 [graus]\).
%
%% References
% [1] GILLESPIE, T. D. Fundamentals of vehicle dynamics. [S.l.]: Society of Automotive Engineers Warrendale, PA, 1992.
%
%% See Also
%
% <../../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function G = egrad(X)
  global P;
  global PA;
  % Same comment here about Xmat.
  Xmat = X.U*X.S*X.V';
  G = P.*Xmat - PA;
end

author: andrewssobral

repo: lrslibrary


text:
function [ a, b ] = p11_lim ( dim_num )

%*****************************************************************************80
%
%% P11_LIM returns the integration limits for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Load NIFTI dataset header. Support both *.nii and *.hdr/*.img file
%  extension. If file extension is not provided, *.hdr/*.img will be
%  used as default.
%
%  Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)
%
%  filename - NIFTI file name.
%
%  Returned values:
%
%  hdr - struct with NIFTI header fields.
%
%  filetype	- 0 for Analyze format (*.hdr/*.img);
%		  1 for NIFTI format in 2 files (*.hdr/*.img);
%		  2 for NIFTI format in 1 file (*.nii).
%
%  fileprefix - NIFTI file name without extension.
%
%  machine    - a string, see below for details. The default here is 'ieee-le'.
%
%    'native'      or 'n' - local machine format - the default
%    'ieee-le'     or 'l' - IEEE floating point with little-endian
%                           byte ordering
%    'ieee-be'     or 'b' - IEEE floating point with big-endian
%                           byte ordering
%    'vaxd'        or 'd' - VAX D floating point and VAX ordering
%    'vaxg'        or 'g' - VAX G floating point and VAX ordering
%    'cray'        or 'c' - Cray floating point with big-endian
%                           byte ordering
%    'ieee-le.l64' or 'a' - IEEE floating point with little-endian
%                           byte ordering and 64 bit long data type
%    'ieee-be.l64' or 's' - IEEE floating point with big-endian byte
%                           ordering and 64 bit long data type.
%
%  Number of scanned images in the file can be obtained by:
%  num_scan = hdr.dime.dim(5)
%
%  Part of this file is copied and modified under GNU license from
%  MRI_TOOLBOX developed by CNSP in Flinders University, Australia
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [hdr, filetype, fileprefix, machine] = load_nii_hdr(fileprefix)

   if ~exist('fileprefix','var'),
      error('Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)');
   end

   machine = 'ieee-le';
   new_ext = 0;

   if findstr('.nii',fileprefix)
      new_ext = 1;
      fileprefix = strrep(fileprefix,'.nii','');
   end

   if findstr('.hdr',fileprefix)
      fileprefix = strrep(fileprefix,'.hdr','');
   end

   if findstr('.img',fileprefix)
      fileprefix = strrep(fileprefix,'.img','');
   end

   if new_ext
      fn = sprintf('%s.nii',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.nii".', fileprefix);
         error(msg);
      end
   else
      fn = sprintf('%s.hdr',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
         error(msg);
      end
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   else
      fseek(fid,0,'bof');

      if fread(fid,1,'int32') == 348
         hdr = read_header(fid);
         fclose(fid);
      else
         fclose(fid);

         %  first try reading the opposite endian to 'machine'
         %
         switch machine,
         case 'ieee-le', machine = 'ieee-be';
         case 'ieee-be', machine = 'ieee-le';
         end

         fid = fopen(fn,'r',machine);

         if fid < 0,
            msg = sprintf('Cannot open file %s.',fn);
            error(msg);
         else
            fseek(fid,0,'bof');

            if fread(fid,1,'int32') ~= 348

               %  Now throw an error
               %
               msg = sprintf('File "%s" is corrupted.',fn);
               error(msg);
            end

            hdr = read_header(fid);
            fclose(fid);
         end
      end
   end

   if strcmp(hdr.hist.magic, 'n+1')
      filetype = 2;
   elseif strcmp(hdr.hist.magic, 'ni1')
      filetype = 1;
   else
      filetype = 0;
   end

   return					% load_nii_hdr

%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)

        %  Original header structures
	%  struct dsr
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

    dsr.hk   = header_key(fid);
    dsr.dime = image_dimension(fid);
    dsr.hist = data_history(fid);

    %  For Analyze data format
    %
    if ~strcmp(dsr.hist.magic, 'n+1') & ~strcmp(dsr.hist.magic, 'ni1')
        dsr.hist.qform_code = 0;
        dsr.hist.sform_code = 0;
    end

    return					% read_header

%---------------------------------------------------------------------
function [ hk ] = header_key(fid)

    fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                     /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */
	%
	% int sizeof_header   Should be 348.
	% char regular        Must be 'r' to indicate that all images and
	%                     volumes are the same size.

    hk.sizeof_hdr    = fread(fid, 1,'int32')';	% should be 348!
    hk.data_type     = deblank(fread(fid,10,'*char')');
    hk.db_name       = deblank(fread(fid,18,'*char')');
    hk.extents       = fread(fid, 1,'int32')';
    hk.session_error = fread(fid, 1,'int16')';
    hk.regular       = fread(fid, 1,'*char')';
    hk.dim_info      = fread(fid, 1,'uchar')';

    return					% header_key

%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
        %       /*
        %           dim[0]      Number of dimensions in database; usually 4.
        %           dim[1]      Image X dimension;  number of *pixels* in an image row.
        %           dim[2]      Image Y dimension;  number of *pixel rows* in slice.
        %           dim[3]      Volume Z dimension; number of *slices* in a volume.
        %           dim[4]      Time points; number of volumes in database
        %       */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%	/*
	%		pixdim[] specifies the voxel dimensions:
	%		pixdim[1] - voxel width, mm
	%		pixdim[2] - voxel height, mm
	%		pixdim[3] - slice thickness, mm
	%		pixdim[4] - volume timing, in msec
	%					..etc
	%	*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

    dime.dim        = fread(fid,8,'int16')';
    dime.intent_p1  = fread(fid,1,'float32')';
    dime.intent_p2  = fread(fid,1,'float32')';
    dime.intent_p3  = fread(fid,1,'float32')';
    dime.intent_code = fread(fid,1,'int16')';
    dime.datatype   = fread(fid,1,'int16')';
    dime.bitpix     = fread(fid,1,'int16')';
    dime.slice_start = fread(fid,1,'int16')';
    dime.pixdim     = fread(fid,8,'float32')';
    dime.vox_offset = fread(fid,1,'float32')';
    dime.scl_slope  = fread(fid,1,'float32')';
    dime.scl_inter  = fread(fid,1,'float32')';
    dime.slice_end  = fread(fid,1,'int16')';
    dime.slice_code = fread(fid,1,'uchar')';
    dime.xyzt_units = fread(fid,1,'uchar')';
    dime.cal_max    = fread(fid,1,'float32')';
    dime.cal_min    = fread(fid,1,'float32')';
    dime.slice_duration = fread(fid,1,'float32')';
    dime.toffset    = fread(fid,1,'float32')';
    dime.glmax      = fread(fid,1,'int32')';
    dime.glmin      = fread(fid,1,'int32')';

    return					% image_dimension

%---------------------------------------------------------------------
function [ hist ] = data_history(fid)

	%  Original header structures
	%  struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

    hist.descrip     = deblank(fread(fid,80,'*char')');
    hist.aux_file    = deblank(fread(fid,24,'*char')');
    hist.qform_code  = fread(fid,1,'int16')';
    hist.sform_code  = fread(fid,1,'int16')';
    hist.quatern_b   = fread(fid,1,'float32')';
    hist.quatern_c   = fread(fid,1,'float32')';
    hist.quatern_d   = fread(fid,1,'float32')';
    hist.qoffset_x   = fread(fid,1,'float32')';
    hist.qoffset_y   = fread(fid,1,'float32')';
    hist.qoffset_z   = fread(fid,1,'float32')';
    hist.srow_x      = fread(fid,4,'float32')';
    hist.srow_y      = fread(fid,4,'float32')';
    hist.srow_z      = fread(fid,4,'float32')';
    hist.intent_name = deblank(fread(fid,16,'*char')');
    hist.magic       = deblank(fread(fid,4,'*char')');

    fseek(fid,253,'bof');
    hist.originator  = fread(fid, 5,'int16')';

    return					% data_history

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
close all;
clear all;
clc;

data_file_path = 'bin/ar_omp_phase_transition_gaussian_dict_gaussian_data.mat';
options.export = false;
options.export_dir = 'bin';
options.export_name = 'gaussian_dict_gaussian_data';
%options.chosen_ks = [2, 4, 8, 16, 32, 64];
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'OMP-AR', options);

data_file_path = 'bin/omp_phase_transition_gaussian_dict_gaussian_data.mat';
options.export = false;
options.export_dir = 'bin';
options.export_name = 'gaussian_dict_gaussian_data';
%options.chosen_ks = [2, 4, 8, 16, 32, 64];
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'OMP', options);


text:
function w = nc_compute ( n, x_min, x_max, x )

%*****************************************************************************80
%
%% NC_COMPUTE computes a Newton-Cotes quadrature rule.
%
%  Discussion:
%
%    For the interval [X_MIN,X_MAX], the Newton-Cotes quadrature rule
%    estimates
%
%      Integral ( X_MIN <= X <= X_MAX ) F(X) dX
%
%    using N abscissas X and weights W:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) ).
%
%    For the CLOSED rule, the equally spaced abscissas include A and B.
%    For the OPEN rule, the equally spaced abscissas do not include A and B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Input, real X_MIN, X_MAX, the endpoints of the interval.
%
%    Input, real X(N), the abscissas.
%
%    Output, real W(N,1), the weights.
%
  d = zeros ( n, 1 );
  w = zeros ( n, 1 );

  for i = 1 : n
%
%  Compute the Lagrange basis polynomial which is 1 at X(I),
%  and zero at the other nodes.
%
    d(1:n) = 0.0;
    d(i) = 1.0;

    for j = 2 : n
      for k = j : n
        d(n+j-k) = ( d(n+j-k-1) - d(n+j-k) ) / ( x(n+1-k) - x(n+j-k) );
      end
    end

    for j = 1 : n - 1
      for k = 1 : n - j
        d(n-k) = d(n-k) - x(n-k-j+1) * d(n-k+1);
      end
    end
%
%  Evaluate the antiderivative of the polynomial at the endpoints.
%
    yvala = d(n) / n;
    for j = n - 1 : -1 : 1
      yvala = yvala * x_min + d(j) / j;
    end
    yvala = yvala * x_min;

    yvalb = d(n) / n;
    for j = n - 1 : -1 : 1
      yvalb = yvalb * x_max + d(j) / j;
    end
    yvalb = yvalb * x_max;

    w(i) = yvalb - yvala;

  end

  return
end

author: johannesgerer