text:
classdef calc_across < ZmapVGridFunction
    % CALC_ACROSS calculate a-value parameters for a cross section

    properties
        mc_auto         McAutoEstimate      = ZmapGlobal.Data.UseAutoEstimate
        aval_method
    end

    properties(Constant)
        PlotTag         = 'calc_across'
        ReturnDetails   = cell2table({ ... VariableNames, VariableDescriptions, VariableUnits
            'a_value','a-value',''...
            }, 'VariableNames', {'Names','Descriptions','Units'})

        CalcFields      = {} % cell array of charstrings, matching into ReturnDetails.Names

        ParameterableProperties = []; % array of strings matching into obj.Properties
        References="";

        aval_methods = {' a(M=0) for fixed b',...
            'a(M=0) for Mc(MaxC) + M(corr)',...
            'a(M=0) for M(EMR)',...
            'a(M=0) by r1 & Mc by r2 '}
    end

    methods
        function obj=calc_across(zap, varargin)
            % CALC_ACROSS
            % obj = CALC_ACROSS() takes catalog, grid, and eventselection from ZmapGlobal.Data
            %
            % obj = CALC_ACROSS(ZAP) where ZAP is a ZmapAnalysisPkg

            obj@ZmapVGridFunction(zap, 'a_value');

            report_this_filefun();
            unimplemented_error()
            obj.parseParameters(varargin);
            obj.StartProcess();
        end

        function InteractiveSetup(obj)
            % create a dialog that allows user to select parameters neccessary for the calculation

            %% make the interface
            zdlg = ZmapDialog();

            zdlg.AddHeader('Choose stuff');
            zdlg.AddMcAutoEstimateCheckbox('mc_auto',obj.mc_auto);
            zdlg.AddPopup('aval_method', 'A-value calculation method', obj.aval_methods,1,...
                'methods used to calculate a value');
            [res,okPressed] = zdlg.Create('Name', 'A-Value Parameters [xsec]');
            if ~okPressed
                return
            end
            obj.SetValuesFromDialog(res);
            obj.doIt()
        end

        function SetValuesFromDialog(obj, res)
            % called when the dialog's OK button is pressed
            obj.mc_auto = res.mc_auto;

            obj.aval_method = res.aval_method;

            tgl1 = tgl1.Value;
            tgl2 = tgl2.Value;
            bRandom = get(chkRandom, 'Value');
            bGridEntireArea = get(chkGridEntireArea, 'Value');
            close,sel ='ca',
            calc_across_orig(sel)

        end

        function results=Calculate(obj)
            % once the properties have been set, either by the constructor or by interactive_setup
            % get the grid-size interactively and calculate the values in the grid by sorting the
            % seismicity and selecting the appropriate neighbors to each grid point

            function out=calculation_function(catalog)
                % calulate values at a single point
            end
        end

        function ModifyGlobals(obj)
        end
    end

    methods(Hidden)
    end
    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin)
            % create a menu item
            label = 'a-value [xsec]';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XZfun.calc_across(zapFcn()),...
                varargin{:});
        end
    end
end

function calc_across_orig(sel)
    % The aValue, Mc and a resolution estimation in each
    % volume around a grid point, or defined by a radius
    % (ra for Mc, ri for aValue) containing ni earthquakes
    % will be calculated
    %
    %   This subroutine provides 4 methods for calculation:
    % 1. calculate a-value for const b-value
    % 2. calculate a-value by maxlikelihood (MaxLikelihoodA.m)
    %    of b-value and Mc defined by MaxC
    % 3. calculate a-value by maxlikelihood (MaxLikelihoodA.m)
    %    of b-value and Mc defined by Mc(EMR)
    % 4. calculate a-value within the radius ri and Mc within ra, where ra > ri
    %
    %   This subrouting is based on bcross.m
    %       by Stefan Wiemer 1/95
    %   and was modified
    %       by Thomas van Stiphout 3/2004

    report_this_filefun();

    % Do we have to create the dialogbox?
    if sel == 'in'
        % Set the grid parameter
        % initial values
        %
        ra=10;
        ri=5;
        dd = 1.00;
        dx = 1.00 ;
        ni = 100;
        fFixbValue=0.9
        bv2 = NaN;
        Nmin = 50;
        stan2 = NaN;
        stan = NaN;
        prf = NaN;
        av = NaN;
        nRandomRuns = 1000;
        bGridEntireArea = false;

        % Create the dialog box
        figure_w_normalized_uicontrolunits(...
            'Name','Grid Input Parameter',...
            'NumberTitle','off', ...
            'units','points',...
            'Visible','on', ...
            'Position',[ ZG.welcome_pos + [200, -200], 550, 300], ...
            'Color', [0.8 0.8 0.8]);
        axis off

        % Dropdown list
        labelList2=[' a(M=0) for fixed b | a(M=0) for Mc(MaxC) + M(corr) | a(M=0) for M(EMR) | a(M=0) by r1 & Mc by r2 '];
        hndl2=uicontrol(...
            'Style','popup',...
            'Position',[ 0.2 0.77  0.6  0.08],...
            'Units','normalized',...
            'String',labelList2,...
            'BackgroundColor','w',...
            'callback',@callbackfun_001);

        % Set selection to 'Radius check'
        set(hndl2,'value',1);

        % Edit fields, radiobuttons, and checkbox
        freq_field=uicontrol('Style','edit',...
            'Position',[.30 .60 .12 .06],...
            'Units','normalized','String',num2str(ni),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_002);

        freq_field0=uicontrol('Style','edit',...
            'Position',[.30 .50 .06 .06],...
            'Units','normalized','String',num2str(ra),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_003);

        freq_field1=uicontrol('Style','edit',...
            'Position',[.36 .50 .06 .06],...
            'Units','normalized','String',num2str(ri),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_004);

        freq_field2=uicontrol('Style','edit',...
            'Position',[.30 .40 .06 .06],...
            'Units','normalized','String',num2str(dx),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_005);

        freq_field3=uicontrol('Style','edit',...
            'Position',[.36 .40 .06 .06],...
            'Units','normalized','String',num2str(dd),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_006);

        tgl1 = uicontrol('BackGroundColor', [0.8 0.8 0.8], ...
            'Style','radiobutton',...
            'string','Number of events:',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'Position',[.02 .60 .28 .06], 'callback',@callbackfun_007,...
            'Units','normalized');

        % Set to constant number of events
        set(tgl1,'value',1);

        tgl2 =  uicontrol('BackGroundColor',[0.8 0.8 0.8],'Style','radiobutton',...
            'string','Constant radius [km]:',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'Position',[.52 .60 .40 .06], 'callback',@callbackfun_008,...
            'Units','normalized');

        chkRandom = uicontrol('BackGroundColor',[0.8 0.8 0.8],'Style','checkbox',...
            'String', 'Additional random simulation',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'Position',[.52 .40 .40 .06],...
            'Units','normalized');
        txtRandomRuns = uicontrol('Style','edit',...
            'Position',[.80 .30 .12 .06],...
            'Units','normalized','String',num2str(nRandomRuns),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_009);

        freq_field4 =  uicontrol('Style','edit',...
            'Position',[.30 .30 .12 .06],...
            'Units','normalized','String',num2str(Nmin),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_010);

        freq_field5 =  uicontrol('Style','edit',...
            'Position',[.30 .20 .12 .06],...
            'Units','normalized','String',num2str(fFixbValue),...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'callback',@callbackfun_011);

        chkGridEntireArea = uicontrol('BackGroundColor', [0.8 0.8 0.8], ...
            'Style','checkbox',...
            'string','Create grid over entire area',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'Position',[.52 .50 .40 .06], 'Units','normalized', 'Value', 0);

        % Buttons
        uicontrol('BackGroundColor', [0.8 0.8 0.8], 'Style', 'pushbutton', ...
            'Units', 'normalized', 'Position', [.80 .05 .15 .10], ...
            'Callback', @(~,~)close, 'String', 'Cancel');

        % Labels
        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [0.2 1 0], 'HorizontalAlignment', 'left',  ...
            'FontSize', fontsz.l, 'FontWeight', 'bold', 'String', 'Please select a Mc estimation option');

        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [0.3 0.75 0], 'HorizontalAlignment', 'left',  ...
            'FontSize', fontsz.l, 'FontWeight', 'bold', 'String', 'Grid parameters');

        text('Color',[0 0 0], 'Units', 'normalized', ...
            'Position', [-.14 .51 0], 'HorizontalAlignment', 'left',  ...
            'FontSize',ZmapGlobal.Data.fontsz.m, 'FontWeight', 'bold', 'String','Radius ra / ri [km]:');

        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [-0.14 .39 0],  'HorizontalAlignment', 'left', ...
            'FontSize',ZmapGlobal.Data.fontsz.m, 'FontWeight', 'bold', 'String', 'Spacing hor / depth [km]:');

        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [-0.14 .27 0],  'HorizontalAlignment', 'left', ...
            'FontSize',ZmapGlobal.Data.fontsz.m, 'FontWeight', 'bold', 'String', 'Min. number of events:');

        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [-0.14 .15 0],  'HorizontalAlignment', 'left', ...
            'FontSize',ZmapGlobal.Data.fontsz.m, 'FontWeight', 'bold', 'String', 'Fixed b-value:');

        text('Color', [0 0 0], 'Units', 'normalized', ...
            'Position', [0.5 0.27 0],  'HorizontalAlignment', 'left', ...
            'FontSize',ZmapGlobal.Data.fontsz.m, 'FontWeight', 'bold', 'String', 'Number of runs:');

        set(gcf,'visible','on');
        watchoff
    end

    % get the grid-size interactively and
    % calculate the b-value in the grid by sorting
    % the seismicity and selecting the ni neighbors
    % to each grid point

    if sel == 'ca'
        figure(xsec_fig());
        set(gca,'NextPlot','add')

        if bGridEntireArea % Use entire area for grid
            vXLim = get(gca, 'XLim');
            vYLim = get(gca, 'YLim');
            x = [vXLim(1); vXLim(1); vXLim(2); vXLim(2)];
            y = [vYLim(2); vYLim(1); vYLim(1); vYLim(2)];
            x = [x ; x(1)];
            y = [y ; y(1)];     %  closes polygon
            clear vXLim vYLim;
        else

            set(gca,'NextPlot','add')
            ax=findobj(gcf,'Tag','mainmap_ax');
            [x,y, mouse_points_overlay] = select_polygon(ax);
        end % of if bGridEntireArea

        plos2 = plot(x,y,'b-');        % plot outline
        sum3 = 0.;
        pause(0.3)

        %create a rectangular grid
        xvect=[min(x):dx:max(x)];
        yvect=[min(y):dd:max(y)];
        gx = xvect;gy = yvect;
        tmpgri=zeros((length(xvect)*length(yvect)),2);
        n=0;
        for i=1:length(xvect)
            for j=1:length(yvect)
                n=n+1;
                tmpgri(n,:)=[xvect(i) yvect(j)];
            end
        end
        %extract all gridpoints in chosen polygon
        XI=tmpgri(:,1);
        YI=tmpgri(:,2);

        ll = polygon_filter(x,y, XI, YI, 'inside');
        %grid points in polygon
        newgri=tmpgri(ll,:);

        % Plot all grid points
        plot(newgri(:,1),newgri(:,2),'+k')

        if length(xvect) < 2  ||  length(yvect) < 2
            errordlg('Selection too small! (not a matrix)');
            return
        end

        itotal = length(newgri(:,1));
        if length(gx) < 4  ||  length(gy) < 4
            errordlg('Selection too small! ');
            return
        end

        %  make grid, calculate start- endtime etc.  ...
        %
        [t0b, teb] = bounds(newa.Date) ;
        n = newa.Count;
        tdiff = round((teb-t0b)/ZG.bin_dur);
        loc = zeros(3, length(gx)*length(gy));

        % loop over  all points
        %
        i2 = 0.;
        i1 = 0.;
        avg = nan(length(newgri),10);
        allcount = 0.;
        wai = waitbar(0,' Please Wait ...  ');
        set(wai,'NumberTitle','off','Name','a-value grid - percent done');
        drawnow
        %
        % loop

        % overall b-value
        [bv, magco, stan, av] =  bvalca3(newa.Magnitude, obj.mc_auto);
        overall_b_value = bv;
        ZG.overall_b_value = bv;
        no1 = newa.Count;
        %
        globalcatalog = ZG.primeCatalog;
        for i= 1:length(newgri(:,1))
            x = newgri(i,1);y = newgri(i,2);
            allcount = allcount + 1.;
            i2 = i2+1;

            % calculate distance from center point and sort wrt distance
            l = sqrt(((xsecx' - x)).^2 + ((xsecy + y)).^2) ;
            [s,is] = sort(l);
            b = newa(is(:,1),:) ;       % re-orders matrix to agree row-wise

            if tgl1 == 0   % take point within r
                l3 = l <= ra;
                l4 = l <= ri;
                b = globalcatalog.subset(l3);        % new data per grid point (b) is sorted in distance
                bri = globalcatalog.subset(l4);
                rd = ra;
            else
                % take first ni points
                b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
                rd = s(ni);

            end

            % Number of earthquakes per node
            [nX,nY] = size(b);

            %estimate the completeness and b-value
            ZG.newt2 = b;

            if length(b) >= Nmin  % enough events?

                if obj.aval_method == 1;   % Calculation of a-value by const b-value, and Mc
                    bv2 = fFixbValue;           % read fixed bValue to the bv2
                    magco=calc_Mc(b, Methods.MaxCurvature, 0.1);
                    l = b.Magnitude >= magco-0.05;
                    if length(b(l,:)) >= Nmin   % calculation of the a-value according to determined Mc (magco)
                        faValue = calc_MaxLikelihoodA(b, bv2);
                        stan2 = NaN;
                        bv = NaN;
                    else
                        bv = NaN; bv2 = NaN; magco = NaN; av = NaN; faValue = NaN; stan2 = NaN; stan = NaN;
                    end

                    % a(0) for Mc(MAxCurv) + Mc(Corr)
                elseif obj.aval_method == 2
                    [bv, magco, stan, av, pr] =  bvalca3(b.Magnitude, McAutoEstimate.auto, overall_b_value);
                    magco = magco + 0.2;    % Add 0.2 to Mc (Tobias)
                    l = b.Magnitude >= magco-0.05;
                    if sum(l) >= Nmin
                        [bv2, stan2,  faValue] = calc_bmemag(b.Magnitude(l), 0.1);
                    else
                        bv = NaN; bv2 = NaN, magco = NaN; av = NaN; faValue = NaN;
                    end

                elseif obj.aval_method == 3; % a(0) for Mc(EMR)
                    [magco, bv2, faValue, stan2, stan] = calc_McEMR(b, 0.1);
                    l = b.Magnitude >= magco-0.05;
                    if length(b(l,:)) >= Nmin
                        faValue = calc_MaxLikelihoodA(b, bv2);
                    else
                        bv = NaN; bv2 = NaN, magco = NaN; av = NaN; faValue = NaN;
                    end

                elseif obj.aval_method == 4
                    % a(0) by r1 and Mc by r2
                    if length(b) >= Nmin
                        [bv, magco, stan, av, pr] =  bvalca3(b.Magnitude, McAutoEstimate.auto, overall_b_value);
                        magco = magco + 0.2;    % Add 0.2 to Mc (Tobias)
                        bv2 = fFixbValue;
                        l = bri(:,6) >= magco-0.05;
                        faValue = calc_MaxLikelihoodA(bri(l,:), bv2);
                    else
                        bv = NaN; bv2 = NaN, magco = NaN; av = NaN; faValue = NaN;
                    end
                end

                % Perform random simulation
                if bRandom
                    nNumberPerNode = length(b);
                    %        [fAverageBValue, fAverageStdDev] = calc_RandomBValue(newa, nNumberPerNode, nRandomRuns);
                    [fStdDevB, fStdDevMc, q1, q2] = calc_BootstrapB(b, nRandomRuns, Nmin);
                else
                    fStdDevB = NaN;
                    fStdDevMc = NaN;
                end

            else % of if length(b) >= Nmin
                bv = NaN; bv2 = NaN; stan = NaN; stan2 = NaN; prf = NaN; magco = NaN; av = NaN; faValue = NaN; fStdDevB = NaN; fStdDevMc = NaN;
                b = [NaN NaN NaN NaN NaN NaN NaN NaN NaN];
                nX = NaN;
            end
            mab = max(b.Magnitude) ;
            if isempty(mab); mab = NaN; end
            avg(allcount,:)  = [bv magco x y rd bv2 stan2 av stan faValue ];
            waitbar(allcount/itotal)
        end  % for  newgri

        if bRandom
            clear nNumberPerNode q1 q2 fStdDevB fStdDevMc;
        end
        clear bRandom;
        % save data
        %
        drawnow
        gx = xvect;gy = yvect;

        catsave3('calc_across_orig');
        close(wai)
        watchoff

        save across2.mat ll a tmpgri newgri lat1 lon1 lat2 lon2 ZG.xsec_defaults.WidthKm  avg ra xvect yvect gx gy dx dd ZG.bin_dur newa ;

        % reshape a few matrices
        %

        normlap2=nan(length(tmpgri(:,1)),1)
        normlap2(ll)= avg(:,1);
        bls=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,5);
        reso=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,6);
        bValueMap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,2);
        MaxCMap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,7);
        MuMap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,8);
        avm=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,9);
        SigmaMap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= avg(:,10);
        aValueMap=reshape(normlap2,length(yvect),length(xvect));

        valueMap = aValueMap;
        kll = ll;
        % View the a-value map
        view_av2

    end

    % Load exist b-grid
    if sel == 'lo'
        [file1,path1] = uigetfile(['*.mat'],'a-value gridfile');
        if length(path1) > 1

            load([path1 file1])
            xsecx = newa(:,end)';
            xsecy = newa(:,7);
            xvect = gx; yvect = gy;
            tmpgri=zeros((length(xvect)*length(yvect)),2);

            normlap2=nan(length(tmpgri(:,1)),1)
            normlap2(ll)= avg(:,1);
            bls=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,5);
            reso=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,6);
            bValueMap=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,2);
            MaxCMap=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,7);
            MuMap=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,8);
            avm=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,9);
            SigmaMap=reshape(normlap2,length(yvect),length(xvect));

            normlap2(ll)= avg(:,10);
            aValueMap=reshape(normlap2,length(yvect),length(xvect));

            valueMap = aValueMap;

            nlammap
            [xsecx xsecy,  inde] =mysect(globalcatalog.Y',globalcatalog.X',globalcatalog.Z,ZG.xsec_defaults.WidthKm,0,lat1,lon1,lat2,lon2);
            % Plot all grid points
            set(gca,'NextPlot','add')
            plot(newgri(:,1),newgri(:,2),'+k')
            view_av2
        else
            return
        end
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2double(freq_field.String);
        freq_field.String=num2str(ni);
        tgl2.Value=0;
        tgl1.Value=1;
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ra=str2double(freq_field0.String);
        freq_field0.String=num2str(ra);
        tgl2.Value=1;
        tgl1.Value=0;
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ri=str2double(freq_field1.String);
        freq_field1.String=num2str(ri);
        tgl2.Value=1;
        tgl1.Value=0;
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dx=str2double(freq_field2.String);
        freq_field2.String=num2str(dx);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dd=str2double(freq_field3.String);
        freq_field3.String=num2str(dd);
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        tgl2.Value=0;
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        tgl1.Value=0;
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        nRandomRuns=str2double(txtRandomRuns.String);
        txtRandomRuns.String=num2str(nRandomRuns);
    end

    function callbackfun_010(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        Nmin=str2double(freq_field4.String);
        freq_field4.String=num2str(Nmin);
    end

    function callbackfun_011(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        fFixbValue=str2double(freq_field4.String);
        freq_field4.String=num2str(fFixbValue);
    end
end

author: CelsoReyes


text:
function euler_m = dcm2euler_m(DCMnb_m)
% dcm2euler_m: transforms a nav-to-body DCM matrix to an Euler angles matrix.
%
% INPUT
%   DCMnb_m: Nx9 matrix with nav-to-body direct cosine matrices (DCM).
% Each row of DCMnb_m contains the 9 elements of a particular DCMnb
% matrix ordered as [a11 a21 a31 a12 a22 a32 a13 a23 a33].
%
% OUTPUT
%   euler_m: Nx3 Euler angles ordered by column, [roll pitch yaw] (rad, rad, rad).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
% 	Titterton, D.H. and Weston, J.L. (2004). Strapdown
% Inertial Navigation Technology (2nd Ed.). Institution
% of Engineering and Technology, USA. Eq. 11.4. Eq. 3.66, p. 46.
%
%   R. Gonzalez, J. Giribet, and H.D. Patiño,. An approach to
% benchmarking of loosely coupled low-cost navigation systems,
% Mathematical and Computer Modelling of Dynamical Systems, vol. 21,
% issue 2, pp. 272-287. Eq. 15.
%
% Version: 001
% Date:    2020/12/17
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

[N,~] = size (DCMnb_m);

euler_m = zeros(N,3);

for i=1:N

    DCMnb = reshape(DCMnb_m (i,:), 3, 3);
    DCMbn = DCMnb';                         % nav-to-body > body-to-nav

    euler_m (i,:) = dcm2euler(DCMbn);       % phi theta psi
end

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [scoremaps, locreg_pred, nextreg_pred] = cnn_process_image_tiled(input, net, sigmoid, im_bg_width, im_bg_height, stride)

    max_size = 1000;
    rf = 224; %receptive field
    cut_off = rf/stride;

    num_tiles_x = get_num_tiles(im_bg_width, max_size);
    num_tiles_y = get_num_tiles(im_bg_height, max_size);

    scoremaps = [];
    locreg_pred = [];
    nextreg_pred = [];

    for j = 1:num_tiles_y
        if j == 1
            start_y = 1;
        else
            start_y = (j-1)*(max_size-2*rf)+1;
        end
        if j == num_tiles_y
            end_y = im_bg_height;
        else
            end_y = start_y+max_size-1;
        end

        scoremaps_line = [];
        locreg_pred_line = [];
        nextreg_pred_line = [];

        for i = 1:num_tiles_x
            if i == 1
                start_x = 1;
            else
                start_x = (i-1)*(max_size-2*rf)+1;
            end
            if i == num_tiles_x
                end_x = im_bg_width;
            else
                end_x = start_x+max_size-1;
            end

            input_tile = input(start_y:end_y, start_x:end_x, :);

            [scoremaps_tile, locreg_pred_tile, nextreg_pred_tile] = cnn_process_image(input_tile, net, sigmoid);

            scoremaps_tile = cutoff_tile(scoremaps_tile, num_tiles_x, i, cut_off, true);
            locreg_pred_tile = cutoff_tile(locreg_pred_tile, num_tiles_x, i, cut_off, true);
            nextreg_pred_tile = cutoff_tile(nextreg_pred_tile, num_tiles_x, i, cut_off, true);

            scoremaps_line = cat(2, scoremaps_line, scoremaps_tile);
            locreg_pred_line = cat(2, locreg_pred_line, locreg_pred_tile);
            nextreg_pred_line = cat(2, nextreg_pred_line, nextreg_pred_tile);
        end

        scoremaps_line = cutoff_tile(scoremaps_line, num_tiles_y, j, cut_off, false);
        locreg_pred_line = cutoff_tile(locreg_pred_line, num_tiles_y, j, cut_off, false);
        nextreg_pred_line = cutoff_tile(nextreg_pred_line, num_tiles_y, j, cut_off, false);

        scoremaps = cat(1, scoremaps, scoremaps_line);
        locreg_pred = cat(1, locreg_pred, locreg_pred_line);
        nextreg_pred = cat(1, nextreg_pred, nextreg_pred_line);
    end

end

function [feat_prob, locreg_pred, next_pred] = cnn_process_image_tile(input, net, sigmoid)

    % switch width and height for Caffe
    input = permute(input, [2 1 3]);

    blob_size = [size(input), 1];
    net.blobs('data').reshape(blob_size);

    blob_size = [size(input), 1];

    batch = zeros(blob_size, 'single');
    batch(:,:,:,1) = input;
    batches = cell(1,1);
    batches{1} = batch;
    zs = net.forward(batches);

    % fc6_data = net.blobs('fc6-conv').get_data();

    locreg_pred = [];
    next_pred = [];

    for k = 1:length(zs);
        fmap = zs{k};
        if size(fmap, 3) == 14
            feat_prob = fmap;
            feat_prob = permute(feat_prob, [2 1 3]);
            continue;
        end
        fmap = zs{k};
        sz = size(fmap);
        fmap = reshape(fmap, sz(1), sz(2), 2, []);
        fmap = permute(fmap, [2 1 4 3]);
        if size(fmap, 3) == 14
            locreg_pred = fmap;
        else
            next_pred = fmap;
        end
    end

    if ~sigmoid
        % remove background class
        feat_prob = feat_prob(:,:,2:end);
    end
end

function num = get_num_tiles(sz, max_size)

rf = 224; %receptive field

if sz <= max_size
    num = 1;
    return;
end

k = 0;
while true
    % first and last tiles + middle tiles
    new_size = (max_size-rf)*2 + (max_size-2*rf)*k;
    if new_size >= sz
        break;
    end
    k = k + 1;
end

num = 2 + k;
end

function sm = cutoff_tile(sm, num_tiles, idx, cut_off, is_x)
    new_sz = [2 1 3 4];
    if is_x
        sm = permute(sm, new_sz);
    end

    if num_tiles == 1
    elseif idx == 1
        sm = sm(1:end-cut_off-1, :, :, :);
    elseif idx == num_tiles
        sm = sm(cut_off+1:end, :, :, :);
    else
        sm = sm(cut_off+1:end-cut_off-1, :, :, :);
    end

    if is_x
        sm = permute(sm, new_sz);
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
%
%   Generates html versions of all MATLAB files
%
opts.outputDir = '.';
opts.evalCode  = false;

d = dir('*.m');
l = size(d);
l = l(1);
for i=1:l
  n = d(i).name;
  publish(n,opts);
  n = n(1:(length(n)-2));
  nhtml = [n '.html'];
  nmhtml = [n '.m.html'];
  movefile(nhtml,nmhtml);
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [matched,only1,only2] = Match(list1,list2,varargin)

%Match - Replace values in one list with closest values in a second list.
%
%  USAGE
%
%    [matched,only1,only2] = Match(list1,list2,<options>)
%
%    list1          first list
%    list2          second list
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'match'       optional matching criterion: 'up', 'down' (default) or
%                   'closest' (see below)
%     'error'       matched values for which the difference exceeds this
%                   threshold are set to NaN (default = inf)
%    =========================================================================
%
%  OUTPUT
%
%    matched        values in list2 that matched those of list1
%    only1          elements in list1 that were not matched by any element in list2
%    only2          elements in list2 that did not match any element in list1
%
%    only1 and only2 are lists of logical indices.
%
%  NOTE
%
%    Pairs are matched in the following way. Assume the lists are sorted in
%    ascending order and let list1 = [a(1)...a(n)] and list2 = [b(1)...b(m)].
%    We are trying to match each a(i) with one of the b(1..m). Suppose
%
%       b(1) <= ... <= b(j) <= a(i) <= b(j+1) <= ... <= b(m)
%
%    Then a(i) will be matched with either b(j) or b(j+1), depending on the
%    'match' option:
%
%       1) 'down' (default) => b(j)
%       2) 'up' => b(j+1)
%       3) 'closest' => the one that minimizes the distance
%          (b(j) if the distances are equal)
%
%    If however the distance between a(i) and b(j) (or b(j+1)) is too large
%    (see option 'error'), then none of the b(1..m) matches a(i).
%
%    Note that this function is generally asymmetrical.

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

% Default values
match = 'down';
err = inf;

% Check parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help Match'' for details).');
end
if ~isdvector(list1) | ~isdvector(list2),
	error('Incorrect sizes: ''list1'' and ''list2'' must be vectors (type ''help Match'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~isa(varargin{i},'char'),
    error(['Parameter ' num2str(i+2) ' is not a property (type ''help Match'' for details).']);
  end
  switch(lower(varargin{i})),

    case 'match',
		match = lower(varargin{i+1});
		if ~isstring_FMAT(match,'up','down','closest'),
			error('Incorrect value for property ''match'' (type ''help Match'' for details).');
      end

	case 'error',
		err = varargin{i+1};
		if ~isdscalar(err,'>=0'),
			error('Incorrect value for property ''error'' (type ''help Match'' for details).');
		end

    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help Match'' for details).']);

  end
end

% Sort elements in both list so that list1 = [a1 a2 a3 ... an] and list2 = [b1 b2 b3 ... bm]
[list1,idx1] = sort(list1(:));
[list2,idx2] = sort(list2(:));

% Find 'up' and 'down' matching indices
% (We need both, because this will allow us to compare matches if the criterion is 'closest')
n = length(list2);
up = MatchUpIndices(list1,list2);
down = up-1;

% List corresponding values in list2
% Special cases: when matching down (resp. up), values in list1 lesser (resp. greater)
% than all values in list2 cannot be matched; set them to NaN
goodUp = up>0 & up<=n;
matchedUp = nan(size(up));
matchedUp(goodUp) = list2(up(goodUp));
goodDown = down>0 & down<=n;
matchedDown = nan(size(down));
matchedDown(goodDown) = list2(down(goodDown));

% Use match criterion
switch(match),
	case 'up',
		matched = matchedUp;
	case 'down',
		matched = matchedDown;
	case 'closest',
		[unused,idx] = min(abs([matchedUp matchedDown]-[list1 list1]),[],2);
		matched = matchedUp;
		matched(idx==2) = matchedDown(idx==2);
end

% Check error threshold
bad = abs(matched-list1) > err+eps; % eps is to compensate for inevitable numerical rounding errors such as 1.1-1.0>0.1
matched(bad) = nan;

% Output
only1 = isnan(matched);
matched = matched(~only1);
only1(idx1) = only1;
only2 = logical(zeros(size(list2)));
m = unique(matched);
[unused,i] = setdiff(list2,m);
only2(i) = 1;
only2(idx2) = only2;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef DagNN < matlab.mixin.Copyable
%DagNN Directed acyclic graph neural network
%   DagNN is a CNN wrapper alternative to SimpleNN. It is object
%   oriented and allows constructing networks with a directed acyclic
%   graph (DAG) topology. It is therefore far more flexible, although
%   a little more complex and slightly slower for small CNNs.
%
%   A DAG object contains the following data members:
%
%   - `layers`: The network layers.
%   - `vars`: The network variables.
%   - `params`: The network parameters.
%   - `meta`: Additional information relative to the CNN (e.g. input
%     image format specification).
%
%   There are additional transient data members:
%
%   `mode`:: `normal`
%      This flag can either be `normal` or `test`. In the latter case,
%      certain blocks switch to a test mode suitable for validation or
%      evaluation as opposed to training. For instance, dropout
%      becomes a pass-through block in `test` mode.
%
%   `accumulateParamDers`:: `false`
%      If this flag is set to `true`, then the derivatives of the
%      network parameters are accumulated rather than rewritten the
%      next time the derivatives are computed.
%
%   `conserveMemory`:: `true`
%      If this flag is set to `true`, the DagNN will discard
%      intermediate variable values as soon as they are not needed
%      anymore in the calculations. This is particularly important to
%      save memory on GPUs.
%
%   `device`:: `cpu`
%      This flag tells whether the DagNN resides in CPU or GPU
%      memory. Use the `DagNN.move()` function to move the DagNN
%      between devices.
%
%   The DagNN is copyable handle, i.e. allows to create a deep copy using
%   `copy` operator `deep_copy = copy(dagnet);`. In all cases the deep copy
%   is located in CPU memory (i.e. is transfered from GPU before copy).
%   Remark: As a side effect the original network is being reset (all
%   variables are cleared) and only the network structure and parameters
%   are copied.
%
%   See Also: matlab.mixin.Copyable

% Copyright (C) 2015-2016 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  properties
    layers
    vars
    params
    meta
  end

  properties (Transient)
    mode = 'normal'
    accumulateParamDers = false
    conserveMemory = true
  end

  properties (Transient, SetAccess = private, GetAccess = public)
    device = 'cpu' ;
  end

  properties (Transient, Access = {?dagnn.DagNN, ?dagnn.Layer}, Hidden = true)
    numPendingVarRefs
    numPendingParamRefs
    computingDerivative = false
    executionOrder
  end

  properties (Transient, Access = private, Hidden = true)
    modifed = false
    varNames = struct()
    paramNames = struct()
    layerNames = struct()
    layerIndexes = {}
  end

  methods
    function obj = DagNN()
    %DAGNN  Initialize an empty DaG
    %   OBJ = DAGNN() initializes an empty DaG.
    %
    %   See Also addLayer(), loadobj(), saveobj().
      obj.vars = struct(...
        'name', {}, ...
        'value', {}, ...
        'der', {}, ...
        'fanin', {}, ...
        'fanout', {}, ...
        'precious', {}) ;
      obj.params = struct(...
        'name', {}, ...
        'value', {}, ...
        'der', {}, ...
        'fanout', {}, ...
        'trainMethod', {}, ...
        'learningRate', {}, ...
        'weightDecay', {}) ;
      obj.layers = struct(...
        'name', {}, ...
        'inputs', {}, ...
        'outputs', {}, ...
        'params', {}, ...
        'inputIndexes', {}, ...
        'outputIndexes', {}, ...
        'paramIndexes', {}, ...
        'forwardTime', {[]}, ...
        'backwardTime', {[]}, ...
        'block', {}) ;
    end

    function set.mode(obj, mode)
      switch lower(mode)
        case {'normal', 'train'}
          obj.mode = 'normal' ;
        case {'test'}
          obj.mode = 'test' ;
      end
    end

    % Manage the DagNN
    reset(obj)
    move(obj, direction)
    s = saveobj(obj)

    % Manipualte the DagNN
    addLayer(obj, name, block, inputs, outputs, params, varargin)
    removeLayer(obj, name)
    setLayerInputs(obj, leyer, inputs)
    setLayerOutput(obj, layer, outputs)
    setLayerParams(obj, layer, params)
    renameVar(obj, oldName, newName, varargin)
    rebuild(obj)

    % Process data with the DagNN
    initParams(obj)
    eval(obj, inputs, derOutputs)

    % Get information about the DagNN
    varSizes = getVarSizes(obj, inputSizes)

    % ---------------------------------------------------------------------
    %                                                           Access data
    % ---------------------------------------------------------------------

    function inputs = getInputs(obj)
    %GETINPUTS Get the names of the input variables
    %   INPUTS = GETINPUTS(obj) returns a cell array containing the name
    %   of the input variables of the DaG obj, i.e. the sources of the
    %   DaG (excluding the network parameters, which can also be
    %   considered sources).
      v = find([obj.vars.fanin]==0) ;
      inputs = {obj.vars(v).name} ;
    end

    function outputs = getOutputs(obj)
    %GETOUTPUTS Get the names of the output variables
    %    OUTPUT = GETOUTPUTS(obj) returns a cell array containing the name
    %    of the output variables of the DaG obj, i.e. the sinks of the
    %    DaG.
      v = find([obj.vars.fanout]==0) ;
      outputs = {obj.vars(v).name} ;
    end

    function l = getLayerIndex(obj, name)
    %GETLAYERINDEX Get the index of a layer
    %   INDEX = GETLAYERINDEX(obj, NAME) returns the index of the layer
    %   NAME. NAME can also be a cell array of strings. If no layer
    %   with such a name is found, the value NaN is returned for the
    %   index.
    %
    %   Layers can then be accessed as the `obj.layers(INDEX)`
    %   property of the DaG.
    %
    %   Indexes are stable unless the DaG is modified (e.g. by adding
    %   or removing layers); hence they can be cached for faster
    %   variable access.
    %
    %   See Also getParamIndex(), getVarIndex().
      if iscell(name)
        l = zeros(1, numel(name)) ;
        for k = 1:numel(name)
          l(k) = obj.getLayerIndex(name{k}) ;
        end
      else
        if isfield(obj.layerNames, name)
          l = obj.layerNames.(name) ;
        else
          l = NaN ;
        end
      end
    end

    function v = getVarIndex(obj, name)
    %GETVARINDEX Get the index of a variable
    %   INDEX = GETVARINDEX(obj, NAME) obtains the index of the variable
    %   with the specified NAME. NAME can also be a cell array of
    %   strings. If no variable with such a name is found, the value
    %   NaN is returned for the index.
    %
    %   Variables can then be accessed as the `obj.vars(INDEX)`
    %   property of the DaG.
    %
    %   Indexes are stable unless the DaG is modified (e.g. by adding
    %   or removing layers); hence they can be cached for faster
    %   variable access.
    %
    %   See Also getParamIndex(), getLayerIndex().
      if iscell(name)
        v = zeros(1, numel(name)) ;
        for k = 1:numel(name)
          v(k) = obj.getVarIndex(name{k}) ;
        end
      else
        if isfield(obj.varNames, name)
          v = obj.varNames.(name) ;
        else
          v = NaN ;
        end
      end
    end

    function p = getParamIndex(obj, name)
    %GETPARAMINDEX Get the index of a parameter
    %   INDEX = GETPARAMINDEX(obj, NAME) obtains the index of the
    %   parameter with the specified NAME. NAME can also be a cell
    %   array of strings. If no parameter with such a name is found,
    %   the value NaN is returned for the index.
    %
    %   Parameters can then be accessed as the `obj.params(INDEX)`
    %   property of the DaG.
    %
    %   Indexes are stable unless the DaG is modified (e.g. by adding
    %   or removing layers); hence they can be cached for faster
    %   parameter access.
    %
    %   See Also getVarIndex(), getLayerIndex().
      if iscell(name)
        p = zeros(1, numel(name)) ;
        for k = 1:numel(name)
          p(k) = obj.getParamIndex(name{k}) ;
        end
      else
        if isfield(obj.paramNames, name)
          p = obj.paramNames.(name) ;
        else
          p = NaN ;
        end
      end
    end

    function layer = getLayer(obj, layerName)
    %GETLAYER Get a copy of a layer definition
    %   LAYER = GETLAYER(obj, NAME) returns a copy of the layer definition
    %   structure with the specified NAME. NAME can also be a cell array
    %   of strings or an array of indexes. If no parameter with a
    %   specified name or index exists, an error is thrown.
    %
    %   See Also getLayerIndex().
      if isnumeric(layerName)
        idxs = layerName;
        if any(idxs > numel(obj.layers) || idxs < 0)
          error('Invalid layer indexes.');
        end
      else
        if ischar(layerName), layerName = {layerName}; end;
        idxs = obj.getLayerIndex(layerName);
        if any(isnan(idxs))
          error('Invalid layer name `%s`', ...
            strjoin(layerName(isnan(idxs)), ', '));
        end
      end
      layer = obj.layers(idxs);
    end

    function var = getVar(obj, varName)
    %GETVAR Get a copy of a layer definition
    %   VAR = GETVAR(obj, NAME) returns a copy of the network variable
    %   with the specified NAME. NAME can also be a cell array
    %   of strings or an array of indexes. If no variable with a
    %   specified name or index exists, an error is thrown.
    %
    %   See Also getVarIndex().
      if isnumeric(varName)
        idxs = varName;
        if any(idxs > numel(obj.vars) || idxs < 0)
          error('Invalid var indexes.');
        end
      else
        if ischar(varName), varName = {varName}; end;
        idxs = obj.getVarIndex(varName);
        if any(isnan(idxs))
          error('Invalid variable name `%s`', ...
            strjoin(varName(isnan(idxs)), ', '));
        end
      end
      var = obj.vars(idxs);
    end

    function param = getParam(obj, paramName)
    %GETPARAM Get a copy of a layer parameter
    %   PARAM = GETPARAM(obj, NAME) returns a copy of the network parameter
    %   with the specified NAME. NAME can also be a cell array
    %   of strings or an array of indexes. If no parameter with a
    %   specified name or index exists, an error is thrown.
    %
    %   See Also getParamIndex().
      if isnumeric(paramName)
        idxs = paramName;
        if any(idxs > numel(obj.params) || idxs < 0)
          error('Invalid param indexes.');
        end
      else
        if ischar(paramName), paramName = {paramName}; end;
        idxs = obj.getParamIndex(paramName);
        if any(isnan(idxs))
          error('Invalid param name `%s`', ...
            strjoin(paramName(isnan(idxs)), ', '));
        end
      end
      param = obj.params(idxs);
    end

    function order = getLayerExecutionOrder(obj)
    %GETLAYEREXECUTIONORDER Get the order in which layers are evaluated
    %   ORDER = GETLAYEREXECUTIONORDER(obj) returns a vector with
    %   the indexes of the layers in the order in which they are
    %   executed. This needs not to be the trivial order 1,2,...,L
    %   as it depends on the graph topology.
      order = obj.executionOrder ;
    end
  end

  methods (Static)
    obj = fromSimpleNN(net, varargin)
    obj = loadobj(s)
  end

  methods (Access = {?dagnn.DagNN, ?dagnn.Layer})
    function v = addVar(obj, name)
    %ADDVAR  Add a variable to the DaG
    %   V = ADDVAR(obj, NAME) adds a varialbe with the specified
    %   NAME to the DaG. This is an internal function; variables
    %   are automatically added when adding layers to the network.
      v = obj.getVarIndex(name) ;
      if ~isnan(v), return ; end
      v = numel(obj.vars) + 1 ;
      obj.vars(v) = struct(...
        'name', {name}, ...
        'value', {[]}, ...
        'der', {[]}, ...
        'fanin', {0}, ...
        'fanout', {0}, ...
        'precious', {false}) ;
      obj.varNames.(name) = v ;
    end

    function p = addParam(obj, name)
    %ADDPARAM  Add a parameter to the DaG
    %   V = ADDPARAM(obj, NAME) adds a parameter with the specified NAME
    %   to the DaG. This is an internal function; parameters are
    %   automatically added when adding layers to the network.
      p = obj.getParamIndex(name) ;
      if ~isnan(p), return ; end
      p = numel(obj.params) + 1 ;
      obj.params(p) = struct(...
        'name', {name}, ...
        'value', {[]}, ...
        'der', {[]}, ...
        'fanout', {0}, ...
        'trainMethod', {'gradient'}, ...
        'learningRate', {1}, ...
        'weightDecay', {1}) ;
      obj.paramNames.(name) = p ;
    end
  end

  methods (Access = protected)
    function cp = copyElement(obj)
      % Create a deep copy of the network
      cp = dagnn.DagNN.loadobj(obj.saveobj());
    end
  end
end

author: jiangqy


text:
function x = p17_root ( i )

%*****************************************************************************80
%
%% P17_ROOT returns a root for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the requested root.
%
%    Output, real X, the value of the root.
%
  if ( i == 1 )
    x = 2.0945514815423265;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P17_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal root index = %d\n', i );
    error ( 'P17_ROOT - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Formatting_step3_align

clear all;close all;clc

%% set manually
M_stimset = GetFishStimset();
M_dir = GetFishDirectories();
save_masterdir = GetNestedDataDir();

%%
range_fish = GetFishRange();

for i_fish = range_fish,
    disp(['i_fish = ', num2str(i_fish)]);

    %% load data
    save_dir = fullfile(save_masterdir,['subject_' num2str(i_fish)]);
    load(fullfile(save_dir,'CoreInfo.mat'));
    load(fullfile(save_dir,'OptionalInfo.mat'));
    load(fullfile(save_dir,'AdditionalInfo.mat'));

    %% index processing
    if M_stimset(i_fish)==1, % fish 1-7
        % 'periods'
        M_period = {480,160,160,320,480,140,150,8,9,10,11,12,13,14,15,200}; %,{120,150,360}};% 8-15 are place holders
        periods = M_period{i_fish};

        % 'timelists'
        period = periods;
        nrep = floor(size(frame_keys,1)/period)-1;
        shift = period;
        timelists_raw = 1+shift:period*nrep+shift;

        % 'timelists_name'
        if i_fish<6,
            timelists_names = {'16 permut: B/W/phototaxis*2'};
        else % i_fish = 6 or 7
            timelists_names = {'phototaxis'};
        end

        % stimuluskey_raw
        stimuluskey_raw = round(frame_keys(:,17)');

    else % fish recorded with multiple protocols
        % 'stimset': struct that contains all info parsed from recorded params
        [stimset,stimuluskey_raw] = StimulusKeyPreprocessing(frame_keys,i_fish); % StimulusKeyPreprocessing(frame_keys,i_fish,'isplotting') to show plot

        % 'periods','timelists_names'
        periods = [stimset.period];
        timelists_names = {stimset.name};

        % very first 100 frames of experiment should be discarded (prone to artifacts)
        if stimset(1).starts(1)<100,
            if stimset(1).nReps(1)>1,
                stimset(1).starts(1) = stimset(1).starts(1)+stimset(1).period;
                stimset(1).nReps(1) = stimset(1).nReps(1)-1;
            elseif stimset(1).nReps(1)==1, % (n/a for current fish, improbable for future fish)
                stimset(1).ij(1,:) = [];
                stimset(1).rawstarts(1) = [];
                stimset(1).rawstops(1) = [];
                stimset(1).starts(1) = [];
                stimset(1).stops(1) = [];
                stimset(1).nReps(1) = [];
            end
        end

        % 'timelists' = list of time-frames numbers for a certain stimulus
        nTypes = length(stimset);
        timelists_raw = cell(1,nTypes); % time-lists, i.e. frame indices
        for i_type = 1:nTypes,
            nSets = length(stimset(i_type).starts);
            IX = [];
            for i_set = 1:nSets,
                IX = horzcat(IX,stimset(i_type).starts(i_set):stimset(i_type).stops(i_set));
            end
            timelists_raw{i_type} = IX;
        end
    end

    %% Behavioral data: fictive motor recordings
    rows = [7,8,9,13,14];
    Behavior_raw = frame_keys(:,rows)';
    Behavior_raw(2,:) = -Behavior_raw(2,:);

    %% save to nested directory
    save_dir_nested = fullfile(save_masterdir,['subject_' num2str(i_fish)]);

    filename = fullfile(save_dir_nested,'CoreInfo.mat');
    varList = {'timelists_raw','timelists_names','periods','stimuluskey_raw'}; % {'CellXYZ','anat_stack','fpsec'};
    save(filename,varList{:},'-append');

    filename = fullfile(save_dir_nested,'OptionalInfo.mat');
    varList = {'Behavior_raw'}; % {'numcell_full','CellXYZ_norm'};
    save(filename,varList{:},'-append');

    if exist('stimset','var'),
        filename = fullfile(save_dir_nested,'AdditionalInfo.mat');
        save(filename,'stimset','-append'); % {'frame_keys','IX_inval_anat'};
    end
end

%     % renameing
%     Behavior_raw = F;
%     stimuluskey_raw = stim_full_raw;
%     timelists = tlists_raw;
%     timelists_names = stimrangenames;

author: xiuyechen

repo: FishExplorer


text:
%--------------------------------------------------------------------------
%
%                    %%%%%  %%%%%  %   %    %    %   %
%                    %        %    %%  %
%                    %  %%%   %    % % %   %%%    % %
%                    %   %    %    %  %%
%                    %%%%%  %%%%%  %   %  %   %    %
%
%
%--------------------------------------------------------------------------
%                            GINav v0.1.0
%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%--------------------------------------------------------------------------
% This script is used to plot and analyze results.It should be noted that
% all file formats must follow GINav's format requirements. For more
% details, please refer to the GINav user manual.
%--------------------------------------------------------------------------

clc
clear all %#ok
close all
addpath(genpath(pwd))
global_variable;

GINavPlot;

%--------------------------------------------------------------------------

repo: GINav


text:
From fessler@eecs.umich.edu Thu Jun  1 16:56 EDT 2000
To: millsk@engin.umich.edu
Subject: ml-em
Cc: fessler@eecs.umich.edu
Content-Type: text
Content-Length: 539

function f = mlem(f, g, h)

if nargin < 3
	ftrue = zeros(11); ftrue(3:5,4:7) = 1; ftrue(4,5) = 0;
	h = ones(3,3);
	g = conv2(ftrue, h, 'same');
	g = g + 0.1 * randn(size(g));
	g = max(g, 0);

	f = g;
	imagesc(f), drawnow
	for ii=1:100
		f = mlem(f, g, h);
		imagesc(f), drawnow
	end
return
end

if any(g < 0), error 'need nonnegative data', end
gp = conv2(f, h, 'same');
if any(g > 0 & ~gp), error 'model mismatch', end
ratio = g ./ (gp + eps * (gp==0));

sens = conv2(ones(size(f)), h, 'same');
f = f .* conv2(ratio, h, 'same') ./ sens;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


clearvars;
close all;
% Add path to stlTools
% You can download the package for free from:
% https://es.mathworks.com/matlabcentral/fileexchange/51200-stltools
addpath('./stlTools');
% Set the name of the mat file containing all the info of the 3D model
MatFileName = '../../3d_models/a320_3d_model.mat';
% Define the list of parts which will be part of the rigid aircraft body
% First part should always be main aircraft body
rigid_body_list   = { 'Body_A320.stl','Engine_L_A320.stl','Engine_R_A320.stl'};
% Define the color of each part
rigid_body_colors = {0.8 * [1, 1, 1], 0.1 * [1, 1, 1], 0.1 * [1, 1, 1]};
% Define the transparency of each part
alphas            = [              1,               1,               1];
% Define the model ofset to center the A/C Center of Gravity
offset   = [0,0,0];
% Define the control surfaces
ControlsFieldNames = {...
'model'                 'label',             'color',                'rot_point',           'rot_vect', 'max_deflection'};
Controls = {
'Aileron_L_A320.stl',     'A_L',   0.3*[0.8, 0.8, 1], [515.1,-1434,80.5]-offset, [ 0.2932,-0.9523, 0.0843], [-30, +30];
'Aileron_R_A320.stl',     'A_R',   0.3*[0.8, 0.8, 1], [515.1,+1434,80.5]-offset, [ 0.2932, 0.9523, 0.0843], [-30, +30];
'Rudder_A320.stl',        'RUD',   0.3*[0.8, 0.8, 1], [ 1954,     0,550.1]-offset,-[ 0.3064, 0, 0.9519], [-30, +30];
'Elevator_L_A320.stl', 'STAB_L',   0.3*[0.8, 0.8, 1], [ 2037,-352.9,168.5]-offset, [-0.2952, 0.9481,-0.1179], [-30, +30];
'Elevator_R_A320.stl', 'STAB_R',   0.3*[0.8, 0.8, 1], [ 2037,+352.9,168.5]-offset, [ 0.2952, 0.9481, 0.1179], [-30, +30];
};
% Definition of the Model3D data structure
% Rigid body parts
for i = 1:length(rigid_body_list)
    Model3D.Aircraft(i).model = rigid_body_list{i};
    Model3D.Aircraft(i).color = rigid_body_colors{i};
    Model3D.Aircraft(i).alpha = alphas(i);
    % Read the *.stl file
   [Model3D.Aircraft(i).stl_data.vertices, Model3D.Aircraft(i).stl_data.faces, ~, Model3D.Aircraft(i).label] = stlRead(rigid_body_list{i});
    Model3D.Aircraft(i).stl_data.vertices  = Model3D.Aircraft(i).stl_data.vertices - offset;
end
% Controls parts
for i = 1:size(Controls, 1)
    for j = 1:size(Controls, 2)
        Model3D.Control(i).(ControlsFieldNames{j}) = Controls{i, j};
    end
    % Read the *.stl file
    [Model3D.Control(i).stl_data.vertices, Model3D.Control(i).stl_data.faces, ~, ~] = stlRead( Model3D.Control(i).model);
    Model3D.Control(i).stl_data.vertices = Model3D.Control(i).stl_data.vertices - offset;
end

%% Save mat file
save(MatFileName, 'Model3D');

%% Check the results
plot3Dmodel(MatFileName)


text:
function [ all_results ] = evaluate_average_precision_pascal( ...
    all_bbox_gt, all_detected_bbox, classes, varargin)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

%************************** OPTIONS *************************************
ip = inputParser;
ip.addParamValue('minoverlap', 0.5,   @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

minoverlap  = opts.minoverlap;
all_results = compute_average_precision_of_detection(all_bbox_gt, all_detected_bbox, minoverlap);
% fprintf('\n~~~~~~~~~~~~~~~~~~~~\n');
% fprintf('Results:\n');
% aps = [all_results(:).ap]' * 100;
% disp(mean(aps));
% fprintf('~~~~~~~~~~~~~~~~~~~~\n');
end

function [false_positives, true_positives] = find_detection_labels(bbox_detections, bbox_gt, is_difficult, minoverlap)
false_positives = zeros(size(bbox_detections,1), 1);
true_positives  = zeros(size(bbox_detections,1), 1);
overlap         = zeros(size(bbox_detections,1), size(bbox_gt,1));
num_bbox_gt     = size(bbox_gt,1);

for j = 1:num_bbox_gt
    overlap(:,j) = boxoverlap(bbox_detections(:,1:4), bbox_gt(j,1:4));
end

overlap(overlap==0) = -inf;
[max_overlap, jmax] = max(overlap,[],2);

does_overlap = max_overlap >= minoverlap;

false_positives(~does_overlap) = 1; % false positive

does_overlap = does_overlap & ~is_difficult(jmax); % dont care about the difficult ones

if any(does_overlap)
    bbox_indices = find(does_overlap);
    jmax         = jmax(does_overlap);

    for j = 1:num_bbox_gt
        overlap_j = bbox_indices(jmax == j);
        if ~isempty(overlap_j)
            true_positives(overlap_j(1)) = 1; % true positive
            if numel(overlap_j) > 1
                false_positives(overlap_j(2:end)) = 1; % false positive - multiple detections
            end
        end
    end
end
end

function all_results = compute_average_precision_of_detection(all_bbox_gt, all_detected_bbox, minoverlap)
num_imgs       = length(all_bbox_gt);
num_classes    = length(all_detected_bbox);
true_positives   = cell(num_classes,1);
false_positives  = cell(num_classes,1);
detection_scores = cell(num_classes,1);

for class_idx = 1:num_classes
    true_positives{class_idx}   = cell(num_imgs,1);
    false_positives{class_idx}  = cell(num_imgs,1);
    detection_scores{class_idx} = cell(num_imgs,1);
end

num_positives = zeros(num_classes, 1);
for img_idx = 1:num_imgs
    for class_idx = 1:num_classes
        ground_truth_idx = all_bbox_gt{img_idx}(:,5) == class_idx;
        bbox_gt          = all_bbox_gt{img_idx}(ground_truth_idx,1:4);
        is_difficult     = all_bbox_gt{img_idx}(ground_truth_idx,6) > 0;
        num_positives(class_idx) = num_positives(class_idx) + sum(~is_difficult);

        bbox_detections = all_detected_bbox{class_idx}{img_idx};
        num_detections  = size(bbox_detections, 1);
        if num_detections > 0
            if isempty(bbox_gt)
                false_positives{class_idx}{img_idx} = ones(num_detections, 1);
                true_positives{class_idx}{img_idx}  = zeros(num_detections, 1);
            else
                [~, order]      = sort(bbox_detections(:,5), 'descend');
                bbox_detections = bbox_detections(order, :);

                [false_positives{class_idx}{img_idx},...
                    true_positives{class_idx}{img_idx}] = ...
                        find_detection_labels(bbox_detections(:,1:4), ...
                        bbox_gt, is_difficult, minoverlap);
            end
            detection_scores{class_idx}{img_idx} = double(bbox_detections(:,5));
        end
    end
end

all_results = compute_average_precision(true_positives, false_positives, ...
    detection_scores, 1:num_classes, num_positives);
end

function res = compute_average_precision(true_positives, false_positives, detection_scores, class_indices, num_positives)
for class_idx = class_indices
    all_true_positives  = cell2mat(true_positives{class_idx}(:));
    all_false_positives = cell2mat(false_positives{class_idx}(:));
    all_scores          = cell2mat(detection_scores{class_idx}(:));
    [all_scores, order] = sort(all_scores, 'descend');
    all_true_positives  = all_true_positives(order);
    all_false_positives = all_false_positives(order);

    % compute precision/recall
    all_false_positives = cumsum(all_false_positives);
    all_true_positives  = cumsum(all_true_positives);

    recall    = all_true_positives  /  num_positives(class_idx);
    precision = all_true_positives ./ (all_false_positives + all_true_positives);
    f1_score  = 2 * (recall.*precision) ./ (recall + precision + eps);
    [f1_score, max_idx] = max(f1_score);
    f1_thresh    = all_scores(max_idx);

    % compute average precision
    average_precion = 0;
    for t=0:0.1:1
        p = max(precision(recall>=t));
        if ~isempty(p)
            average_precion = average_precion + p/11;
        end
    end
    ap_auc = xVOCap(recall, precision);
    res(class_idx).recall    = recall;
    res(class_idx).precision = precision;
    res(class_idx).ap        = average_precion;
    res(class_idx).ap_auc    = ap_auc;
    res(class_idx).f1_score  = f1_score;
    res(class_idx).f1_thresh = f1_thresh;
%     fprintf('!!! %s : %.4f %.4f\n', classes{class_idx}, average_precion, ap_auc);
end
end

function ap = xVOCap(rec,prec)
mrec=[0 ; rec ; 1];
mpre=[0 ; prec ; 0];
for i=numel(mpre)-1:-1:1
    mpre(i)=max(mpre(i),mpre(i+1));
end
i=find(mrec(2:end)~=mrec(1:end-1))+1;
ap=sum((mrec(i)-mrec(i-1)).*mpre(i));
end

author: gidariss

repo: mrcnn-object-detection


text:
function [err, g] = metric_logssd(vol_fix, voldef, dvol)
    vol_diff = voldef - vol_fix;
    vdsq = vol_diff.^2;
    err = log(1 + sum(sum(sum(vdsq)))) * dvol;
    g = 2 * vol_diff ./ (1 + vdsq) * dvol;
end

repo: pTVreg


text:
function markdown2matlab(infile,outfile,varargin)

% MARKDOWN2MATLAB converts a Markdown file to a MATLAB script or function. All text
% is converted to comments, headings are converted to comment lines starting with %%
% sections with code are properly formatted, and words that appear in bold, italic or
% monospace are converted.
%
% Use as
%   markdown2matlab(infile, outfile)
%
% If no outfile is specified, it will write it to a .m file with the same name as
% the infile. In case the file exists, it will be written with a numeric suffix.
%
% The best is to provide the full filepath, otherwise it will look for the file within
% the current path.
%
% Optional input arguments can be specified as key-value pairs and can include
%   ...
%
% See also MATLAB2MARKDOWN

% Copyright (C) 2018 Sophie Arana and Robert Oostenveld
%
% This file is part of FieldTrip.
%
% megconnectome is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% megconnectome 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 megconnectome.  If not, see <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following links are useful during development
% https://regexr.com
% https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin>2 && ismember(outfile, {})
  varargin = [outfile varargin];
  outfile = [];
end

% parse the optional input arguments
% val = ft_getopt(varargin, 'key', default);

% check input
[inpath, inname, inext] = fileparts(infile);
if isempty(inpath), [inpath, inname, inext] = fileparts(which([inname inext])); end
if ~strcmp(inext,'.md')
  error('please specify a Markdown file')
end

if nargin < 2 || isempty(outfile)
  outfile = infile;
end

% check output
[outpath, outname,outext] = fileparts(outfile);
if isempty(outpath), outpath = inpath; end
if ~strcmp(outext,'.m')
  outext = '.m';
end

% add a suffix to avoid overwriting files
suffix = 1;
newname = outname;
while exist(fullfile(outpath,[newname,outext]),'file') == 2
  newname = [outname sprintf('-%0d',suffix)];
  suffix = suffix+1;
end
outname = newname;

infile = fullfile(inpath,[inname,inext]);
outfile = fullfile(outpath,[outname,outext]);

% read & convert file line by line
infid = fopen(infile,'r');
outfid = fopen(outfile,'w');

state = 'comment';    % keeps track of multi-line formatting decisions
linenumber = 0;       % keep track of the line number

while ~feof(infid)

  % reset the state, unless explicitly specified
  reset_state = true;

  line = fgetl(infid);
  linenumber = linenumber + 1;

  if ~ischar(line), break, end

  % replace each tab by two spaces
  line = strrep(line, sprintf('\t'), '  ');

  if match(line, '^---$') && linenumber==1
    state = 'jekyllheader';
    reset_state = false;

  elseif strcmp(state, 'jekyllheader') && ~match(line, '^---$')
    % do not output the header
    reset_state = false;

  elseif strcmp(state, 'jekyllheader') && match(line, '^---$')
    % do not output the end of the header

  elseif strcmp(state, 'codeblock') && match(line, '^```$')
    % do not output the end of the code block
    reset_state = true;

  elseif match(line, '^```')
    state = 'codeblock';
    reset_state = false;

  elseif strcmp(state, 'codeblock') && ~match(line, '^```')
    % output the code as-is
    fprintf(outfid, '%s\n', line);
    reset_state = false;

  elseif match(line, '{%[ -~]*%}')
    % jekyll code should not be converted into the MATLAB script

  elseif match(line, '!\[[ -~]*\]([ -~]*)')
    % inline images should not be converted into the MATLAB script

  elseif match(line, '^ *[-+*] ')
    % unordered list
    [dum,endIndex] = regexp(line, '^ *[-+*] ');
    [dum,level] = regexp(line, '^ *');
    if isempty(level)
      level = 0;
    end
    level = repmat(' ', 1, level);
    remainder = reformat(line((endIndex+1):end));
    fprintf(outfid, '%%%s* %s\n', level, remainder);

  elseif match(line, '^    ')
    % normal code
    [dum,endIndex] = regexp(line, '^    ');
    remainder = line((endIndex+1):end);
    fprintf(outfid, '%s\n', remainder);
    state = 'code';
    % assume that the next line will be in the same format
    reset_state = false;

  elseif match(line, '^[A-Za-z]')
    % normal text
    remainder = reformat(line);
    fprintf(outfid, '%% %s\n', remainder);

  elseif match(line, '^ *[0-9]*\.')
    % ordered list
    [dum,endIndex] = regexp(line, '^ *[0-9]*\. ');
    remainder = reformat(line((endIndex+1):end));
    fprintf(outfid, '%% # %s\n', remainder);

  elseif match(line, '^#')
    % heading
    [dum,endIndex] = regexp(line, '^# *');
    remainder = reformat(line((endIndex+1):end));
    fprintf(outfid, '%%%% %s\n', remainder);

  elseif isempty(line)
    % keep the complete line as it is
    if strcmp(state, 'comment')
      fprintf(outfid, '%%\n');
    else
      fprintf(outfid, '\n');
    end

  end % if the line matches some pattern

  if reset_state
    state = 'comment';
  end

end % while not feof

fclose(infid);
fclose(outfid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function status = match(str, pattern)
status = ~isempty(regexp(str, pattern, 'once'));

function str = reformat(str)
str = formatbold(str);
str = formatitalic(str);
str = formatunderline(str);
str = formatmonospace(str);

function str = formatbold(str)
[startIndex,endIndex] = regexp(str, '\*\*[a-zA-Z0-9]*\*\*');
for i=1:length(startIndex)
  part1 = str(1:startIndex(i)-1);
  part2 = str(startIndex(i):endIndex(i));
  part3 = str((endIndex(i)+1):end);
  str = [part1 '*' part2(3:end-2) '*' part3];
  % correct the subsequent indices
  startIndex(i+1:end) = startIndex(i+1:end)-2;
  endIndex(i+1:end) = endIndex(i+1:end)-2;
end

function str = formatitalic(str)
% the format is actually the same
[startIndex,endIndex] = regexp(str, '_[a-zA-Z0-9]*_');
for i=1:length(startIndex)
  part1 = str(1:startIndex(i)-1);
  part2 = str(startIndex(i):endIndex(i));
  part3 = str((endIndex(i)+1):end);
  str = [part1 '_' part2(2:end-1) '_' part3];
end

function str = formatunderline(str)
% this is not supported in the conversion from MATLAB live scripts to normal MATLAB code.

function str = formatmonospace(str)
[startIndex,endIndex] = regexp(str, '`[a-zA-Z0-9]*`');
for i=1:length(startIndex)
  part1 = str(1:startIndex(i)-1);
  part2 = str(startIndex(i):endIndex(i));
  part3 = str((endIndex(i)+1):end);
  str = [part1 '|' part2(2:end-1) '|' part3];
  % correct the subsequent indices
  startIndex(i+1:end) = startIndex(i+1:end);
  endIndex(i+1:end) = endIndex(i+1:end);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = r8vec_scalar_triple_product ( v1, v2, v3 )

%*****************************************************************************80
%
%% R8VEC_SCALAR_TRIPLE_PRODUCT finds the scalar triple product in 3D.
%
%  Discussion:
%
%    [A,B,C] = A dot ( B cross C )
%            = B dot ( C cross A )
%            = C dot ( A cross B )
%
%    The volume of a parallelepiped, whose sides are given by
%    vectors A, B, and C, is abs ( A dot ( B cross C ) ).
%
%    Three vectors are coplanar if and only if their scalar triple
%    product vanishes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eric Weisstein,
%    "Scalar Triple Product",
%    CRC Concise Encyclopedia of Mathematics, 1999
%
%  Parameters:
%
%    Input, real V1(3), V2(3), V3(3), the vectors.
%
%    Output, real VALUE, the scalar triple product.
%
  v4 = r8vec_cross_product_3d ( v2, v3 );

  value = v1' * v4;

  return
end


text:
datasetPath='/home/bernard/HD1/ABC/sun';
addpath('..');
lambdas = 10.^[-3:3];
sigmas=10.^[-3:3];
tasks=[0];

validationClasses_generalscript(datasetPath, lambdas, sigmas, tasks);

repo: Embarrassingly-simple-ZSL


text:
function [v] = tapas_mh_mc3_arc(ollh, olpp, nllh, nlpp, ...
    ratio, T)
%% Acceptance rejection criterion for metropolis hastings in the context of
% population mcmc.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

v = bsxfun(@times, nllh, T) + nlpp - (bsxfun(@times, ollh, T) + olpp) ...
    + ratio;

nansv = isnan(v);
v(nansv) = -inf;

v = rand(size(v)) < exp(v);

assert(all(-inf < nllh(v) + nlpp(v)), 'tapas:mh', ...
    '-inf value in the new samples');

end

author: translationalneuromodeling

repo: tapas


text:
% This is a simple example of how to fit
% a cubic spline contour
% to data using active geometric shape model.

% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Quan Wang, Kim L. Boyer,
% The active geometric shape model: A new robust deformable shape model and its applications,
% Computer Vision and Image Understanding, Volume 116, Issue 12, December 2012, Pages 1178-1194,
% ISSN 1077-3142, 10.1016/j.cviu.2012.08.004.
%
% For commercial use, please contact the authors.

clear;clc;close all;
addpath('../force field');
addpath('../math');

%% 1. experiment set up

% image size = 400*500
rows=400;
cols=500;

% number of trials
Num_Of_Trials=10;

% number of data points and outliers
num_data=100;
num_outlier=5;
noise=5;

% number of landmarks for data/model
Nlm=6; % data
NN=6; % model

% ground truth parameters
x0=250+(rand(1)-0.5)*100;
y0=200+(rand(1)-0.5)*100;
D0=60+rand(1,Nlm)*100;

% image blur parameters
sigma1=20;
sigma2=5;

%% 2. generate data

% generate data points and outliers
theta=rand(num_data,1)*2*pi;
lm_theta=linspace(0,2*pi,Nlm+1);
DD=myspline(lm_theta,D0,theta);
x=x0+DD.*cos(theta)+noise*randn(size(theta));
y=y0+DD.*sin(theta)+noise*randn(size(theta));
x=[x;rand(num_outlier,1)*cols];
y=[y;rand(num_outlier,1)*rows];
x=round(x);
y=round(y);
x(x<1)=1;
x(x>cols)=cols;
y(y<1)=1;
y(y>rows)=rows;
xp=x;
yp=y;

% display data points and outliers
plot(x,y,'.');
axis equal;
axis([1 cols 1 rows]);
title('data points and outliers');
xlabel('x');
ylabel('y');
drawnow;

% generate image
I=zeros(rows,cols);
for i=1:max(size(y))
    I(y(i),x(i))=100;
end
[m n]=size(I);

%% 3. GVF field
I2=gaussianBlur(I,sigma1);
[u,v] = GVF(I2, 1 , 0.1, 50);
dx2=u;dy2=v;

I3=gaussianBlur(I,sigma2);
[dx3 dy3]=gradient(I3);

best_fit=Inf;

%% 4. many trials
disp('Fitting the cubic spline contour...');
for trial=1:Num_Of_Trials
    close all;
    fprintf('Trial: %d\n',trial);
    init=[250+(rand(1)-0.5)*100,200+(rand(1)-0.5)*100,60+rand(1,NN)*40];
    increment=[0.2,0.2,0.2];
    threshold=[0.000001,0.000001,0.000001];
    bound=[40,150];
    %% fit a spline contour
    % stage 1: big sigma
    [xc yc D fit_save]=fit_spline_contour_force(...
        init,increment,threshold,bound,dx2,dy2,500,0);

    % stage 2: small sigma
    [xc yc D]=fit_spline_contour_force(...
        [xc yc D],increment,threshold,bound,dx3,dy3,100,0);

    %% correction
    clear D2;
    for k=1:NN
        dt=2*pi/100000;
        angle=2*pi/NN*(k-1);

        lm_theta=linspace(0,2*pi,NN+1);
        distances=myspline(lm_theta,D,[angle-dt angle angle+dt]);
        r=D(k);
        r1=(distances(3)-distances(1)) / (2*dt);
        r2=(distances(3)+distances(1)-2*distances(2)) / (dt^2);
        D2(k)=correctCurve_polar(r,r1,r2,sigma2,100);
    end
    D=D2;

    %% fitness function
    beta=0.9;
    fit0=0;
    fit1=0;
    fit2=0;
    [x,y,theta]=spline_contour_in_image(m,n,xc,yc,D);
    [x1,y1,theta1]=spline_contour_in_image(m,n,xc,yc,D*beta);
    [x2,y2,theta2]=spline_contour_in_image(m,n,xc,yc,D/beta);
    for i=1:max(size(theta))
        fit0=fit0+norm([dx2(y(i),x(i)),dy2(y(i),x(i))]);
    end
    for i=1:max(size(theta1))
        fit1=fit1+norm([dx2(y1(i),x1(i)),dy2(y1(i),x1(i))]);
    end
    for i=1:max(size(theta2))
        fit2=fit2+norm([dx2(y2(i),x2(i)),dy2(y2(i),x2(i))]);
    end
    fit0=fit0/max(size(theta));
    fit1=fit1/max(size(theta1));
    fit2=fit2/max(size(theta2));
    current_fit=fit0-fit1/2-fit2/2;

    if current_fit<best_fit
        best_fit=current_fit;
        best_xc=xc;
        best_yc=yc;
        best_D=D;
        best_fit_save=fit_save;
    end
end

% display fitting results
close all;

figure;
plot(xp,yp,'.');
axis equal;
axis([1 cols 1 rows]);
hold on;
theta=0:0.01:2*pi;
[x2,y2,temp]=spline_contour_in_image(rows,cols,best_xc,best_yc,best_D);
plot(x2,y2,'-.r','LineWidth',2);
plot(xp,yp,'.');
legend('noisy data points and outliers','active geometric shape model fit');
title('fitting a cubic spline contour using active geometric shape model');
xlabel('x');
ylabel('y');

%% 5. show numerical results
fprintf('num_data=%-4d    num_outlier=%-4d    \n',num_data,num_outlier);
fprintf('Nlm(data)=%-4d   Nlm(model)=%-4d    \n',Nlm,NN);
fprintf('------------------------------------------------------------\n');
fprintf('         xc        yc        D   \n');
fprintf('True:  ');
disp([x0 y0 D0]);
fprintf('Fit:   ');
disp([best_xc best_yc best_D]);

%% 6. show fitness function
figure;hold on;
plot(1:500,best_fit_save,'b','LineWidth',2);
legend('fitness function');
grid on;
xlabel('iteration');
ylabel('fitness function');
title('fitness function in each iteration');

author: Sable

repo: mcbench-benchmarks


text:
function err = get_cross_entropy(model, data, l)
% cross-entropy(reconstruction error, I use L2) of each layer. This is used
% to monitor the pretraining progress.
% data: can be real data(for rbm) or data filelist(crbm).

global kConv_backward  kConv_backward_c kConv_forward2 kConv_forward_c;

fraction = 5;
if l == 2
    n = length(data);
    batch_size = 32; batch_num = n / batch_size;
    assert(batch_num == floor(batch_num));
    stride = model.layers{l}.stride;

    batch_num = floor(batch_num / fraction); % evaluate 1/fraction% of data for speed
    shuffle_index = randperm(n);
    err = 0;
    for b = 1 : batch_num
        batch_index = shuffle_index((b-1)*batch_size + 1 : b * batch_size);
        batch = read_batch(model, data(batch_index), false);

        hidden_presigmoid = myConvolve2(kConv_forward2, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        hidden_prob = sigmoid(hidden_presigmoid);
        hidden_sample = single(hidden_prob > rand(size(hidden_prob)));

        visible_presigmoid = myConvolve(kConv_backward, hidden_sample, model.layers{l}.w, stride, 'backward');
        visible_presigmoid = bsxfun(@plus, visible_presigmoid, permute(model.layers{l}.b, [5,1,2,3,4]));
        visible_prob = sigmoid(visible_presigmoid);

        this_err = batch - visible_prob;
        err = err + sum(this_err(:).^2);
    end
    err = err / (batch_num * batch_size);
elseif strcmp(model.layers{l}.type, 'convolution')
    n = length(data);
    batch_size = 32; batch_num = n / batch_size;
    assert(batch_num == floor(batch_num));
    stride = model.layers{l}.stride;

    batch_num = floor(batch_num / fraction); % evaluate 1/fraction% of data for speed
    shuffle_index = randperm(n);
    err = 0;
    for b = 1 : batch_num
        batch_index = shuffle_index((b-1)*batch_size + 1 : b * batch_size);
        batch = read_batch(model, data(batch_index), false);
        batch = propagate_batch(model, batch, l);

        hidden_presigmoid = myConvolve(kConv_forward_c, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        hidden_prob = sigmoid(hidden_presigmoid);
        hidden_sample = single(hidden_prob > rand(size(hidden_prob)));

        visible_presigmoid = myConvolve(kConv_backward_c, hidden_sample, model.layers{l}.w, stride, 'backward');
        visible_presigmoid = bsxfun(@plus, visible_presigmoid, permute(model.layers{l}.b, [5,1,2,3,4]));
        visible_prob = sigmoid(visible_presigmoid);

        this_err = batch - visible_prob;
        err = err + sum(this_err(:).^2);
    end
    err = err / (batch_num * batch_size);
elseif model.classes > 0 && l == length(model.layers) % for the last layer
    n = size(data,1);
    batch_size = 32; batch_num = n / batch_size;
    err = 0;
    temp_w = model.layers{l}.w;
    temp_w(1:model.classes,:) = temp_w(1:model.classes,:) * model.duplicate;
    for b = 1 : batch_num
        batch = data((b-1)*batch_size + 1 : b * batch_size,:,:,:,:);
        hidden_presigmoid = bsxfun(@plus, ...
			batch * temp_w, model.layers{l}.c);
		hidden_prob = 1 ./ ( 1 + exp(- hidden_presigmoid) );
		hidden_sample = single(hidden_prob > rand(size(hidden_prob)));

        visible_presigmoid = bsxfun(@plus, ...
					hidden_sample * model.layers{l}.w', model.layers{l}.b);
        visible_prob_post = sigmoid(visible_presigmoid(:,model.classes+1:end));
        temp_exponential = exp(bsxfun(@minus,visible_presigmoid(:,1:model.classes),max(visible_presigmoid(:,1:model.classes),[],2)));
        visible_prob_pre = bsxfun(@rdivide, temp_exponential, sum(temp_exponential,2));

        data_err = batch(:,model.classes+1:end) - visible_prob_post;
        label_err = batch(:, 1:model.classes) - visible_prob_pre;
        err = err + model.duplicate * sum(label_err(:).^2) + sum(data_err(:).^2);
    end
    err = err / n;
else % for original rbm
    n = size(data,1);
    batch_size = 32; batch_num = n / batch_size;
    err = 0;
    for b = 1 : batch_num
        batch = data((b-1)*batch_size + 1 : b * batch_size,:,:,:,:);
        hidden_presigmoid = bsxfun(@plus, ...
			batch * model.layers{l}.w, model.layers{l}.c);
		hidden_prob = 1 ./ ( 1 + exp(- hidden_presigmoid) );
		hidden_sample = single(hidden_prob > rand(size(hidden_prob)));

        visible_presigmoid = bsxfun(@plus, ...
					hidden_sample * model.layers{l}.w', model.layers{l}.b);
        visible_prob = 1 ./ ( 1 + exp(-visible_presigmoid) );

        this_err = batch - visible_prob;
        err = err + sum(this_err(:).^2);
    end
    err = err / n;
end

function y = sigmoid(x)
    y = 1 ./ ( 1 + exp(-x) );

author: zhirongw


text:
%% Floyd algorithm
function [S, P]=FloydSPR(AdjMax)
% *INPUT:*
% AdjMax: Adjacent matrix that represents a weighted, directed graph
%
% *OUTPUT:*
% S: distance to destination node
% P: next hop node
%
% *DESCRIPTION*
% Given a input adjacent matrix (AdjMax) that represents a weighted, directed graph.
% The function finds the shorest path from one vertex 'i' to another 'j'.
% The return values includes a matrix (S) that denotes the shortest distance
% between vertices 'i' and 'j', and a matrix (P) that denotes the next vertex 'k'
% on the path from vertex 'i' to vertex 'j'

N=min(length(AdjMax(:,1)),length(AdjMax(1,:)));
P=-1*ones(N,N);
S=AdjMax;

for k=1:N
    for i=1:N
        for j=1:N
            if S(i,k)==inf
                continue;
            end
            if S(k,j)==inf
                continue;
            end
            if S(i,j)>S(i,k)+S(k,j)
                if P(i,k)==-1
                    P(i,j)=k;
                else
                    P(i,j)=P(i,k);
                end
                S(i,j)=S(i,k)+S(k,j);
            end
        end
    end
end

%% END


text:
function [Xmat,idxisa] = workvolume(cam,room,imres,idxcams)

STEP = 0.1;

if nargin < 4
	idxcams = [1:size(cam,2)];
end

if nargin < 3
	imres = repmat([640 480],size(idxcams,2),1);
end
imres = imres(idxcams,:);

if nargin < 2
	% room [x_min, x_max, y_min, y_max, z_min, z_max]
	room = [-3 3 -3 3 0 3];
end

% compose Pmat containing all P matrices
Pmat = [];
for i=idxcams,
	Pmat = [Pmat; cam(i).P];
end

% create points

zcoor = room(5):STEP:room(6);
znum = size(zcoor,2);
ycoor = room(3):STEP:room(4);
ynum = size(ycoor,2);
xcoor = room(1):STEP:room(2);
xnum = size(xcoor,2);

buff  = repmat(ycoor,znum,1);
yvec  = reshape(buff,prod(size(buff)),1);
zvec  = repmat(zcoor',ynum,1);

buff = repmat(xcoor,znum*ynum,1);
xvec = reshape(buff,prod(size(buff)),1);

Xmat = [xvec,repmat([yvec,zvec],xnum,1)];
Xmat = [Xmat,ones(size(Xmat(:,1)))];

clear buff

size(Xmat)
umat = Pmat*Xmat';

% normalize projected points in umat
scalemat = [];
for i=1:size(idxcams,2),
	scalemat = [scalemat; repmat(umat(3*i,:),3,1)];
end
umat = umat./scalemat;
clear scalemat;
mask = zeros(size(umat));
mask(1:3:end,:) = umat(1:3:end,:)<repmat(imres(:,1),1,size(mask,2));
mask(2:3:end,:) = umat(2:3:end,:)<repmat(imres(:,2),1,size(mask,2));
mask(3:3:end,:) = 1;
mask = mask.*(umat>0);
idxisa = find(sum(mask)==size(mask,1));

return;

author: strawlab

repo: MultiCamSelfCal


text:
function fp = fmf_write_header( filename, header_size, version, ...
                                f_height, f_width, bytes_per_chunk, ...
                                max_n_frames, data_format)
% fp = fmf_write_header( filename, header_size, version,
%                        f_height, f_width, bytes_per_chunk,
%                        max_n_frames, data_format )
%
% writes FlyMovieFormat header data to FILENAME and returns a file
% pointer FP.
%
% VERSION is the fmf version
% HEADER_SIZE is in bytes for the appropriate file format (fmf VERSION)
% F_HEIGHT and F_WIDTH are the number of pixels in a frame
% BYTES_PER_CHUNK is the number of bytes per frame in the file
% MAX_N_FRAMES is the number of frames in the file
% DATA_FORMAT is the format in which each pixel is stored.
%
% KMB 04/07

fp = fopen( filename, 'w' );
if fp <= 0,
  error(['Could not open file ',filename,' for writing.']);
end;

% write header
if version ~= 1 & version ~= 3,
  error( 'version not supported -- FMF versions 1 and 3 only' );
end;

fwrite(fp,version,'uint32');
if version == 1,
  data_format = 'MONO8';
  bits_per_pixel = 8;
end
if version == 3,
  fwrite(fp,length(data_format),'uint32');
  fwrite(fp,data_format);
  bits_per_pixel = fmf_get_bits_per_pixel(data_format);
  fwrite(fp,bits_per_pixel,'uint32');
end

fwrite(fp,f_height,'uint32');
fwrite(fp,f_width,'uint32');
fwrite(fp,bytes_per_chunk,'integer*8');
fwrite(fp,max_n_frames,'integer*8');
curr = ftell(fp);
if header_size > curr,
  fwrite(fp,zeros(header_size-curr,1), 'uint8');
end;

author: kristinbranson

repo: JAABA


text:
%kkhypot 'Output = sqrt[(Input 1)**2 + (Input 2 or Constant)**2]'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros khypot.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'First input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% InputFile: i2 'Input 2', optional: 'Second input data object'
%
% Example: o = kkhypot({i1, i2}, {'i1','';'o','';'i2',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% khypot - Output = sqrt[(Input 1)**2 + (Input 2 or Constant)**2]
%
%  DESCRIPTION
% The "Hypotenuse" operator returns the square root of the result
% produced by adding the square of each point in \fBInput 1" and the square
% of either the corresponding point in \fBInput 2" or the \fBConstant\fP
% value, which ever is specified by the user.
%
% Executing "Hypotenuse" runs the program \fIkarith2\fP with the -hypot flag.
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
%
%  "Map Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "Map Data - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_2input
%
%  "Validity Mask" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_2input
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The values used to pad the data when input files are not the same size are
% (0.0, 0.0).
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_2input
%
%  "Failure Modes - Single Input"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%  "Failure Modes - Two Input Objects"
% .cI $DATAMANIP/repos/shared/man/sections/fail_2input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith2, DATAMANIP::kcmplx2real, DATAMANIP::kmag
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.  If the
% input object is complex, only the real component of the data is
% processed, and the output object is real.  To calculate the hypotenuse,
% or magnitude of a complex number, use DATAMANIP::kcmplx2real.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkhypot(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkhypot(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input'};
maxval={0,0,1};
minval={0,0,1};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith2"  -hypot'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function c = power(a,b)
% POWER implements a.^b, where either a or bor both is an adiff object.

switch [class(a),class(b)]

case 'adiffdouble'
   c = adiff( a.x.^b, rowmult(b.*a.x.^(b-1), a.dx), a.root);

case 'doubleadiff'
   c = adiff( a.^b.x, rowmult(a.^b.x.*log(a), b.dx), b.root);

case 'adiffadiff'
   c = exp(log(a).*b);
   % Commented out code does it the hard way:
   % cx  = a.x.^b.x;
   % cdx = rowmult(cx.*log(a.x), b.dx) + rowmult(a.x.^(b.x-1), a.dx);
   % c   = class(struct('x',cx,'dx',cdx),'adiff');

otherwise
   error(['Can''t do ',class(a),'.^',class(b)]);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [f] = spm_mc_fx_1(x,v,P)
% equations of motion for the mountain car problem using basis functions
% problem
% FORMAT [f] = spm_mc_fx_1(x,v,P)
%
% x   - hidden states
% v   - exogenous inputs
% P.p - parameters for gradient function:     G(x(1),P.p)
% P.q - parameters for cost or loss function: C(x(1),P.q)
%
% returns f = dx/dt = f  = [x(2);
%                           G - x(2)*C]*dt;
%
% where C determines divergence of flow x(2) at any position x(1).
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mc_fx_1.m 3140 2009-05-21 18:38:17Z karl $

% gradient (G)
%--------------------------------------------------------------------------
G  = spm_DEM_basis(x(1),v,P.p);

% cost function (C)
%--------------------------------------------------------------------------
C  = spm_DEM_basis(x(1),v,P.q);

% flow
%--------------------------------------------------------------------------
dt = 1/8;
f  = [x(2);
      G - x(2)*C]*dt;

% true scalar potential gradient (see spm_moutaincar_fx)
%--------------------------------------------------------------------------
% if x(1) < 0
%     G  = 2*x(1) + 1;
% else
%     xx = x(1)^2;
%     G  = (1 + 5*xx)^(-1/2) - 5*xx/(1 + 5*xx)^(3/2) + (x(1)/2)^4;
% end

author: spm


text:
function XYZ = quat2cube(q)
%
% transforms unit quaternions (representing rotations) into cubochoric
%   representation
% the transformation is achieved by going from representation as unit
%   quaternions to homochoric representation (Lambert) and then further to
%   cubochoric representation (ball2cube)
%
% Input
%  q - @quaternion
%
% Output
%  XYZ - cubochoric coordinates

xyz = Lambert([q.a(:) q.b(:) q.c(:) q.d(:)]);
XYZ = ball2cube(xyz);

end

author: mtex-toolbox

repo: mtex

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


text:
function data = noZero(data)
% Return the array without zeros
%
% SYNTAX:
%   data = zeros(data)
%
% SEE ALSO
%   noZero nanRep

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    data = data(data ~= 0);
end

author: goGPS-Project


text:
function view = cmapRedgreenblue(view, displayMode, hemi);
%
% view = cmapRedgreenblue(view, [displayMode='ph'], [hemi=1]);
%
% Set the color map for a view's display mode (default phase mode) to
% be a 'redgreenblue' cmap, a custom map Rory has found is useful in
% lookint at retinotopy data.
%
% The 'hemi' flag can be set to 0, 1 or 2: if 1, will make the second half of
% the colors be uniform gray, as if it was representing a non-mapped range
% of phases. If 0, will fill it with a blue->yellow->red submap, for
% looking at bilateral activations (e.g., in two hemispheres). If 2, will
% fill the entire cmap with the r->g->b cmap. (You can later reduce the range
% using cmapExtended. [Default is 1: only show 1 hemifield.]
%
% The redgreenblue cmap shows simple transititions between three colors
%
%
% ras, 10/2006.
% ras, 04/2009 -- got the 'both' option to work well.
if notDefined('view'),              view = getCurView;      	end
if notDefined('displayMode'),       displayMode = 'ph';         end
if notDefined('hemi'),              hemi = 1;                   end
modeInfo = view.ui.([displayMode 'Mode']);

nG = modeInfo.numGrays;
nC = modeInfo.numColors;

% build two colormaps: a red->green->blue map,
% and a blue->yellow->red map, each using half
% the available colors:
rgb = mrvColorMaps('redgreenblue', round(nC/2));
switch hemi
	case 0
		byr = repmat(.6, [nC/2 3]);
		subCmap = [rgb; byr];

	case 1
			byr = mrvColorMaps('redyellowblue', round(nC/2));
			byr = flipud(byr);

% 		% alternate mnemonic:
% 		% red is RIGHT, green is LEFT, blue is UP (like sky), yellow is DOWN
% 		% (like sand)
% 		gbr = mrvColorMaps('greenbluered', round(nC/2), 2);
% 		ryg = flipud(gbr);
% 		ryg(:,2) = ryg(:,2) + ryg(:,3);
% 		ryg(:,1) = ryg(:,1) + ryg(:,3);
% 		ryg(:,3) = 0;
% 		ryg(ryg > 1) = 1;  ryg(ryg < 0) = 0;
% 		subCmap = flipud([gbr; ryg]);
% 		subCmap = circshift(subCmap, round(nC/4));
		subCmap = [rgb; byr];

	case 2
		byr = flipud( mrvColorMaps('redgreenblue', nC-size(rgb,1)) );
		subCmap = [rgb; byr];

	otherwise
		error('Invalid hemifield flag.');
end

modeInfo.cmap = [gray(nG); subCmap];
modeInfo.name = 'redgreenblue cmap';

% override: let you set a the rgb cmap as the whole deal:
if hemi==2, modeInfo.cmap = [gray(nG); mrvColorMaps('redgreenblue', nC)]; end
% if hemi==2, modeInfo.cmap = [gray(nG); mrvColorMaps('gbr', nC, 4)]; end

view.ui.([displayMode 'Mode']) = modeInfo;

% view = refreshScreen(view);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft