text:
classdef (Sealed) Toolbar < wt.abstract.BaseWidget & wt.mixin.TitleColorable ...
        & wt.mixin.FontStyled
    % A configurable toolbar

    % Copyright 2020-2021 The MathWorks Inc.

    %% Events
    events (HasCallbackProperty, NotifyAccess = protected)

        % Triggered when a button is pushed
        ButtonPushed

    end %events

    %% Public properties
    properties (AbortSet)

        % Sections that are part of the toolbar
        Section (:,1) wt.toolbar.HorizontalSection

    end %properties

    properties (Dependent, AbortSet, UsedInUpdate = false)

        % Color of the dividers between horizontal sections
        DividerColor

    end %properties

    %% Internal Properties
    properties ( Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % The listbox control
        ListBox (1,1) matlab.ui.control.ListBox

        % The label for each section
        SectionLabel (:,1) matlab.ui.control.Label

        % Buttons used when space is limited
        SectionButton (:,1) matlab.ui.control.Button

        % The dummy section filling up remaining space
        DummySection (:,1) matlab.ui.container.Container

        % Listen to section changes
        SectionChangedListener event.listener

        % Listen to size changes
        SizeChangedListener event.listener

        % Listen to button pushes in sections
        ButtonPushedListener event.listener

    end %properties

    properties (Dependent, NonCopyable, Access = {?wt.test.BaseWidgetTest})

        % Indicates which sections are open
        SectionIsOpen (:,1) logical

    end %properties

    properties (Constant, Access = private)

        % The down arrow mask for the button icons
        BUTTON_MASK (:,:) logical = sectionButtonIconMask()

    end %properties

    %% Protected methods
    methods (Access = protected)

        function setup(obj)

            % Call superclass setup first to establish the grid
            obj.setup@wt.abstract.BaseWidget();

            % Set default size
            obj.Position(3:4) = [500 90];

            % Configure style defaults
            obj.FontColor = [1 1 1] * 0.3;
            obj.TitleColor = [1 1 1] * 0.5;

            % Configure grid
            obj.Grid.Padding = [0 0 0 0]; %If changed, modify updateLayout!
            obj.Grid.RowHeight = {'1x'};
            obj.Grid.ColumnWidth = {};
            obj.Grid.ColumnSpacing = 1; %If changed, modify updateLayout!

            obj.Grid.RowHeight = {'1x',15};
            obj.Grid.ColumnWidth = {'1x'};
            obj.Grid.RowSpacing = 0;
            obj.Grid.BackgroundColor = [1 1 1]*0.5;

            % Add a dummy section to color the empty space
            obj.DummySection = uicontainer(obj.Grid);
            obj.DummySection.Layout.Row = [1 2];
            obj.BackgroundColorableComponents = obj.DummySection;

            % Listen to size changes
            obj.SizeChangedListener = event.listener(obj,'SizeChanged',...
                @(h,e)obj.updateLayout());

        end %function

        function update(obj)
            % Update the toolbar

            % Update each section
            numSections = numel(obj.Section);
            for sIdx = 1:numSections

                % Do we need to add more labels/buttons?
                if sIdx > numel(obj.SectionLabel)
                    obj.SectionLabel(sIdx) = uilabel(obj.Grid);
                    obj.SectionButton(sIdx) = uibutton(obj.Grid);
                end

                % Update each label
                obj.SectionLabel(sIdx).Parent = obj.Grid;
                obj.SectionLabel(sIdx).Layout.Row = 2;
                obj.SectionLabel(sIdx).Layout.Column = sIdx;
                obj.SectionLabel(sIdx).BackgroundColor = obj.BackgroundColor;
                obj.SectionLabel(sIdx).FontSize = 10;
                obj.SectionLabel(sIdx).HorizontalAlignment = 'center';
                obj.SectionLabel(sIdx).Text = upper(obj.Section(sIdx).Title);

                % Update each button
                obj.SectionButton(sIdx).Parent = obj.Grid;
                obj.SectionButton(sIdx).Layout.Row = [1 2];
                obj.SectionButton(sIdx).Layout.Column = sIdx;
                obj.SectionButton(sIdx).BackgroundColor = obj.BackgroundColor;
                obj.SectionButton(sIdx).FontSize = 10;
                obj.SectionButton(sIdx).IconAlignment = 'bottom';
                obj.SectionButton(sIdx).WordWrap = 'on';
                obj.SectionButton(sIdx).Visible = 'off';
                obj.SectionButton(sIdx).ButtonPushedFcn = @(h,e)obj.onPanelButtonPushed(e);
                obj.SectionButton(sIdx).Text = upper(obj.Section(sIdx).Title);

                % Update each section
                if ~isequal(obj.Section(sIdx).Parent, obj.Grid)
                    obj.Section(sIdx).Parent = obj.Grid;
                end
                obj.Section(sIdx).Layout.Row = 1;
                obj.Section(sIdx).Layout.Column = sIdx;

            end %for

            % Remove any extra labels/buttons
            numLabels = numel(obj.SectionLabel);
            if numLabels > numSections
                idxRemove = (numSections+1):numLabels;
                delete(obj.SectionLabel(idxRemove));
                obj.SectionLabel(idxRemove) = [];
                delete(obj.SectionButton(idxRemove));
                obj.SectionButton(idxRemove) = [];
            end

            % Unparent removed components
            isHorizontalSection = get(obj.Grid.Children,'Type') == ...
                lower("wt.toolbar.HorizontalSection");
            oldSections = obj.Grid.Children(isHorizontalSection);
            removedSections = setdiff(oldSections, obj.Section);
            set(removedSections,'Parent',[]);

            % Update component style lists
            obj.TitleColorableComponents = obj.SectionLabel;
            obj.BackgroundColorableComponents = [
                obj.Section
                obj.SectionButton
                obj.SectionLabel
                obj.DummySection
                ];
            obj.FontStyledComponents = [
                obj.Section
                obj.SectionButton
                ];

            % Update listeners
            obj.ButtonPushedListener = event.listener(obj.Section,...
                'ButtonPushed',@(h,e)obj.onButtonPushed(e));
            obj.SectionChangedListener = event.listener(obj.Section,...
                'PropertyChanged',@(h,e)obj.update());

            % Update the layout
            obj.updateLayout();

        end %function

        function updateLayout(obj)
            % Dynamically configure the toolbar based on space

            % Return if no sections to show
            if isempty(obj.Section)
                return
            end

            % How many sections
            numSections = numel(obj.Section);

            % Total available width to place sections
            pos = getpixelposition(obj);
            wAvail = pos(3);

            % Subtract off needed spacing between components
            spacingNeeded = numSections - 1;
            wAvail = wAvail - spacingNeeded;

            % Calculate the width options
            sectionWidth = [obj.Section.TotalWidth];
            minSectionWidth = [obj.Section.MinimizedWidth];
            widthOptions = cumsum(sectionWidth) + ...
                flip( cumsum([0 minSectionWidth(1:end-1)]) ) ;

            % Which panels can be fully shown?
            isFullyShown = widthOptions <= wAvail;

            % Loop on each obj.Section
            for idx = 1:numel(obj.Section)

                % If unparented, fix it
                if ~isequal(obj.Section(idx).Parent, obj.Grid)
                    obj.Section(idx).Parent = obj.Grid;
                    obj.Section(idx).Layout.Row = 1;
                    obj.Section(idx).Layout.Column = idx;
                end

                % Toggle visibilities
                obj.Section(idx).Visible = isFullyShown(idx);
                if idx <= numel(obj.SectionLabel)
                    obj.SectionLabel(idx).Visible = isFullyShown(idx);
                    obj.SectionButton(idx).Visible = ~isFullyShown(idx);
                end

            end %for idx = 1:numel(section)

            % Update the grid column widths
            sectionWidth(~isFullyShown) = minSectionWidth(~isFullyShown);
            obj.Grid.ColumnWidth = [num2cell(sectionWidth),'1x'];

            % Place the dummy section
            obj.DummySection.Layout.Column = numSections + 1;

        end %function

        function updateBackgroundColorableComponents(obj)

            % Update button icons
            obj.updateButtonIcons();

            % Override the default, not setting the Grid background
            hasProp = isprop(obj.BackgroundColorableComponents,'BackgroundColor');
            set(obj.BackgroundColorableComponents(hasProp),...
                "BackgroundColor",obj.BackgroundColor);

        end %function

        function updateFontStyledComponents(obj,varargin)

            % Update button icons
            obj.updateButtonIcons();

            % Call the superclass method
            obj.updateFontStyledComponents@wt.mixin.FontStyled(varargin{:});

        end %function

        function updateButtonIcons(obj)
            % Color the down arrow on buttons the same as font color

            % Has foreground or background color changed?
            if ~isempty(obj.SectionButton) && ( ...
                    ~isequal(obj.SectionButton(end).FontColor, obj.FontColor) || ...
                    ~isequal(obj.DummySection.BackgroundColor, obj.BackgroundColor) )

                % Create the button icon
                icon = cell(1,3);
                for cIdx = 1:3
                    icon{cIdx} = ~obj.BUTTON_MASK * obj.BackgroundColor(cIdx);
                    icon{cIdx}(obj.BUTTON_MASK) = obj.FontColor(cIdx);
                end
                icon = cat(3,icon{:});

                % Update the icon on each button
                set(obj.SectionButton,'Icon',icon);

            end %if

        end %function

        function onButtonPushed(obj,evt)

            % Close any open sections
            obj.updateLayout();

            % Trigger event
            notify(obj,"ButtonPushed",evt);

        end %function

    end %methods

    %% Private methods
    methods (Access = private)

        function onPanelButtonPushed(obj,e)

            % Which panel?
            sectionButton = e.Source;
            isActivatedSection = obj.SectionButton == sectionButton;
            section = obj.Section(isActivatedSection);

            % Get the status of section panels to change
            isOpenSection = obj.SectionIsOpen;
            isDeactivating = isOpenSection & ~isActivatedSection;
            isActivating = ~isOpenSection & isActivatedSection;

            % Is the pushed section already open?
            if any(isOpenSection(isActivatedSection))
                % Instead update the layout to closse it
                obj.updateLayout();
                return
            end

            % Close any deactivated sections
            if any(isDeactivating)
                set(obj.Section(isDeactivating), 'Visible', 'off');
            end

            % Activate the specified section
            if any(isActivating)

                % Where are things located now?
                bPos = getpixelposition(sectionButton, true);
                fig = ancestor(obj,'figure');
                wt.utility.fastSet(fig,"Units","pixels");
                figureWidth = fig.Position(3);

                % Where should the panel go?
                panelX = bPos(1);
                panelWidth = section.TotalWidth;
                panelHeight = bPos(4) - obj.Grid.RowHeight{2} - obj.Grid.RowSpacing;
                panelY = bPos(2) - panelHeight - 1;

                % Adjust panel X position if needed
                panelRightEdge = panelX + panelWidth;
                if panelRightEdge > figureWidth
                    panelX = figureWidth - panelWidth;
                end

                % Now, position and show the panel as a dropdown
                panelPos = [panelX panelY panelWidth panelHeight];

                % Turn it on
                set(obj.Section(isActivating), 'Parent', fig);
                set(obj.Section(isActivating), 'Position', panelPos);
                set(obj.Section(isActivating), 'Visible', 'on');

            end %if any(isActivating)

        end %function

    end %methods

    %% Accessors
    methods

        function value = get.DividerColor(obj)
            value = obj.Grid.BackgroundColor;
        end
        function set.DividerColor(obj,value)
            obj.Grid.BackgroundColor = value;
        end

        function value = get.SectionIsOpen(obj)
            value = false(size(obj.Section));
            for idx = 1:numel(obj.Section)
                value(idx) = ~isequal(obj.Section(idx).Parent, obj.Grid);
            end
        end

    end %methods

end % classdef

%% Helper Functions

function mask = sectionButtonIconMask()

mask = logical([
    0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0
    1 1 1 1 1 1 1 1 1
    0 1 1 1 1 1 1 1 0
    0 0 1 1 1 1 1 0 0
    0 0 0 1 1 1 0 0 0
    0 0 0 0 1 0 0 0 0
    ]);

end %function

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%sculling (essentially sculling+rotation) correction for each minor
%interval
%gyro & acc must be increment type (3,length) outputs. alg is the desired
%sculling algorithm
function [vinc, ainc, corr]=sculling_minor(gyro, acc, alg)
inlen=size(acc,2);
switch (alg)
    case (0)
        %ignagni(1998):algorithm 1 -> (1996):algo2
        outlen=floor((inlen-3)/2)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:2:inlen
            vinc(:,ind)=sum(acc(:,i-1:i),2);
            ainc(:,ind)=sum(gyro(:,i-1:i),2);

            rot=0.5*cross(ainc(:,ind), vinc(:,ind));
            vr_a=cross(((-1/30)*gyro(:,i-2)+(11/15)*gyro(:,i-1)), acc(:,i));
            vr_b=cross(((-1/30)*acc(:,i-2)+(11/15)*acc(:,i-1)), gyro(:,i));
            corr(:,ind)=rot+vr_a+vr_b;
            ind=ind+1;
        end
    case (1)
        %ignagni(1998):algorithm 2 -> (1996):algo3
        outlen=floor((inlen-3)/3)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:3:inlen
            vinc(:,ind)=sum(acc(:,i-2:i),2);
            ainc(:,ind)=sum(gyro(:,i-2:i),2);

            rot=0.5*cross(ainc(:,ind), vinc(:,ind));
            vr_a=cross(((9/20)*gyro(:,i-2)+(27/20)*gyro(:,i-1)), acc(:,i));
            vr_b=cross(((9/20)*acc(:,i-2)+(27/20)*acc(:,i-1)), gyro(:,i));
            corr(:,ind)=rot+vr_a+vr_b;
            ind=ind+1;
        end
    case(2)
        %ignagni(1998):algorithm 3 -> (1996):algo5
        outlen=floor((inlen-4)/3)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=4:3:inlen
            vinc(:,ind)=sum(acc(:,i-2:i),2);
            ainc(:,ind)=sum(gyro(:,i-2:i),2);

            rot=0.5*cross(ainc(:,ind), vinc(:,ind));
            vr_a=cross(((3/280)*gyro(:,i-3)+(57/140)*gyro(:,i-2)+(393/280)*gyro(:,i-1)), acc(:,i));
            vr_b=cross(((3/280)*acc(:,i-3)+(57/140)*acc(:,i-2)+(393/280)*acc(:,i-1)), gyro(:,i));
            corr(:,ind)=rot+vr_a+vr_b;
            ind=ind+1;
        end
    case(3)
        %ignagni(1998):algorithm 4 -> (1996):algo6
        outlen=floor((inlen-5)/3)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=5:3:inlen
            vinc(:,ind)=sum(acc(:,i-2:i),2);
            ainc(:,ind)=sum(gyro(:,i-2:i),2);

            rot=0.5*cross(ainc(:,ind), vinc(:,ind));
            vr_a=cross(((-1/420)*gyro(:,i-4)+(1/40)*gyro(:,i-3)+(157/420)*gyro(:,i-2)+(1207/840)*gyro(:,i-1)), acc(:,i));
            vr_b=cross(((-1/420)*acc(:,i-4)+(1/40)*acc(:,i-3)+(157/420)*acc(:,i-2)+(1207/840)*acc(:,i-1)), gyro(:,i));
            corr(:,ind)=rot+vr_a+vr_b;
            ind=ind+1;
        end
    case(4)
        %ignagni(1998):algorithm 5 -> (1996):algo7
        outlen=floor((inlen-4)/4)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=4:4:inlen
            vinc(:,ind)=sum(acc(:,i-3:i),2);
            ainc(:,ind)=sum(gyro(:,i-3:i),2);

            rot=0.5*cross(ainc(:,ind), vinc(:,ind));
            vr_a=cross(((54/105)*gyro(:,i-3)+(92/105)*gyro(:,i-2)+(214/105)*gyro(:,i-1)), acc(:,i));
            vr_b=cross(((54/105)*acc(:,i-3)+(92/105)*acc(:,i-2)+(214/105)*acc(:,i-1)), gyro(:,i));
            corr(:,ind)=rot+vr_a+vr_b;
            ind=ind+1;
        end
    case(5) %%quadratic fit to both gyro and accel outputs (similar to the coning algorithm in ignagni(1990):algoD
        outlen=floor((inlen-3)/3)+1;
        vinc=zeros(3,outlen);
        ainc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:3:inlen
            vinc(:,ind)=sum(acc(:,i-2:i),2);
            ainc(:,ind)=sum(gyro(:,i-2:i),2);

            %compute the polynomial coefficients
            b=(11/6)*gyro(:,i-2)+(-7/6)*gyro(:,i-1)+(1/3)*gyro(:,i);
            c=(-2)*gyro(:,i-2)+(3)*gyro(:,i-1)+(-1)*gyro(:,i);
            d=(1/2)*gyro(:,i-2)+(-1)*gyro(:,i-1)+(1/2)*gyro(:,i);

            e=(11/6)*acc(:,i-2)+(-7/6)*acc(:,i-1)+(1/3)*acc(:,i);
            f=(-2)*acc(:,i-2)+(3)*acc(:,i-1)+(-1)*acc(:,i);
            g=(1/2)*acc(:,i-2)+(-1)*acc(:,i-1)+(1/2)*acc(:,i);

            %Compute sculling(+rotation)
            corr(:,ind)=(3^2/2)*cross(b,e)+(3^3/3)*(cross(b,f)+cross(c,e)/2)+(3^4/4)*(cross(b,g)+cross(c,f)/2+cross(d,e)/3)+...
                (3^5/5)*(cross(c,g)/2+cross(d,f)/3)+(3^6/6)*cross(d,g)/3;

            ind=ind+1;
        end

    otherwise
        disp('undefined algorithm');
end

repo: ekfmonoslam


text:
function [ A ] = opt_catStruct( varargin )
% CATSTRUCT   Concatenate or merge structures with different fieldnames
%   X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ...
%   into one new structure X. X contains all fields present in the various
%   structures. An example:
%
%     A.name = 'Me' ;
%     B.income = 99999 ;
%     X = catstruct(A,B)
%     % -> X.name = 'Me' ;
%     %    X.income = 99999 ;
%
%   If a fieldname is not unique among structures (i.e., a fieldname is
%   present in more than one structure), only the value from the last
%   structure with this field is used. In this case, the fields are
%   alphabetically sorted. A warning is issued as well. An axample:
%
%     S1.name = 'Me' ;
%     S2.age  = 20 ; S3.age  = 30 ; S4.age  = 40 ;
%     S5.honest = false ;
%     Y = catstruct(S1,S2,S3,S4,S5) % use value from S4
%
%   The inputs can be array of structures. All structures should have the
%   same size. An example:
%
%     C(1).bb = 1 ; C(2).bb = 2 ;
%     D(1).aa = 3 ; D(2).aa = 4 ;
%     CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa
%
%   The last input can be the string 'sorted'. In this case,
%   CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically.
%   To sort the fieldnames of a structure A, you could use
%   CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that.
%
%   When there is nothing to concatenate, the result will be an empty
%   struct (0x0 struct array with no fields).
%
%   NOTE: To concatenate similar arrays of structs, you can use simple
%   concatenation:
%     A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;

%   NOTE: This function relies on unique. Matlab changed the behavior of
%   its set functions since 2013a, so this might cause some backward
%   compatibility issues when dulpicated fieldnames are found.
%
%   See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS

% version 4.1 (feb 2015), tested in R2014a
% Copyright (c) 2015, Jos van der Geest
% All rights reserved.
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% Created in 2005
% Revisions
%   2.0 (sep 2007) removed bug when dealing with fields containing cell
%                  arrays (Thanks to Rene Willemink)
%   2.1 (sep 2008) added warning and error identifiers
%   2.2 (oct 2008) fixed error when dealing with empty structs (thanks to
%                  Lars Barring)
%   3.0 (mar 2013) fixed problem when the inputs were array of structures
%                  (thanks to Tor Inge Birkenes).
%                  Rephrased the help section as well.
%   4.0 (dec 2013) fixed problem with unique due to version differences in
%                  ML. Unique(...,'last') is no longer the deafult.
%                  (thanks to Isabel P)
%   4.1 (feb 2015) fixed warning with narginchk

narginchk(1,Inf) ;
N = nargin ;

if ~isstruct(varargin{end}),
    if isequal(varargin{end},'sorted'),
        narginchk(2,Inf) ;
        sorted = 1 ;
        N = N-1 ;
    else
        error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
    end
else
    sorted = 0 ;
end

sz0 = [] ; % used to check that all inputs have the same size

% used to check for a few trivial cases
NonEmptyInputs = false(N,1) ;
NonEmptyInputsN = 0 ;

% used to collect the fieldnames and the inputs
FN = cell(N,1) ;
VAL = cell(N,1) ;

% parse the inputs
for ii=1:N,
    X = varargin{ii} ;
    if ~isstruct(X),
        error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
    end

    if ~isempty(X),
        % empty structs are ignored
        if ii > 1 && ~isempty(sz0)
            if ~isequal(size(X), sz0)
                error('catstruct:UnequalSizes','All structures should have the same size.') ;
            end
        else
            sz0 = size(X) ;
        end
        NonEmptyInputsN = NonEmptyInputsN + 1 ;
        NonEmptyInputs(ii) = true ;
        FN{ii} = fieldnames(X) ;
        VAL{ii} = struct2cell(X) ;
    end
end

if NonEmptyInputsN == 0
    % all structures were empty
    A = struct([]) ;
elseif NonEmptyInputsN == 1,
    % there was only one non-empty structure
    A = varargin{NonEmptyInputs} ;
    if sorted,
        A = orderfields(A) ;
    end
else
    % there is actually something to concatenate
    FN = cat(1,FN{:}) ;
    VAL = cat(1,VAL{:}) ;
    FN = squeeze(FN) ;
    VAL = squeeze(VAL) ;

    [UFN,ind] = unique(FN, 'last') ;
    % If this line errors, due to your matlab version not having UNIQUE
    % accept the 'last' input, use the following line instead
    % [UFN,ind] = unique(FN) ; % earlier ML versions, like 6.5

    if numel(UFN) ~= numel(FN),
%         warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ;
        sorted = 1 ;
    end

    if sorted,
        VAL = VAL(ind,:) ;
        FN = FN(ind,:) ;
    end

    A = cell2struct(VAL, FN);
    A = reshape(A, sz0) ; % reshape into original format

end

author: PatternRecognition

repo: OpenBMI


text:
function S = APlusAtranspose( PSF, Acenter, boundary)
%
% This function exists specifically to handle the addition
% A + Atranspose, where A is a psfMatrix.  It's called by
% the overloaded method plus.m whenever the input matrices
% are transposes of each other.
%
%  INPUT: PSF = the psf.image field of the psfMatrix that's being added
%                to its own transpose
%         Acenter = the center coordinates of the input PSF(s)
%         boundary = the boundary condition of A
%
%  OUTPUT is a new psfMatrix with fields as follows:
%        .psf is created from the summation of the original psf fields.
%             Centers are aligned before adding.
%             Care is taken so that the resulting psf is no larger than
%             the original psf (if possible).
%        .matdata as created as usual (by constructMatrix)
%        .boundary is the same as A
%        .type is determined by number of images in the psf
%        .transpose = 0
%

%  L. Perrone  4/9/02

%disp('You are entering the transpose zone')
sa = size(PSF);

% Make everything look like a 3D image in order to accommodate
% actual 3D images in the loop (i,j,k)
if length(sa)==1
  sa = [sa 1 1];
elseif length(sa) == 2
  sa = [sa 1];
end

newPsfs=cell(sa);
for k=1:sa(3);
 for j=1:sa(2);
  for i=1:sa(1);

    summand1 = PSF{i,j,k};
    summand1 = summand1/sum(summand1(:));
    summand2 = summand1;
    summand2 = flipdim(summand2,1);
    summand2 = flipdim(summand2,2);
    summand2 = flipdim(summand2,3);
    Bcenter{i,j,k} = size(summand1) - Acenter{i,j,k};

    % Pad around each image, however necessary,
    % so that the resulting images are equal in size and have
    % the same center coordinates.
    Asize = size(summand1);
    centerA = Acenter{i,j,k};
    Bsize = size(summand2);
    centerB = Bcenter{i,j,k};
    Aprepad = centerB - centerA;
    Aprepadsize = (Aprepad + abs(Aprepad))/2;
    summand1 = padarray(summand1,Aprepadsize,'pre');
    Bprepad = -Aprepad;
    Bprepadsize = (Bprepad + abs(Bprepad))/2;
    summand2 = padarray(summand2,Bprepadsize,'pre');
    Apost = Asize - centerA;
    Bpost = Bsize - centerB;
    Apostpad = Bpost - Apost;
    Apostpadsize = (Apostpad + abs(Apostpad))/2;
    summand1 = padarray(summand1,Apostpadsize,'post');
    Bpostpad = -Apostpad;
    Bpostpadsize = (Bpostpad + abs(Bpostpad))/2;
    summand2 = padarray(summand2,Bpostpadsize,'post');

    % Add the psfs together to get the sum psf.
    newPsf = summand1 + summand2;
    newPsf = newPsf/sum(newPsf(:));

    % Now for the messy stuff.
    % Cut off excess rows &/or columns so the resultant PSF is
    % no bigger than the input PSFs... but check the size of
    % what you want to cut off before chopping.  If it's too large
    % a piece, leave it in.
    maxSize = max( Asize, Bsize);
    sumSize = size(newPsf);
    chopSize = sumSize - maxSize;
    if length(chopSize)==1
       chopSize=[chopSize 0 0];
    elseif length(chopSize)==2
       chopSize=[chopSize 0];
    end
    if (chopSize(1) > 0)
      newPsf = circshift(newPsf,[chopSize(1) 0 0]);
      rowNorms = zeros(2*chopSize(1),1);
      for ii=1:2*chopSize(1)
       rowNorms(ii) = norm(squeeze(newPsf(ii,:,:)));
      end
      blooble = zeros(chopSize(1)+1,1);
      for ii=1:chopSize(1)+1
        blooble(ii) = norm(rowNorms(ii:ii+chopSize(1)-1));
      end
      cutStart = find(min(blooble)==blooble);
      if length(cutStart) ~= 1 %% (if more than one min is identified)
         cutStart = cutStart(fix(length(cutStart)/2)+1);
      end
      if blooble(cutStart) < .01*norm(newPsf(:))
         newPsf = cat(1, newPsf(1:cutStart-1,:,:), newPsf(cutStart+chopSize(1):sumSize(1),:,:));
         newPsf = circshift(newPsf, [cutStart-1 0 0]);
     else
        newPsf = circshift(newPsf,[-chopSize(1) 0 0]);
        disp('Row norms are too large for chopping')
      end
    end %%%   if (chopSize(1) > 0)

    if (chopSize(2) > 0)
      newPsf = circshift(newPsf,[0 chopSize(2) 0]);
      colNorms = zeros(2*chopSize(2));
      for ii=1:2*chopSize(2)
      colNorms(ii) = norm(squeeze(newPsf(:,ii,:)));
      end
      blooble2 = zeros(chopSize(2)+1,1);
      for ii=1:chopSize(2)+1
        blooble2(ii) = norm(colNorms(ii:ii+chopSize(2)-1));
      end
      cutStart2 = find(min(blooble2)==blooble2);
      if length(cutStart2) ~= 1 %% (if more than one min is identified)
         cutStart2 = cutStart2(fix(length(cutStart2)/2)+1);
      end
      if blooble2(cutStart2) < .01*norm(newPsf(:))
         newPsf = cat(2,newPsf(:,1:cutStart2-1,:), newPsf(:,cutStart2+chopSize(2):sumSize(2),:));
         newPsf = circshift(newPsf,[0 cutStart2-1 0]);
     else
        newPsf = circshift(newPsf, [0 chopSize(2) 0]);
        disp('Column norms are too large for chopping')
      end
    end %%%   if (chopSize(2) > 0)

    if (chopSize(3) > 0)
      newPsf = circshift(newPsf,[0 0 chopSize(3)]);
      stackNorms = zeros(2*chopSize(3));
      for ii=1:2*chopSize(3)
      stackNorms(ii) = norm(squeeze(newPsf(:,:,ii)));
      end
      blooble3 = zeros(chopSize(3)+1,1);
      for ii=1:chopSize(3)+1
        blooble3(ii) = norm(stackNorms(ii:ii+chopSize(3)-1));
      end
      cutStart3 = find(min(blooble3)==blooble3);
      if length(cutStart3) ~= 1 %% (if more than one min is identified)
         cutStart3 = cutStart3(fix(length(cutStart3)/2)+1);
      end
      if blooble3(cutStart3) < .01*norm(newPsf(:))
         newPsf = cat(3,newPsf(:,:,1:cutStart3-1), newPsf(:,:,cutStart3+chopSize(3):sumSize(3)));
         [ci,cj,ck]=find(max(newPsf(:))==newPsf);
         newPsf = circshift(newPsf,[0 0 cutStart3-1]);
      else
         newPsf = circshift(newPsf,[0 0 chopSize(3)]);
         disp('Stack norms are too large for chopping')
      end
    end %%%   if (chopSize(3) > 0)

    % Now that the chopping is complete, normalize the sum psf
    % and assign it to its compartment in newPsfs
    newPsfs{i,j,k} = newPsf/sum(newPsf(:));
  end
 end
end

% Now build the new psfMatrix from the summed psf image
NewPSF = psf(newPsfs);
S = psfMatrix(NewPSF,boundary);
S=S*2;

repo: IRtools


text:
function g = diffy(f, n)
%DIFFY   Differentiate a CHEBFUN2 with respect to its second argument.
%
%   G = DIFFY(F) returns a CHEBFUN2 representing the derivative of F in its
%   second argument. This is the same as DIFF(F,1,1).
%
%   G = DIFFY(F,N) returns a CHEBFUN2 representing the Nth derivative of F in
%   its second argument. This is the same as DIFF(F,N,1).
%
%   This command is for convenience as the syntax for DIFF, inherited from the
%   DIFF command for matrices, can be confusing.
%
% See also DIFFX, DIFF.

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

% Default to first derivative:
if ( nargin == 1 )
    n = 1;
end

% Call diff:
g = diff(f, n, 1);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Data Quality specification snippet to populate default values
% and verify times are in range
% axs June 2019, ERPLAB
% ams Dec 2022, ERPLAB update: incl'd smaller window for short epoch(<100ms)
function DQ_spec_out = make_DQ_spec(timelimits_ms)

% Check and populate missing args
if exist('timelimits_ms','var') == 0 || isempty(timelimits_ms)

    try
        evalin('base','timelimits = [EEG.xmin EEG.xmax];');
        timelimits_ms = timelimits * 1000;

    catch
        %warning('Couldn''t find time limits for DQ range, trying defaults')
        timelimits_ms = [-200 500];
    end
end

% Prepare time range and time-windows
win_size = 100; %ms
%     If crossing zero, ensure windows are centered on zero
if (min(timelimits_ms) < 0) && (max(timelimits_ms) > 0)
    % if crossing zero, grow 0:-win_size:min and 0:win_size:max

    tw_neg = sort([0:-win_size:timelimits_ms(1)]);
    tw_pos = [0:win_size:timelimits_ms(2)];

    % concatenate
    tw_starts = unique([tw_neg tw_pos]);
    tw_ends = tw_starts + win_size;

else
    % else, grow min:win_size:term

    tw_starts = min(timelimits_ms):win_size:max(timelimits_ms);
    tw_ends = tw_starts + win_size;
end

% remove time-windows out of range
if tw_ends(end) > max(timelimits_ms)
    tw_starts(end) = [];
    tw_ends(end) = [];
end

% for tw less than 100ms
if isempty(tw_starts) & isempty(tw_ends)
    win_size = 2;  %considers typical auditory brainstem responses (ABRs) (eg, 25ms epochs)
                   %thanks Dr. Kelsey Mankel!
    if (min(timelimits_ms) < 0) && (max(timelimits_ms) > 0)
        % if crossing zero, grow 0:-win_size:min and 0:win_size:max

        tw_neg = sort([0:-win_size:timelimits_ms(1)]);
        tw_pos = [0:win_size:timelimits_ms(2)];

        % concatenate
        tw_starts = unique([tw_neg tw_pos]);
        tw_ends = tw_starts + win_size;

    else
        % else, grow min:win_size:term

        tw_starts = min(timelimits_ms):win_size:max(timelimits_ms);
        tw_ends = tw_starts + win_size;
    end

    % remove time-windows out of range
    if tw_ends(end) > max(timelimits_ms)
        tw_starts(end) = [];
        tw_ends(end) = [];
    end

end

% make labels
root_str = 'aSME at ';
for i = 1:numel(tw_starts)
    tw_labels{i} = [root_str num2str(tw_starts(i)) ' to ' num2str(tw_ends(i)) ];
end

root_str = 'aSD Across Trials at ';
for i = 1:numel(tw_starts)
    tw_labels_sd{i} = [root_str num2str(tw_starts(i)) ' to ' num2str(tw_ends(i)) ];
end

tw = [tw_starts; tw_ends]';

% Build spec with these times
DQ_defaults(1).type = 'Baseline Measure - SD';

DQ_defaults(2).type = 'Point-wise SEM';

DQ_defaults(3).type = 'aSME';
DQ_defaults(3).times = tw;
DQ_defaults(3).time_window_labels = tw_labels;
DQ_defaults(3).comments = [];

%add SD across trials (01/05/22)
DQ_defaults(4).type = 'SD Across Trials';
DQ_defaults(4).times = tw;
DQ_defaults(4).time_window_labels = tw_labels_sd;
DQ_defaults(4).comments=[];

DQ_spec_out = DQ_defaults;

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function min_max = bounds2(someVector)
    % bounds2 is get min and max limits of vector.
    % like bounds, but returns a nx2 instead of two nx1 vectors.
    % min_max = bounds2([1 2 3 -1]); % returns [-1 3];
    %
    % see also bounds
    [min_max(1), min_max(2)] = bounds(someVector);
end

author: CelsoReyes


text:
function [Sbus, dSbus_dVm] = makeSbus(baseMVA, bus, gen, mpopt, Vm, Sg)
%MAKESBUS   Builds the vector of complex bus power injections.
%   SBUS = MAKESBUS(BASEMVA, BUS, GEN)
%   SBUS = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM)
%   SBUS = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM, SG)
%   returns the vector of complex bus power injections, that is, generation
%   minus load. Power is expressed in per unit. If the MPOPT and VM arguments
%   are present it evaluates any ZIP loads based on the provided voltage
%   magnitude vector. If VM is empty, it assumes nominal voltage. If SG is
%   provided, it is a complex ng x 1 vector of generator power injections in
%   p.u., and overrides the PG and QG columns in GEN, using GEN only for
%   connectivity information.
%
%   [SBUS, DSBUS_DVM] = MAKESBUS(BASEMVA, BUS, GEN, MPOPT, VM)
%   With two output arguments, it computes the partial derivative of the
%   bus injections with respect to voltage magnitude, leaving the first
%   return value SBUS empty. If VM is empty, it assumes no voltage dependence
%   and returns a sparse zero matrix.
%
%   See also MAKEYBUS.

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into bus, gen matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% default inputs
if nargin < 5
    Vm = [];
    if nargin < 4
        mpopt = [];
    end
end
nb = size(bus, 1);

%% get load parameters
Sd = makeSdzip(baseMVA, bus, mpopt);

if nargout == 2
    Sbus = [];
    if isempty(Vm)
        dSbus_dVm = sparse(nb, nb);
    else
        dSbus_dVm = -(spdiags(Sd.i + 2 * Vm .* Sd.z, 0, nb, nb));
    end
else
    %% compute per-bus generation in p.u.
    on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
    gbus = gen(on, GEN_BUS);                %% what buses are they at?
    ngon = size(on, 1);
    Cg = sparse(gbus, (1:ngon)', 1, nb, ngon);  %% connection matrix
                                                %% element i, j is 1 if
                                                %% gen on(j) at bus i is ON
    if nargin > 5 && ~isempty(Sg)
        Sbusg = Cg * Sg(on);
    else
        Sbusg = Cg * (gen(on, PG) + 1j * gen(on, QG)) / baseMVA;
    end

    %% compute per-bus loads in p.u.
    if isempty(Vm)
        Vm = ones(nb, 1);
    end
    Sbusd = Sd.p + Sd.i .* Vm + Sd.z .* Vm.^2;

    %% form net complex bus power injection vector
    %% (power injected by generators + power injected by loads)
    Sbus = Sbusg - Sbusd;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function f = restrict(f, dom)
% RESTRICT    Restrict the domain of a CHEBFUN3.
%
%   F = RESTRICT(F, DOM) returns a
%   - CHEBFUN3 on the domain DOM that approximates F on that domain if DOM
%     is a nondegenarate cuboid,
%   - CHEBFUN2 that approximates F on the domain if DOM is a plane, and
%   - CHEBFUN that approximates F on the domain if DOM is a line.
%   DOM should be a vector of length 6 specifying the underlying cuboid.

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

if ( isa(dom, 'double') )
    if ( numel(dom) == 6 )                   % Restrict to DOM.
        xlen = diff(dom(1:2));
        ylen = diff(dom(3:4));
        zlen = diff(dom(5:6));

        if ( ( xlen == 0 ) && ( ylen == 0) && ( zlen == 0 ) )
            % DOM is a point.
            f = feval(f, dom(1), dom(3), dom(5));

        elseif ( ( xlen ~= 0 ) && ( ylen == 0 ) && ( zlen == 0 ) )
            % DOM is a vertical line (corresponding to a column).
            cols = restrict(f.cols, dom(1:2));
            rows = feval(f.rows, dom(3));
            tubes = feval(f.tubes, dom(5));
            core = chebfun3.txm(chebfun3.txm(f.core, rows, 2), tubes, 3);
            f = chebfun(cols*core);

        elseif ( ( xlen == 0 ) && ( ylen ~= 0 ) && ( zlen == 0 ) )
            % DOM is a horizontal line (corresponding to a row).
            cols = feval(f.cols, dom(1));
            rows = restrict(f.rows, dom(3:4));
            tubes = feval(f.tubes, dom(5));
            core = chebfun3.txm(chebfun3.txm(f.core, cols, 1), tubes, 3);
            f = chebfun(rows*core.');

        elseif ( ( xlen == 0 ) && ( ylen == 0 ) && ( zlen ~= 0 ) )
            % DOM is an oblique line (corresponding to a tube).
            cols = feval(f.cols, dom(1));
            rows = feval(f.rows, dom(3));
            tubes = restrict(f.tubes, dom(5:6));
            core = squeeze(chebfun3.txm(chebfun3.txm(f.core, cols, 1), ...
                rows, 2));
            f = chebfun(tubes*core);

        elseif ( ( xlen == 0 ) && ( ylen ~= 0 ) && ( zlen ~= 0 ) )
            % DOM is a horizontal plane (corresponding to a horizontal slice).
            cols = feval(f.cols, dom(1));
            rows = restrict(f.rows, dom(3:4));
            tubes = restrict(f.tubes, dom(5:6));
            core = squeeze(chebfun3.txm(f.core, cols, 1));
            f = chebfun2(tubes*core.'*rows.');

        elseif ( ( xlen ~= 0 ) && ( ylen == 0 ) && ( zlen ~= 0 ) )
            % DOM is a lateral plane (corresponding to a lateral slice).
            cols = restrict(f.cols, dom(1:2));
            rows = feval(f.rows, dom(3));
            tubes = restrict(f.tubes, dom(5:6));
            core = squeeze(chebfun3.txm(f.core, rows, 2));
            f = chebfun2(tubes*core.'*cols.');

        elseif ( ( xlen ~= 0 ) && ( ylen ~= 0 ) && ( zlen == 0 ) )
            % DOM is a frontal plane (corresponding to a frontal slice).
            cols = restrict(f.cols, dom(1:2));
            rows = restrict(f.rows, dom(3:4));
            tubes = feval(f.tubes,dom(5));
            core = chebfun3.txm(f.core, tubes, 3);
            f = chebfun2(rows*core.'*cols.');

        else
            % DOM is not degenerate
            f.cols = restrict(f.cols, dom(1:2));
            f.rows = restrict(f.rows, dom(3:4));
            f.tubes = restrict(f.tubes, dom(5:6));
            f.domain = dom;
        end
    else
        error('CHEBFUN:CHEBFUN3:restrict:domain', 'Domain not determined.');
    end

else
    error('CHEBFUN:CHEBFUN3:restrict:domain', 'Unrecognizable domain.');
end

end


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

function y = sprice_2(a, b, r, n, f, k, t,cp)
% sabr prices using an admissible region kl <= k <= ku where sabr is used
% for 0 <= k <= kl
% we use a put pricing function f(x) = K^mu exp(a + bx + cx^2)
% for ku < k < +infty
% we use a call pricing function f(x) = K^(-nu) exp(a + bx^(-1) + cx(-2)
eps = 1e-004;

%find index where density is positive
s = @(x) psabr(a, b, r, n, f, x, t);
index = find(s(k)<0,1,'last');

if isempty(index)
   kl = .25 *f;
else
   kl = max(.25 * f,k(index+1));% lower strike level (free parameter)
end

%kl = .5 * f;                               % lower strike level

s = @(x) sprice(a, b, r, n, f, x, t,-1);    % sabr price (-1) for puts

s1 = s(kl-eps);                             % for calc of derivatives
s2 = s(kl);                                 % for calc of derivatives
s3 = s(kl+eps);                             % for calc of derivatives

V1 = log(s2);                               % log put price

U2 = (s3-s1)/(2*eps);                           % derivative of the price
V2 = U2/s2;                                 % derivative of the log price

U3 = (s3-2*s2+s1)/eps^2;                    % second derivative of the price
V3 = U3/s2 - V2^2;                     % second derivative of the log price

% fix mu and solve
%       V1 = mu log kl + a + b kl + c kl^2
%       V2 = mu / kl + b + 2c kl
%       V3 = - mu / kl^2 + 2c

mu = -V3*kl^2; bl = V2 - mu / kl; al = V1 - mu *log(kl) - bl * kl;

ku = 15 * f;                               % upper strike level
s = @(x) sprice(a, b, r, n, f, x, t,1);     % sabr price (1) for calls
s1 = s(ku-eps);                             % for calc of derivatives
s2 = s(ku);                                 % for calc of derivatives
s3 = s(ku+eps);                             % for calc of derivatives

V1 = log(s2);                               % log call price

U2 = (s3-s1)/(2*eps);                           % derivative of the price
V2 = U2/s2;                                 % derivative of the log price

U3 = (s3-2*s2+s1)/eps^2;                    % second derivative of the price
V3 = U3/s2 - (U2/s2)^2;                     % second derivative f the price

% fix nu and solve
%       V1 = -nu log ku + a + b/ku + c/ku^2
%       V2 = -nu / ku - b / ku^2 - 2c/ku^3
%       V3 = nu / ku^2 + 2 b / ku^3 - 6 c / ku^4

nu = 2; cu = (-1.5*nu / ku + .5*V3 * ku - V2)*ku^3/5;
bu = -ku^2*(V2 + nu/ku +2*cu/ku^3);
au = V1 + nu * log(ku) - bu / ku - cu / ku^2;

% the sabr volatility for the admissible region kl <= k <= ku
sigma = svol(a,b,r,n,f,k((kl<=k)&(k<=ku)),t);

    d1= (log(f./k((kl<=k)&(k<=ku)))+(0.5*t*sigma.*sigma))./(sqrt(t)*sigma); % BS d1 with sabr vol
    d2= (log(f./k((kl<=k)&(k<=ku)))-(0.5*t*sigma.*sigma))./(sqrt(t)*sigma); % BS d2 with sabr vol

if cp==1
    yu = k(k>ku).^(-nu) .* exp(au+bu./k(k>ku)+cu./k(k>ku).^2);              % call price ku < k < +infty
    ym = f.*normcdf( d1,0,1)-k((kl<=k)&(k<=ku)).*normcdf( d2,0,1);          % admissible region kl <= k <= ku
    yl = k(k<kl).^mu .* exp(al + bl.*k(k<kl)) + f - k(k<kl); % call price using call-put parity 0<= k < kl
else
    yl = k(k<kl).^mu .* exp(al + bl .* k(k<kl));               % put price 0<= k < kl
    ym = k((kl<=k)&(k<=ku)).*normcdf(-d2,0,1)-f.*normcdf(-d1,0,1);          % admissible region kl <= k <= ku
    yu = k(k>ku).^(-nu) .* exp(au+bu./k(k>ku)+cu./k(k>ku).^2) - f + k(k>ku);% put prices using call-out parity ku < k < +infty
end

y = [yl ym yu];                                                             % output
end

author: Sable

repo: mcbench-benchmarks


text:
% Visualization of the classifier neurons. Modify the targetclass variable
% to see other classes.
% Please feel free to modify all the switches in line 40 - 46 and width,
% height, scale, rotate variables to see what happens.
% The best visualization can be produced by setting use_image_blur = true
% and use_image_deblur = true, with all other switches are turned to false.
% However, with no constraints, you can still get recognizable images.
% Feel the magic by your self!
%
% This code is messay, I will refactor it in the future.
% Based on paper:
% Feng Wang, Haijun Liu, Jian Cheng,
% Visualizing Deep Neural Network by Alternately Image Blurring and Deblurring
caffe.reset_all();
caffe.set_mode_gpu();
gpu_id = 0;  % we will use the first gpu in this demo
caffe.set_device(gpu_id);

% net_model = 'D:\deeplearning\caffe-windows\examples\GoogLeNet\googlenet_deploy_prob3.prototxt';
% net_weights = 'D:\deeplearning\caffe-windows\examples\GoogLeNet\imagenet_googlenet.caffemodel';
net_model ='D:\deepLearning\caffe-windows\models\bvlc_googlenet\deploy.prototxt';
net_weights = 'D:\deepLearning\caffe-windows\models\bvlc_googlenet\thinned.caffemodel';
% mean_file = 'D:\deeplearning\caffe-windows\examples\GoogLeNet\imagenet_mean.binaryproto';

% net_model = 'D:\deeplearning\caffe-windows\examples\PlaceCNN\hybridCNN_deploy_upgraded.prototxt';
% net_weights = 'D:\deeplearning\caffe-windows\examples\PlaceCNN\hybridCNN_iter_700000_upgraded.caffemodel';
% mean_file = 'D:\deeplearning\caffe-windows\examples\PlaceCNN\hybridCNN_mean.binaryproto';

% net_model = 'D:\deeplearning\caffe-windows\examples\VGG\VGG_ILSVRC_16_layers_deploy.prototxt';
% net_weights = 'D:\deeplearning\caffe-windows\examples\VGG\VGG_ILSVRC_16_layers.caffemodel';
% mean_file = [];
% vgg_mean =  [103.939, 116.779, 123.68];
% net_model = 'D:\deeplearning\caffe-windows\models\bvlc_googlenet\deploy.prototxt';
% net_weights = 'D:\deeplearning\caffe-windows\models\bvlc_googlenet\bvlc_googlenet.caffemodel';
mean_file = [];
vgg_mean =  [103.939, 116.779, 123.68];

%%%%%%%%%extract the train features
train_net = caffe.Net(net_model,net_weights,'test');
use_clip = false;
use_cv_norm = false;
use_weight_decay = false;
use_image_blur = true;
use_image_deblur = true;
use_gradient_blur = false;
use_dropout = false;

for targetclass = 13:13
%     if exist(['gallery\' num2str(targetclass) '.png'],'file')
%         continue;
%     end;
if ~isempty(mean_file)
    mean_image = caffe.read_mean(mean_file);
    mean_image = mean_image(16:242,16:242,:);
else
    mean_image = permute(repmat(vgg_mean',[1,224,224]),[2,3,1]);
end;
% mean_image = mean_image + randn(size(mean_image));
input_data = zeros(size(mean_image,1), size(mean_image,2), 3, 1, 'single');

width = 300;
height = 300;
meanV = mean(mean_image(:));
BigImage = permute(repmat(vgg_mean',[1,height,width]),[2,3,1]);%this image will be transposed
H = fspecial('gaussian',[7 7],1.2);

if length(targetclass) > 1
    sample_density = 2;
    class_center = zeros(sample_density*length(targetclass),2);
    for i = 1:sample_density*length(targetclass)
        class_center(i,1) = randi(height - size(mean_image,1),1,1) + size(mean_image,1) / 2;
        class_center(i,2) = randi(width - size(mean_image,2),1,1) + size(mean_image,2) / 2;
        if i > 1
            min_distance = min(pdist2(class_center(i,:),class_center(1:i-1,:)));
            while min_distance < size(mean_image,1) / 2
                class_center(i,1) = randi(height - size(mean_image,1),1,1) + size(mean_image,1) / 2;
                class_center(i,2) = randi(width - size(mean_image,2),1,1) + size(mean_image,2) / 2;
                min_distance = min(pdist2(class_center(i,:),class_center(1:i-1,:)));
            end;
        end;
    end;

    class_label = repmat(targetclass,1,sample_density);
end;

blur_data = zeros(size(input_data));
base_lr = 5000;
max_lr = 5000;
lambda1 = 0.01;
lambda2 = 0.1;
lambda3 = 300;
last_prob = -999;
momentum = 0.8;
lastgrad = zeros(size(mean_image));
mask = ones(size(mean_image,1), size(mean_image,2));
iter = 1;
dropout = 0.5;
sep_iter = zeros(ceil((height - size(mean_image,1)) / size(mean_image,1) *2),ceil((width - size(mean_image,2)) / size(mean_image,2) * 2));
erode = 5;
start_sigma = 1.2;
end_sigma = 0.5;
this_sigma = start_sigma;
sigma_decay = 0.99;
deblur_sigma_ratio = 0.8;
scale = [0.8 1.2];
rotate = [0 0] / 180 * pi;
func_loss = [];
BigGrad = zeros(size(BigImage));

while 1
%     for i = 0:size(mean_image,1) / 2:height - size(mean_image,1)
%         for j = 0:size(mean_image,2) / 2:width - size(mean_image,2)
            this_sigma = this_sigma * sigma_decay;
            if this_sigma < end_sigma
                this_sigma = end_sigma;
            end;
            s = rand(1) * (scale(2)-scale(1)) + scale(1);
            r = rand(1) * (rotate(2)-rotate(1)) + rotate(1);
            if length(targetclass)>1 && iter < 500
                idx = randi(size(class_center,1),1,1);
                x = class_center(idx,1) - size(mean_image,1) / 2;
                y = class_center(idx,2) - size(mean_image,2) / 2;
            else
                if rand() > size(mean_image,1) / height
                    x = randi(height - floor(size(mean_image,1) * s),1,1);% + size(mean_image,1) / 2;
                else
                    pos = randi(2,1,1)-1;
                    x = pos*(height - floor(size(mean_image,1) * s));
                end;
                if rand() > size(mean_image,2) / width
                    y = randi(width - floor(size(mean_image,2) * s),1,1);% - size(mean_image,2) / 2;
                else
                    pos = randi(2,1,1)-1;
                    y = pos*(width - floor(size(mean_image,2) * s));
                end;
            end;

            II = zeros(floor(size(mean_image,1) * s),floor(size(mean_image,2) * s),3);
            II(:,:,1) = BigImage(x + 1:x + floor(size(mean_image,1) * s),y + 1:y + floor(size(mean_image,2) * s),1) - vgg_mean(1);
            II(:,:,2) = BigImage(x + 1:x + floor(size(mean_image,1) * s),y + 1:y + floor(size(mean_image,2) * s),2) - vgg_mean(2);
            II(:,:,3) = BigImage(x + 1:x + floor(size(mean_image,1) * s),y + 1:y + floor(size(mean_image,2) * s),3) - vgg_mean(3);
            input_data(:,:,:,1) = imresize(imrotate(II,r),[size(mean_image,1), size(mean_image,2)]);
            GG = BigGrad(x + 1:x + floor(size(mean_image,1) * s),y + 1:y + floor(size(mean_image,2) * s),:);
            last_grad = imresize(GG,[size(mean_image,1), size(mean_image,2)]);
            for k = 1 : 10

                lr = base_lr;% * sqrt(this_prob / (1 - this_prob));

                prob = train_net.forward({input_data});

                    this_prob = prob{end}(targetclass);
                [max_prob,max_idx] = max(this_prob);
                fprintf('iter=%d,lr=%f,max_idx=%d,prob1=%f,last_prob=%f,this_sigma=%f\n',iter,lr,max_idx,this_prob(max_idx),last_prob,this_sigma);
                func_loss = [func_loss; this_prob(max_idx)];

                if k==1
                    back_data = ones(size(prob{1}),'single') * -1;
%                     back_data = zeros(size(prob{end}),'single');
                    if length(targetclass) >1
                        center_x = x + size(mean_image,1) / 2;
                        center_y = y + size(mean_image,2) / 2;
                        distances = pdist2([center_x center_y],class_center);
                        distances = exp(-(distances / size(mean_image,1) * 4).^2);
                        distances = distances / sum(distances);
                        distances = reshape(distances,length(targetclass),sample_density);
                        distances = sum(distances,2);
                        distances = distances * 2 - 1;
                        back_data(targetclass) = distances';
                    else
                        back_data(targetclass(max_idx)) = 1;
                    end;

                    back_cell = cell(length(prob),1);
%                     back_cell{1} = zeros(size(back_data));
%                     back_cell{2} = zeros(size(back_data));
%                     back_cell{1} = back_data;
%                     back_cell{2} = back_data;
                    back_cell{1} = back_data;
                end;

                iter = iter + 1;

%                 sep_iter(floor(i / size(mean_image,1) * 2)+1,floor(j/size(mean_image,1)*2) + 1) = sep_iter(floor(i / size(mean_image,1) * 2)+1,floor(j/size(mean_image,1)*2) + 1) + 1;

                if max_prob<last_prob
                    base_lr = base_lr * 0.99;
            %         input_data = bak_data;
                end;
                if max_prob>last_prob&&base_lr<max_lr%&& (this_prob-last_prob) / this_prob < 0.001
                    base_lr = base_lr * 1.01;
            %         input_data = bak_data;
                end;
            %     if this_prob>last_prob
                    last_prob = max_prob;
            %     end;

                res = train_net.backward(back_cell);

                bak_data = input_data;

                if use_clip
                    app_gradient = sum(abs(res{1} .* input_data(:,:,:,1)),3);
                    app_gradient = app_gradient < mean(app_gradient(:));
                    grad = reshape(res{1},[size(mean_image,1)*size(mean_image,2) 3]);
                    grad(app_gradient==1,:) = 0;
                    grad = reshape(grad,size(input_data));
                    res{1} = grad;
                end;

                input_data(:,:,:,1) = input_data(:,:,:,1) + res{1} * 1.5 / mean(abs(res{1}(:)));

                if use_cv_norm
                    I = input_data(:,:,:,1);
                    Gx = smoothL1(I(2:end-1,:,:) - I(1:end-2,:,:)) - smoothL1(I(3:end,:,:) - I(2:end-1,:,:));
                    Gx = [smoothL1(I(1,:,:) - I(2,:,:)); Gx; smoothL1(I(end,:,:) - I(end-1,:,:))];
                    Gy = smoothL1(I(:,2:end-1,:) - I(:,1:end-2,:)) - smoothL1(I(:,3:end,:) - I(:,2:end-1,:));
                    Gy = [smoothL1(I(:,1,:) - I(:,2,:)) Gy smoothL1(I(:,end,:) - I(:,end-1,:))];
                    input_data(:,:,:,1) = input_data(:,:,:,1) - lr * lambda2 * (Gx + Gy);
                end;
                if use_weight_decay
                    input_data(:,:,:,1) = input_data(:,:,:,1) - lambda1 * input_data(:,:,:,1);
                end;

                if k==1
                    H = fspecial('gaussian',[randi(6)+4 7],this_sigma);
                    if use_image_blur
                        input_data(:,:,:,1) = imfilter(input_data(:,:,:,1),H,'same');
                    end
                end;
                if k==6
                    H = fspecial('gaussian',[randi(6)+4 7],this_sigma * deblur_sigma_ratio);
                    if use_image_deblur
                        input_data(:,:,:,1) = deconvlucy(input_data(:,:,:,1), H);
                    end;
                end;

                if lr<0.000001
                    break;
                end;
            end;
            II = imresize(imrotate(input_data(:,:,:,1),-r),floor([size(mean_image,1), size(mean_image,2)] * s));
            BigImage(x + 1+erode:x + floor(size(mean_image,1) * s)-erode,y + 1+erode:y + floor(size(mean_image,2) * s)-erode,1) = II(erode+1:end-erode,erode+1:end-erode,1) + vgg_mean(1);
            BigImage(x + 1+erode:x + floor(size(mean_image,1) * s)-erode,y + 1+erode:y + floor(size(mean_image,2) * s)-erode,2) = II(erode+1:end-erode,erode+1:end-erode,2) + vgg_mean(2);
            BigImage(x + 1+erode:x + floor(size(mean_image,1) * s)-erode,y + 1+erode:y + floor(size(mean_image,2) * s)-erode,3) = II(erode+1:end-erode,erode+1:end-erode,3) + vgg_mean(3);
            GG = imresize(last_grad,floor([size(mean_image,1), size(mean_image,2)] * s));
            BigGrad(x + 1+erode:x + floor(size(mean_image,1) * s)-erode,y + 1+erode:y + floor(size(mean_image,2) * s)-erode,:) = GG(erode+1:end-erode,erode+1:end-erode,:);
%             sort_data = sort(BigImage(:));
%             BigImage = BigImage .* (255 / sort_data(floor(length(sort_data) * 0.999)));
            if mod(iter,100)==1
                figure(2);
                % imshow(uint8(mean_image + input_data));
                output = BigImage;
                output = output(:, :, [3, 2, 1]);
                output = permute(output, [2 1 3]);
                output = output(erode+1:end-erode,erode+1:end-erode,:);
                sort_output = sort(output(:));
                output = output .* (255 / sort_output(floor(length(sort_output) * 0.99)));
                imshow(uint8(output));
                figure(3);
%                 hist(output(:),1000);
                plot(1:iter-1, func_loss);
            end;
            if iter > 1500
                break;
            end;
%         end;
%     end;
end;
LastBigImage = BigImage;
imwrite(uint8(output),['gallery\' num2str(targetclass) '.png']);
end;

author: happynear

repo: DeepVisualization


text:
function Cpos = ShiftBaselinePosSignal(C,dim)
% for 2D arrays only.
% 'dim' (dimension along which the signal is processed) can be 1 or 2

if ~exist('dim','var')
    [s1,s2] = size(C);
    if s1>s2
        dim = 1;
    else
        dim = 2;
    end
end

Cpos = zeros(size(C));
if dim==2
    C = C';
end
for i = 1:size(C,2)
    x = C(:,i);
    shift = prctile(x,5);
    Cpos(:,i) = C(:,i)-shift;
end
if dim==2
    Cpos = Cpos';
end

end

author: xiuyechen

repo: FishExplorer


text:
function Offspring = Operator(Problem,Particle,Pbest,Gbest,W)
% The particle swarm optimization in dMOPSO

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

    %% Parameter setting
    if nargin < 5
        W = 0.4;
    end
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    r1     = repmat(rand(N,1),1,D);
    r2     = repmat(rand(N,1),1,D);
    OffVel = W.*ParticleVel + r1.*(PbestDec-ParticleDec) + r2.*(GbestDec-ParticleDec);
    OffDec = ParticleDec + OffVel;

    %% Deterministic back
    repair = OffDec < repmat(Problem.lower,N,1) | OffDec > repmat(Problem.upper,N,1);
    OffVel(repair) = -OffVel(repair);
    Offspring = Problem.Evaluation(OffDec,OffVel);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ y, m, d, f ] = jed_to_ymdf_gregorian2 ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_GREGORIAN2 converts a JED to a Gregorian YMDF date.
%
%  Discussion:
%
%    The theory behind this routine is very clean.  The Gregorian
%    calendar has cycles of 1, 4, 100 and 400 years, and we can
%    analyze a date by determining where it lies within these cycles.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Reingold, Nachum Dershowitz, Stewart Clamen,
%    Calendrical Calculations, II: Three Historical Calendars,
%    Software - Practice and Experience,
%    Volume 23, Number 4, pages 383-404, April 1993.
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F, the YMDF date.
%
  g1 = 365;
  g4 = 1461;
  g100 = 36524;
  g400 = 146097;

  jed_epoch = epoch_to_jed_gregorian ( );

  j = floor ( jed - jed_epoch );
  f1 = ( jed - jed_epoch ) - j;

  d0 = j;
  n400 = 0;

  while ( d0 < 0 )
    d0 = d0 + g400;
    n400 = n400 - 1;
  end

  n400 = n400 + floor ( d0 / g400 );
  d1 = i4_modp ( d0, g400 );

  n100 = floor ( d1 / g100 );
  d2 = i4_modp ( d1, g100 );

  n4 = floor ( d2 / g4 );
  d3 = i4_modp ( d2, g4 );

  n1 = floor ( d3 / g1 );
  d4 = i4_modp ( d3, g1 );

  if ( n100 == 4 || n1 == 4 )
    j1 = 366;
    y1 = 400 * n400 + 100 * n100 + 4 * n4 + n1;
  else
    j1 = d4 + 1;
    y1 = 400 * n400 + 100 * n100 + 4 * n4 + n1 + 1;
  end
%
%  Any year before 1 AD must be moved one year further back, since
%  this calendar does not include a year 0.
%
  y1 = y_astronomical_to_common ( y1 );

  [ y, m, d, f ] = yjf_to_ymdf_gregorian ( y1, j1, f1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [tfr,t,f] = tfrbud(x,t,N,g,h,sigma,trace);
%TFRBUD	Butterworth time-frequency distribution.
%	[TFR,T,F]=TFRBUD(X,T,N,G,H,SIGMA,TRACE) computes the Butterworth
%	distribution of a discrete-time signal X, or the
%	cross Butterworth representation between two signals.
%
%	X     : signal if auto-BUD, or [X1,X2] if cross-BUD.
%	T     : time instant(s)          (default : 1:length(X)).
%	N     : number of frequency bins (default : length(X)).
%	G     : time smoothing window, G(0) being forced to 1.
%	                                 (default : Hamming(N/4)).
%	H     : frequency smoothing window, H(0) being forced to 1.
%	                                 (default : Hamming(N/4)).
%	SIGMA : kernel width             (default : 1).
%	TRACE : if nonzero, the progression of the algorithm is shown
%                                        (default : 0).
%	TFR   : time-frequency representation. When called without
%               output arguments, TFRBUD runs TFRQVIEW.
%	F     : vector of normalized frequencies.
%
%	Example :
%	 sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4);
%	 g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser');
%	 t=1:128; tfrbud(sig,t,128,g,h,3.6,1);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

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

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
end

if (nargin <= 2),
 N=xrow;
elseif (N<0),
 error('N must be greater than zero');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
end;

hlength=floor(N/4); hlength=hlength+1-rem(hlength,2);
glength=floor(N/10);glength=glength+1-rem(glength,2);

if (nargin == 1),
 t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 2)|(nargin == 3),
 g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 4),
 h = tftb_window(hlength); sigma = 1.0; trace = 0;
elseif (nargin == 5),
 sigma = 1.0; trace = 0;
elseif (nargin == 6),
 trace = 0;
end;

[trow,tcol] = size(t);
if (trow~=1),
 error('t must only have one row');
end;

[grow,gcol]=size(g); Lg=(grow-1)/2;
if (gcol~=1)|(rem(grow,2)==0),
 error('G must be a smoothing window with odd length');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1);
if (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
end;

if (sigma<=0.0),
 error('SIGMA must be strictly positive');
end;

taumax = min([round(N/2),Lh]); tau = 1:taumax; points = -Lg:Lg;
BudKer = exp(-kron( abs(points.'), 1.0 ./ (2.0*tau/sqrt(sigma))));
BudKer = diag(g) * BudKer;

tfr= zeros (N,tcol) ;
if trace, disp('Butterworth distribution'); end;
for icol=1:tcol,
 ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]);
 if trace, disprog(icol,tcol,10); end;
 tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol));

 for tau=1:taumax,
  points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]);
  g2 = BudKer(Lg+1+points,tau); g2=g2/sum(g2);
  R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)));
  tfr(  1+tau,icol)=h(Lh+tau+1)*R;
  R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)));
  tfr(N+1-tau,icol)=h(Lh-tau+1)*R;
 end;

 tau=round(N/2);
 if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh),
  points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]);
  g2 = BudKer(Lg+1+points,tau); g2=g2/sum(g2);
  tfr(tau+1,icol) = 0.5 * ...
   (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+...
    h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))));
 end;
end;

clear BudKer;

if trace, fprintf('\n'); end;
tfr= fft(tfr);
if (xcol==1), tfr=real(tfr); end ;

if (nargout==0),
 tfrqview(tfr,x,t,'tfrbud',g,h,sigma);
elseif (nargout==3),
 f=(0.5*(0:N-1)/N)';
end;

author: HeLiangHIT


text:
function f = times(f, g, varargin)
%.*   CLASSICFUN multiplication.
%   F.*G multiplies CLASSICFUN objects F and G or a CLASSICFUN by a scalar if
%   either F or G is a scalar.
%
%   If F is an array-valued CLASSICFUN, then F.*C is supported if C is a row
%   vector of doubles with the same number of columns as F. Similarly if F if a
%   row vector of doubles and C is an array-valued CLASSICFUN
%
%   If F and G are both CLASSICFUN objects, they are assumed to have the same
%   domain and, if they are array-valued, the same number of columns. The
%   method gives no warning if their domains don't agree, but the output of the
%   method will be meaningless.
%
% See also MTIMES, RDIVIDE.

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

if ( isempty(f) || isempty(g) ) % CLASSICFUN * [] = []

    % Return empty:
    f = [];

elseif ( ~isa(f, 'classicfun') )       % double * CLASSICFUN

    % If f is not a CLASSICFUN, g must be. Call ONEFUN/TIMES():
    g.onefun = times(f, g.onefun, varargin{:});
    % Swap arguments for output variable:
    f = g;

elseif ( isa(g, 'classicfun'))         % CLASSICFUN * CLASSICFUN

    % Multiply the ONEFUN objects of f and g together.
    f.onefun = times(f.onefun, g.onefun, varargin{:});

elseif ( isa(g, 'double') )     % CLASSICFUN * double

    % Multiply the ONEFUN of f with g.
    f.onefun = times(f.onefun, g, varargin{:});

end

end


text:
% DEMOILLLE4 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('oil');

options = lleOptions(32, 2);
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model);

lvmScatterPlot(model, lbls);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, 'Oil', 4);
end

errors = lvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function out = isinf(f)
%ISINF   Test if a CHEBFUN is infinite.
%   ISINF(F) returns TRUE if F has any infinite values and FALSE otherwise.

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

out = ~isfinite(f);

end


text:
function [hr] = us2hr(us)
% Convert time from microseconds to hours.
% Chad Greene 2012
hr = us*2.777777777778e-10 ;


text:
function ye = fL(t,D,nalpha,nbeta,u,h)		% Left boundary value
alpha = nalpha*pi; beta = nbeta*pi;
ye = h*cos(beta*u*t) + exp(-D*alpha*alpha*t)*cos(alpha*u*t);

repo: Aero-matlab


text:
function errors = test_jtv_filter()
%TEST_JTV_FILTER This function test all the time-vertex filters
close all
errors = 0;
warning on

errors = errors + test_jtv_graph();
errors = errors + test_jtv_graph_diagonalization();

errors = errors + test_diffusion();
errors = errors + test_wave();
errors = errors + test_jft_ijft();
errors = errors + test_swap_transform();

errors = errors + test_theory()
errors = errors + test_swap_localization();

errors = errors + test_jtv_cheb_coeff();
errors = errors + test_jtv_cheb_op();

errors = errors + test_jtv_filter_evaluate();
errors = errors + test_jtv_filter_array();

errors = errors + test_jtv_filter_analysis();
errors = errors + test_jtv_filter_synthesis();

errors = errors + test_jtv_matrix_analysis_op();
errors = errors + test_jtv_matrix_synthesis_op();

errors = errors + test_jtv_filter_dual();
errors = errors + test_jtv_filter_inverse();

try  %#ok<TRYNC>
    close(100)
end

end

function errors = test_jtv_graph()
errors=0;

try

    N = 100;
    T = 200;
    fs = 1;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T,fs);

    fprintf('Test JTV - timevertex graph:  OK\n');
catch
    errors = errors + 1;
    warning('Test JTV - Error in timevertex graph test')
end

end

function errors = test_jtv_graph_diagonalization()
errors  = 0;
T=30;
N=20;
G = gsp_sensor(N);
G = gsp_compute_fourier_basis(G);
DFT = dftmtx(T)';

Gt = gsp_ring(T);
Gp = gsp_graph_product(G,Gt);
Gp = gsp_create_laplacian(Gp,'combinatorial');
Gp = gsp_compute_fourier_basis(Gp);

U = kron(G.U,DFT);

errors = errors + gsp_assert_test(U'*Gp.L*U,diag(diag(U'*Gp.L*U)),1e-10,'JTV - diagonalization test');

end

function errors = test_diffusion()

errors = 0;
try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);
    [g,ft] = gsp_jtv_design_diffusion(G);
    gsp_plot_jtv_filter(G,g,ft);
    close

    fprintf('JTV: diffusion kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error diffusion kernel 1 test')
end

try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);
    tau=.1;
    [g,ft] = gsp_jtv_design_diffusion(G,tau);
    gsp_plot_jtv_filter(G,g,ft);
    close

    fprintf('JTV: diffusion kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error diffusion kernel 2 test')
end

try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);

    tau=[.01 0.1 1 2];
    param.normalize=1;
    param.show_sum=0;
    param.title=num2cell(tau);
    [g,ft] = gsp_jtv_design_diffusion(G,tau,param);

    gsp_plot_jtv_filter(G,g,ft);
    close

    fprintf('JTV: diffusion kernel 3 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error diffusion kernel 3 test')
end

end

function errors = test_wave()

errors = 0;
try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);

    [g,ft] = gsp_jtv_design_wave(G);
    gsp_plot_jtv_filter(G,g,ft);

    close

    fprintf('JTV: wave kernel 1 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error wave kernel 1 test')
end

try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);

    alpha=1;
    [g,ft] = gsp_jtv_design_wave(G,alpha);
    gsp_plot_jtv_filter(G,g,ft);

    close

    fprintf('JTV: wave kernel 2 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error wave kernel 2 test')
end

try

    N = 100;
    T = 200;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);

    alpha=1;
    param.normalize=1;
    [g,ft] = gsp_jtv_design_wave(G,alpha,param);
    gsp_plot_jtv_filter(G,g,ft);

    close

    fprintf('JTV: wave kernel 3 ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error wave kernel 3 test')
end

end

function errors = test_jft_ijft()
errors = 0;

try

    N=50;
    T=100;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_compute_fourier_basis(G);

    x = randn(N,T);

    xhat = gsp_jft(G,x);

catch
    errors = errors + 1;
    warning('JTV: Error jtgft test')
end

try

    N=50;
    T=100;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_compute_fourier_basis(G);

    xhat = randn(N,T);
    x = gsp_ijft(G,x);

catch
    errors = errors + 1;
    warning('JTV: Error ijtgft test')
end

N=50;
T=100;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

x = randn(N,T);

x2 = gsp_ijft(G,gsp_jft(G,x));

errors = errors + gsp_assert_test(x,x2,eps(1000),'JTV - inverse jtv Fourier transform');

x = rand(N,T,1);

s1 = gsp_jft(G,x);
s2 = gsp_jft_simple(G,x);

errors = errors + gsp_assert_test(s1,s2,eps(1000),'JTV - JFT 1 dim');

x = rand(N,T,10);

s1 = gsp_jft(G,x);
s2 = gsp_jft_simple(G,x);

errors = errors + gsp_assert_test(s1,s2,eps(1000),'JTV - JFT 2 dim');

end

function errors = test_swap_transform()
errors = 0;

N=50;
T=100;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

x = randn(N,T);

x1 = fft(gsp_gft(G,x),[],2)/sqrt(T);
x2 = gsp_gft(G,fft(x,[],2))/sqrt(T);
x3 = gsp_jft(G,x);

e1 = norm(x1-x2,'fro')/norm(x1,'fro');
e2 = norm(x1-x3,'fro')/norm(x1,'fro');
e3 = norm(x2-x3,'fro')/norm(x2,'fro');

errors = errors + gsp_assert_test(0,mean([e1 e2 e3]),eps(1000),'JTV - commutative fft gft transform');

%

xhat = gsp_jft(G,x);

x1 = ifft(gsp_igft(G,xhat),[],2)*sqrt(T);
x2 = gsp_igft(G,ifft(xhat,[],2))*sqrt(T);
x3 = gsp_ijft(G,xhat);

e1 = norm(x1-x2,'fro')/norm(x1,'fro');
e2 = norm(x1-x3,'fro')/norm(x1,'fro');
e3 = norm(x2-x3,'fro')/norm(x2,'fro');

errors = errors + gsp_assert_test(0,mean([e1 e2 e3]),eps(1000),'JTV - commutative ifft igft transform');

end

function errors = test_swap_localization()
errors = 0;

N = 100;
T = 50;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

x = rand(N,T);
vertex=randi(N);
time=randi(T);

% to be updated (with gsp_localize?)
error('Does the boundary still make sense?')
param.boundary='periodic';
x1p = gsp_time_translate(G,gsp_translate_old(G,x,vertex),time,param);
x2p = gsp_translate_old(G,gsp_time_translate(G,x,time,param),vertex);
param.boundary='symmetric';
x1s = gsp_time_translate(G,gsp_translate_old(G,x,vertex),time,param);
x2s = gsp_translate_old(G,gsp_time_translate(G,x,time,param),vertex);
param.boundary='absorbing';
x1a = gsp_time_translate(G,gsp_translate_old(G,x,vertex),time,param);
x2a = gsp_translate_old(G,gsp_time_translate(G,x,time,param),vertex);

n1=norm(x1p-x2p)/norm(x1p);
n2=norm(x1s-x2s)/norm(x1s);
n3=norm(x1a-x2a)/norm(x1a);

if  (n1+n2+n3)/3< eps(1000)
    fprintf('JTV: commutative translation operator ok\n');

else
    errors = errors + 1;
    warning('JTV: error commutative translation operator test')
end

end

function errors = test_jtv_cheb_coeff()
errors = 0;

try
    N = 100;
    T = 50;
    alpha = 1;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);
    [g,ft] = gsp_jtv_design_wave(G, alpha);
    c = gsp_jtv_cheby_coeff(G, g,ft);

    fprintf('FILTER: jtv cheb coeff 1 ok\n');

catch
    errors = 1;
    warning('FILTER: Error in jtv cheb coeff 1 test')
end

try
    N = 100;
    T = 50;
    alpha = [ 0.1 0.5 1 ];
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);
    [g,ft] = gsp_jtv_design_wave(G, alpha);
    c = gsp_jtv_cheby_coeff(G, g,ft);

    fprintf('FILTER: jtv cheb coeff 2 ok\n');

catch
    errors = 1;
    warning('FILTER: Error in jtv cheb coeff 2 test')
end

end

function errors = test_jtv_cheb_op()

errors = 0;
try
    N = 100;
    T = 50;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_estimate_lmax(G);
    alpha = 1;
    [g,ft] = gsp_jtv_design_wave(G, alpha);
    c = gsp_jtv_cheby_coeff(G, g,ft);
    f = randn(N,T);
    r = gsp_jtv_cheby_op(G, c, f);

    fprintf('FILTER: cheb op ok\n');
catch
    errors = 1;
    warning('FILTER: Error in cheb op test')
end

end

function errors = test_jtv_filter_evaluate()
errors = 0;

N = 100;
T = 100;
G = gsp_sensor(N);
G = gsp_compute_fourier_basis(G);
G = gsp_jtv_graph(G,T);

alpha=1;
[g,ft] = gsp_jtv_design_wave(G,alpha);
x = G.e;
param.filtertype = ft;
t = gsp_jtv_ta(G);
f = gsp_jtv_fa(G);

x1 = gsp_jtv_filter_evaluate_simple(G,g,x,param);

x2 = gsp_jtv_filter_evaluate(g,ft,x,t,param);

errors = errors + gsp_assert_test(x1,x2,1e-12,'JTV filter evaluate 1');

g = @(x,y) x+y;
ft = 'js';
param.filtertype = ft;
x1 = gsp_jtv_filter_evaluate_simple(G,g,x,param);

x2 = gsp_jtv_filter_evaluate(g,ft,x,f);

errors = errors + gsp_assert_test(x1,x2,1e-12,'JTV filter evaluate 2');

end

function errors = test_jtv_filter_array()
errors = 0;

N = 500;
T = 100;
G = gsp_sensor(N);
G = gsp_compute_fourier_basis(G);
G = gsp_jtv_graph(G,T);

alpha=[0 0.5 0.7 1];
[g,ftg] = gsp_jtv_design_wave(G,alpha);
x = G.e;
t = gsp_jtv_ta(G);

[h,fth] = gsp_jtv_filter_array(G,g,ftg);

x1 = gsp_jtv_filter_evaluate(g,ftg,x,t);
x2 = gsp_jtv_filter_evaluate(h,fth,x,t);

errors = errors + gsp_assert_test(x1,x2,1e-12,'JTV filter array 1');

g = @(x,y) exp((x-y.^2)/0.1);
ftg = 'js';
f = gsp_jtv_fa(G);

[h,fth] = gsp_jtv_filter_array(G,g,ftg);

x1 = gsp_jtv_filter_evaluate(g,ftg,x,f);
x2 = gsp_jtv_filter_evaluate(h,fth,x,f);

errors = errors + gsp_assert_test(x1,x2,1e-12,'JTV filter array 2');

end

function errors = test_jtv_filter_analysis()
%%
errors = 0;

try
    N = 100;
    T = 500;
    G = gsp_sensor(N);
    G = gsp_jtv_graph(G,T);
    G = gsp_compute_fourier_basis(G);

    alpha=[0.1:0.1:0.7];
    [g,ft] = gsp_jtv_design_wave(G,alpha);

    x = rand(N,T);

    param.method = 'exact';
    coeff = gsp_jtv_filter_analysis(G,g,ft,x,param);

    fprintf('JTV: analysis ok\n');
catch
    errors = errors + 1;
    warning('JTV: Error analysis test')
end

N = 1000;
T = 100;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_estimate_lmax(G);

alpha=[0.1:0.1:0.7];
[g,ft] = gsp_jtv_design_wave(G,alpha);

x = randn(N,T);

param.method = 'exact';
tic
G = gsp_compute_fourier_basis(G);
c_exact = gsp_jtv_filter_analysis(G,g,ft,x,param);
param.filtertype = ft;
c_exact2 = gsp_jtv_filter_analysis_simple(G,g,x,param);
time=toc;
fprintf(['Eigendecomposition + Exact filtering time: ' num2str(time) '\n'])

errors = errors + gsp_assert_test(c_exact,c_exact2,1e-10,'JTV - exact analysis 1');

param.method = 'cheby';
param.order=40;
tic
c_cheby = gsp_jtv_filter_analysis(G,g,ft,x,param);
time=toc;
fprintf(['Cheby filtering time: ' num2str(time) '\n'])

errors = errors + gsp_assert_test(c_exact,c_cheby,1e-4,'JTV - cheby analysis 1');

param.method = 'exact';
g = @(x,y) exp(-x.*(y.^2)/0.1);
ft = 'js';
c_exact = gsp_jtv_filter_analysis(G,g,ft,x,param);

param.filtertype = ft;
c_exact2 = gsp_jtv_filter_analysis_simple(G,g,x,param);
errors = errors + gsp_assert_test(c_exact,c_exact2,1e-10,'JTV - exact analysis 2');

param.method = 'cheby';
param.order=20;
c_cheby = gsp_jtv_filter_analysis(G,g,ft,x,param);

errors = errors + gsp_assert_test(c_exact,c_cheby,1e-4,'JTV - cheby analysis 2');

end

function errors = test_jtv_filter_synthesis()
errors = 0;

N = 500;
T = 300;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

alpha=[0.1:0.1:0.7];
[g,ft] = gsp_jtv_design_wave(G,alpha);

param.lag=1;
x = rand(N,T,numel(alpha));
param.method = 'exact';
s = gsp_jtv_filter_synthesis(G,g,ft,x,param);

param.filtertype = ft;
s2 = gsp_jtv_filter_synthesis_simple(G,g,x,param);

errors = errors + gsp_assert_test(s,s2,1e-10,'JTV  - exact synthesis 1');

param.method = 'cheby';
param.order=200;
s2 = gsp_jtv_filter_synthesis(G,g,ft,x,param);

errors = errors + gsp_assert_test(s,s2,1e-10,'JTV - cheby synthesis 1');

param.method = 'exact';
g = @(x,y) exp(-x.*(y.^2)/0.1);
ft = 'js';
x = rand(N,T,1,3);
s = gsp_jtv_filter_synthesis(G,g,ft,x,param);

param.filtertype = ft;
s2 = gsp_jtv_filter_synthesis_simple(G,g,x,param);
errors = errors + gsp_assert_test(s,s2,1e-10,'JTV  - exact synthesis 2');

param.method = 'cheby';
param.order=20;
s2 = gsp_jtv_filter_synthesis(G,g,ft,x,param);
time=toc;

errors = errors + gsp_assert_test(s,s2,1e-4,'JTV - cheby analysis 2');

end

function errors = test_jtv_matrix_analysis_op()
errors = 0;

N = 20;
T = 20;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

alpha = 1;
[g,ft] = gsp_jtv_design_wave(G,alpha);

F = gsp_jtv_compute_frame(G,g,ft);

x = rand(N,T);

c1 = gsp_jtv_frame_analysis(F,x);

param.method = 'exact';
c2 = gsp_jtv_filter_analysis(G,g,ft,x,param);

errors = errors + gsp_assert_test(c1,c2,eps(1000),'JTV - analysis operator matrix');

end

function errors = test_jtv_matrix_synthesis_op()
errors = 0;

N = 20;
T = 20;
G = gsp_sensor(N);
param.extension=1;
G = gsp_jtv_graph(G,T,[],param);
G = gsp_compute_fourier_basis(G);

[g,ft] = gsp_jtv_design_wave(G);

F = gsp_jtv_compute_frame(G,g,ft);

x = randn(N,T);

param.method = 'exact';

coeff = gsp_jtv_filter_analysis(G,g,ft,x,param);

x1 = gsp_jtv_frame_synthesis(F,coeff);

x2 = gsp_jtv_filter_synthesis(G,g,ft,coeff,param);

errors = errors + gsp_assert_test(x1,x2,eps(1000),'JTV - synthesis operator matrix');

end

function errors = test_jtv_filter_dual()
errors = 0;

N = 50;
T = 50;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

alpha=linspace(0,2,5);
[g,ft] = gsp_jtv_design_damped_wave(G,alpha,[0.01 0.1]);

try
    gd = gsp_jtv_design_can_dual(g,ft);
    gsp_plot_jtv_filter(G,gd,ft);
    close
    fprintf('Test JTV - Canonical 1: OK \n');
catch
    errors = errors +1;
    warning('Test JTV - Error canonical dual')
end

N = 50;
T = 50;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

Nf = 4;
g = gsp_jtv_design_meyer(G,Nf);

try
    gd = gsp_jtv_design_can_dual(g,ft);
    gsp_plot_jtv_filter(G,gd,ft);
    close
    fprintf('Test JTV - Canonical 2: OK \n');
catch
    errors = errors +1;
    warning('Test JTV - Error canonical dual')
end

end

function errors = test_jtv_filter_inverse()
%%
errors = 0;

N = 50;
T = 50;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);

alpha=linspace(0,2,5);
[g,ft] = gsp_jtv_design_damped_wave(G,alpha);
f = randn(N,T);

f1 = gsp_jtv_filter_analysis(G,g,ft,f);
f2 = gsp_jtv_filter_inverse(G,g,ft,f1);

errors = errors + gsp_assert_test(f,f2,eps(1000),'JTV - filter inverse 1');

[h,ft] = gsp_jtv_filter_array(G,g,ft);
f1 = gsp_jtv_filter_analysis(G,h,ft,f);
f2 = gsp_jtv_filter_inverse(G,h,ft,f1);
errors = errors + gsp_assert_test(f,f2,eps(1000),'JTV - filter inverse 2');

Nf = 4;
[g,ft] = gsp_jtv_design_meyer(G,Nf);
f1 = gsp_jtv_filter_analysis(G,g,ft,f);
f2 = gsp_jtv_filter_inverse(G,g,ft,f1);

errors = errors + gsp_assert_test(f,f2,eps(1000),'JTV - filter inverse 3');

[h,ft] = gsp_jtv_filter_array(G,g,ft);
f1 = gsp_jtv_filter_analysis(G,h,ft,f);
f2 = gsp_jtv_filter_inverse(G,h,ft,f1);
errors = errors + gsp_assert_test(f,f2,eps(1000),'JTV - filter inverse 4');

end

function errors = test_theory()

errors = 0;

N = 100;
T = 50;
G = gsp_sensor(N);
G = gsp_jtv_graph(G,T);
G = gsp_compute_fourier_basis(G);
g = @(x,w) abs(exp(-w.*x)).^2;
ft = 'js';
garray = gsp_jtv_filter_array(G,g,ft);
x1p = gsp_itft(G, gsp_vec2mat(gsp_filter_analysis(G,garray,gsp_delta(G,1) ),T))/sqrt(T);
x2p = gsp_jtv_filter_analysis(G,g,ft,gsp_jtv_delta(G,1,1));
% norm(x1p(:)/norm(x1p(:))-x2p(:)/norm(x2p(:)))

errors = errors +gsp_assert_test(x1p,x2p,1e-10,'TEST THEORY')

end

repo: gspbox


text:
%% Basic Data Structure representing a Mesh
%
%% Data structure: node and elem
%
% Two matrices |node(1:N,1:d)| and |elem(1:NT,1:d+1)| are used to represent
% a d-dimensional triangulation embedded in $\mathbf R^d$, where |N| is the number
% of vertices and |NT| is the number of elements.
%
% |node(k,1)| and |node(k,2)| are the x- and y-coordinates of the k-th node
% for points in 2-D. In 3-D, |node(k,3)| gives the additional z-coordinates
% of the k-th node.
%
% |elem(t,1:3)| are indices of 3 vertices of triangle t. |elem(t,1:4)| are
% indices of 4 vertices of tetrahedron t. By convention, the vertices are
% ordered such that the signed area/volume is positive. Therefore in 2-D,
% three vertices of a triangle is ordered counterclockwise and in 3-D, the
% four vertices of a tetrahedron follows the right-hand rule.

clear all; close all;
%% Example: L-shape domain in 2-D
node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0];  % coordinates
elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];     % connectivity
showmesh(node,elem)
axis on
findnode(node)       % plot indices of all vertices
findelem(node,elem)  % plot indices of all triangles
%%
% Apply uniform refinement three times to obtain a fine mesh.
for i = 1:3
    [node,elem] = uniformrefine(node,elem);
end
showmesh(node,elem)

%% Example: Cube in 3-D
node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
clf; showmesh3(node,elem,[],'FaceAlpha',0.25);
view([-53,8]);
axis on
findnode3(node)
findelem3(node,elem)
%%
% Apply uniform refinement twice to obtain a fine mesh.
for i = 1:2
    [node,elem] = uniformrefine3(node,elem);
end
showmesh3(node,elem,[],'FaceAlpha',0.25); view([-24 10]);

%% Example: Sphere in 3-D
% a simple triangular mesh in the space
node = [1,0,0; 0,1,0; -1,0,0; 0,-1,0; 0,0,1; 0,0,-1];
elem = [6,1,2; 6,2,3; 6,3,4; 6,4,1; 5,1,4; 5,3,4; 5,3,2; 5,2,1];
showmesh(node,elem,'FaceAlpha',0.5);
axis on;
findnode3(node);
%%
% uniformly refined mesh
for i = 1:3
    [node,elem] = uniformrefine(node,elem);
end
showmesh(node,elem)
%%
% project vertices onto the unit sphere
r = sqrt(node(:,1).^2 + node(:,2).^2 + node(:,3).^2);
node = node./[r r r];
showmesh(node,elem)

%% Order of vertices: orientation
% Any permutation of vertices of an element will represent the same
% abstract simplex. By convention, the vertices of a simplex is ordered
% such that the signed volume is positive. Therefore in 2-D, three vertices
% of a triangle is ordered counterclockwise and in 3-D, the ordering of
% vertices follows the right-hand rule. The function |fixorientation| will
% compute the signed area or volume and permute vertices if necessary.

node = [1,0; 1,1; 0,1];
elem = [1 3 2];
figure(4);
subplot(1,2,1)
showmesh(node,elem)
findnode(node,elem)
display('Clockwise'); display(elem)
elem = fixorientation(node,elem);
display('Counter-Clockwise'); display(elem)

%% Order of vertices: longest edge
% An even permutation of vertices is still allowed to represent the same
% simplex with the same orientation.
%
% * The function |label| will permute the vertices such that |elem(t,2:3)|
% is the longest edge of |t|.
% * The function |label3| will permute the vertices such that |elem(t,1:2)|
% is the longest edge of |t|.
%
% These functions are important for bisection methods of triangulations.

node = [1,0; 1,1; 0,1];
elem = [1 2 3];
subplot(1,2,1)
showmesh(node,elem)
findnode(node,elem)
display('Before labeling'); display(elem)
elem = label(node,elem);
display('After labeling'); display(elem)
display('elem(t,2:3) is the longest edge')

%%
node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
elem = [1,2,3,7; 1,6,2,7; 1,5,6,7];
figure(2)
showmesh3(node,elem,[],'FaceAlpha',0.25); view([-53,8]);
findnode3(node,[1 2 3 5 6 7]);
findelem3(node,elem);
display('Before labeling'); display(elem)
elem = label3(node,elem);
display('After labeling'); display(elem)
display('elem(t,1:2)=[1 7] is the longest edge')


text:
% Demo script for granulometry applied on grayscale image of coins.
%
%   output = granulo_coins(input)
%
%   Example
%   granulo_coins
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2023-01-17,    using Matlab 9.13.0.2049777 (R2022b)
% Copyright 2023 INRAE.

%% Sample image

% read image
img = imread('coins.png');
imwrite(img, 'coins.png');

% compute several steps (before, between and after the two peaks)
radiusList = [15 27 35];
for iPlot = 1:3
    se = strel('disk', radiusList(iPlot), 0);
    imgOp = imopen(img, se);

    imwrite(imgOp, sprintf('coins_OpDisk%02d.tif', radiusList(iPlot)));
end

%% compute granulo.

% granulometry analysis setup
xi = 1:60;
vol0 = sum(img(:));
[gr, diams, vol] = imGranulo(img, 'opening', 'disk', xi);

% display granulo
figure; plot(diams, gr, 'color', 'b', 'linewidth', 2);
xlim([0 100]);
xlabel('Diameter of structuring element (pixels)');
ylabel('Variation of gray levels (%)');
title('Gray level granulometry by opening', 'Interpreter', 'none');
print(gcf, 'coins_grOpDk30.png', '-dpng');

% display volume curve
figure; plot([0 diams], vol, 'color', 'b', 'linewidth', 2);
xlim([0 100]);
xlabel('Diameter of structuring element (pixels)');
ylabel('Sum of gray levels');
title('Gray level granulometry by opening', 'Interpreter', 'none');
print(gcf, 'coins_grOpDk40_sumOfGrays.png', '-dpng');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function Deff = electrolyteDiffusionCoefficients(ce,T,param,batterySection)
% electrolyteDiffusionCoefficients  Evaluates the diffusion coefficients for the electrolyte phase [m^2/s].
%
%   [Deff_p, Deff_s, Deff_n] = electrolyteDiffusionCoefficients(ce,T,param) evaluates
%   the diffusion coefficients for the anode, separator and cathode of the
%   battery.
%
%   Note that this is an interface for the main program. The authors
%   suggest to maintain the name of the function and its signature, while
%   modifying only the body of the script.
%
%   The diffusion coefficients can be evaluated in isothermal case
%   (param.TemperatureEnabled=0) or adiabatic case
%   (param.TemperatureEnabled=1 or 2).
%
%   You can modify the way that the diffusion coefficients are computed, as
%   function of electrolyte concentration and temperature. The main script
%   will pass also the param array.
%
%   The user can modify this script to meet particular requirements.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

if(param.TemperatureEnabled>=1)
    switch(batterySection)
        case 'p'
            Deff = param.eps_p^param.brugg_p*1e-4*10.^((-4.43-54./(T-229-5e-3*ce)-0.22e-3*ce));
        case 's'
            Deff = param.eps_s^param.brugg_s*1e-4*10.^((-4.43-54./(T-229-5e-3*ce)-0.22e-3*ce));
        case 'n'
            Deff = param.eps_n^param.brugg_n*1e-4*10.^((-4.43-54./(T-229-5e-3*ce)-0.22e-3*ce));
    end
else
    switch(batterySection)
        case 'p'
            Deff = repmat(param.Dp*param.eps_p^param.brugg_p,param.Np,1);
        case 's'
            Deff = repmat(param.Ds*param.eps_s^param.brugg_s,param.Ns,1);
        case 'n'
            Deff = repmat(param.Dn*param.eps_n^param.brugg_n,param.Nn,1);
    end
end

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function test_suite=SimTest_amico_test
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;

function TestSetup
setenv('ISDISPLAY','0') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','1')

function test_Sim
disp('===========================================')
disp('Running simulation test for amico');
disp('testing Simulation Single Voxel Curve...');

Model = str2func('amico'); Model = Model();
savedModel_fname = fullfile(fileparts(which('qMRLab')),'Test','MoxUnitCompatible','static_savedModelsforRetrocompatibility',['amico.qmrlab.mat']);
if ~exist(savedModel_fname,'file')
Model.saveObj(savedModel_fname);
else
Model = Model.loadObj(savedModel_fname);
end

disp(class(Model))
try Opt = button2opts(Model.Sim_Single_Voxel_Curve_buttons,1); end
try st = Model.st; catch, try st = mean([Model.lb(:),Model.ub(:)],2); catch, st = ones(length(Model.xnames),1); end; end
if exist('Opt','var') && length(Opt)>1
[Opt(:).SNR] = deal(1000);
else
Opt.SNR=1000;
end
for iopt=1:length(Opt) % Test all simulation options
disp(['Testing ' class(Model) ' simulation option:'])
disp(Opt(iopt))
FitResults = Model.Sim_Single_Voxel_Curve(st,Opt(iopt));
% Compare inputs and outputs
fnm=fieldnames(FitResults);
FitResults = rmfield(FitResults,fnm(~ismember(fnm,Model.xnames))); fnm=fieldnames(FitResults);
[~,FitResults,GroundTruth]=comp_struct(FitResults,mat2struct(st,Model.xnames),[],[],.30);
assertTrue(isempty(FitResults) & isempty(GroundTruth),evalc('FitResults, GroundTruth'))
end
disp ..ok

function TestTeardown
setenv('ISDISPLAY','') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
clc;
close all;
clearvars;
rand('seed', 0);
randn('seed', 0);

m = 1000; % number of examples
n = 100;  % number of features
% A random matrix
A = randn(m,n);
% Original signal
x0 = 10*randn(n,1);
% Measurements
b = A*x0;
% Introduce sparse noise in the signal
% Select a subset of indices [2% of m indices]
idx = randsample(m,ceil(m/50));
% Introduce noise in these indices
b(idx) = b(idx) + 1e2*randn(size(idx));
% Perform Least absolute deviation optimization
tstart = tic;
options.verbose  = 1;
x = spx.opt.admm.lad(A, b, options);
toc(tstart)
% measure the maximum difference
recovery_error = x - x0;
max_diff = max(abs(x - x0));
fprintf('Maximum difference: %e\n', max_diff);
fprintf('Recovery error norm: %e\n', norm(recovery_error));

author: indigits

repo: sparse-plex