text:
classdef(Abstract) ZmapFunction < matlab.mixin.Copyable
    %ZmapFunction base class for functions in Zmap, providing a common interface
    %   expects that functions can be called from command prompt or menus
    %
    %  To use, first create a class that inherits from ZmapFunction and provides
    %  the following:
    %
    %
    %     define variables here that would be used in calculation
    %
    %   ZmapFunction Properties(Constant):
    %      PlotTag - string name used to access all items from the plot this creates.
    %
    %
    %   ZmapFunction Methods:
    %
    %       The following methods are RECOMMENDED, but do not need to be defined in your derrived
    %       class:
    %
    %       ModifyGlobals - change ZmapGlobal.Data items here. If possible: ONLY here.
    %
    %       The following methods MUST BE DEFINED.  They provide a
    %       coherent framework for this function's behaviors and to make it easy to
    %       control when this routine does things (IO, for example)
    %
    %       function obj=<name of this class/function>(varargin)
    %           % this is your function constructor that should behave so:
    %           % if varargin is empty, then do the interactive setup.
    %           % otherwise, check pronditions, calculate values, plot, and then modify globals.
    %
    %       InteractiveSetup - create a ZmapDialog that provides the user a way to manipulate parameters used by the Calculate() method.
    %
    %       Calculate - do the calculations, and put all important results into fields of the obj.Results variable.
    %         optionally return the Results so that the caller can directly manipulate them
    %
    %       plot - plot the results of this function
    %
    %  METHODS(Static)
    %       AddMenuItem - create a menu item that calls the constructor
    %
    %
    % TROUBLESHOOTING HELP.
    %    If not all methods have been created, then MATLAB will fail to create
    %    your function with a message "Abstract classes cannot be instantiated."
    %    that includes a clickable link.  Click it go get a list of methodes/properties
    %    that still need to be created.
    %
    %
    % see also ZmapDialog, sample_ZmapFunction, blank_ZmapFunction, ZmapQuickResultPcolor

    properties
        % THESE ARE ACCESSIBLE BY ALL DERRIVED CLASSES

         % holds complete catalog to be analyzed. defaults to the primary catalog
        RawCatalog      {mustBeZmapCatalog}  = ZmapGlobal.Data.primeCatalog
        Result          struct            % results of the calculation, stored in a struct

        ZG                           = ZmapGlobal.Data; % provides access to the ZMAP global variables.
        hPlot                               % tracks the plot(s) for each function
        ax                           = [];       % axis where plotting will go

        AutoShowPlots   logical      = true; % will plots be generated upon calculation?
        InteractiveMode logical      = true;
        DelayProcessing logical      = false;
        CalcFinishedFcn function_handle = @do_nothing ; %call this once calculation ("do it") has finished.
    end

    properties(Constant,Abstract)
        PlotTag % string used for tracking the plot for each function (assign in your function)
        ParameterableProperties % properties from each class that can be used as constructor parameters
        References % paper references for calculations contained within this fn
    end

    properties(Transient)
        ResultDisplayer ResultsDisplay.ZmapResultsPlugin = ResultsDisplay.Hdisplay% chosen from +ResultsDisplay
    end

    properties(Dependent)
        FunctionCall    char; % text representation of the function call.
    end

    methods
        function obj=ZmapFunction(catalog)
            if exist('catalog','var')
                obj.RawCatalog=catalog;
                obj.verify_catalog(catalog);
            end
            obj.ResultDisplayer.Parent = obj;
        end

        function fcall = get.FunctionCall(obj)
            % FUNCTIONCALL provides probable function call for CURRENT STATE of object
            % obj.SETFUNCTIONCALL(varargin)

            fcall = class(obj) + "(ZAP";
            for j=1:numel(obj.ParameterableProperties)
                fld = obj.ParameterableProperties(j);
                val=obj.(fld);
                switch class(val)
                    case 'char'
                        fcall = fcall + sprintf(", '%s', '%s'", fld, val);
                    case 'string'
                        fcall = fcall + sprintf(', ''%s'', "%s"', fld, val);
                    case 'datetime'
                        fcall = fcall + sprintf(", '%s', datetime('%s')", fld, string(val));
                    case 'duration'
                        val.Format = 'dd:hh:mm:ss.SSS';
                        fcall = fcall + sprintf(", '%s', duration('%s')", fld,string(val));
                    case 'McMethods'
                        fcall = fcall + sprintf(", '%s', McMethods.%s", fld,string(val));

                    otherwise
                        if isnumeric(val) && isempty(val)
                            val = '[]';
                        elseif (isnumeric(val) || islogical(val)) && numel(val)>1
                            v="";
                            for r=1:size(val,1)
                                v(r)=strjoin(string(val(r,:)),',');
                            end
                            val="[" + strjoin(v,';') + "]";
                        else
                            try
                            val=string(val);
                            if ismissing(val)
                                val = "<missing>";
                            end
                            catch
                                val = "<some " + class(val) + ">";
                            end
                        end
                        fcall = fcall + sprintf(", '%s', %s", fld, val);
                end
            end
            fcall = fcall + ", 'AutoShowPlots', " + string(obj.AutoShowPlots);
            fcall = fcall + ", 'InteractiveMode', " + string(obj.InteractiveMode);
            % fcall = fcall + ", 'DelayProcessing', " + string(obj.DelayProcessing);
            fcall = fcall + ")";
        end

        function p = parseParameters(obj, varginstuff)

            % property list is a list of this objects properties that
            % can be set via the parameter list.
            p=inputParser();
            for i=1:numel(obj.ParameterableProperties)
                % use this object's default properties as the... uh. default.
                p.addParameter(obj.ParameterableProperties(i), obj.(obj.ParameterableProperties(i)));
            end

            % add parameters inheritable by all functions
            p.addParameter('AutoShowPlots',   obj.AutoShowPlots)
            p.addParameter('InteractiveMode', obj.InteractiveMode);
            p.addParameter('DelayProcessing', obj.DelayProcessing);
            p.addParameter('CalcFinishedFcn', obj.CalcFinishedFcn);
            p.parse(varginstuff{:});

            % assign values from parameter to this object
            for i=1:numel(obj.ParameterableProperties)
                obj.(obj.ParameterableProperties(i)) = p.Results.(obj.ParameterableProperties(i));
            end
            obj.AutoShowPlots = p.Results.AutoShowPlots;
            obj.InteractiveMode = p.Results.InteractiveMode;
            obj.DelayProcessing = p.Results.DelayProcessing;
            obj.CalcFinishedFcn = p.Results.CalcFinishedFcn;
        end

        function StartProcess(obj)
            if obj.DelayProcessing
                msg.dbdisp('Delaying processing',class(obj))
                return
            end

            if obj.InteractiveMode
                obj.InteractiveSetup();
            else
                obj.doIt();
            end
        end

        function clearPlot(obj)
            % CLEARPLOT remove all plot items added by this function
            delete(obj.hPlot);
            obj.hPlot=[];
        end

        function doIt(obj)
            % DOIT called by the interactive dialog box once OK is pressed.
            % Calls, in turn: Calculate, plot, ModifyGlobals, saveToDesktop
            %
            % each of these functions is defined in the subclass.
            % change DOIT behavior by redefining DOIT in the sublcass.
            obj.Calculate();
            if obj.AutoShowPlots
                obj.plot();
            end
            obj.ModifyGlobals();
            obj.saveToDesktop();
        end

        function  ModifyGlobals(obj)
            % MODIFYGLOBALS modify zmapglobals (define in subclass)
            % do any changing of zmap globals here
            do_nothing(obj)
        end

        function saveToDesktop(obj)
            % SAVETODESKTOP saves results to desktop
            vname = matlab.lang.makeValidName([class(obj),'_result']);
            %assert(~isfield(obj.Result,'FunctionCall'), 'FunctionCall is a reserved field in the results');
            obj.Result.FunctionCall = obj.FunctionCall;

            obj.Result.InCatalogName = obj.RawCatalog.Name; %was OperatingCatalog
            assignin('base', vname, obj.Result);
            fprintf('%s  %% called by Zmap : %s\n',obj.FunctionCall, char(datetime));
            fprintf('%% results in: %s\n',vname);
        end

        function plot(obj, varargin)
            obj.ResultDisplayer.plot(varargin{:})
        end
    end % PUBLIC METHODS

    methods(Abstract)  % these functions MUST be defined in every class derived from ZmapFunction

        % INTERACTIVESETUP ask the user for additional parameters necessary for the calculation
        InteractiveSetup(obj);

        % CALCULATE perform calculations
        % obj.CALCULATE() stores results in obj.Result  (obj.Result.values = results)
        % RESULTS = OBJ.CALCULATE() also, return results
        Results=Calculate(obj);

        % PLOT visually communicate results to the user
        %
        % obj.PLOT(varargin) plot results, passing along all parameters

    end %ABSTRACT METHODS

    methods % will be protected
        function t = helptext(obj, subject)
            % provides help text, based on the subject.
            % -all lists everything
            % -choice gives list of topics one can ask about
            switch subject
                case 'PlotTag'
                    t = "Title for this plot will be '" + obj.PlotTag+"'";
                case 'ParameterableProperties'
                    m=metaclass(obj);
                    t = "Name-Value pairs that can be specified when creating this function:"+newline+...
                        "";
                    plist = [string(obj.ParameterableProperties),"AutoShowPlots","InteractiveMode","DelayProcessing"];
                    for j=1:numel(plist)
                        pidx=string({m.PropertyList.Name})==plist(j);
                        p = m.PropertyList(pidx);
                        t=t+newline+"   "+ pad(plist(j),20,'right') + " : " + class(obj.(plist(j)));
                        if p.HasDefault
                            if numel(p.DefaultValue)>0 && (isnumeric(p.DefaultValue) || islogical(p.DefaultValue))
                                dv=string(mat2str(p.DefaultValue));
                            else
                                dv = string(p.DefaultValue);
                            end
                            if ismissing(dv)
                                if isnumeric(p.DefaultValue)
                                    dv="nan";
                                else
                                    dv="<missing>";
                                end
                            end

                        if dv.strlength > 33
                            dv=dv.extractBefore(15) + "..."+dv.extractAfter(dv.strlength-14);
                        end
                            t=t+" [default: "+dv+"]";
                        end
                        try
                        if ~isempty(p.Validation) && ~isempty(p.Validation.ValidatorFunctions)
                            fn = cellfun(@func2str,p.Validation.ValidatorFunctions,'UniformOutput',false);
                            t=t+ " that " + strjoin(fn,',');
                        end
                        catch ME
                            disp("error wile parsing validator functions:"+ ME.identifier)
                        end
                    end

                    if ~isempty(plist)
                        pv = obj.(obj.ParameterableProperties(1));
                        if ischar(pv)
                            pp = sprintf("'%s'",pv);
                        elseif isstring(pv)
                            pp = sprintf('"%s"',pv);
                        elseif isa(pv,'function_handle')
                            pp = sprintf('@%s', func2str(pv)');
                        elseif (isnumeric(pv) || islogical(pv)) && length(pv)>1
                            pp=string(mat2str(pv));
                        else
                            try
                                pp = string(pv);
                            catch
                                pp="some_"+ class(pv);
                            end
                        end

                        t = t+ newline+newline+ "For example:"+newline+...
                            "    obj = "+class(obj) +"(... , '" + ...
                            obj.ParameterableProperties(1) + "', " + ...
                            pp + ")";
                    end

                    assert(numel(t)==1)

                case "ReturnDetails"
                    t = "This function calculates :" + newline + ...
                        strjoin("  "+pad(obj.ReturnDetails.Names,20,'right') + " "+ obj.ReturnDetails.Descriptions,newline);

                case "References"
                    t = "References for this function:"+newline;
                    padValue = "  ";
                    r = strjoin(padValue + string(obj.References),newline);
                    if r~=padValue
                        t = t + r;
                    else
                        t = t + "   <missing>";
                    end

                case "-choices"
                    t = ["PlotTag","ParameterableProperties","ReturnDetails", "References"];
                    t = [t,  obj.ResultDisplayer.helptext(subject)];
                    t = t(~ismissing(t)) ;
                case "-all"
                    ch = obj.helptext("-choices");
                    t="";
                    for i=1:numel(ch)
                        t = t+ newline + newline + obj.helptext(ch(i));
                    end

                otherwise
                    placesToLook = "ResultDisplayer"; %comma separated string list
                    for i = placesToLook
                        t = obj.ResultDisplayer.helptext(subject);
                        if ~ismissing(t)
                            return
                        end
                    end
                    error('help for item %s not found', subject)
             end
        end
    end

    methods(Static, Abstract)

        % create a menu item that will be used to call this function/class (Define this as a static member of your class)
        %    sample implementation here:
        %
        %    function h = AddMenuItem(parent, zapFcn, varargin)
        %        % ADDMENUITEM create a menu item that will be used to call this function/class
        %        % h = AddMenuItem(parent, zapFcn) creates a menu item under the container PARENT
        %        %     that will call this class/fuction.  ZAPFCN is a function_handle that, when
        %        %     called returns a ZmapAnalysisPkg containing all the state details used to
        %        %     process the data.
        %
        %       label = 'testmenuitem';                        % CHANGE THIS TO YOUR MENUNAME
        %       cb = @(~,~)sample_ZmapFunction(zapFcn());      % CHANGE THIS TO YOUR CALLBACK
        %       h = uimenu(parent, 'Label', label, 'MenuSelectedFcn', cb, varargin{:});
        %    end
        %
        % you can copy/paste the above into your Static methods section
        AddMenuItem(parent, zapFcn, varargin)

    end % STATIC ABSTRACT METHODS

    methods(Static, Access=protected)
        function verify_catalog(c)
            % VERIFY_CATALOG makes sure the catalog is valid and has events
            % verify_catalog( CATALOG )
            assert(~isempty(c),'Catalog is empty');
            assert(isa(c,'ZmapCatalog'),'Please provide a Zmap catalog, not a ',class(c));
        end
    end
end

author: CelsoReyes


text:
classdef TestAreImagesInSameGroup < CoreTest
    % TestAreImagesInSameGroup. Tests for the DMAreImagesInSameGroup function.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods
        function obj = TestAreImagesInSameGroup
            this_metadata = [];
            other_metadata = [];

            % For images to be in the same group, the following tags must either
            % be present in both metadata, or not present in both metadata. If
            % present, they must be equal. We check that the test passes in
            % these cases and fails otherwise

            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PatientName', 'Name1', 'Name2');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PatientID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PatientBirthDate');

            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'StudyInstanceUID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SeriesInstanceUID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'StudyID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'StudyDescription');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SeriesNumber');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SeriesDescription', 'Series Description 1', 'Series Description 2');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'StudyDate');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SeriesDate');

            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'Rows', 256, 512);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'Columns', 256, 512);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PixelSpacing', [1 2], [3 4]);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PatientPosition');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'ImageOrientationPatient', [1;0;0;0;1;0], [0;1;0;0;0;1]);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'FrameOfReferenceUID');

            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'Modality', 'CT', 'MR');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'MediaStorageSOPClassUID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'ImageType', 'ORIGINAL\PRIMARY\AXIAL', 'ORIGINAL\PRIMARY');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SOPClassUID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'ImplementationClassUID');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'ImagesInAcquisition', 100, 102);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'SamplesPerPixel', 1, 2);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PhotometricInterpretation');
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'BitsAllocated', 1, 2);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'BitsStored', 1, 2);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'HighBit', 1, 2);
            [this_metadata, other_metadata] = obj.CheckMatchingOfTags(this_metadata, other_metadata, 'PixelRepresentation');

            % Slice location tags: existence must match but values do not have
            % to match
            this_metadata.SliceLocation = 1;
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should fail');
            other_metadata.SliceLocation = 2;
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should pass');

            % ImagePositionPatient tag: check consistency of image locations
            % Image positions should lie approximately on a line; there is a
            % tolerance of 10mm
            extra_metadata = other_metadata;
            this_metadata.ImagePositionPatient = [10;20;30];
            extra_metadata.ImagePositionPatient = [20;40;60];
            other_metadata.ImagePositionPatient = [30;60;110];
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should fail');
            other_metadata.ImagePositionPatient = [30;60;90];
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should pass');
            extra_metadata.ImagePositionPatient = [20;-40;60];
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should fail');

            this_metadata.ImagePositionPatient = [0;30;40];
            extra_metadata.ImagePositionPatient = [0;90;120];
            other_metadata.ImagePositionPatient = [0;60;80];
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should pass');
            other_metadata.ImagePositionPatient = [0;55;80];
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should pass');
            this_metadata.ImagePositionPatient = [0;30;-40];
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should fail');
            this_metadata.ImagePositionPatient = [0;15;40];
            extra_metadata.ImagePositionPatient = [0;90;120];
            other_metadata.ImagePositionPatient = [0;60;80];
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata, extra_metadata), 'Should fail');
        end
    end

    methods (Access = private)
        function [this_metadata, other_metadata] = CheckMatchingOfTags(obj, this_metadata, other_metadata, tag_name, value_1, value_2)
            if nargin < 6
                value_2 = CoreSystemUtilities.GenerateUid;
            end
            if nargin < 5
                value_1 = CoreSystemUtilities.GenerateUid;
            end

            % Images should match if the tag is not present in either
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should pass');

            % Images should not match if one contains a tag
            this_metadata.(tag_name) = value_1;
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should fail');

            % Images should not match if the tags don't match
            other_metadata.(tag_name) = value_2;
            obj.Assert(~DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should fail');

            % Images should match if the tags match
            other_metadata.(tag_name) = value_1;
            obj.Assert(DMAreImagesInSameGroup(this_metadata, other_metadata), 'Should pass');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [k, sk, n2] = rbfperiodicKernCompute(kern, x, x2)

% RBFPERIODICKERNCOMPUTE Compute the RBFPERIODIC kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the RBF derived periodic
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the RBF derived periodic
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : rbfperiodicKernParamInit, kernCompute, kernCreate, rbfperiodicKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2007, 2009

% KERN
if isfield(kern, 'period')
  factor = 2*pi/kern.period;
else
  factor = 1;
end
if nargin < 3
  n2 = sin(0.5*factor*(repmat(x, 1, size(x, 1)) - repmat(x', size(x, 1), 1)));
  n2 = n2.*n2;
  wi2 = (2 .* kern.inverseWidth);
  sk = exp(-n2*wi2);
else
  n2 = sin(0.5*factor*(repmat(x, 1, size(x2, 1)) - repmat(x2', size(x, 1), 1)));
  n2 = n2.*n2;
  wi2 = (2 .* kern.inverseWidth);
  sk = exp(-n2*wi2);
end
k = kern.variance*sk;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function ExecFlag=DoSARWatchExec(handles)

global VObj;
global VMag;
global VCoi;
global VCtl;
global VSig;
global VSeq;
global VVar;

%preserve VObj
VTmpObj=VObj;

handles.Simuh=guidata(handles.Simuh.SimuPanel_figure);
%% Do spin execution
try
    % Read tab parameters
    fieldname=fieldnames(handles.Attrh1);
    for i=1:length(fieldname)/2
        try
            eval(['SW.' fieldname{i*2} '=[' get(handles.Attrh1.(fieldname{i*2}),'String') '];']);
        catch me
            TAttr=get(handles.Attrh1.(fieldname{i*2}),'String');
            eval(['SW.' fieldname{i*2} '=''' TAttr{get(handles.Attrh1.(fieldname{i*2}),'Value')}  ''';']);
        end
    end

    if SW.N_Second<VCtl.TR
        error('''N_Second'' needs to be no less than one TR.');
    end

    % Prescan config
    DoPreScan(handles.Simuh);
    DoUpdateBar(handles.TimeBar_axes,10,40);

    if ~isequal(size(VCoi.TxE1x(:,:,:,1)),[max(max(sum(VMag.FRange,1))),max(max(sum(VMag.FRange,2))),max(max(sum(VMag.FRange,3)))])
        error('E1+ map is not available or invalid.');
    end

    % Create Executing Virtual Structure VOsw
    VOsw=VObj;
    VOsw.MassDen(VMag.FRange==0)=[];
    VOsw.ECon(repmat(VMag.FRange,[1,1,1,3])==0)=[];

    % Kernel uses Mz to determine SpinMx size
    VOsw.MassDen=reshape(VOsw.MassDen,[max(max(sum(VMag.FRange,1))),max(max(sum(VMag.FRange,2))),max(max(sum(VMag.FRange,3)))]);
    VOsw.ECon=reshape(VOsw.ECon,[max(max(sum(VMag.FRange,1))),max(max(sum(VMag.FRange,2))),max(max(sum(VMag.FRange,3))),3]);

    % Spin execution
    VObj=VOsw;

    if sum(VObj.MassDen(VObj.MassDen~=0))*(VObj.XDimRes*VObj.YDimRes*VObj.ZDimRes)*0.01<SW.N_Gram/1000
        error(['''N_Gram'' needs to be less than ' num2str(sum(VObj.MassDen(VObj.MassDen~=0))*(VObj.XDimRes*VObj.YDimRes*VObj.ZDimRes)*10) 'g (1% of the total mass).']);
    end

    % Generate Pulse line
    DoPulseGen(handles.Simuh);
    DoUpdateBar(handles.TimeBar_axes,20,40);

    % Initialize SAR, Muts in VSig & VVar
    seq=VSeq.utsLine;
    dseq=diff(seq);
    dseq(dseq<0)=0;
    tSeqEnd=sum(dseq); % sequence end time point
    if isinf(SW.N_Second)
        tStart=0;
        tEnd=tSeqEnd;
        tSample=unique([tStart, tEnd]);
    else
        tEnd=max(0,SW.tStart):SW.dt:min(tSeqEnd,SW.tEnd); % prepare SAR sampling time point
        tStart=tEnd-SW.N_Second;
        tStart(tStart<0)=0;
        tStart(tEnd==0)=[];
        tEnd(tEnd==0)=[];
        if isempty(tEnd) | isempty(tStart)
            error('Input sample time is invalid.');
        end
        tSample=unique([tStart, tEnd]);
    end

    [row,col,layer]=size(VObj.MassDen);
    VSig.tSample=tSample;
    VSig.tRealSample=zeros(size(tSample));
    VSig.SAR=zeros([[row,col,layer] length(tSample)]);
    VSig.Muts=zeros(max(size(VSeq.utsLine)),1);
    VVar.SAR=zeros([[row,col,layer] 3]);

    % Simulation Process
    fprintf('Calculating unaveraged local SAR...\n');
    try
        % convert data type
        DoDataTypeConv(handles.Simuh);
        VObj.MassDen=single(VObj.MassDen);
        VObj.ECon=single(VObj.ECon);
        VCoi.TxE1x=single(VCoi.TxE1x);
        VCoi.TxE1y=single(VCoi.TxE1y);
        VCoi.TxE1z=single(VCoi.TxE1z);
        VSig.SAR=single(VSig.SAR);
        VVar.SAR=single(VVar.SAR);

        % calculate dissipated power
        VCtl.MaxThreadNum=int32(handles.Simuh.CPUInfo.NumThreads);
        VVar.SARi=int32(0);
        DoCalSARAtCPU;

        % calculate unaveraged SAR
        SAR=zeros([size(VObj.MassDen) length(tEnd)]);
        Power=zeros([size(VObj.MassDen) length(tEnd)]);
        tSARSample=zeros(size(tEnd));
        tSARSecond=zeros(size(tEnd));
        for i=1:length(tEnd)
            if tEnd(i)>max(VSig.tRealSample)
                SAR(:,:,:,i:end)=[];
                tSARSample(i:end)=[];
                break;
            end
            idx1 = find(tSample==tStart(i));
            idx2 = find(tSample==tEnd(i));
            tSARSecond(i) = VSig.tRealSample(idx2)-VSig.tRealSample(idx1);
            % Collins et.al. MRM 2001 and Tang et.al. PIERS 2007
            SAR(:,:,:,i)=((VSig.SAR(:,:,:,idx2) - VSig.SAR(:,:,:,idx1))/(VSig.tRealSample(idx2)-VSig.tRealSample(idx1)))./(2*VObj.MassDen);
            Power(:,:,:,i)=((VSig.SAR(:,:,:,idx2) - VSig.SAR(:,:,:,idx1))/(VSig.tRealSample(idx2)-VSig.tRealSample(idx1)))*(VObj.XDimRes*VObj.YDimRes*VObj.ZDimRes/2);
            tSARSample(i)=VSig.tRealSample(idx2);
        end
        SAR(isnan(SAR))=0;
        fprintf('Calculating unaveraged local SAR completed.\n');
        DoUpdateBar(handles.TimeBar_axes,30,40);

    catch me
        error_msg{1,1}='ERROR!!! SAR calculation process aborted.';
        error_msg{2,1}=me.message;
        errordlg(error_msg);
        ExecFlag=0;

        %recover VObj
        VObj=VTmpObj;
        return;
    end

    %N-gram average SAR
    if SW.N_Gram~=0
        fprintf('Calculating %f-gram local SAR...\n', SW.N_Gram);
%         aveSAR=zeros(size(SAR));
%         avePower=zeros(size(Power));
%         aveKGram=zeros(size(SAR));
        % Carluccio et.al. TBME 2013
        [aveSAR,avePower]=DoSARAverageAtCPU(SAR,Power,double(VObj.MassDen),VObj.XDimRes,VObj.YDimRes,VObj.ZDimRes,SW.N_Gram);
%         for i=1:size(SAR,4)
%             [aveSAR(:,:,:,i),avePower(:,:,:,i),aveKGram(:,:,:,i)]=DoSARAverage(SAR(:,:,:,i),Power(:,:,:,i),VObj.MassDen*(VObj.XDimRes*VObj.YDimRes*VObj.ZDimRes),SW.N_Gram/1000);
%             fprintf('(%d/%d) %f-gram local SAR at time point %fs completed.\n', i, size(SAR,4), SW.N_Gram, tSARSample(i));
%             DoUpdateBar(handles.TimeBar_axes,i,size(SAR,4));
%             pause(0.001);
%         end
    else
        aveSAR=SAR;
        avePower=Power;
%         aveKGram=zeros(size(SAR));
    end

    handles.SW = SW;
    handles.aveSAR = aveSAR;
    handles.avePower = avePower;
%     handles.aveKGram = aveKGram;
    handles.tSAR = tEnd; % prescribed sample point
    handles.tSARSample = tSARSample; % actual sample point
    handles.tSARSecond = tSARSecond; % actual N_Second
    handles.MassDen = VObj.MassDen;

    DoUpdateBar(handles.TimeBar_axes,40,40);
catch me
    error_msg{1,1}='ERROR!!! SAR calculation process aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    ExecFlag=0;

    %recover VObj
    VObj=VTmpObj;
    return;
end

%recover VObj
VObj=VTmpObj;

guidata(handles.SARWatcherPanel_figure, handles);
ExecFlag=1;

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function value = p14_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P14_F evaluates the integrand for problem 14.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    sum ( 1 <= i <= dim_num ) (-1)**i * product ( 1 <= j <= i ) x(j)
%
%  Exact Integral:
%
%    -1/3 ( 1 - (-1/2)**DIM_NUM )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Harald Niederreiter,
%    Implementation and Tests of Low-Discrepancy Sequences,
%    ACM Transactions on Modeling and Computer Simulation,
%    Volume 2, Number 3, July 1992, pages 195-213.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 0.0;

  for point = 1 : point_num

    factor = 1.0;

    for dim = 1 : dim_num

      factor = - factor * x(dim,point);

      value(point) = value(point) + factor;

    end

  end

  p14_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [matrix_data] = Dvec2Dmat(vector_data,WIDTH,LENGTH,IND,nodata)
% function which converts a vector to a matrix:
% option 1: WIDTH only where it is assumed that the original grid was
%                  made using vector = reshape(matrix,[],1);
% option 2: WIDTH, LENGTH, IND where the IND has the same length as the vector_data
%                  and where IND relates to the positon in a matrix. Note
%                  that this method allows for a mask file to be included.
%                  You can use the matlab function IND = sub2ind([WIDTH LENGTH], I, J)
%                  to generate the IND variable
% nodata argument is optional and us used to fill the nodata. By default
% this is a NaN.
%
%     Copyright (C) 2016  Bekaert David
%     davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% Bekaert David
% modifications:
% 8/10/2016     DB  Add nodata option

%  checking of the input arguments
method =1;          % this is the width method
if nargin <2
    error(['Expecting two inputs: vector_data and WIDTH'])
end
if nargin>2 && nargin<4
    fprintf('WARNING: you need to specify 4 inputs when planning on using IND option, will fall back on width only \n')
end
if length(WIDTH)~=1
    error('WIDTH is given incorrect')
end
if nargin >=4
   if isempty(LENGTH) || isempty(IND)
       error('You need to specify LENGTH and IND')
   end
   if length(LENGTH)~=1
       error('LENGTH is given incorrect')
   end
   if length(IND)~=length(vector_data)
       error('Length of vector_data and IND needs to be consistent')
   end
   method = 2;      % change to the width, length and ind method
end
if nargin<5
    nodata=NaN;
end

%% storing the data again
if method ==1
    matrix_data = reshape(vector_data,WIDTH,[]);
elseif method ==2
    if isnan(nodata)
        matrix_data = NaN([WIDTH LENGTH]);
    elseif nodata==0
        matrix_data = zeros([WIDTH LENGTH]);
    else
        matrix_data = ones([WIDTH LENGTH]).*nodata;
    end
    matrix_data(IND)= vector_data;
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
% MANIPULATOR TRAJECTORY GENERATION
% Generates Joint space trajectories by performing inverse kinematics on
% each waypoint and interpolating between the joint angles.
%
% Copyright 2019 The MathWorks, Inc.

%% Setup
clear, clc, close all

% Define waypoint information
createWaypointData;

% Define IK
ik = inverseKinematics('RigidBodyTree',gen3);
ikWeights = [1 1 1 1 1 1];
ikInitGuess = jointAnglesHome';
ikInitGuess(ikInitGuess > pi) = ikInitGuess(ikInitGuess > pi) - 2*pi;
ikInitGuess(ikInitGuess < -pi) = ikInitGuess(ikInitGuess < -pi) + 2*pi;

% Set up plot
plotMode = 1; % 0 = None, 1 = Trajectory, 2 = Coordinate Frames
show(gen3,gen3.homeConfiguration,'Frames','off','PreservePlot',false);
xlim([-1 1]), ylim([-1 1]), zlim([0 1.2])
hold on
if plotMode == 1
    hTraj = plot3(waypoints(1,1),waypoints(2,1),waypoints(3,1),'b.-');
end
plot3(waypoints(1,:),waypoints(2,:),waypoints(3,:),'ro','LineWidth',2);

%% Solve IK for all waypoints
includeOrientation = false; % Set this to use zero vs. nonzero orientations

numWaypoints = size(waypoints,2);
numJoints = numel(gen3.homeConfiguration);
jointWaypoints = zeros(numJoints,numWaypoints);

for idx = 1:numWaypoints
    if includeOrientation
        tgtPose = trvec2tform(waypoints(:,idx)') * eul2tform(orientations(:,idx)');
    else
        tgtPose =  trvec2tform(waypoints(:,idx)');
    end
    [config,info] = ik(eeName,tgtPose,ikWeights,ikInitGuess);
    jointWaypoints(:,idx) = config';
end

%% Generate trajectory on joint space
trajType = 'trap';
switch trajType
    case 'trap'
        [q,qd,qdd] = trapveltraj(jointWaypoints,numel(trajTimes), ...
            'AccelTime',repmat(waypointAccelTimes,[numJoints 1]), ...
            'EndTime',repmat(diff(waypointTimes),[numJoints 1]));

    case 'cubic'
        [q,qd,qdd] = cubicpolytraj(jointWaypoints,waypointTimes,trajTimes, ...
            'VelocityBoundaryCondition',zeros(numJoints,numWaypoints));

    case 'quintic'
        [q,qd,qdd] = quinticpolytraj(jointWaypoints,waypointTimes,trajTimes, ...
            'VelocityBoundaryCondition',zeros(numJoints,numWaypoints), ...
            'AccelerationBoundaryCondition',zeros(numJoints,numWaypoints));

    case 'bspline'
        ctrlpoints = jointWaypoints; % Can adapt this as needed
        [q,qd,qdd] = bsplinepolytraj(ctrlpoints,waypointTimes([1 end]),trajTimes);

    otherwise
        error('Invalid trajectory type! Use ''trap'', ''cubic'', ''quintic'', or ''bspline''');
end

% To visualize the trajectory, run the following line
% plotTrajectory(trajTimes,q,qd,qdd,'Names',"Joint " + string(1:numJoints),'WaypointTimes',waypointTimes)

%% Trajectory following loop
for idx = 1:numel(trajTimes)

    config = q(:,idx)';

    % Find Cartesian points for visualization
    eeTform = getTransform(gen3,config,eeName);
    if plotMode == 1
        eePos = tform2trvec(eeTform);
        set(hTraj,'xdata',[hTraj.XData eePos(1)], ...
                  'ydata',[hTraj.YData eePos(2)], ...
                  'zdata',[hTraj.ZData eePos(3)]);
    elseif plotMode == 2
        plotTransforms(tform2trvec(eeTform),tform2quat(eeTform),'FrameSize',0.05);
    end

    % Show the robot
    show(gen3,config,'Frames','off','PreservePlot',false);
    title(['Trajectory at t = ' num2str(trajTimes(idx))])
    drawnow

end

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function [y, z, p] = mixLinPred(model, X, t)
% Prediction function for mxiture of linear regression
% input:
%   model: trained model structure
%   X: d x n data matrix
%   t:(optional) 1 x n responding vector
% output:
%   y: 1 x n prediction
%   z: 1 x n cluster label
%   p: 1 x n predict probability for t
% Written by Mo Chen (sth4nth@gmail.com).
W = model.W;
alpha = model.alpha;
beta = model.beta;

X = [X;ones(1,size(X,2))]; % adding the bias term
y = W'*X;
D = bsxfun(@minus,y,t).^2;
logRho = (-0.5)*beta*D;
logRho = bsxfun(@plus,logRho,log(alpha));
T = logsumexp(logRho,1);
p = exp(T);
logR = bsxfun(@minus,logRho,T);
R = exp(logR);
z = max(R,[],1);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function result = rectangle_sub_2d ( func, xval, yval, nsub, norder, xtab, ...
  ytab, weight )

%*****************************************************************************80
%
%% RECTANGLE_SUB_2D carries out a composite quadrature over a rectangle in 2D.
%
%  Integration region:
%
%    XVAL(1) <= X <= XVAL(2),
%    YVAL(1) <= Y <= YVAL(2).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, external FUNC, the name of the function to be
%    integrated.  The user must declare the name an EXTERNAL
%    parameter in the calling program, pass the name of the
%    function in FUNC, and write a function of the form
%      function value = func ( x, y )
%    which evaluates the function at the point (X,Y).
%
%    Input, real XVAL(2), the left and right X coordinates.
%
%    Input, real YVAL(2), the lower and upper Y coordinates.
%
%    Input, integer NSUB(2).
%    NSUB(1) is the number of subintervals to use in the X direction,
%    and NSUB(2) is the same thing for Y.
%
%    Input, integer NORDER, the order of the rule.
%
%    Input, real XTAB(NORDER), YTAB(NORDER), the abscissas.
%
%    Input, real WEIGHT(NORDER), the weights of the rule.
%
%    Output, real RESULT, the approximate integral of the function.
%
  a(1) = xval(1);
  a(2) = yval(1);
  b(1) = xval(2);
  b(2) = yval(2);

  for i = 1 : 2
    if ( a(i) == b(i) )
      result = 0.0E+00;
      return
    end
  end

  for i = 1 : 2
    if ( nsub(i) < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'RECTANGLE_SUB_2D - Fatal error!\n' );
      fprintf ( 1, '  Nonpositive value of NSUB(I) = %d\n', nsub(i) );
      fprintf ( 1, '  for index I = %d\n', i );
      error ( 'RECTANGLE_SUB_2D - Fatal error!' );
    end
  end
%
%  Break up the X interval into NSUB(1) subintervals.
%
  volume = 0.0E+00;
  result = 0.0E+00;

  for i = 1 : nsub(1)

    xlo = r8vec_even_select ( a(1), b(1), nsub(1)+1, i );
    xhi = r8vec_even_select ( a(1), b(1), nsub(1)+1, i+1 );
%
%  Break up the Y interval into NSUB(2) subintervals.
%
    for j = 1 : nsub(2)

      ylo = r8vec_even_select ( a(2), b(2), nsub(2)+1, j );
      yhi = r8vec_even_select ( a(2), b(2), nsub(2)+1, j+1 );

      quad_sub = 0.0E+00;
      for k = 1 : norder

        x = xlo + 0.5E+00 * ( xtab(k) + 1.0E+00 ) * ( xhi - xlo );
        y = ylo + 0.5E+00 * ( ytab(k) + 1.0E+00 ) * ( yhi - ylo );

        quad_sub = quad_sub + weight(k) * feval ( func, x, y ) / 4.0E+00;

      end

      volume_sub = ( xhi - xlo ) * ( yhi - ylo );
      result_sub = quad_sub * volume_sub;

      volume = volume + volume_sub;
      result = result + result_sub;

    end

  end

  return
end


text:
function [sCodebook, settings] = gmlvq(sCodebook, sData, varargin)

%GMLVQ trains a codebook using the Generalized Matrix LVQ (GMLVQ)
%algorithm.
%
% sM = gmlvq(sM, sD, [argID, value, ...])
%
%  sM = gmlvq(sM, sD)
%  sM = gmlvq(sM, sD, 'learningRatePrototypes', [0.01 0.001], 'relevanceStart', 10);
%
% Required Input Arguments:
%   sM                     (struct) map struct containing the labels as the
%                                   first column of .labels
%   sD                     (struct) data struct containing the labels as
%                                   the first column of .labels
%
% Optional Input Arguments:
%   PrototypesPerClass     (vector) number of prototypes per class used;
%                                   a number or a vector with a number
%                                   for each class (default=1)
%   initialMatrix          (matrix) initial Matrix omega
%   regularization         (scalar) regularization for relevances
%   dim                    (scalar) maximum rank or projection dimension
%   comparable             (scalar) a flag which resets the random
%                                   generator to produce comparable results
%                                   if set to 1
%   optimization           (string) choice for the optimization technique:
%                                   sgd or fminlbfgs (default=fminlbfgs)
%  Parameter for the stochastic gradient descent sgd:
%   nb_epochs              (scalar) the number of epochs (default=100)
%   learningRatePrototypes (vector) learning rate for the prototypes; could
%                                   be the start and end value or a vector
%                                   of length nb_epochs
%   learningRateMatrix     (vector) learning rate for the relevance matrix;
%                                   could be the start and end value or a
%                                   vector of length nb_epochs
%   MatrixStart            (scalar) epoch to start the matrix learning
%                                   (default=1)
%  Parameter for the build-in function fminlbfgs:
%   threshstop             (scalar) the training error for early stopping
%                                   (default=0)
%   nb_reiterations        (scalar) the number of optimization reiterations
%                                   performed (default=100)
%   useEarlyStopping       (scalar) use early stopping based on threshstop
%                                   (default=1)
%   Display                (string) the optimization output 'iter' or 'off'
%                                   (default='off')
%   GradObj                (string) turn the usage of gradient information
%                                   on or off (default='on')
%   HessUpdate             (string) the update can be 'lbfgs', 'bfgs' or
%                                   'steepdesc' (default='lbfgs')
%   TolFun                 (scalar) the tolerance (default=1e-6)
%   MaxIter                (scalar) the maximal number of iterations
%                                   (default=2500)
%   MaxFunEvals                     the maximal number of function
%                                   evaluations (default=1000000)
%   TolX                   (scalar) tolerance on the minimum
%                                   (default=1e-10)
%   DiffMinChange          (scalar) minimal change (default=1e-10)
%
% Output Arguments:
%   sM                     (struct) map struct containing the optimized
%                                   prototypes and the decomposed
%                                   relevances omega
%   settings               (struct) information on the settings of the
%                                   algorithm
%
% NOTE: does not take the vector mask into account but trains a relevance
%       matrix.
%
% For more help, try 'type gmlvq', or check out the online documentation.
% See also GRLVQ, GMLVQ_CORE, LVQ3, LVQ1.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% gmlvq
%
% PURPOSE
%
% Trains a codebook and a global metric with the Generalized Matrix LVQ
% (GMLVQ) algorithm.
%
% SYNTAX
%
%   sM = gmlvq(sM, sD)
%   sM = gmlvq(sM, sD, 'learningRatePrototypes', [0.01 0.001], 'MatrixStart', 10);
%   sM = gmlvq(sM, sD, argID, value, ...);
%
% DESCRIPTION
%
% GMLVQ constitutes an enhancement of the LVQ algorithm by minimizing an
% explicit error function and optimizing the metric of the space where the
% computation is performed. Distances are calculated by
%   d(x,y) = (x - y) * lambda * (x - y)', where lambda = omega' * omega.
% The matrix omega is optimized by the algorithm.
%
% Two optimization techniques are implemented:
%     - stochastic gradient descent (sgd)
%     - limited memory Quasi Newton Broyden-Fletcher-Goldfarb-Shanno
%       (L-BFGS)
%
% Classification is performed similarly to standard LVQ, with the only
% difference that the above specified formula is applied for the
% computation of distances.
%
% This function provides an interface in the style of the som-toolbox and
% is basically a wrapper for the method gmlvq_core. It deals with structs
% rather then directly with data matrices.
%
% REFERENCES
%
%   Petra Schneider, Michael Biehl, Barbara Hammer: Adaptive Relevance
%   Matrices in Learning Vector Quantization. Neural Computation 21(12):
%   3532-3561 (2009).
%
%   K. Bunte, P. Schneider, B. Hammer, F.-M. Schleif, T. Villmann and M.
%   Biehl: Limited Rank Matrix Learning - Discriminative Dimension
%   Reduction and Visualization, Neural Networks, vol. 26, nb. 4,
%   pp. 159-173, 2012.
%
%   P. Schneider, K. Bunte, B. Hammer and M. Biehl: Regularization in
%   Matrix Relevance Learning, IEEE Transactions on Neural Networks, vol.
%   21, nb. 5, pp. 831-840, 2010.
%
% REQUIRED INPUT ARGUMENTS
%   sM                     (struct) map struct containing the initialized
%                                   prototype vectors as rows of the
%                                   .codebooks matrix and the labels as the
%                                   first column of .labels
%   sD                     (struct) data struct containing the data vectors
%                                   as rows of .data and the labels as the
%                                   first column of .labels
%
% OPTIONAL INPUT ARGUMENTS
%   PrototypesPerClass     (vector) number of prototypes per class used;
%                                   a number or a vector with a number
%                                   for each class (default=1)
%   initialMatrix          (matrix) initial Matrix omega
%                                   (default: identity matrix)
%   regularization         (scalar) regularization for relevances
%                                   (detault=0)
%   dim                    (scalar) maximum rank or projection dimension,
%                                   i.e. of the matrix omega
%                                   (default=nb of features for training)
%   comparable             (scalar) a flag which resets the random
%                                   generator to produce comparable results
%                                   if set to 1 (default=0)
%   optimization           (string) choice for the optimization technique:
%                                   sgd (stochastic gradient descent) or
%                                   fminlbfgs (Limited memory Quasi Newton
%                                   Broyden-Fletcher-Goldfarb-Shanno)
%                                   (default=fminlbfgs)
%  Parameter for the stochastic gradient descent sgd:
%   nb_epochs              (scalar) the number of epochs (default=100)
%   learningRatePrototypes (vector) learning rate for the prototypes; could
%                                   be the value for the first and the last
%                                   epoch or a vector of length nb_epochs
%   learningRateMatrix     (vector) learning rate for the relevance matrix;
%                                   could be the value for the first and
%                                   the last epoch or a vector of length
%                                   nb_epochs
%   MatrixStart            (scalar) epoch to start the matrix learning
%                                   (default=1)
%  Parameter for the build-in function fminlbfgs:
%   threshstop             (scalar) the training error for early stopping
%                                   (default=0)
%   nb_reiterations        (scalar) the number of optimization reiterations
%                                   performed (default=100)
%   useEarlyStopping       (scalar) use early stopping based on threshstop
%                                   (default=1)
%   Display                (string) the optimization output 'iter' or 'off'
%                                   (default='off')
%   GradObj                (string) turn the usage of gradient information
%                                   on or off (default='on')
%   HessUpdate             (string) the update can be 'lbfgs', 'bfgs' or
%                                   'steepdesc' (default=lbfgs)
%   TolFun                 (scalar) the tolerance (default=1e-6)
%   MaxIter                (scalar) the maximal number of iterations
%                                   (default=2500)
%   MaxFunEvals                     the maximal number of function
%                                   evaluations (default=1000000)
%   TolX                   (scalar) tolerance on the minimum
%                                   (default=1e-10)
%   DiffMinChange          (scalar) minimal change (default=1e-10)
%
% OUTPUT ARGUMENTS
%   sM                     (struct) map struct containing the optimized
%                                   prototypes and the decomposed
%                                   relevances omega
%   settings               (struct) information on the settings of the
%                                   algorithm
%
% EXAMPLES
%
%   lab   = unique(sD.labels(:,1));           % different classes
%   nProt = length(lab)*5;                    % 5 prototypes for each
%   sM = som_randinit(sD,'msize',[nProt 1]);  % initial prototypes
%   sM.labels = [lab;lab;lab;lab;lab];        % and their classes
%   sM = gmlvq(sM,sD);                        % use GMLVQ to adjust the
%                                             % prototypes and the metric
%
% SEE ALSO
%
%   grlvq         Use the GRLVQ algorithm for training.
%   gmlvq_core    Access the GMLVQ functionality without using structs.
%   lvq3          Use LVQ3 algorithm for training.
%   lvq1          Use LVQ1 algorithm for training.

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% Copyright (c) Alexander Schulz
% http://www.cis.hut.fi/projects/somtoolbox/

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

cods      = sCodebook.codebook;
codLabels = sCodebook.labels;
data      = sData.data;
labels    = sData.labels;

% convert labels from cell to a vector
labelsV    = som_label2num(labels);
codLabelsV = som_label2num(codLabels);

% calculate the number of prototypes per class
n_prot_per_class = zeros(1, max(codLabelsV));
for i=1:max(codLabelsV)
    n_prot_per_class(i) = sum(codLabelsV == i);
end

% call the main function
[model, settings]  = gmlvq_core(data, labelsV, 'initialPrototypes', ...
    [cods, codLabelsV], 'PrototypesPerClass', n_prot_per_class, varargin{:});

% write the results in the output struct
sCodebook.codebook = model.w;
sCodebook.omega    = model.omega;
%sCodebook.lambda   = sCodebook.omega'*sCodebook.omega;

if strcmp(settings.optimization, 'sgd')
    trainlen   = settings.nb_epochs;
    alpha_ini  = settings.learningRatePrototypes(1);
    alpha_type = 'power';
else
    trainlen   = NaN;
    alpha_ini  = NaN;
    alpha_type = '';
end

sTrain = som_set('som_train','algorithm','GMLVQ',...
		 'data_name',sData.name,...
		 'neigh','',...
		 'mask',ones(size(model.w,2),1),...
		 'radius_ini',NaN,...
		 'radius_fin',NaN,...
		 'alpha_ini',alpha_ini,...
		 'alpha_type',alpha_type,...
		 'trainlen',trainlen,...
		 'time',datestr(now,0));
sCodebook.trainhist(end+1) = sTrain;

author: ilarinieminen

repo: SOM-Toolbox


text:
function modelEFM = efmSubmodelExtractionAsSBML(model, rxnsForSubmodel, outputFileName, remFlag, ubiquitousMets)
% This function takes an input array containing reaction indices (in an EFM) and extracts as a sub-model (SBML file) which can be visualised in Cytoscape
%
% USAGE:
%    modelEFM = extractSBMLFromEFM(model, data, outputFileName); % without optional inputs
%    modelEFM = extractSBMLFromEFM(model, data, outputFileName, remFlag, ubiquitousMets); % with optional inputs
%
% INPUTS:
%    model:              COBRA model that was used for EFM calculation
%    rxnsForSubmodel:    array of reaction indices (such as those in an EFM)
%    outputFileName:     file name of the output sbml file e.g., <name>.xml
%
% OPTIONAL INPUTS:
%    remFlag:           boolean indicating whether ubiquitous metabolites should be removed
%    ubiquitousMets:    list of metabolites to remove.
%                       An example list of ubiquitous metabolites:
%                       ubiquitousMets = {'h';'h2';'fe3';'fe2';'h2o';'na1';'atp';'datp';'hco3';'pi';'adp';'dadp';'nadp';'nadph';'coa';'o2';'nad';'nadh';'ppi';'pppi';'amp';'so4';'fad';'fadh2';'udp';'dudp';'co2';'h2o2';'nh4';'ctp';'utp';'gtp';'cdp';'gdp';'dcdp';'dgdp';'dtdp';'dctp';'dttp';'dutp';'dgtp';'cmp';'gmp';'ump';'dcmp';'dtmp';'dgmp';'dump';'damp';'cl'};
%
% OUTPUTS:
%    modelEFM:    submodel containing reactions, metabolites and genes in the EFM of interest
%
% EXAMPLE:
%    data = [1 2 3 4 5]; % select first 5 reactions
%    modelEFM = extractSBMLFromEFM(model, data) ;
%
% .. Author: Last modified: Chaitra Sarathy, 1 Oct 2019

if nargin < 4
    remFlag = 0;
    ubiquitousMets = [];
end

%remove reactions other than those in EFM
rxnRemoveList = model.rxns(setdiff(1:length(model.rxns), rxnsForSubmodel));

% write a model with the reactions (and corresponding metabolites) from the
% EFM
modelEFM = removeRxns(model, rxnRemoveList);

if (remFlag)
    % Define ubiquitous metabolites
%     ubiquitousMets = {'h';'h2';'fe3';'fe2';'h2o';'na1';'atp';'datp';'hco3';'pi';'adp';'dadp';'nadp';'nadph';'coa';'o2';'nad';'nadh';'ppi';'pppi';'amp';'so4';'fad';'fadh2';'udp';'dudp';'co2';'h2o2';'nh4';'ctp';'utp';'gtp';'cdp';'gdp';'dcdp';'dgdp';'dtdp';'dctp';'dttp';'dutp';'dgtp';'cmp';'gmp';'ump';'dcmp';'dtmp';'dgmp';'dump';'damp';'cl'};

    % Get compartments in the model so that ubiquitous molecules from all compartments can be removed
    [~, uniqComps] = getCompartment(model.mets);
    temp = cellfun(@(x) strcat(x, '[', uniqComps, ']'), ubiquitousMets,'UniformOutput' ,0);
    listOfAbundantMets = vertcat(temp{:});

    % Remove metabolites and unused genes, if any
    modelEFM = removeMetabolites(modelEFM, listOfAbundantMets, true);
end

% the unused genes in the model do not get removed, so remove mannually
modelEFM = removeUnusedGenes(modelEFM);

writeCbModel(modelEFM, 'format','sbml', 'fileName', outputFileName);

end

author: opencobra

repo: cobratoolbox


text:
function varargout = drawEdgeLabels(p, e, value)
%DRAWEDGELABELS Draw values associated to graph edges
%
%   usage:
%   drawEdgeLabels(NODES, EDGES, VALUES);
%   NODES: array of double, containing x and y values of nodes
%   EDGES: array of int, containing indices of in and out nodes
%   VALUES is an array the same length of EDGES, containing values
%   associated to each edges of the graph.
%
%   The function computes the center of each edge, and puts the text with
%   associated value.
%
%   H = drawEdgeLabels(...) return array of handles to each text structure,
%   making possible to change font, color, size
%
%
%   -----
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2003.
%

%   HISTORY
%   10/03/2004 included into lib/graph library

if length(p) > 1 && length(e) > 1
    h = zeros(length(e), 1);
    hold on;
    for l=1:length(e)
        % indices of source and target nodes
        n1 = e(l, 1);
        n2 = e(l, 2);

        % node coordinates
        x1 = p(n1, 1);
        y1 = p(n1, 2);
        x2 = p(n2, 1);
        y2 = p(n2, 2);

        % display the edge
        line([x1 x2], [y1 y2]);

        % coordinates of edge label
        xm = (x1 + x2)/2;
        ym = (y1 + y2)/2;

        % display label
        h(l) = text(xm, ym, sprintf('%3d', floor(value(l))));
    end
end

if nargout == 1
    varargout = {h};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function iccvWriteSupplementalResults(fn)

gtdir1 = './iccvGroundTruth/gtsave/';
ourdir1 = './data4/smallsegs/3/';
ncutdir1 = './data4/ncutsegs/';
geomdir1 = './data4/geomsegs/';

gtdir2 = './results/supp/gt/';
ourdir2 = './results/supp/result/';
ncutdir2 = './results/supp/ncut/';
geomdir2 = './results/supp/geom/';

imdir = './iccvGroundTruth/images/';

for f = 1:numel(fn)

    disp(num2str(f))

    bn = strtok(fn{f}, '.');

    im = im2double(imread([imdir bn '.jpg']));

    %% gt
    if 0
    load([gtdir1 bn '_gt.mat']);
    outname = [gtdir2 bn '_gt.jpg'];

    contact = {}; %cell(bndinfo.nseg, 1);

    bndinfo2 = updateBoundaryInfo2(bndinfo, bndinfo.labels);

    bndinfo2.names = bndinfo.names;
    bndinfo2.type = bndinfo.type;
    bndinfo2 = orderfields(bndinfo2);

    bndinfo2 = processGtBoundaryLabels(bndinfo2);

    lab = bndinfo2.edges.boundaryType;
    lab = (lab(1:end/2)>0) + 2*(lab(end/2+1:end)>0);

    printOcclusionResultWithContact(im, bndinfo2, lab, outname, contact, 1);
    end

    %% results
    resdirs = {ourdir1, geomdir1, ncutdir1};
    outdirs = {ourdir2, geomdir2, ncutdir2};

    for k = 1 %:numel(resdirs)

        load([resdirs{k} bn '_seg.mat']);
        %try
            if k==1
                writeDepthImages(resdirs{k}, fn(f), outdirs{k});
                load([outdirs{k} bn '_contact.mat']);
            end
%         catch
%             lasterr
%             contact = {};
%         end
if 0
        try
            lab = bndinfo.edges.boundaryType;
            lab = lab(1:end/2) + 2*lab(end/2+1:end);
        catch
            lab = [];
        end

        outname = [outdirs{k} bn '_res.jpg'];
        printOcclusionResultWithContact(im, bndinfo, lab, outname, contact, 1);
end
    end
end

author: Cloud-CV

repo: object-proposals


text:
function [u,v,w]=v_sphrharm(m,a,b,c,d)
%V_SPHRHARM  forward and inverse spherical harmonic transform
%
% Usage: (1) y=('f',n,x)      % Calculate complex transform of spatial data x up to order n
%
%        (2) y=('fr',n,x)     % Calculate real transform of spatial data x(ne,na) up to order n
%                             % x is specified on a uniform grid with inclination e=(0.5:ne-0.5)*pi/ne
%                             % and azimuth a=(0:na-1)*2*pi/na. The North pole has inclination 0 and
%                             % azimuths increase going East.
%
%        (3) y=('fd',n,e,a,v) % Calculate transform of impulse at (e(i),a(i)) of height v(i) up to order n
%                             % e(i), a(i) and v(i) should be the same dimension; v defaults to 1.
%
%        (4) x=('i',y,ne,na)  % Calculate spatial data from spherical harmonics y on
%                             % a uniform grid with ne inclinations and na azimuths
%
%        (5) [e,a,w]=('cg',ne,na)   % Calculate the inclinations, azimuths and
%                                   % quadrature weights of a Gaussian sampling grid
%
%        (6) n=2;             % illustrate real basis functions upto order 2
%            for i=0:n
%              for j=-i:i
%                subplot(n+1,2*n+1,i*(2*n+1)+j+n+1);
%                v_sphrharm('irp',[zeros(1,i^2+i+j) 1],25,25);
%              end
%            end
%
%  Inputs:  m string specifying various options:
%               'f' forward transform
%               'i' inverse transform
%               'c' calculate the coordinates of the sampling grid [default if f or i omitted]
%               'r' real spherical harmonics instead of complex
%               'u' uniform inclination grid:  e=(0.5:ne-0.5)*pi/ne (includes neither pole) [default]
%                      for invertibility, ne>=2N+1 for order N
%               'U' uniform inclination grid:  e=(0:ne-1)*pi/ne (includes North pole only)
%                      for invertibility, ne>=2N+1 for order N
%               'g' gaussian inclination grid (non-uniform but fewer samples needed)
%                      for invertibility, ne>=N+1 for order N
%               'a' arbitrary (specified by user - inverse transform only)
%               'd' delta function [forward transform only]
%               'p' plot result
%               'P' plot with colourbar
%
%           The remaining inputs depend on the mode specified:
%               'f'  a         order of transform
%                    b(ne,na)  input data array on the chosen inclination-azimuth sampling grid
%               'fd' a         order of transform
%                    b(k)      inclinations of k delta functions
%                    c(k)      azimuths of k delta functions
%                    d(k)      amplitudes of k delta functions [default=1]
%               'i'  a()       spherical harmonics as a single row vector in the order:
%                                 (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
%                              To access as a 2-D harmonic: Y(n,m)=a(n^2+n+m+1) where n=0:N, m=-n:n
%                    b         number of inclination values in output grid
%                    c         number of azimuth values in output grid
%               'c'  a         number of inclination values in output grid
%                    b         number of azimuth values in output grid
%
% Outputs:  u output data according to transform direction:
%                'f': spherical harmonics as a single row vector in the order:
%                        (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
%                     To access as a 2-D harmonic: Y(n,m)=u(n^2+n+m+1) where n=0:N, m=-n:n
%                'i': u(ne,na) is spatial data sampled on an azimuth-inclination grid
%                     with ne inclination points (in 0:pi) and na azimuth points (in 0:2*pi)
%                'c': u(ne) gives inclination grid positions with 0 being the North pole
%           v(na) gives azimuth grid positions
%           w(ne) gives the quadrature weights used in the transform
%
% Suppose f(e,a) is a complex-valued function defined on the surface of the
% sphere (0<=e<=pi, 0<=a<2pi) where e=inclination=pi/2-elevation and
% a=azimuth. (0,*) is the North pole, (pi/2,0) is on the equator near
% Ghana and (pi/2,pi/2) is on the equator near India.
%
% We can approximate f(e,a) using complex spherical harmonics as
% f(e,a) = sum(F(n,m)*u(n,m,e,a)) where the sum
% is over n=0:N,m=-n:n giving a total of (N+1)^2 coefficients F(n,m).
%
% If f(e,a) happens to be real-valued, then we can transform the
% coefficients using G(n,0)=F(n,0) and, for m>0,
% [G(n,m); G(n,-m)]=sqrt(0.5)*[1 1; 1i -1i]*[F(n,m); F(n,-m)]
% to give the real spherical harmonic coefficients G(n,m).
%
% The basis functions u(n,m,e,a) are orthonormal under the inner product
% <u(e,a),v(e,a)> = integral(u(e,a)*v'(e,a)*sin(e)*de*da).
%
%  Minimum spatial grid for invertibility:
%     Gaussian grid: ne >= N+1,  na >= 2N+1
%     Uniform grid:  ne >= 2N+1, na >= 2N+1
%     An inverse transform followed by a forward transform will restore the
%     original coefficient values provided that the sampling grid is large
%     enough. The advantage of the Gaussian grid is that it can be smaller.
%
%   Data formats:
%     (1) Spatial Data: x=v_sphrharm('i',y,ne,na)
%            x(1:ne,1:na) is spatial data sampled on an azimuth-inclination grid
%                         with ne inclination points (in 0:pi) and
%                         na azimuth points (in 0:2*pi)
%     (2) Spherical harmonics: y=v_sphrharm('f',n,x)
%            y(1:(n+1)^2)  spherical harmonics as a single row vector
%                          in the order: (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
%                          To access as a 2-D harmonic, use
%                          Y(n,m)=y(n^2+n+m+1)   where n=0:N, m=-i:i
%     (3) Sampling Grid:  [e,a,w]=v_sphrharm('c',ne,na)
%            e(1:ne)   monotonically increasing inclination values (0<=e<=pi) where 0 is the North pole
%            a(1:na)   monotonically increasing azimuth values (0<=a<2pi)
%            w(1:ne)   Quadrature weights
%

% future options
%    direction: [m=v_rotation matrix]
%    transform: [h=complex harmonics but only the positive ones specified]
%    grid       d=delta function [forward transform only]
%               [s=sparse azimuth]
%               [z=include both north and south pole]
%
% bugs:
%   (1) we could save space and time by taking advantage of symmetry of legendre polynomials
%   (2) the number of points in the inclination (elevation) grid must, for now, be even if the 'U' option is chosen
%   (3) should ensure that the negative nyquist azimuth frequency is not used
%   (4) save time by manipulating only the necessary 2*m columns of the da matrix
%   (5) should make non-existant coefficients black in plots
%   (6) using 'surf' for plots adds an offset in azimuth and elevation
%       because colours correspond to vertices not faces
%   (7) the normalization for mode 'fd' seems incorrect
%   (8) mode 'fd' should allow multiple impulses to be summed

% errors:

% tests:
% check for inverse transform n=4; m=4; [ve,va]=spvals(8,8,n,m); ve*va-v_sphrharm('iur',[zeros(1,n^2+n+m) 1],8,8)
% check inverse followed by forward: no=4; h=rand(1,(no+1)^2); max(abs(v_sphrharm('fur',v_sphrharm('iur',h,10,10),no)-h))
% same but complex: no=4; h=rand(1,(no+1)^2); max(abs(v_sphrharm('fu',v_sphrharm('iu',h,10,10),no)-h))
% same but gaussian grid: no=4; h=rand(1,(no+1)^2); max(abs(v_sphrharm('fg',v_sphrharm('ig',h,6,10),no)-h))

%      Copyright (C) Mike Brookes 2009
%      Version: $Id: v_sphrharm.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% decode option string
mv='c u ';       % mv(1)=transform, mv(2)=real/complex, mv(3)=grid, mv(4)=plot
if ~nargout
    mv(4)='P';
end
mc='ficruUgadpP';
mi=[1 1 1 2 3 3 3 3 3 4 4];
for i=1:length(m)
    j=find(m(i)==mc);
    if ~isempty(j)
        mv(mi(j))=m(i);
    end
end

switch mv(1)
    case 'f' % forward transform [sp]=('f',order,data)
        if mv(3)~='d'
            % input data has elevations down the columns and azimuths along the rows
            % output data is in the form:
            [ne,na]=size(b);
            da=fft(b,[],2)*sqrt(pi)/na;
            if mv(2)=='r' % only actually need to do this for min(a,floor(na/2)) values (but nyquist is tricky)
                ix=2:ceil(na/2);
                iy=na:-1:na+2-ix(end);
                da(:,ix)=da(:,ix)+da(:,iy);
                da(:,iy)=1i*(da(:,ix)-2*da(:,iy));  % note
                da(:,1)=da(:,1)*sqrt(2);
            else
                da=da*sqrt(2);
            end
            [ue,we,lgp]=sphrharp(mv(3),ne,a);
            da=da.*repmat(we',1,na);
            u=zeros(1,(a+1)^2);
            i=0;
            for nn=0:a
                % we could vectorize this to avoid the inner loop
                % we should ensure the the negative nyquist value is not used
                for mm=-nn:nn
                    i=i+1;
                    u(i)=lgp(1+nn*(nn+1)/2+abs(mm),:)*da(:,1+mod(mm,na));
                end
            end
        else       % forward transform of impulses [sp]=('fd',order,e,a,value)
            if nargin<5
                d=ones(size(b));
            end
            [ue,we,lgp]=sphrharp('a',b(:),a);
            uu=zeros(1,(a+1)^2);     % reserve space for spherical harmonic coefficients
            u=uu;
            for j=1:numel(b)
                i=0;
                if mv(2)=='r'
                    for nn=0:a
                        % we could vectorize this to avoid the inner loop
                        % we should ensure the the negative nyquist value is not used
                        i=i+1;
                        uu(i+nn)=lgp(1+nn*(nn+1)/2,1);
                        for mm=1:nn
                            uu(i+nn-mm)=lgp(1+nn*(nn+1)/2+abs(mm),j)*sin(mm*c(j))*sqrt(2);
                            uu(i+nn+mm)=lgp(1+nn*(nn+1)/2+abs(mm),j)*cos(mm*c(j))*sqrt(2);
                        end
                        i=i+2*nn;
                    end
                else
                    for nn=0:a
                        % we could vectorize this to avoid the inner loop
                        % we should ensure the the negative nyquist value is not used
                        for mm=-nn:nn
                            i=i+1;
                            uu(i)=lgp(1+nn*(nn+1)/2+abs(mm),j)*exp(-1i*mm*c(j));
                        end
                    end
                end
                u=u+d(j)*uu/sqrt(2*pi);
            end
        end
    case 'i' % inverse transform [data]=('i',sp,nincl,nazim)
        %or [data]=('a',sp,e,a) where e is a list of elevations and a is either a corresponding list of azimuths
        % or else a cell array contining several azimuths for each inclination
        % input data is sp=[(0,0) (1,-1) (1,0) (1,1) (2,-2) (2,-1) (2,0) ... ]
        % length is (n+1)^2 where n is the order
        % output data is an array [ne,na]
        nsp=ceil(sqrt(length(a)))-1;
        [ue,we,lgp]=sphrharp(mv(3),b,nsp);
        if mv(3)=='a'
            na=nsp*2+1;
            ne=length(b);
        else
            na=c;
            ne=b;
        end
        i=0;
        da=zeros(ne,na);
        for nn=0:nsp
            % this could be vectorized somewhat to speed it up
            for mm=-nn:nn
                i=i+1;
                if i>length(a)
                    break;
                end
                da(:,1+mod(mm,na))=da(:,1+mod(mm,na))+a(i)*lgp(1+nn*(nn+1)/2+abs(mm),:)';
            end
        end
        if na>1 && mv(2)=='r' % convert real to complex only actually need to do this for min(b,floor(na/2)) values (but nyquist is a bit tricky)
            ix=2:ceil(na/2);
            iy=na:-1:na+2-ix(end);
            da(:,iy)=(da(:,ix)+1i*da(:,iy))*0.5;
            da(:,ix)=da(:,ix)-da(:,iy);
            da(:,1)=da(:,1)/sqrt(2);
        else
            da=da/sqrt(2);
        end
        if mv(3)=='a' % do a slow fft
            if iscell(c)
                u{ne,1}=[]; % reserve space for output cell array
                for i=1:ne
                    ai=c{i};
                    ui=repmat(da(i,1),1,length(ai));
                    for j=1:nsp
                        exj=exp(1i*j*ai(:).'); % we could vectorize this more
                        ui=ui+da(i,j+1).'.*exj+da(i,na+1-j).'.*conj(exj);
                    end
                    u{i}=ui/sqrt(pi);
                end
            else
                u=da(:,1);
                for j=1:nsp
                    exj=exp(1i*j*c(:)); % we could vectorize this more
                    u=u+da(:,j+1).*exj+da(:,na+1-j).*conj(exj);
                end
                u=u/sqrt(pi);
                if mv(2)=='r' && isreal(a)
                    u=real(u);
                end
            end
        else
            u=ifft(da,[],2)*na/sqrt(pi); % could put the scale factor 1/sqrt(pi) earlier
            if mv(2)=='r' && isreal(a)
                u=real(u);
            end
        end
    case 'c' % just output coordinates [inclination,azim,weights]=('c',nincl,nazim)
        % if m!='g', then order, a, must be even
        [u,w]=sphrharp(mv(3),a,0);
        if nargin<3
            b=ceil(a/1+(mv(3)=='g'));
        end
        v=(0:b-1)*2*pi/b;
end
if mv(4)~=' '
    switch mv(1)
        case 'f'
            if mv(2)=='r'
                ua=u;
                tit='Real Coefficients';
            else
                ua=abs(u);
                tit='Complex Coefficient Magnitudes';
            end
            nu=length(ua);
            no=ceil(sqrt(nu))-1;
            yi=zeros(no,2*no+1);
            for i=0:no
                for j=-i:i
                    iy=i^2+i+j+1;
                    if iy<=nu
                        yi(i+1,j+no+1)=ua(iy);
                    end
                end
            end
            imagesc(-no:no,0:no,yi);
            axis 'xy';
            if mv(4)=='P'
                colorbar;
            end
            xlabel('Harmonic');
            ylabel('Order');
            title(tit);
        case 'i'            % [data]=('i',sp,nincl,nazim)
            vv=(0:na)*2*pi/na;  % azimuth array
            gr=sin(ue)';
            surf(gr*cos(vv),gr*sin(vv),repmat(cos(ue)',1,na+1),real(u(:,[1:na 1])));
            axis equal;
            xlabel('X');
            ylabel('Y');
            zlabel('Z');
            if mv(4)=='P'
                colorbar;
            end
            %             v_cblabel('Legendre Weight');
            %             title(sprintf('Sampling Grid: %s, %d, %d',mv(3),length(u),na-1));
        case 'c'            % [inclination,azim,weights]=('c',nincl,nazim)
            vv=[v v(1)];    % replicate initial azimuth point
            na=length(vv);
            gr=sin(u)';
            mesh(gr*cos(vv),gr*sin(vv),repmat(cos(u)',1,na),repmat(w',1,na));
            axis equal;
            xlabel('X');
            ylabel('Y');
            zlabel('Z');
            if mv(4)=='P'
                colorbar;
                v_cblabel('Quadrature Weight');
            end
            title(sprintf('Sampling Grid: %s, %d, %d',mv(3),length(u),na-1));
    end
end

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

function [ueo,weo,lgpo]=sphrharp(gr,ne,nsp)
% calculate persistent variables for grid points and legendre polynomials
% we recalculate if transform order or inclination grid size or type are changed
% gr = grid type:
%       'u'=uniform inclination starting at pi/2n (default)
%       'U'=uniform but starting at north pole
%       'g'=gaussian inclination
%       'a'=arbitrary (specified by user - inverse transform only)
% ne = number of inclination values
% nsp = maximum order
% Outputs:
%    ueo(ne)     vector containing the inclination values
%    weo(ne)     vector containing the quadrature weights
%    lgpo((nsp+1)*(nsp+2)/2,ne)  evaluates the Schmidt seminormalized
%                                associated Legendre function at each value
%                                of cos(ue) for n=0:nsp, m=0:n.
persistent gr0 lgp ne0 ue we nsp0
if isempty(ne0)
    ne0=-1;
    nsp0=-1;
    gr0='a';
end
if gr=='a'
    ue=ne;
    ne=length(ue);
    ne0=-1;
    gr0=gr;
    nsp0=-1; % delete previous legendre polynomials
    lgp=[];
    we=[];
elseif gr~=gr0 || ne~=ne0
    if gr=='g'
        r = 1:ne-1;
        r = r ./ sqrt(4*r.^2 - 1);
        p = zeros( ne, ne );
        p( 2 : ne+1 : ne*(ne-1) ) = r;
        p( ne+1 : ne+1 : ne^2-1 ) = r;
        [q, ue] = eig(p);
        [ue, k] = sort(diag(ue));
        ue=acos(-ue)';
        we = 2*q(1,k).^2;
    elseif gr=='U'
        if rem(ne,2)
            error('inclination grid size must be even when using ''U'' option');
        end
        ue=(0:ne-1)*pi/ne;
        xx=zeros(1,ne);
        ah=ne/2;
        xx(1:ah)=(1:2:ne).^(-1);
        we=-4*sin(ue).*imag(fft(xx).*exp(-1i*ue))/ne;
    else % default is m='u'
        ue=(1:2:2*ne)*pi*0.5/ne; % vector of elevations
        vq=(ne-abs(ne+1-2*(1:ne))).^(-1).*exp(-1i*(ue+0.5*pi));
        we=(-2*sin(ue).*real(fft(vq).*exp(-1i*(0:ne-1)*pi/ne))/ne);
    end
    gr0=gr;
    ne0=ne;
    nsp0=-1; % delete previous legendre polynomials
    lgp=[];
end
if nsp>nsp0
    lgp((nsp+1)*(nsp+2)/2,ne)=0; % reserve space
    for i=nsp0+1:nsp
        lgp(1+i*(i+1)/2:(i+1)*(i+2)/2,:)=legendre(i,cos(ue),'sch')*sqrt(0.5*i+0.25);
        lgp(1+i*(i+1)/2,:)=lgp(1+i*(i+1)/2,:)*sqrt(2);
    end
    nsp0=nsp;
end
lgpo=lgp;
weo=we;
ueo=ue;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [pE,pC] = spm_null_priors(A,B,C)
% prior moments for null (Jacobian) model
% FORMAT [pE,pC] = spm_null_priors(A,B,C)
%
% A{1},B{m},C  - binary constraints on extrinsic connections
%
% pE - prior expectation
% pC - prior covariance
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_null_priors.m 5908 2014-03-05 20:31:57Z karl $

% default: a single source model
%--------------------------------------------------------------------------
if nargin < 3
    A   = {1};
    B   = {};
    C   = 1;
end

% orders
%--------------------------------------------------------------------------
n     = size(C,1);                                % number of sources
u     = size(C,2);                                % number of inputs

% parameters for Jacobian
%==========================================================================

% canonical source
%--------------------------------------------------------------------------
a     = 8;
Hz    = [8 16 48];
for i = 1:length(Hz)
   b      = 2*pi*Hz(i);
   J{i,i} = [a b;-b a];
end

% Jacobian
%--------------------------------------------------------------------------
J     = full(spm_cat(J));
nx    = length(J);
pE.A  = logm(kron(eye(n,n),J));
pC.A  = kron(A{1},ones(nx,nx));

% Bilinear terms
%--------------------------------------------------------------------------
for i = 1:u
    pE.B{i} = kron(B{i},zeros(nx,nx));
    pC.B{i} = kron(B{i},ones(nx,nx));
end

% input coeficicents
%--------------------------------------------------------------------------
pE.C  = kron(C,0);
pC.C  = C;

% input coeficicents
%--------------------------------------------------------------------------
pE.D  = kron(zeros(n,n),ones(nx,1));
pC.D  = kron(speye(n,n),ones(nx,1));

author: spm


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: test_prop_diff.m
%
% use: run "test_prop_diff" and then examine "test_prop_diff_results.csv" in the work folder and
%      compare it to the "test_prop_diff_results_reference.csv" distributed with the code.
%      Failure to run or significant differences in results may indicate installation problems.
%
% change log:
%   030106 tdr created
%
% function y = prop_diff(x1,n1,x2,n2,delta)
% x - number of positive events
% n - number of trials
% delta - the difference of interest
% y - Pr (p1 - p2 >= delta), where p1_hat = x1/n1 and p2_hat = x2/n2

function test_prop_diff

n1s = [1 30 100 1e5];
n2s = n1s;
deltas = [-0.5 0 1e-2 0.9];

fid = fopen('test_prop_diff_results.csv','w');
fprintf(fid,'%s\n','delta, n1, n2, x1, x2, p1_hat, p2_hat, Pr(p1 - p2 >= delta)');

for delta = deltas
   for n1 = n1s
      if n1 < 10
         x1s = [0 1 round(n1/2) n1];
      else
         x1s = [0 1 round(n1/3) round(n1/2) n1-1 n1];
      end;
      for n2 = n2s
         if n2 < 10
            x2s = [0 1 round(n2/2) n2];
         else
            x2s = [0 1 round(n2/3) round(n2/2) n2-1 n2];
         end;
         for x1 = x1s
            for x2 = x2s
               result_array = [delta n1 n2 x1 x2 x1/n1 x2/n2 prop_diff(x1,n1,x2,n2,delta)];
               fprintf(fid,'%8.6f, %d, %d, %d, %d, %8.6f, %8.6f, %8.6f\n',result_array);
		    end;
         end;
      end;
   end;
end;

fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
clear all;
clc;

addpath util;

load('profile.mat');

count=length(profile);

pts_num = 39;

% for i=1:length(profile)
%     profile{i}.tmp=LoadPTS(sprintf('./profile/%s/%s.pts', profile{i}.status, profile{i}.name));
% end

gt=zeros(pts_num,2,count);
for i=1:length(profile)
    gt(:,:,i)=profile{i}.gt;
end

% num = 5;
num = 4;

spacing=0.0005;
sampling = 0:spacing:0.35;
x = sampling;
y = repmat(x',1,num);
profile_errs=cell(num,1);
det=zeros(pts_num,2,count);

% participants = {'yang','he','wu','deng','tmp'};
participants = {'yang','he','wu','deng'};

for iter = 1:num
    for i=1:length(profile)
        det(:,:,i)=profile{i}.(participants{iter});
    end
    profile_errs{iter} = compute_diag_error( gt, det );
    for i=1:numel(sampling)
        y(i,iter) = sum(profile_errs{iter} < sampling(i)) / numel(profile_errs{iter});
    end
end

title = 'Menpo Profile Test Set';

% methods={'J. Yang et al','Z. He et al', 'W. Wu et al', 'J. Deng et al','tmp'};
methods={'J. Yang et al','Z. He et al', 'W. Wu et al', 'J. Deng et al'};

colors = lines(num);

x_limit = 0.03;

linewidth = 3;

fontsize = 12;

plot_ced(x, y, methods, strrep(title,'_',' '), ...
         x_limit, colors, linewidth, fontsize);

ax = gca;
outerpos = ax.OuterPosition;
ti = ax.TightInset;
left = outerpos(1) + ti(1);
bottom = outerpos(2) + ti(2);
ax_width = outerpos(3) - ti(1) - ti(3);
ax_height = outerpos(4) - ti(2) - ti(4);
ax.Position = [left bottom ax_width ax_height];

fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)]* 1.05;

print(fig,'-dpdf',[ title '.pdf']);

author: jiankangdeng

repo: MenpoBenchmark


text:
function model = createToyModelForLooplessFVA()
% Create a toy model that has a loop and will give different solutions when
% minimizing the 0-, 1- and 2-norms respectively
%
%    <=> B    10  10   10  10
%         \    <===> F <===>
% <=> A ----> D -----------> E <=>
%      2 2/
%       2/
%   <=> C

model = createModel();
Reactions = {'R1', 'A + B -> D'; ...
             'R2', '2 A + 2 C -> 2 D';...
             'R3', 'D -> E';...
             'R4', '10 D <=> 10 F';...
             'R5', '10 F <=> 10 E';...
             'Ex_A', 'A <=>'; ...
             'Ex_B', 'B <=>'; ...
             'Ex_C', 'C <=>'; ...
             'Ex_E', 'E <=>'};

%Add Reactions
for i = 1:size(Reactions,1)
    %All reactions are irreversible
    model = addReaction(model, Reactions{i,1}, 'reactionFormula', Reactions{i,2}, 'printLevel', -1);
end

% uptake bound for A = 1
model = changeRxnBounds(model, 'Ex_A', -1, 'l');
% objective: max production of E
model = changeObjective(model, 'Ex_E', 1);
end


text:
function [ pn, dist, t ] = segment_point_near ( p1, p2, p )

%*****************************************************************************80
%
%% SEGMENT_POINT_NEAR finds the line segment point nearest a point.
%
%  Discussion:
%
%    A line segment is the finite portion of a line that lies between
%    two points.
%
%    The nearest point will satisfy the condition
%
%      PN = (1-T) * P1 + T * P2.
%
%    T will always be between 0 and 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), the endpoints of the line segment.
%
%    Input, real P(2,1), the point whose nearest neighbor
%    on the line segment is to be determined.
%
%    Output, real PN(2,1), the point on the line segment which is
%    nearest the point (X,Y).
%
%    Output, real DIST, the distance from the point to the
%    nearest point on the line segment.
%
%    Output, real T, the relative position of the point (XN,YN)
%    to the points (X1,Y1) and (X2,Y2).
%

%
%  If the line segment is actually a point, then the answer is easy.
%
  if ( p1(1:2,1) == p2(1:2,1) )

    t = 0.0;

  else

    bot = sum ( ( p2(1:2,1) - p1(1:2,1) ).^2 );

    t = ( p(1:2,1) - p1(1:2,1) )' * ( p2(1:2,1) - p1(1:2,1) ) / bot;

    t = max ( t, 0.0 );
    t = min ( t, 1.0 );

  end

  pn(1:2,1) = p1(1:2,1) + t * ( p2(1:2,1) - p1(1:2,1) );

  dist = sqrt ( sum ( ( pn(1:2,1) - p(1:2,1) ).^2 ) );

  return
end


text:
function drawRobotStopAction(q,p)
% drawRobotStopAction(q,p)
%
% This function draws the robot with configuration q and parameters p over
% a series of configurations, leaving hold on so that the motion can be
% seen in a single (static) frame.
%
% INPUTS:
%   q = [5, 1] = column vector of a single robot configuration
%   p = parameter struct
%

nTime = size(q,2);

% Compute the points that will be used for plotting
[P, G] = getPoints(q,p);

x = [zeros(1,nTime); P(1:2:end,:)];
y = [zeros(1,nTime); P(2:2:end,:)];

P1 = P(1:2,:);
P2 = P(3:4,:);
P3 = P(5:6,:);
P4 = P(7:8,:);
P5 = P(9:10,:);

% G1 = G(1:2,:);
% G2 = G(3:4,:);
% G3 = G(5:6,:);
% G4 = G(7:8,:);
% G5 = G(9:10,:);

% Heuristics:
L = (p.l1 + p.l2);  % Maximum extended leg length
xBnd = L*[-1.2,1.2];
yBnd = [-0.2*L, L + p.l3];

% Colors:
colorGround = [118,62,12]/255;
colorStance = [200,60,60]/255;
colorSwing = [60,60,200]/255;
colorTorso = [160, 80, 160]/255;

% Plot parameters:
legWidth = 2;
jointSize = 20;

% Set up the figure
hold off;

% Plot the ground:
plot(xBnd,[0,0],'LineWidth',6,'Color',colorGround);

hold on;

% Plot the links:
for i=1:nTime
plot(x(1:3,i),y(1:3,i),'--','LineWidth',legWidth,'Color',colorStance);
plot(0, 0,'k.','MarkerSize',jointSize,'Color',colorStance);
plot(P1(1,i), P1(2,i),'k.','MarkerSize',jointSize,'Color',colorStance);
plot(P2(1,i), P2(2,i),'k.','MarkerSize',jointSize,'Color',colorStance);
% plot(G1(1,i), G1(2,i),'ko','MarkerSize',8,'LineWidth',2);
% plot(G2(1,i), G2(2,i),'ko','MarkerSize',8,'LineWidth',2);
end
for i=1:nTime
plot(x(3:4,i),y(3:4,i),'LineWidth',legWidth+1,'Color',colorTorso);
plot(P3(1,i), P3(2,i),'k.','MarkerSize',jointSize,'Color',colorTorso);
% plot(G3(1,i), G3(2,i),'ko','MarkerSize',8,'LineWidth',2,'Color',colorTorso);
end
for i=1:nTime
plot(x([3,5,6],i),y([3,5,6],i),'LineWidth',legWidth,'Color',colorSwing);
plot(P4(1,i), P4(2,i),'k.','MarkerSize',jointSize,'Color',colorSwing);
plot(P5(1,i), P5(2,i),'k.','MarkerSize',jointSize,'Color',colorSwing);
% plot(G4(1,i), G4(2,i),'ko','MarkerSize',8,'LineWidth',2,'Color',colorSwing);
% plot(G5(1,i), G5(2,i),'ko','MarkerSize',8,'LineWidth',2,'Color',colorSwing);
end

% Format the axis:
axis([xBnd,yBnd]); axis equal; axis off;

end

author: MatthewPeterKelly

repo: OptimTraj


%  Script file: plotsinc.m
%
%  Purpose:
%    This program illustrates the use of handle graphics
%    commands by creating a plot of sinc(x) from -3*pi to
%    3*pi, and modifying the characteristics of the figure,
%    axes, and line using the "set" function.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    04/02/07    S. J. Chapman        Original code
%
% Define variables:
%   hndl         -- Handle of line
%   x            -- Independent variable
%   y            -- sinc(x)

% Calculate sinc(x)
x = -3*pi:pi/10:3*pi;
y = sin(x) ./ x;

% Find the zero value and fix it up.  The zero is
% located in the middle of the x array.
index = fix(length(y)/2) + 1;
y(index) = 1;

% Plot the function.
hndl = plot(x,y);

% Now modify the figure to create a pink background,
% modify the axis to turn on y-axis grid lines, and
% modify the line to be a 2-point wide orange line.
set(gcf,'Color',[1 0.8 0.8]);
set(gca,'YGrid','on');
set(hndl,'Color',[1 0.5 0],'LineWidth',3);


text:
function [spectrogramPosition, tracePosition] = calculatePanelPositions(numchannels, channelNum, fractionalSpectrogramHeight, frameLeft, frameBottom, totalWidth, totalHeight);
debug.printfunctionstack('>');
channelHeight 		= totalHeight/numchannels;
spectrogramHeight 	= fractionalSpectrogramHeight * channelHeight;
traceHeight 		= channelHeight - spectrogramHeight;
spectrogramBottom   = frameBottom + (numchannels - channelNum) * channelHeight;
traceBottom         = spectrogramBottom + spectrogramHeight;
spectrogramPosition = [frameLeft, spectrogramBottom, totalWidth, spectrogramHeight];
tracePosition 		= [frameLeft, traceBottom      , totalWidth, traceHeight];
debug.printfunctionstack('<');

author: geoscience-community-codes

repo: GISMO


text:
function [a,b] = fdct_usfft_pos2idx(NX,NY,SX,SY,s,w,x,y)

%fdct_usfft_pos2idx.m - For a fixed scale and fixed direction, returns
%		the curvelet which is closest to a certain point on the image
%
% Inputs
%   NX,NY,SX,SY     Values returned by fdct_usfft_param
%   s               scale index
%   w               wedge (angular) index
%   x,y             position in image
%
% Outputs
%   a,b             Index of the curvelet at scale s and angle w which is nearest to (x,y)
%

  nx = NX{s}{w};  ny = NY{s}{w};
  bx = SX{s}{w}(1,1);  by = SY{s}{w}(1,1);
  sx = [SX{s}{w}(2,1)-bx, SY{s}{w}(2,1)-by];
  sy = [SX{s}{w}(1,2)-bx, SY{s}{w}(1,2)-by];
  tmp = 1 + [x-bx,y-by] / [sx; sy];
  a = round(tmp(1));  b = round(tmp(2));
  a = mod(a,nx);  b = mod(b,ny);
  if(a<1)    a = a+nx;  end
  if(a>nx)    a = a-nx;  end
  if(b<1)    b = b+ny;  end
  if(b>ny)    b = b-ny;  end

author: thomaskuestner


text:
function H_normalized = spd_normalization(H)

    n = 1 ./ sqrt(diag(H));
    H_normalized = diag(n) * H * diag(n);

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function str = fillStr(str,len,varargin)

if check_option(varargin,'left')
  str = [repmat(' ',1,max(0,len-length(str))),str];
else
  str = [str, repmat(' ',1,max(0,len-length(str)))];
end

author: mtex-toolbox

repo: mtex

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


text:
% angle to the closest roi
function [data,units] = compute_angle2closestroi2(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);

for i = 1:nflies,
  fly = flies(i);

  % roi closest to fly
  closestroi = trx(fly).closestroi2;

  % position of fly
  xnose_mm1 = trx(fly).xnose_mm;
  ynose_mm1 = trx(fly).ynose_mm;
  theta_mm1 = trx(fly).theta_mm;

  ROIdata=trx.roi2{n}.data{fly};

  % loop over all rois
  for j = 1:size(ROIdata,1),

    % frames where this roi is closest
    idx = find(closestroi == j);
    if isempty(idx), continue; end

    % angle to roi
    dx2 = ROIdata(j,1)-xnose_mm1(idx);
    dy2 = ROIdata(j,2)-ynose_mm1(idx);
    theta2 = atan2(dy2,dx2);

    % angle relative to fly's orientation
    data{i}(idx) = modrange(theta2 - theta_mm1(idx),-pi,pi);

  end
end

units = parseunits('rad');

author: kristinbranson

repo: JAABA


text:
 function y = ir_fftshift2(x)
%function y = ir_fftshift2(x)
%|
%| 2D fft shift of 2D+T object

if nargin < 1, ir_usage, end

y = fftshift(fftshift(x, 1), 2);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function h = compute_sparse_spike_filter(type,n,options)

% compute_sparse_spike_filter - compute a discrete seismic filter
%
%   h = compute_sparse_spike_filter(type,n,options);
%
%   type can be 'dergauss', 'seisfreq' or 'bump'
%   the width of the gaussian derivative is options.sigma, expressed in [-1,1]
%   so that the width in pixel is sigma*n/2.
%
%   Copyright (c) 2008 Gabriel Peyre

switch type
    case 'dergauss'
        x = linspace(-1,1,n+1)'; x(end)= [];
        sigma = getoptions(options,'sigma',.02);
        h = x.*exp( -(x.^2)/(2*sigma^2) );
        h = (1-x.^2/sigma^2).*exp( -(x.^2)/(2*sigma^2) );
        h = h-mean(h);
        %sigma = .005;
        %h = exp( -(x.^2)/(2*sigma^2) );
        h0 = h/max(h);
        h0 = h/norm(h,'fro');
        h0 = .1*h0;
        % switch for fft
        h = [h0(end/2+1:end); h0(1:end/2)];
    case {'seisfreq' 'wavelet'}
        rho = getoptions(options, 'rho', .1 * 1024/n);
        % compute by frequency
        x = linspace(-1,1,n+1);
        hf = sin(pi*x/rho).^2;
        hf(abs(x)>rho) = 0;
        hf = [hf(n/2+1:end) hf(1:n/2-1)];
        h = real(ifft(hf));
        h = h'/max(h);
        h = h/sum(abs(h));

    case 'bump'
        sigma = getoptions(options,'sigma',.07);
        x = linspace(-1,1,n+1)'; x(end)=[];
        h = exp( -x.^2 / (2*(sigma)^2) );
        h = h; % /norm(h);
        h = [h(end/2+1:end); h(1:end/2)];
end

author: gpeyre

repo: matlab-toolboxes


text:
function y = residuals(fit,type)

% Residuals (or a few other things) from a locfit() fit.
%
% Input arguments:
%   fit - the locfit() fit.
%   type (optional) type of residuals. Valid types are
%     'dev'    (deviance, the default)
%     'd2'     (deviance squared)
%     'pearson'(Pearson)
%     'raw'    (observed - fitted)
%     'ldot'   (derivative of log-likelihood)
%     'lddot'  (second derivative)
%     'fit'    (fitted values - no transformation)
%     'mean'   (fitted values - with back transformation)
%
%  Author: Catherine Loader.

if (nargin<2) type = 'dev'; end;

y = predict(fit,'d','restyp',type);

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function ranInt = myrandint(outputRow,outputCol,outputRange,varargin)
% MYRANDINT(M,N,RANGE) is an M-by-N matrix with random integer entries
% drawn with replacement from elements of vector RANGE.  The elements in
% vector RANGE do not need to be contiguous or unique.  (Actually, they do
% not even need to be integers: The function works the exact same way with
% noninteger elements, but a warning is generated to alert the user that
% noninteger elements are being sampled.)
%
% To specify a contiguous integer range from Xlow to Xhi, use RANGE = [Xlow:Xhi].
%
% MYRANDINT(M,N,RANGE,'noreplace') is an M-by-N matrix with random integers
% drawn without replacement.
%
% This function is based around RAND and RANDPERM, and is intended as a
% modest imitation of Comm Toolbox's RANDINT.

if isequal(size(outputRange),[1 2]) && ~isequal(outputRange(1),outputRange(2)-1),
    warning('To specify a range [low high] use [low:high].')
end
if ~isequal(round(outputRange),outputRange),
    warning('Specified RANGE contains noninteger values.')
end
if ~isequal(length(outputRange),length(outputRange(:))),
    error('Range must be a vector of integer values.')
end

numElements = outputRow*outputCol;

if isempty(varargin),

    ranInt = zeros(outputRow,outputCol);
    randIx = floor((length(outputRange))*rand(size(ranInt))) + 1;
    ranInt = outputRange(randIx);
    if ~isequal(size(randIx),size(ranInt)),
        ranInt = reshape(ranInt,size(randIx));
    end

elseif isequal(varargin{1},'noreplace'),

    if numElements > length(outputRange),
        error('Not enough elements in range to sample without replacement.')
    else
        % Generate full range of integers
        XfullShuffle = outputRange(randperm(length(outputRange)));
        % Select the first bunch:
        ranInt = reshape(XfullShuffle(1:numElements),outputRow,outputCol);
    end

else
    error('Valid argument is ''noreplace''.')
end