text:
function fh = GNSSspectrumgen
%==========================================================================
%GNSSspectrumgen Functions to generate GNSS analytical spectra.
%   fh = GNSSspectrumgen  returns the handlers of the local functions:
%
%   spectrumgen_call
%   [SI SQ] = spectrumgen_call(signal,F) returns the analytical complex
%   spectra of a GNSS signal with normalized complex power.
%
%       Inputs
%           signal --> GPS L1: L1CA, L1P, L1M, L1C, L1Cd, L1Cp, L1Cp1,L1Cp2
%                              L1, L1_new.
%                      GPS L2: L2C, L2P, L2M, L2, L2_new.
%                      GPS L5: L5I, L5Q, L5, L5_new
%                      Galileo E1: E1PRS/EA, E1OS.
%                      Galileo E6: E6PRS/E6OS.
%                      Galileo E5: E5, E5A, E5B.
%                      BeiDou-2 Current:
%                      	B1: B11, B12, B1
%                       B2: B2I, B2Q, B2
%                       B3: B3
%                      BeiDou-2 Future:
%                       B1: B1Cd, B1Cp, B1C, B1_new
%                       B2: B2_new
%                       B3: B3_new, B3A, B3composite
%           F --> Baseband frequency points.
%
%       Ouputs
%           SI/SQ --> Normalized complex spectrum.
%
%   spectrum_BPSK
%   S = spectrum_BPSK(fc,F) returns spectrum of a BPSK modulation with a
%   chipping rate fc.
%
%   spectrum_BOCs
%   S = spectrum_BOCs(n,m,F) returns spectrum of a sine-phased even BOC
%   modulation with a chipping rate fc=m*1.023e6 and sub-carrier rate
%   fs = n*1.023e6.
%
%   spectrum_BOCc
%   S = spectrum_BOCc(n,m,F) returns spectrum of a cosine-phased even BOC
%   modulation with a chipping  rate fc=m*1.023e6 and sub-carrier rate
%   fs = n*1.023e6.
%
%   spectrum_AltBOC returns spectrum of a sine-phased
%   modified even AltBOC modulation with a chipping rate fc=m*1.023e6 and
%   sub-carrier rate fs = n*1.023e6.
%
%   Observations
%       # I think that L2 can also transmit a C/A component, but that they
%       actually never do it.
%       # I am not sure if the future BeiDou signals will replace or
%       complement the former signals. Now I assume they will replace them.
%       # I cound't find the official transmitted powers of the BeiDou
%       signals anywhere. I have assumed they are -163dBm per component.
%
%   References
%       # L1CA/L2: GPS Interface Control Document IS-GPS-200
%       # L1C:  GPS Interface Control Document IS-GPS-800
%       # L5: GPS Interface Control Document IS-GPS-705
%       # Galileo: Galileo Open Service Signal In Space Interface Control
%       Document (OS SIS ICD)
%       # BeiDou: BeiDou Navigation Satellite System Signal In Space
%       Interface Control Document 2.0
%--------------------------------------------------------------------------
% Version log (main changes)
%   02/03/2017 --> Log started
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

% Copyright 2017 Daniel Pascual
%
% 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/>.

    fh.spectrumgen_call = @spectrumgen_call;
    fh.spectrum_BPSK = @spectrum_BPSK;
    fh.spectrum_BOCs = @spectrum_BOCs;
    fh.spectrum_BOCc = @spectrum_BOCc;
    fh.spectrum_AltBOC = @spectrum_AltBOC;
end

function [SI, SQ] = spectrumgen_call(signal,F)
%==========================================================================
%
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

    switch(signal)
        %-------------------- All GNSS individual -------------------------
        % BPSK
        case {'L1CA', 'L2C'}
            SI = spectrum_BPSK(1.023e6,F);
            SQ = zeros(size(SI));

        case {'E6CS'}
            SI = spectrum_BPSK(5*1.023e6,F);
            SQ = zeros(size(SI));

        case {'L1P', 'L2P', 'L5I', 'L5Q'}
            SI = spectrum_BPSK(10.23e6,F);
            SQ = zeros(size(SI));

        % BOCs
        case{'L1Cd','L1Cp1'}
            SI = spectrum_BOCs(1,1,F);
            SQ = zeros(size(SI));

        case{'L1Cp2'}
            SI = spectrum_BOCs(6,1,F);
            SQ = zeros(size(SI));

        case {'L1M','L2M'}
            SI = spectrum_BOCs(10,5,F);
            SQ = zeros(size(SI));

        case{'L1Cp'}
            P_Cp1 = GNSS_POWERS.GPS_L1_Cp1/GNSS_POWERS.GPS_L1_Cp;
            P_Cp2 = GNSS_POWERS.GPS_L1_Cp2/GNSS_POWERS.GPS_L1_Cp;

            S_Cp1 = P_Cp1*spectrumgen_call('L1Cp1',F);
            S_Cp2 = P_Cp2*spectrumgen_call('L1Cp2',F);

            SI = S_Cp1 + S_Cp2;
            SQ = zeros(size(SI));

        case{'L1C'}
            P_Cd = GNSS_POWERS.GPS_L1_Cd/GNSS_POWERS.GPS_L1_C;
            P_Cp = GNSS_POWERS.GPS_L1_Cp/GNSS_POWERS.GPS_L1_C;

            S_Cd = P_Cd*spectrumgen_call('L1Cd',F);
            S_Cp = P_Cp*spectrumgen_call('L1Cp',F);

            SI = S_Cd;
            SQ = S_Cp;

        case 'E1OS'
            SI = (10/11)*spectrum_BOCs(1,1,F)+(1/11)*spectrum_BOCs(6,1,F);
            SQ = zeros(size(SI));

        % BOCc
        case {'E6A','E6PRS'}
            SI = spectrum_BOCc(10,5,F);
            SQ = zeros(size(SI));

        case {'E1A','E1PRS'}
            SI = spectrum_BOCc(15,2.5,F);
            SQ = zeros(size(SI));

        %-------------------- GPS L1  -------------------------------------
        % Prior to Block III
        case 'L1'
            P_M = GNSS_POWERS.GPS_L1_M/GNSS_POWERS.GPS_L1;
            P_P = GNSS_POWERS.GPS_L1_P/GNSS_POWERS.GPS_L1;
            P_CA =  GNSS_POWERS.GPS_L1_CA/GNSS_POWERS.GPS_L1;

            S_CA = P_CA*spectrumgen_call('L1CA',F);
            S_P = P_P*spectrumgen_call('L1P',F);
            S_M = P_M*spectrumgen_call('L1M',F);

            SI = S_P+S_M;
            SQ = S_CA;

        % Block III (includes L1C)
        case 'L1_new'

            P_M = GNSS_POWERS.GPS_L1_M/GNSS_POWERS.GPS_L1_new;
            P_P = GNSS_POWERS.GPS_L1_P/GNSS_POWERS.GPS_L1_new;
            P_CA = GNSS_POWERS.GPS_L1_CA/GNSS_POWERS.GPS_L1_new;
            P_C = GNSS_POWERS.GPS_L1_C/GNSS_POWERS.GPS_L1_new;

            S_CA = P_CA*spectrumgen_call('L1CA',F);
            S_P = P_P*spectrumgen_call('L1P',F);
            S_M = P_M*spectrumgen_call('L1M',F);
            [S_Cd, S_Cp] = spectrumgen_call('L1C',F);
            S_Cd = S_Cd*P_C;
            S_Cp = S_Cp*P_C;

            SI = S_P+S_M+S_Cd;
            SQ = S_CA+S_Cp;

        %-------------------- GPS L2  -------------------------------------
        case 'L2_new'

            P_P = GNSS_POWERS.GPS_L2_P_new/GNSS_POWERS.GPS_L2_new;
            P_C = GNSS_POWERS.GPS_L2_C/GNSS_POWERS.GPS_L2_new;
            P_M = GNSS_POWERS.GPS_L2_M/GNSS_POWERS.GPS_L2_new;

            S_P = P_P*spectrumgen_call('L2P',F);
            S_C = P_C*spectrumgen_call('L2C',F);
            S_M = P_M*spectrumgen_call('L2M',F);

            SI = S_P+S_M;
            SQ = S_C;

        case 'L2_new_2'

            P_P = GNSS_POWERS.GPS_L2_P_new/GNSS_POWERS.GPS_L2_new_2;
            P_C =  GNSS_POWERS.GPS_L2_C_new/GNSS_POWERS.GPS_L2_new_2;
            P_M =  GNSS_POWERS.GPS_L2_M/GNSS_POWERS.GPS_L2_new_2;

            S_P =   P_P*spectrumgen_call('L2P',F);
            S_C =   P_C*spectrumgen_call('L2C',F);
            S_M =   P_M*spectrumgen_call('L2M',F);

            SI = S_P+S_M;
            SQ = S_C;

        %-------------------- GPS L5  -------------------------------------
        case {'L5', 'L5_new'}
            SI = 0.5*spectrumgen_call('L5I',F);
            SQ = SI;

        %-------------------- Galileo E1 ----------------------------------
        case {'E1'}
            P_E1PRS = 0.5;
            P_E1OS = 0.5;

            S_E1PRS = P_E1PRS*spectrumgen_call('E1PRS',F);
            S_E1OS = P_E1OS*spectrumgen_call('E1OS',F);

            SI = S_E1OS;
            SQ = S_E1PRS;

        %-------------------- Galileo E6 ----------------------------------
        case 'E6'
            P_E6A = 0.5;
            P_E6B = 0.5;

            S_E6PRS = P_E6A*spectrumgen_call('E6PRS',F);
            S_E6CS = P_E6B*spectrumgen_call('E6CS',F);

            SI = S_E6CS;
            SQ = S_E6PRS;

        %-------------------- Galileo E5 ----------------------------------
       case 'E5'
            SI = 0.5*spectrum_AltBOC(15,10,F);
            SQ = SI;

        case {'E5A','E5B'}
            aux = spectrum_AltBOC(15,10,F);
            len = floor(length(aux)/2);
            aux2 = [aux(1:len+1) zeros(1,len)];
            len2 = floor(length(aux2)/2);
            [~, aux3] =  max(aux2);
            SI =  0.5*2*circshift(aux2',len2-aux3)';  % actually is not exactly multply by 2..
            SQ = SI;

        %-------------------- BeiDou-2 current ----------------------------
        % B1
        case {'B1','B11','B12'}
            SI = 0.5*spectrum_BPSK(2*1.023e6,F);
            SQ = 0.5*spectrum_BPSK(2*1.023e6,F);

        % B2
        case {'B2Q'}
            SI = spectrum_BPSK(10*1.023e6,F);
            SQ = zeros(size(SI));
        case {'B2I'}
            SI = spectrum_BPSK(2*1.023e6,F);
            SQ = zeros(size(SI));
        case{'B2'}
            P_B2Q = 0.5;
            P_B2I = 0.5;

            SQ = P_B2Q*spectrumgen_call('B2Q',F);
            SI = P_B2I*spectrumgen_call('B2I',F);

        % B3
        case{'B3'}
             SI = 0.5*spectrum_BPSK(10*1.023e6,F);
             SQ = SI;

        %-------------------- BeiDou-2 future -----------------------------
        % B1
        case {'B1Cd'}
            SI = spectrum_BOCs(1,1,F);
            SQ = zeros(size(SI));

        case {'B1Cp'}
            SI = spectrum_BOCs(6,1,F);
            SQ = zeros(size(SI));

        case {'B1C'}
            P_B1Cd = (10/11);
            P_B1Cp = (1/11);

            S_B1Cd = P_B1Cd*spectrumgen_call('B1Cd',F);
            S_B1Cp = P_B1Cp*spectrumgen_call('B1Cp',F);

            SI = S_B1Cd + S_B1Cp;
            SQ = zeros(size(SI));

        case {'B1_new'}
            SI = spectrum_BOCs(14,2,F);
            SQ = zeros(size(SI));

        case {'B1composite'}
            P_B1C = GNSS_POWERS.BEIDOU_B1C/GNSS_POWERS.BEIDOU_B1_composite;
            P_B1_new = GNSS_POWERS.BEIDOU_B1_new/GNSS_POWERS.BEIDOU_B1_composite;

            S_B1C = P_B1C*spectrumgen_call('B1C',F);
            S_B1_new = P_B1_new*spectrumgen_call('B1_new',F);

            SI = S_B1C+S_B1_new;
            SQ = zeros(size(SI));

        % B2
        case{'B2_new'}
            SI = 0.5*spectrum_AltBOC(15,10,F);
            SQ = 0.5*spectrum_AltBOC(15,10,F);

        % B3
        case{'B3_new'}
            SI = 0.5*spectrum_BPSK(10*1.023e6,F);
            SQ = 0.5*SI;
        case{'B3A'}
            S = spectrum_BOCs(15,2.5,F);
            SQ = zeros(size(SI));

        case {'B3composite'}
            P_B3_new = GNSS_POWERS.BEIDOU_B3_new/GNSS_POWERS.BEIDOU_B3_composite;
            P_B3A = GNSS_POWERS.BEIDOU_B3A/GNSS_POWERS.BEIDOU_B3_composite;

            S_B3_new = P_B3_new*spectrumgen_call('B3_new',F);
            S_B3A = P_B3A*spectrumgen_call('B3A',F);

            SI = S_B3_new;
            SQ = S_B3A;
    end
    end

function S = spectrum_BPSK(fc,F)
%==========================================================================
% BPSK of chipping rate fc.
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

    S = sinc((1/fc)*F);
    aux = sum(abs(S).^2)/length(S);
    S = S/sqrt(aux);
    S = abs(S).^2;
end

function S = spectrum_BOCs(n,m,F)
%==========================================================================
% sine-phased BOC even of a chipping rate fc=m*1.023e6 and sub-carrier rate
% of fs = n*1.023e6.
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

    fs = n*1.023e6;
    fc = m*1.023e6;

    S = sinc((1/fc)*F).*tan(pi*F/(2*fs));
    aux = sum(abs(S).^2)/length(S);
    S = S/sqrt(aux);
    S = abs(S).^2;
end

function S = spectrum_BOCc(n,m,F)
%==========================================================================
% cosine-phased BOC even of a chipping rate fc=m*1.023e6 and sub-carrier
% rate of fs = n*1.023e6.
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

    fs = n*1.023e6;
    fc = m*1.023e6;

    S = 2*sinc((1/fc)*F).*((sin(pi*F/(4*fs)).^2)./cos(pi*F/(2*fs)));
    aux = sum(abs(S).^2)/length(S);
    S = S/sqrt(aux);
    S = abs(S).^2;
end

function S = spectrum_AltBOC(n,m,F)
%==========================================================================
% sine-phased AltBOC even of a chipping rate fc=m*1.023e6 and sub-carrier
% rate of fs = n*1.023e6.
%--------------------------------------------------------------------------
% Author: Daniel Pascual (daniel.pascual [at] protonmail.com)
% Copyright 2017 Daniel Pascual
% License: GNU GPLv3
%==========================================================================

    fs = n*1.023e6;
    fc = m*1.023e6;

    S = (4*fc./((pi*F).^2)).* ((cos(pi*F/(fc))).^2).* (((cos(pi*F/(2*fs))).^2) - cos(pi*F/(2*fs)) - 2*cos(pi*F/(2*fs)).*cos(pi*F/(4*fs)) + 2 ) ./ ((cos(pi*F/(2*fs))).^2);
    S_ = S;
    S_(find(isnan(S))) = 0;
    aux = sum(abs(S_))/length(S);
    S = S/aux;
    S = abs(S);
end

author: danipascual

repo: GNSS-matlab


text:
load demo.mat
pitch =  -120/180*pi;
yaw = -30/180*pi;
pitchR = [1           0           0           ;...
          0           cos(pitch)  -sin(pitch) ;...
          0           sin(pitch)  cos(pitch)  ];
yawR =   [cos(yaw)    0           sin(yaw)  ;...
          0           1           0           ;...
          -sin(yaw)   0           cos(yaw)  ];

R = pitchR * yawR ;

Rt = [R [0.8;0;12]];

width = 640*1;
height= 480*1;
K = [  1100*1    0  width/2 ;...
         0  1100*1  height/2 ;...
         0    0    1];
[render,depth] = RenderPCcameraMatlab(double(pt),color, ptCamera,colorCamera ,K, Rt,width,height);
imshow(render);
imwrite(render,'demo.png');

author: jianxiongxiao

repo: ProfXkit


text:
function visualizeCellsRadial2D(phi)
%VISUALIZECELLS plots the values of cell variable phi
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

L = phi.domain.cellcenters.x(end);
x = [phi.domain.facecenters.x(1); phi.domain.cellcenters.x; phi.domain.facecenters.x(end)];
y = [phi.domain.facecenters.y(1); phi.domain.cellcenters.y; phi.domain.facecenters.y(end)];
[TH,R] = meshgrid(y, x);
[X,Y] = pol2cart(TH,R);
h = polar([0 2*pi], [0 L]);
delete(h);
hold on
pcolor(X,Y,phi.value)
colorbar
hold off

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
classdef TestPSNR
    %TestPSNR

    methods (Static)
        function test_1
            src1 = randi([0 255], [50 50], 'uint8');
            src2 = randi([0 255], [50 50], 'uint8');
            x = cv.PSNR(src1, src2);
        end

        function test_2
            src1 = randi([0 255], [50 50 3], 'uint8');
            src2 = randi([0 255], [50 50 3], 'uint8');
            x = cv.PSNR(src1, src2);
        end

        function test_error_argnum
            try
                cv.PSNR();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function opt=fmin_fvalfind_ppform(pp,x,y)

yf=ppval(pp,x);
opt=abs(yf-y)^2;

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function test06 ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_GL_TEST06 creates a sparse Gauss-Legendre grid and writes it to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06:\n' );
  fprintf ( 1, '  SPARSE_GRID_GL makes a sparse Gauss-Legendre grid.\n' );
  fprintf ( 1, '  Write the data to a set of quadrature files.\n' );

  level_min = max ( 0, level_max + 1 - dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MIN = %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
%
%  Determine the number of points.
%
  point_num = sparse_grid_gl_size ( dim_num, level_max );

  r(1:dim_num,1) = -1.0;
  r(1:dim_num,2) = +1.0;
%
%  Compute the weights and points.
%
  [ w, x ] = sparse_grid_gl ( dim_num, level_max, point_num );
%
%  Write the data out.
%
  r_filename = sprintf ( 'gl_d%d_level%d_r.txt', dim_num, level_max );
  w_filename = sprintf ( 'gl_d%d_level%d_w.txt', dim_num, level_max );
  x_filename = sprintf ( 'gl_d%d_level%d_x.txt', dim_num, level_max );

  r8mat_write ( r_filename, dim_num, 2,         r );
  r8mat_write ( w_filename, 1,       point_num, w );
  r8mat_write ( x_filename, dim_num, point_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R data written to "%s".\n', r_filename );
  fprintf ( 1, '  W data written to "%s".\n', w_filename );
  fprintf ( 1, '  X data written to "%s",\n', x_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [params] = sv_calcB(vResults)
% function [vResults] = sv_calcB(vResults)
% -------------------------------------
% Calculation of b-value uncertainty from given mValueGrid_ determined with
% sv_calcMc
%
%
% Input parameters:
% vResults
%
% Output parameters:
%
%
% J. Woessner; woessner@seismo.ifg.ethz.ch
% last update: 14.04.03

% Initialize

% Create Indices to catalog
[vResults.caNodeIndices] = ex_CreateIndexCatalog(vResults.mCatalog, vResults.mPolygon, vResults.bMap, vResults.nGriddingMode,...
    vResults.nNumberEvents, vResults.fRadius, vResults.fSizeRectHorizontal, vResults.fSizeRectDepth);

% Determine time period of catalog
vResults.fTminCat = min(vResults.mCatalog(:,3));
vResults.fTmaxCat = max(vResults.mCatalog(:,3));
% Adjust to decimal years
fTimePeriod =vResults.fTimePeriod/365;

% Init result matrix
mValueGrid_ = [];

% Loop over all grid nodes
for nNode_ = 1:length(vResults.mPolygon(:,1))
    nNode_
    % Create node catalog
    mNodeCatalog_ = vResults.mCatalog(vResults.caNodeIndices{nNode_}, :);
    % Select time period for mNodeCatalog_
    vSel = (vResults.fTstart <= mNodeCatalog_(:,3) & mNodeCatalog_(:,3) < vResults.fTstart+fTimePeriod);
    mNodeCatalog_ = mNodeCatalog_(vSel,:);
    % Check for constant number of events calculations
    if (vResults.nGriddingMode == 0)
        [mNodeCatalog_] = ex_MaxRadius(vResults.mCatalog, vResults.mPolygon, nNode_, vResults.caNodeIndices, vResults.fMaxRadius, vResults.nNumberEvents, vResults.bMap);
    end
    [nX,nY] = size(mNodeCatalog_);
    if (nX < vResults.nMinimumNumber)
        mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
    else
        try
            % b-value for median Mc(EMR)
            vSel1=(mNodeCatalog_(:,6) >= vResults.mValueGrid(nNode_,6));
            mCat1=mNodeCatalog_(vSel1,:);
            [fMeanMag1, fBValue1, fStdDev1, fAValue1] =  calc_bmemag(mCat1, vResults.fBinning);
            % b-value for 16 percentile Mc
            vSel2=(mNodeCatalog_(:,6) >= vResults.mValueGrid(nNode_,16));
            mCat2=mNodeCatalog_(vSel2,:);
            [fMeanMag2, fBValue2, fStdDev2, fAValue2] =  calc_bmemag(mCat2, vResults.fBinning);
            % b-value for 84 percentile Mc
            vSel3=(mNodeCatalog_(:,6) >= vResults.mValueGrid(nNode_,17));
            mCat3=mNodeCatalog_(vSel3,:);
            [fMeanMag3, fBValue3, fStdDev3, fAValue3] =  calc_bmemag(mCat3, vResults.fBinning);
            % b-value check for Mc(EMR)
            vSel4=(mNodeCatalog_(:,6) >= vResults.mValueGrid(nNode_,5));
            mCat4=mNodeCatalog_(vSel4,:);
            [fMeanMag4, fBValue4, fStdDev4, fAValue4] =  calc_bmemag(mCat4, vResults.fBinning);
            mValueGrid_= [mValueGrid_; fBValue1 fStdDev1 fAValue1 fBValue2 fStdDev2 fAValue2...
                    fBValue3 fStdDev3 fAValue3 fMeanMag4 fBValue4 fStdDev4 fAValue4];
        catch
            mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
        end
    end; % End of if on length(mNodeCatalog_(:,1)
end; % for nNode
save(['Bval_result_Time' num2str(vResults.fTstart) '_Nmin_' num2str(vResults.nMinimumNumber) '.mat'], 'mValueGrid_');

author: CelsoReyes


text:
% S = vgg_vec_swap(x,y)  For square matrix A and vectors x, y it is x'*A*y = vgg_vec(A)'*vgg_vec_swap(x,y).
%
% x ... matrix N-by-K
% y ... matrix N-by-K
% S ... matrix K^2-by-N
%
% Examples :-
%
%  - Estimating fundamental matrix F that should satisfy x(:,k)'*F*y(:,k) = 0 from
%    given points x(:,k), y(:,k). We solve homogeneous system vec(F)*vgg_vec_swap(x,y) = 0.
%
%  - Evaluating values v(k) of x(:,k)'*A{k}*y(:,k) for all k, A{k} are square matrices.
%    If column vectors vec(A{k}) are stacked to a single matrix AA, then
%        v = sum(AA'.*vgg_vec_swap(x,y)).
%
% See also vgg_vech_swap

% Tomas Werner, Oct 2001

function M = vgg_vec_swap(x,y)

if nargin==1
  y = x;
end

[i j] = find(ones(size(x,1)));
M = x(i,:).*y(j,:);

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function value = i4_choose ( n, k )

%*****************************************************************************80
%
%% I4_CHOOSE computes the binomial coefficient C(N,K).
%
%  Discussion:
%
%    The value is calculated in such a way as to avoid overflow and
%    roundoff.  The calculation is done in integer arithmetic.
%
%    The formula used is:
%
%      C(N,K) = N! / ( K! * (N-K)! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    ML Wolfson, HV Wright,
%    Algorithm 160:
%    Combinatorial of M Things Taken N at a Time,
%    Communications of the ACM,
%    Volume 6, Number 4, April 1963, page 161.
%
%  Parameters:
%
%    Input, integer N, K, are the values of N and K.
%
%    Output, integer VALUE, the number of combinations of N
%    things taken K at a time.
%
  mn = min ( k, n - k );

  if ( mn < 0 )

    value = 0;

  elseif ( mn == 0 )

    value = 1;

  else

    mx = max ( k, n - k );
    value = mx + 1;

    for i = 2 : mn
      value = ( value * ( mx + i ) ) / i;
    end

  end

  return
end


text:
function owens_test01 ( )

%*****************************************************************************80
%
%% TEST01 demonstrates the use of T.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  T evaluates Owen''s T function.\n' );
  fprintf ( 1, '  Compare to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          H            A        ' );
  fprintf ( 1, 'T                         T\n' );
  fprintf ( 1, '                                ' );
  fprintf ( 1, '(Tabulated)               (T)                     DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, h, a, t1 ] = owen_values ( n_data );

    if ( n_data == 0 )
      break
    end

    t2 = t ( h, a );

    fprintf ( 1, '  %12.8f  %12.8f  %24.16e  %24.16e  %10.4e\n', ...
    h, a, t1, t2, abs ( t1 - t2 ) );

  end

  return
end


text:
close all
clear all
path(path,'..\..\..\FUZZCLUST')
%the data
load motorcycle.txt
data.X = motorcycle(:,[1 2]);

[N,n]=size(data.X);

%data normalization
data = clust_normalize(data,'range');
plot(data.X(:,1),data.X(:,2),'.')
hold on
%parameters
param.c=4;
param.vis=1;
param.val=1;
%clustering
result=kmeans(data,param);
%validation
result = validity(result,data,param);
result.validity

author: Sable

repo: mcbench-benchmarks


text:
% INVWISHRND - Inverse Wishart Distribution - Random Matrix Value
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
%   [ IW ] = invwishrnd(S,d)
%
% S = p x p symmetric, postitive definite "scale" matrix
% d = "degrees of freedom" parameter (integer)
%   = "precision" parameter (d may be non-integer)
%
% IW = random matrix from the inverse Wishart distribution
%
% Note:
%   Different sources use different parameterizations.
%   This routine uses that of Press and Shigemasu (1989):
%   density (IW) is proportional to
%     exp[-.5*trace(S*inv(IW))] / [det(IW) ^ (d/2)].
%
%   With this density definition:
%   mean of IW = S/(d-2p-2) for d>2p+2,
%   mode of IW = S/d.
%
% See also: INVWISHIRND, WISHRND

function [IW] = invwishrnd(S,d)

[p,p2] = size(S) ;

W = wishrnd(inv(S),d-p-1) ;

IW = inv(W) ;


text:
function [traj, infStates] = tapas_hgf_ar1_mab(r, p, varargin)
% Calculates the trajectories of the agent's representations under the AR(1)-HGF in a multi-armed
% bandit task
%
% This function can be called in two ways:
%
% (1) tapas_hgf_ar1_mab(r, p)
%
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_hgf_ar1_mab(r, ptrans, 'trans')
%
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_hgf_ar1_mab_transp(r, p);
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = length(p)/6;

    if l ~= floor(l)
        error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
    end
end

% Number of bandits
try
    b = r.c_prc.n_bandits;
catch
    error('tapas:hgf:NumOfBanditsConfig', 'Number of bandits has to be configured in r.c_prc.n_bandits.');
end

% Unpack parameters
mu_0 = p(1:l);
sa_0 = p(l+1:2*l);
phi  = p(2*l+1:3*l);
m    = p(3*l+1:4*l);
ka   = p(4*l+1:5*l-1);
om   = p(5*l:6*l-2);
th   = exp(p(6*l-1));
al   = p(6*l);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];
y = [0; r.u(:,2)];

% Number of trials (including prior)
n = size(u,1);

% Construct time axis
if r.c_prc.irregular_intervals
    if size(u,2) > 1
        t = [0; r.u(:,end)];
    else
        error('tapas:hgf:InputSingleColumn', 'Input matrix must contain more than one column if irregular_intervals is set to true.');
    end
else
    t = ones(n,1);
end

% Initialize updated quantities

% Representations
mu = NaN(n,l,b);
pi = NaN(n,l,b);

% Other quantities
muhat = NaN(n,l,b);
pihat = NaN(n,l,b);
v     = NaN(n,l);
w     = NaN(n,l-1);
da    = NaN(n,l);
dau   = NaN(n,1);

% Representation priors
% Note: first entries of the other quantities remain
% NaN because they are undefined and are thrown away
% at the end; their presence simply leads to consistent
% trial indices.
mu(1,:,:) = repmat(mu_0,[1 1 b]);
pi(1,:,:) = repmat(1./sa_0,[1 1 b]);

% Representation update loop
% Pass through trials
for k = 2:1:n
    if not(ismember(k-1, r.ign))

        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%

        % 1st level
        % ~~~~~~~~~
        % Prediction
        muhat(k,1,:) = mu(k-1,1,:) +t(k) *phi(1) *(m(1) -mu(k-1,1,:));

        % Precision of prediction
        pihat(k,1,:) = 1/(1/pi(k-1,1,:) +t(k) *exp(ka(1) *mu(k-1,2,:) +om(1)));

        % Input prediction error
        dau(k) = u(k) -muhat(k,1,y(k));

        % Updates
        pi(k,1,:) = pihat(k,1,:);
        pi(k,1,y(k)) = pi(k,1,y(k)) +1/al;

        mu(k,1,:) = muhat(k,1,:);
        mu(k,1,y(k)) = mu(k,1,y(k)) +1/pihat(k,1,y(k)) *1/(1/pihat(k,1,y(k)) +al) *dau(k);

        % Volatility prediction error
        da(k,1) = (1/pi(k,1,y(k)) +(mu(k,1,y(k)) -muhat(k,1,y(k)))^2) *pihat(k,1,y(k)) -1;

        if l > 2
            % Pass through higher levels
            % ~~~~~~~~~~~~~~~~~~~~~~~~~~
            for j = 2:l-1
                % Prediction
                muhat(k,j,:) = mu(k-1,j,:) +t(k) *phi(j) *(m(j) -mu(k-1,j,:));

                % Precision of prediction
                pihat(k,j,:) = 1/(1/pi(k-1,j,:) +t(k) *exp(ka(j) *mu(k-1,j+1,:) +om(j)));

                % Weighting factor
                v(k,j-1) = t(k) *exp(ka(j-1) *mu(k-1,j,y(k)) +om(j-1));
                w(k,j-1) = v(k,j-1) *pihat(k,j-1,y(k));

                % Updates
                pi(k,j,:) = pihat(k,j,:) +1/2 *ka(j-1)^2 *w(k,j-1) *(w(k,j-1) +(2 *w(k,j-1) -1) *da(k,j-1));

                if pi(k,j,1) <= 0
                    error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
                end

                mu(k,j,:) = muhat(k,j,:) +1/2 *1/pi(k,j,:) *ka(j-1) *w(k,j-1) *da(k,j-1);

                % Volatility prediction error
                da(k,j) = (1/pi(k,j,y(k)) +(mu(k,j,y(k)) -muhat(k,j,y(k)))^2) *pihat(k,j,y(k)) -1;
            end
        end

        % Last level
        % ~~~~~~~~~~
        % Prediction
        muhat(k,l,:) = mu(k-1,l,:) +t(k) *phi(l) *(m(l) -mu(k-1,l,:));

        % Precision of prediction
        pihat(k,l,:) = 1/(1/pi(k-1,l,:) +t(k) *th);

        % Weighting factor
        v(k,l)   = t(k) *th;
        v(k,l-1) = t(k) *exp(ka(l-1) *mu(k-1,l,y(k)) +om(l-1));
        w(k,l-1) = v(k,l-1) *pihat(k,l-1,y(k));

        % Updates
        pi(k,l,:) = pihat(k,l,:) +1/2 *ka(l-1)^2 *w(k,l-1) *(w(k,l-1) +(2 *w(k,l-1) -1) *da(k,l-1));

        if pi(k,l,1) <= 0
            error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
        end

        mu(k,l,:) = muhat(k,l,:) +1/2 *1/pi(k,l,:) *ka(l-1) *w(k,l-1) *da(k,l-1);

        % Volatility prediction error
        da(k,l) = (1/pi(k,l,y(k)) +(mu(k,l,y(k)) -muhat(k,l,y(k)))^2) *pihat(k,l,y(k)) -1;
    else

        mu(k,:,:) = mu(k-1,:,:);
        pi(k,:,:) = pi(k-1,:,:);

        muhat(k,:,:) = muhat(k-1,:,:);
        pihat(k,:,:) = pihat(k-1,:,:);

        v(k,:)  = v(k-1,:);
        w(k,:)  = w(k-1,:);
        da(k,:) = da(k-1,:);

    end
end

% Remove representation priors
mu(1,:,:)  = [];
pi(1,:,:)  = [];

% Check validity of trajectories
if any(isnan(mu(:))) || any(isnan(pi(:)))
    error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
else
    % Check for implausible jumps in trajectories
    dmu = diff(mu);
    dpi = diff(pi);
    rmdmu = repmat(sqrt(mean(dmu.^2)),length(dmu),1);
    rmdpi = repmat(sqrt(mean(dpi.^2)),length(dpi),1);

    jumpTol = 256;
    if any(abs(dmu(:)) > jumpTol*rmdmu(:)) || any(abs(dpi(:)) > jumpTol*rmdpi(:))
        error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
    end
end

% Remove other dummy initial values
muhat(1,:,:) = [];
pihat(1,:,:) = [];
v(1,:)       = [];
w(1,:)       = [];
da(1,:)      = [];
dau(1)       = [];

% Create result data structure
traj = struct;

traj.mu     = mu;
traj.sa     = 1./pi;

traj.muhat  = muhat;
traj.sahat  = 1./pihat;

traj.v      = v;
traj.w      = w;
traj.da     = da;
traj.dau    = dau;

% Updates with respect to prediction
traj.ud = mu -muhat;

% Psi (precision weights on prediction errors)
psi      = NaN(n-1,l);
pi1      = squeeze(pi(:,1,:));
pi1obs   = pi1(sub2ind(size(pi1), (1:size(pi1,1))', y));
psi(:,1) = 1./(al*pi1obs);
for i=2:l
    pihati    = squeeze(pihat(:,i-1,:));
    pihatiobs = pihati(sub2ind(size(pihati), (1:size(pihati,1))', y));
    pii       = squeeze(pi(:,i,:));
    piiobs    = pii(sub2ind(size(pii), (1:size(pii,1))', y));
    psi(:,i)  = pihatiobs./piiobs;
end
traj.psi = psi;

% Epsilons (precision-weighted prediction errors)
epsi        = NaN(n-1,l);
epsi(:,1)   = psi(:,1) .*dau;
epsi(:,2:l) = psi(:,2:l) .*da(:,1:l-1);
traj.epsi   = epsi;

% Full learning rate (full weights on prediction errors)
wt        = NaN(n-1,l);
wt(:,1)   = psi(:,1);
wt(:,2:l) = 1/2 *(v(:,1:l-1) *diag(ka(1:l-1))) .*psi(:,2:l);
traj.wt   = wt;

% Create matrices for use by the observation model
infStates = NaN(n-1,l,b,4);
infStates(:,:,:,1) = traj.muhat;
infStates(:,:,:,2) = traj.sahat;
infStates(:,:,:,3) = traj.mu;
infStates(:,:,:,4) = traj.sa;

return;

author: translationalneuromodeling

repo: tapas


text:
function [loopcheck, pc_, pp_, pk_, sdc, sdp, sdk] = ploop_c_and_p_calcs(MIN_CSTEP, MIN_PSTEP, trackLoop, stdevcall)
    % loop for calculating c and p values
    % if cstep is empty, then neither cstep or err1 will be calculated
    % sdevcall will contain 'kpc', 'kp', or 'kp~' to determine how kcp_stdevs behaves
    % -  'kpc' will calculate the k, p, and c errors
    % - 'kp' calculates only the k & p errors
    % - 'kp~' also calculates the k & p errors, but in a 3x matrix
    %   where the results for c are ignored
    % all these options exist because that is how they were
    % implemented by the ploop functions

    % modified omori law:
    % n(t) = K /(t+c)^p
    % n(t) : frequency of aftershocks per unit time interval,
    % K : the productivity of the sequence,
    % c : adjusts for missing earthquakes in the catalog
    % p : how quickly the activity falls off to the constant background intensity.
    % typically normaized to days
    %
    % ploop parts attributed to A.Allmann and B. Enescu
    % routines deconstructed & merged by C. Reyes 2017

    % TODO: get rid of all these globals!!@!!Q#RQR@!#R$@!
    global pc pp pk
    global tt loop nn
    global nit t
    global err1x err2x  % previous errors
    global ieflag isflag %error small enough, step small enough
    global qp
    global err1 err2
    global cstep pstep
    global ts eps1 eps2
    global pcheck
    %global sdk sdp sdc

    loopcheck=0;

    % pp = initial p = 1.1
    % nn = length of catalog (# of events)
    % tt = end time (last event)
    % ts = start time (1st event)
    % pc = initial c???
    % t  = time of each event (a vector)
    while loopcheck < 500
        if ~isnumeric(t)
            %fprintf('t is a %s',class(t));
            % normalize all times to days after first event
            t = days(t - ts);
            tt = days(tt - ts);
            ts = days(ts - ts);
            %fprintf(' now t is a %s\n',class(t) );
        end

        % preverr was never assigned anywhere
        preverr = 0.0;

        if pp==1.0
            pp=1.001;
        end

        if trackLoop
            loop=loop+1;
        end
        loopcheck=loopcheck+1;
        nit=nit+1;

        qp=1-pp;

        pk=(qp*nn)/((tt+pc)^qp-(ts+pc)^qp);

        update_cof_and_cog(pk, qp, ts, pc)

        if ~isempty(MIN_CSTEP) % kludgy! but still reduces lines tremendously
            err1 = c_err(pk, pc, t, ts, tt, pp);
            ieflag = ieflag || abs(err1)< eps1 ;
            isflag = isflag || cstep <=MIN_CSTEP;
        end
        err2 = p_err(t, pc, pk, qp, tt, ts);

        ieflag = ieflag || ( abs(err2) <  eps2); % errors small enough
        isflag = isflag || ( pstep <= MIN_PSTEP); % steps small enough

        %stop searching if errors or steps are small enough
        if ieflag || isflag
            switch stdevcall
                case 'kpc'
                    [sdk, sdp, sdc] = kcp_stdevs();
                case 'kp'
                    [sdk, sdp] = kcp_stdevs();
                case 'kp~'
                    [sdk, sdp, ~] = kcp_stdevs();
            end
            pcheck = 1;
            break
        else
            if trackLoop
                [cstep, pc, loop] = take_cstep(MIN_CSTEP, cstep, preverr, ts, pc, nit, err1, loop);
            else
                [cstep, pc, ~] = take_cstep(MIN_CSTEP, cstep, preverr, ts, pc, nit, err1, []);
            end
            [pstep, pp] = take_pstep(MIN_PSTEP, pp, pstep, preverr, nit, err2);
        end

        if ~isempty(MIN_CSTEP)
            err1x=err1;
        end
        err2x=err2;

    end

    %p= pp;
    pp_=pp;
    pc_=pc;
    pk_=pk;

end

%%
function cerr = c_err(pk, pc, t, ts, tt, pp)
    %c_err calculate c error
    qsum=pk*((1/(tt+pc)^pp)-(1/(ts+pc)^pp));
    psum=sum(1./(t+pc));
    cerr=qsum+pp*psum;
end

function perr = p_err(t, pc, pk, qp, tt, ts)
    %p_err calculate p error
    sumln=sum(log(t+pc));
    qsumln=pk/qp^2;
    qsumln=qsumln*(((tt+pc)^qp)*(1-qp*log(tt+pc))-((ts+pc)^qp)*(1-qp*log(ts+pc)));
    esumln=qsumln+sumln;
    perr=esumln;
end

function [cstep, pc, loop] = take_cstep(MIN_CSTEP, cstep, preverr, ts, pc, nit, err1, loop)
    if isempty(MIN_CSTEP); return; end

    STEP_REDUCTIONFACTOR = 0.9;

    if nit>1
        %if error has changed sign,reduce the step size
        if ~same_sign(preverr,err1)  &&  cstep >= MIN_CSTEP
            cstep=cstep * STEP_REDUCTIONFACTOR;
        end
    end

    pc = pc - cstep * sign(err1); % move closer to zero

    dt=ts+pc;
    if dt <= 0
        pc=ts+cstep;
    end
    if ~isempty(loop) && loop == 30
        cstep=cstep * STEP_REDUCTIONFACTOR;
        loop=0;
    end

end

function [pstep, pp] = take_pstep(MIN_PSTEP, pp, pstep, preverr, nit, err2)
    % calculate the parameters of p-value
    if isempty(MIN_PSTEP); return; end
    STEP_REDUCTIONFACTOR = 0.9;

    if nit>1
        %if error has changed sign,reduce the step size
        if ~same_sign(preverr,err2) &&  pstep>= MIN_PSTEP
            pstep=pstep * STEP_REDUCTIONFACTOR;
        end
    end

    pp = pp - pstep * sign(err2); % move closer to zero
end

function [sdk, sdp, sdc] = kcp_stdevs()
    % calculate standard deviations for k, c, and p values
    %kcp_stdevs.m                          A.Allmann
    %
    % calculat the parameters of p-value
    %calls itself with different parameters for different loops in programm

    % TODO: get rid of all these globals!!@!!Q#RQR@!#R$@!

    global p c tt pc pp
    global pk ts dk

    te=tt;
    p=pp;
    c=pc;
    dk=pk;

    %case1
    f1=((te+c)^(-p+1))/(-p+1);
    h1=((ts+c)^(-p+1))/(-p+1);
    s(1)=(1/dk)*(f1-h1);

    %case2
    f2=((te+c)^(-p));
    h2=((ts+c)^(-p));
    s(2)=f2-h2;

    %case3

    f3=(-(te+c)^(-p+1))*(((log(te+c))/(-p+1))-(1/((-p+1)^2)));
    h3=(-(ts+c)^(-p+1))*(((log(ts+c))/(-p+1))-(1/((-p+1)^2)));
    s(3)=f3-h3;

    %case4

    s(4)=s(2);

    %case5

    f5=((te+c)^(-p-1))/(p+1);
    h5=((ts+c)^(-p-1))/(p+1);
    s(5)=(-dk)*(p^2)*(f5-h5);

    %case6

    f6=((te+c)^(-p))*(((log(te+c))/(-p))-(1/(p^2)));
    h6=((ts+c)^(-p))*(((log(ts+c))/(-p))-(1/(p^2)));
    s(6)=(dk*p)*(f6-h6);

    %case7

    s(7)=s(3);

    %case8

    s(8)=s(6);

    %case9

    f10=((te+c)^(-p+1))*((log(te+c))^2)/(-p+1);
    f11=(2*((te+c)^(-p+1)))/((-p+1)^2);
    f12=(log(te+c))-(1/(-p+1));
    f9=f10-(f11*f12);

    h10=((ts+c)^(-p+1))*((log(ts+c))^2)/(-p+1);
    h11=(2*((ts+c)^(-p+1)))/((-p+1)^2);
    h12=(log(ts+c))-(1/(-p+1));
    h9=h10-(h11*h12);
    s(9)=(dk)*(f9-h9);

    %assign the values of s to the matrix A(i,j)
    %invert the matrix to calculate the standard deviation
    %for k,c,p .

    if nargout == 3
        A=[s(1) s(2) s(3); s(4) s(5) s(6); s(7) s(8) s(9)];

        A=inv(A);

        sdk=sqrt(A(1,1));
        sdc=sqrt(A(2,2));
        sdp=sqrt(A(3,3));

    elseif nargout == 2
        A=[s(1) s(3); s(3) s(9)];
        A=inv(A);
        sdk=sqrt(A(1,1));
        sdp=sqrt(A(2,2));
    else
        error('wrong number of output arguments');
    end

end


text:
% DEFINE RAY
function [rayObj] = ray(origin,direction,magnitude)
% Define a ray with a given a origin and direction.
assert(numel(origin) == 3 && numel(direction) == 3,'Parameters must be 3D');
if nargin < 3
    magnitude = 1;
end
% NORMALISE THE DIRECTION VECTOR
direction = direction/norm(direction);
% DEFINE THE RAY STRUCTURE
rayObj = struct('origin',origin,...
                'direction',direction,...
                'magnitude',magnitude);
end

repo: OpenMAS


text:
function CODEWORD = RS_ENC(INFO,n,k,g,field)

%CODEWORD = RS_ENC(INFO,n,k,g,field)
%
% m  is the number of bits of each symbol
% n = 2^m-1 => the number of symbols transmitted
% k = the number of code symbols that is going to be codes to a n symbol message
% t = the number of errors that can be found + corrected

%Tripple-error-correcting Reed-Solomon code with symbols from GF(2^4)
% Lin & Costello p.175 and article: Reed_Solomon Codes by Joel Sylvester

%generator polynomial

%field = gftuple([-1:2^m-2]', m, 2);

%p = 2; m = 4;
%primpoly = [0 0 -Inf -Inf 0];
%field = gftuple([-1:p^m-2]',primpoly,p);

%Lin + Costello, p.171

%Encoder (Article)
%shift codeword by X^(n-k)
for ii = 1:n-k
    shiftpol(ii) = -Inf;
end
%shiftpol(n-k+1) = 0;
shiftcode = [shiftpol INFO];

%divide shifted codeword by g(x)
[Q, R] = GFDECONV(shiftcode, g, field);

while length(R) < n-k
    R = [R -inf];
end

CODEWORD = [R INFO];

%CODWORD = gfconv(CODEWORD,0,field);

for i =1:n
    if CODEWORD(i) == -1
        CODEWORD(i) = -Inf;
    end
end


text:
function [Sd] = zdeblank (S)
%
%   [Sd] = zdeblank (s)
%
%Purpose:
%   removes blanks surrounding a word in S
%
%
%Input Parameters:
%   S : a string
%
%
%Output Parameters:
%  Sd : S without the surrounding blanks
%
%
%Key Terms:
%
%More Info :
%   deblank, zdeblankall
%
%
%
%     Author : Ziad Saad
%     Date : Tue Sep 12 12:15:24 PDT 2000
%     LBC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'zdeblank';

blank_chars=sprintf(' \t\r\n\f\v%s',char(0));

non_blank_mask=~any(bsxfun(@eq,blank_chars,S(:)),2);
first_pos=find(non_blank_mask,1,'first');
last_pos=find(non_blank_mask,1,'last');

if isempty(first_pos)
    assert(isempty(last_pos));
    Sd='';
else
    Sd=S(first_pos:last_pos);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%SF_LINE_P5 1D Fifth order Lagrange shape functions for lines (P5).
%
%   [ VBASE, NLDOF, XLDOF, SFUN ] = SF_LINE_P5( I_EVAL, N_SDIM, N_VERT, I_DOF, XI, AINVJAC, VBASE )
%   Evaluates conforming fifth order P5 Lagrange shape functions on 1D line elements
%   with values defined in the nodes and center. XI are Barycentric coordinates.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       i_eval      scalar:  1             Evaluate function values
%                           >1             Evaluate values of derivatives
%       n_sdim      scalar: 1              Number of space dimensions
%       n_vert      scalar: 2              Number of vertices per cell
%       i_dof       scalar: 1-6            Local basis function to evaluate
%       xi          array [2,1]            Local coordinates of evaluation point
%       aInvJac     [n,3]                  Inverse of transformation Jacobian
%       vBase       [n]                    Preallocated output vector
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       vBase       [n]                    Evaluated function values
%       nLDof       [4]                    Number of local degrees of freedom on
%                                          vertices, edges, faces, and cell interiors
%       xLDof       [2,n_ldof]             Local coordinates of local dofs
%       sfun        string                 Function name of called shape function
%
%   See also SF_LINE_P1

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [len, labels] = imLength(img, varargin)
% Compute total length of a binary 1D structure.
%
%   LEN = imLength(IMG);
%   When IMG is a binary image, returns the total length of the structure,
%   equal to the number of pixels.
%   This function is intended to be used for debugging purpose.
%
%   LEN = imLength(IMG, RESOL);
%   Compute the length in user unit, by multiplying length in pixel by the
%   resolution.
%
%   [LEN, LABELS] = imLength(LBL, ...);
%   When LBL is a label image, returns the total length of each label
%   different from 0. Returns also the set of unique values in LBL.
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-10-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY

% check image dimension
if ndims(img)>2 || min(size(img))>1 %#ok<ISMAT>
    error('first argument should be a 1D image');
end

% in case of a label image, return a vector with a set of results
labels = 1;
if ~islogical(img)
    labels = unique(img);
    labels(labels==0) = [];

    len = zeros(length(labels), 1);
    for i=1:length(labels)
        len(i) = imLength(img==labels(i), varargin{:});
    end
    return;
end

% extract resolution if present
resol = 1;
if ~isempty(varargin)
    resol = varargin{1};
end

% Total length of stucture is equal to the number of vertices multiplied by
% the resolution
len = sum(img(:)) * resol;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function cmap = winterCmap(numGrays,numColors)
%
% cmap = winterCmap(numGrays,numColors)
%
% Makes colormap array with:
%   gray scale - 1:numGrays
%   winter colors - numGrays+1:numGrays+numColors
%
% djh 1/98
if ~exist('numGrays','var')
  numGrays=128;
end
if ~exist('numColors','var')
  numColors=96;
end

cmap = zeros(numGrays+numColors,3);
cmap(1:numGrays+numColors,:) = [gray(numGrays); winter(numColors)];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
clear classes

W = diag([0.1, 1*pi/180].^2);
P0 = diag([0.005, 0.005, 0.001].^2);
V = diag([0.005, 0.5*pi/180].^2);

randinit
map = Map(20, 10);
veh = Vehicle(V);
RandomPath(veh, map.dim);
sensor = RangeBearingSensor(veh, map, W);
sensor.interval = 5;
ekf = EKF(veh, W, P0, sensor, W);

randinit
ekf.run(1000);

f1
clf
map.visualize()
veh.plot_xy('b');
ekf.plot_xy('r');
ekf.plot_ellipse([], 'k')
grid on
xyzlabel

f2
clf
ekf.plot_P()

author: petercorke

repo: robotics-toolbox-matlab


text:
function [dfdx, fx] = VBA_numericDiff(fName, idxArg2Diff, varargin)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [dfdx, fx] = numericDiff(fname, idxArg2Diff, fArg1, fArg2, ...)
% numerical evaluation of derivatives
%
% This function evaluates numerically the derivatives of the function
% 'fname' with respect to its 'idxArg2Diff' input argument at the point
% defined by the input arguments [fArg1, fArg2, ..., fArgN]
%
% IN:
%   - fname: handle or name of the function to be differentiated.
%     ie either its name or a function handle. The latter can be very useful
%     if the function is itself a subfunction of a function!
%   - idxArg2Diff: the index of the argument to differentiate the function
%     with (1 <= idxArg2Diff <= numel(varargin))
%   - fArg1, fArg2, ..., fArgN: list of arguments which is required to call
%     the function 'fname', and which defines the ordinate at which the
%     derivative will be numerically evaluated
%
% OUT:
%   - dfdx: an m x p matrix containing the numerical differentiation
%     of the function evaluated at {fArg1, ..., fArgN}, where p is the
%     number of elements of the function's first output and m is the number
%     of elements of the idxArg2Diff input argument. If the function output
%     or the differentiated argument are in matrix form, they will be
%     vectorized first.
%   - f: the function evaluation at x
%
% NB: Mixed partials can be obtained by recursive call of this routine,
% e.g. consider:
%
% dfdx_idx_j = numericDiff(@numericDiff,i+2,fname,j,arg1,arg2,...,argn)
% dfdx_idx_j = reshape(dfdx_idx_j,mi,mj,p)
%
% The first line evaluates the numerical derivative of the function
% numericDiff(fname,j,arg1,arg2,...,argn) wrt to its (i+2)th entry, which
% is arg_i.
% The second line reshapes the output such that the first dimension is the
% dimension of argi (mi), the second is of the dimension of argj (mj), and
% the last dimension is the one of the function output itself (p).
%
% /////////////////////////////////////////////////////////////////////////

% perturbation scale
epsilon = 1e-4;

% shortcut
fArgs = varargin;

% evaluate function at the specified argument
% =========================================================================
try
    fx = VBA_vec(fName(fArgs{:}));
catch
    message = sprintf(...
        '*** VBA_numericDiff: Can not call %s with the provided arguments (nArgs = %d).', ...
        func2str(fName), numel(fArgs));
    error(message);
end

% evaluate the perturbation df of the function fname in the neighbourhood
% of the specified argument (ie at x + dx) and calculate the derivative wrt x
% =========================================================================

% pre-allocate the variables
m = numel (fArgs{idxArg2Diff});
p = numel (fx);
dfdx = zeros (m, p);

% loop over the dimension of idxArg2Diff-th argument
dx = epsilon * fArgs{idxArg2Diff};
dx (abs(dx) <= eps) = epsilon;

% compute effect of pertubations
for i = 1 : m
    xpdx = fArgs;
    xpdx{idxArg2Diff}(i) = xpdx{idxArg2Diff}(i) + dx(i);
    dfdx(i,:) = (VBA_vec (fName (xpdx{:})) - fx)' / dx(i);
end

author: MBB-team

repo: VBA-toolbox


text:
filename='Arch_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function dcmFolder2MATobject(varargin)

% function dcmFolder2MATobject(PathName,MaxVarSize,reOrderOpt,dicomDictFactory,fileExtension)
% ------------------------------------------------------------------------
% This function converts the DICOM files in the folder PathName to a MAT
% object which is stored as a file called IMDAT in a new subfolder called
% IMDAT.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% Change log
% 2018/06/11 Suppressed/removed triggertime related warning (this appears
% to be a custom field which is nearly always missing).
%------------------------------------------------------------------------

%% PARSE INPUT

switch nargin
    case 1
        PathName=varargin{1};
        MaxVarSize=[];
        reOrderOpt=0;
        dicomDictFactory=0;
        fileExtension='.dcm';
        optionStruct=[];
    case 2
        PathName=varargin{1};
        MaxVarSize=varargin{2};
        reOrderOpt=0;
        dicomDictFactory=0;
        fileExtension='.dcm';
        optionStruct=[];
    case 3
        PathName=varargin{1};
        MaxVarSize=varargin{2};
        reOrderOpt=varargin{3};
        dicomDictFactory=0;
        fileExtension='.dcm';
        optionStruct=[];
    case 4
        PathName=varargin{1};
        MaxVarSize=varargin{2};
        reOrderOpt=varargin{3};
        dicomDictFactory=varargin{4};
        fileExtension='.dcm';
    case 5
        PathName=varargin{1};
        MaxVarSize=varargin{2};
        reOrderOpt=varargin{3};
        dicomDictFactory=varargin{4};
        fileExtension=varargin{5};
        optionStruct=[];
    case 6
        PathName=varargin{1};
        MaxVarSize=varargin{2};
        reOrderOpt=varargin{3};
        dicomDictFactory=varargin{4};
        fileExtension=varargin{5};
        optionStruct=varargin{6};
end

defaultOptionStruct.ignoreDynamic=0;
defaultOptionStruct.nDownSample=[1 1 1];
defaultOptionStruct.skipInfo=0;
%Fix option structure, complete and remove empty values
[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1);
nDownSample=optionStruct.nDownSample;
skipInfo=optionStruct.skipInfo;

%%
%Get/set maximum variable size (sets save steps in case the variable size
%is exceeded)
if isempty(MaxVarSize)
    MaxVarSize=1e12; %Default
end

if isempty(reOrderOpt)
    reOrderOpt=0; %Default
end

if isempty(dicomDictFactory)
    dicomDictFactory=0; %Default
end

if isempty(fileExtension)
%     fileExtension='.dcm'; %Default
end

%% Getting DICOM file names and image dimension parameters
if ~isempty(fileExtension)
    files = dir(fullfile(PathName,['*',fileExtension]));
else %No file extension
    files = dir(fullfile(PathName));

    %Remove files and hidden folders containing dot in path
    logicKeep=false(1,numel(files));
    for q=1:1:numel(files)
        if isempty(strfind(files(q).name,'.')) && isempty(strfind(files(q).name,'IMDAT'))
            logicKeep(q)=1;
        end
    end
    files=files(logicKeep);
end
files={files(1:end).name};
NumberOfFiles=numel(files);

%% Attempt to sort files
% This sorting assumes the files names make sense in terms of slice
% numbers.

try %Try to get numbers from image file names
    nameNumbers=cellfun(@str2double,regexp(files,'\d+\d*','match'));
    [~,indSort]=sort(nameNumbers); %Get sort order
    files=files(indSort); %Sort file list based on numeric data in names
catch %Normal sort instead
    files=sort(files(:)); %Attempt sort (might cause dcm1,dcm10,...)
end

%%

firstWarning_TriggerTime=0;
firstWarning_Rotation=0;
if NumberOfFiles>0
    %% SPECIFY Tags to collect

    collectTags = { 'Basename', ...
        'PatientID', ...
        'StudyInstanceUID', ...
        'SeriesInstanceUID', ...
        'StudyID', ...
        'InstanceNumber', ...
        'SeriesNumber', ...
        'SeriesDescription', ...
        'StudyDescription', ...
        'ImageType', ...
        'Rows', ...
        'Columns', ...
        'PixelSpacing',...
        'ImagePositionPatient', ...
        'ImageOrientationPatient',...
        'SliceLocation', ...
        'SliceThickness', ...
        'SpacingBetweenSlices', ...
        'SliceGap', ...
        'NumberOfTemporalPositions', ...
        'TemporalPositionIdentifier', ...
        'HeartRate', ...
        'TriggerTime', ...
        'RepetitionTime', ...
        'EchoTime', ...
        'FlipAngle', ...
        'RescaleSlope', ...
        'ScaleSlope',...
        'RescaleIntercept',...
        'RescaleType',...
        'Height',...
        'Width',...
        'MRImagePhaseNumber'...
        'SliceNumberMR'...
        'NumberOfPCDirections',...
        'NumberOfPhasesMR',...
        'NumberOfSlicesMR',...
        'MRSeriesNrOfSlices',...
        'MRSeriesNrOfPhases'...
        'MRSeriesNrOfDynamicScans',...
        'MRImagePhaseNumber',...
        'ImagePlaneNumber'...
        'MRScaleSlope'...
        'MRScaleIntercept'...
        'AcquisitionMatrix'...
        'PixelBandwidth'...
        'ImagingFrequency'...
        'MagneticFieldStrength'...
        'InPlanePhaseEncodingDirection'...
        'PhaseEncodingDirection'...
        };
    % 'EchoNumbers' is private and just a echotime identifier or sequence number
    % 'SliceNumberMR' is also private (2001,100a) and can be usefull for identifing slices that are on same position
    % N.B. Filename is added by matlab including the path, Basename is one we add ourself

    %% Creating output file name and folder
    fName='IMDAT';
    foldername_out=fullfile(PathName,fName,filesep);
    if ~exist(foldername_out,'dir') %create output folder if it does not exist already
        mkdir(foldername_out); %Create new dir
    end
    savename_out1=fullfile(foldername_out,'IMDAT.mat');
    if exist(savename_out1,'file') %Check if file exist already
        delete(savename_out1); %Remove file
    end
    matObj = matfile(savename_out1,'Writable',true);

    %% SETTING DICOM DICTIONARY
    fName=fullfile(PathName,files{1});

    %First import using factory settings
    dicomdict('factory');
    dcmInfo_full=dicominfo(fName);
    m=dicomread(fName); %Test upload first
    sizSlice=size(m);
    dictSetting=3;

    if dicomDictFactory==0
        try
            if ~isempty(strfind(lower(dcmInfo_full.Manufacturer),'philips'))
                disp(['Detected ',dcmInfo_full.Manufacturer,' files']);
                dicomdict('set','gibbon_dict.txt');
                disp('DICOM dictionary set to: gibbon_dict.txt');
                dictSetting=1;
            elseif ~isempty(strfind(lower(dcmInfo_full.Manufacturer),'pms'))
                disp(['Detected ',dcmInfo_full.Manufacturer,' files']);
                dicomdict('set','gibbon_dict.txt');
                disp('DICOM dictionary set to: gibbon_dict.txt');
                dictSetting=1;
            elseif ~isempty(strfind(lower(dcmInfo_full.Manufacturer),'siemens'))
                disp(['Detected ',dcmInfo_full.Manufacturer,' files']);
                dicomdict('set','dicom-dict-siemens.txt');
                disp('DICOM dictionary set to: dicom-dict-siemens.txt');
                dictSetting=2;
            elseif ~isempty(strfind(lower(dcmInfo_full.Manufacturer),'GE MEDICAL SYSTEMS'))
                disp(['Detected ',dcmInfo_full.Manufacturer,' files']);
                warning('No settings for this vendor using DICOM dictionary factory settings');
                dictSetting=3;
            else
                dicomdict('factory');
                warning('Unknown vendor, using DICOM dictionary factory settings');
                dictSetting=3;
            end
            %Test info import with new dictionary (if this fails we resume with
            % factory
            dcmInfo_full=dicominfo(fName);
        catch %e.g. if the manufacturer field is missing
            dicomdict('factory');
            warning('Unknown vendor, using DICOM dictionary factory settings');
            dictSetting=3;
            dcmInfo_full=dicominfo(fName);
        end
    elseif dicomDictFactory==1
        dicomdict('factory');
        dcmInfo_full=dicominfo(fName);
        dictSetting=3;
    end

    %% LOADING DICOM INFO

    if skipInfo==1
        numInfoLoad=1; %Only first will be loaded
    else
        numInfoLoad=numel(files); %All will be loaded
    end

    hw = waitbar(0,'Loading DICOM info...');
    for c=1:1:numInfoLoad
        fName=fullfile(PathName,files{c});
        %dicomdict('get');
        dcmInfo_full=dicominfo(fName);
        dcmInfo_full.Basename=fName; %Add custom field
        iFields = find(isfield(dcmInfo_full,collectTags)); %indices of existing fields
        for iField=1:numel(iFields)
            fieldName = collectTags{iFields(iField)};
            dcmInfo(c).(fieldName) = dcmInfo_full.(fieldName); %allocated after first iteration
        end
        if c==1
            dcmInfo=repmat(dcmInfo,1,numel(files)); %Allocate after size of first is known
        end
        waitbar(c/numInfoLoad,hw,['Loading DICIM info...',num2str(round(100.*c/numel(files))),'%']);
    end
    waitbar(c/numInfoLoad,hw,['Saving DICOM info to MAT-file...',num2str(round(100.*c/numel(files))),'%']);
    % matObj.dcmInfo=dcmInfo;
    close(hw);

    %% Reorder files based on InstanceNumber (file naming may deviate from slice order)
    if reOrderOpt
        [~,indSort]=sort([dcmInfo(:).InstanceNumber]);
        files=files(indSort);
        dcmInfo=dcmInfo(indSort);
    end

    %%

    %Geometry information based on first info file
    [G]=dicom3Dpar(dcmInfo(1));
    G.v=G.v(:).*nDownSample(:);

    matObj.G=G; %Saving geometry information

    %TODO this part is not memory protected yet, all DICOM info is loaded into
    %one variable

    % hw = waitbar(0,'Loading DICOM info...');
    % step_size=100;
    % q=1;
    % ind_save=[];
    % for c=1:1:numel(files)
    %     ind_save=[ind_save c];
    %     dcmInfo(q)=dicominfo(fullfile(foldername,files{c}));
    %     waitbar(c/numel(files),hw,['Loading DICIM info...',num2str(round(100.*c/numel(files))),'%']);
    %     if numel(ind_save)>=step_size
    %         waitbar(c/numel(files),hw,['Saving to MAT-file...',num2str(round(100.*c/numel(files))),'%']);
    %         matObj.dcmInfo(1,ind_save)=dcmInfo;
    %         ind_save=[];
    %         dcmInfo=[];
    %         q=1;
    %     end
    %     q=q+1;
    % end
    % matObj.dcmInfo(1,ind_save)=dcmInfo;
    % close(hw);

    %% LOADING IMAGE DATA

    if isfield(dcmInfo,'EchoTime')
        EchoTimesAll=[dcmInfo(:).EchoTime];
    else
        EchoTimesAll=nan;
    end
    EchoTimesUni=unique(EchoTimesAll);

    matObj.EchoTimesUni=EchoTimesUni;
    NumEchoTimes=numel(EchoTimesUni);

    %Get image types
    if isfield(dcmInfo,'ImageType')
        ImageTypesAll={dcmInfo.ImageType};
        checkType=1;
    else
        ImageTypesAll=repmat({'Unknown'},1,numel(files));
        checkType=0;
        warning('ImageType information missing, using Unknown as type, assuming single type');
    end

    ImageTypesUni=unique(ImageTypesAll);
    matObj.ImageTypesUni=ImageTypesUni;
    NumImageTypes=numel(ImageTypesUni);

    if isfield(dcmInfo,'TriggerTime') %Multiple dynamics
        TriggerTimesAll=[dcmInfo(:).TriggerTime];
        TriggerTimesUni=unique(TriggerTimesAll);
        NumberOfTemporalPositions=numel(TriggerTimesUni); %NumberOfPhasesMR
    else
        TriggerTimesAll=zeros(1,NumberOfFiles);
        TriggerTimesUni=unique(TriggerTimesAll);
        if isfield(dcmInfo,'NumberOfTemporalPositions') %Multiple dynamics
            NumberOfTemporalPositions=double(dcmInfo(1).NumberOfTemporalPositions);
        else
            NumberOfTemporalPositions=1;
        end
    end

    if optionStruct.ignoreDynamic==1
        NumberOfTemporalPositions=1;
    end

    NumberOfSlices=NumberOfFiles/NumImageTypes/NumEchoTimes/NumberOfTemporalPositions; %NumberOfSlicesMR
%     NumberOfFilesPerType=NumberOfSlices*NumberOfTemporalPositions;

    switch dictSetting
        case 1 %PHILIPS
            %             NumberOfRows=double(dcmInfo(1).Width);
            %             NumberOfColumns=double(dcmInfo(1).Height);
            NumberOfRows=double(dcmInfo(1).Rows);
            NumberOfColumns=double(dcmInfo(1).Columns);
        case 2 % SIEMENS
            NumberOfRows=double(dcmInfo(1).Rows);
            NumberOfColumns=double(dcmInfo(1).Columns);
        case 3 %FACTORY
            NumberOfRows=double(dcmInfo(1).Rows);
            NumberOfColumns=double(dcmInfo(1).Columns);
            %             NumberOfRows=double(dcmInfo(1).Width);
            %             NumberOfColumns=double(dcmInfo(1).Height);
    end

    %Determine image size
    ImageSize=[NumberOfRows NumberOfColumns NumberOfSlices NumberOfTemporalPositions];

    rowSelectSet=1:nDownSample(1):NumberOfRows;
    NumberOfRowsKeep=numel(rowSelectSet);
    columnSelectSet=1:nDownSample(2):NumberOfColumns;
    NumberOfColumnsKeep=numel(columnSelectSet);
    sliceSelectSet=1:nDownSample(3):NumberOfSlices;
    NumberOfSlicesKeep=numel(sliceSelectSet);

    ImageSizeKeep=[NumberOfRowsKeep NumberOfColumnsKeep NumberOfSlicesKeep NumberOfTemporalPositions];

%     %Check if slice size matching DICOM parameters
%     if sizSlice(1)~=ImageSize(1)
%         ImageSize(1)=sizSlice(1);
%         warning('Image size 2 does not match size specified in DICOM info. Using image size from image instead')
%     end
%     if sizSlice(2)~=ImageSize(2)
%         ImageSize(2)=sizSlice(2);
%         warning('Image size 2 does not match size specified in DICOM info. Using image size from image instead')
%     end

    matObj.ImageSize=ImageSizeKeep;

    NumClass=['uint',num2str(dcmInfo_full.BitsAllocated)];

    c=0;
    hw = waitbar(0,'Loading DICOM info...');

    for iEcho=1:NumEchoTimes
        %Finding files for current EchoTime
        EchoTimeNow=EchoTimesUni(iEcho); %The current echo time

        L_Echo=EchoTimesAll==EchoTimeNow;

        %String to add to type spec
        if NumEchoTimes==1
            echoNameAppend=[];
        else
            echoNameAppend=['_EchoTime_',num2str(iEcho)];
        end

        for iType=1:NumImageTypes

            %Finding files for current type
            if checkType==1
                ImageTypeNow=ImageTypesUni(iType); %The current image type
                L_Type=strcmpi(ImageTypesAll,ImageTypeNow);
            else
                L_Type=true(size(ImageTypesAll));
            end

            %Fix L_Echo in case EchoTIme is not defined
            if isnan(EchoTimesAll)
                L_Echo=true(size(L_Type));
            end

            L_now=L_Echo&L_Type;

            TypeFiles=files(L_now);
%             iStart=((iType-1)*NumberOfFilesPerType)+1; %     iEnd=iStart-1+NumberOfFilesPerType;

            TypeName=['type_',num2str(iType),echoNameAppend];
            TypeNameDcmInfo=[TypeName,'_info'];

            VarSize=(prod(ImageSizeKeep))*dcmInfo_full.BitsAllocated/8;
            if VarSize>MaxVarSize %Save to MAT-file for each temporal step
                if ImageSizeKeep(4)>1
                    matObj.(TypeName)=zeros([ImageSizeKeep(1:3) 2],NumClass);
                else
                    matObj.(TypeName)=zeros(ImageSizeKeep,NumClass);
                end
                for iTemp=1:1:NumberOfTemporalPositions
                    %Finding files for current temporal position
                    TriggerTimesType=TriggerTimesAll(L_now);%[dcmInfo(L_now).TriggerTime];
                    TriggerTime=TriggerTimesUni(iTemp);
                    L_Temp=(TriggerTimesType==TriggerTime);
                    %Creating volume by assembling slices
                    TempTypeFiles=TypeFiles(L_Temp);
                    M=zeros(ImageSizeKeep(1:3),NumClass);
                    sliceIndexResample=1;
                    for iSlice=1:nDownSample(3):numel(TempTypeFiles)

                        waitbar(c/numel(files),hw,['Loading DICIM image data...',num2str(round(100.*c/numel(files))),'%']);
                        load_name=fullfile(PathName,TempTypeFiles{iSlice});
                        m=dicomread(load_name);
                        m=m(1:nDownSample(1):end,1:nDownSample(2):end,:); %Resample
                        if size(m,3)>1
                            m=mean(double(m),3);
                            warning('Multi-dimensional (e.g. RGB) slices where converted to 2D grayscale');
                        end
                        M(:,:,sliceIndexResample)=m;

                        sliceIndexResample=sliceIndexResample+1;
                        if iSlice==1
                            c=c+1;
                        else
                            c=c+nDownSample(3);
                        end
                    end
                    waitbar(c/numel(files),hw,['Saving DICOM image data to MAT-file...',num2str(round(100.*c/numel(files))),'%']);
                    eval(['matObj.',TypeName,'(:,:,:,iTemp)=M;']);
                end
            else %create 4D array and save whole array to MAT-file
                N=zeros(ImageSizeKeep,NumClass);
                for iTemp=1:1:NumberOfTemporalPositions
                    %Finding files for current temporal position
                    TriggerTimesType=TriggerTimesAll(L_now);%[dcmInfo(L_now).TriggerTime];
                    TriggerTime=TriggerTimesUni(iTemp);
                    L_Temp=(TriggerTimesType==TriggerTime);

                    %Creating volume by assembling slices
                    TempTypeFiles=TypeFiles(L_Temp);
                    M=zeros(ImageSizeKeep(1:3),NumClass);
                    sliceIndexResample=1;
                    for iSlice=1:nDownSample(3):numel(TempTypeFiles)
                        waitbar(c/numel(files),hw,['Loading DICIM image data...',num2str(round(100.*c/numel(files))),'%']);
                        load_name=fullfile(PathName,TempTypeFiles{iSlice});
                        m=dicomread(load_name); %Current slice
                        m=m(1:nDownSample(1):end,1:nDownSample(2):end,:); %Resample
                        if size(m,3)>1
                            m=mean(double(m),3);
                            warning('Multi-dimensional (e.g. RGB) slices where converted to 2D grayscale');
                        end
                        if size(M,1)~=size(m,1)
                            M(:,:,sliceIndexResample)=m'; %Try transpose
                            if firstWarning_Rotation==0
                                warning('Image data was rotated to match expected image size!');
                                firstWarning_Rotation=1;
                            end
                        else
                            M(:,:,sliceIndexResample)=m;
                        end
                        sliceIndexResample=sliceIndexResample+1;
                        if iSlice==1
                            c=c+1;
                        else
                            c=c+nDownSample(3);
                        end
                    end
                    N(:,:,:,iTemp)=M;
                end
                waitbar(c/numel(files),hw,['Saving DICOM image data to MAT-file...',num2str(round(100.*c/numel(files))),'%']);
                eval(['matObj.',TypeName,'=N;']);
                eval(['matObj.',TypeNameDcmInfo,'=dcmInfo(L_now);']);

            end
        end
    end
    close(hw);
else
    warning(['Folder ',PathName,' skipped since it does not contain .dcm files']);
end

dicomdict('factory');

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.


text:
function [relamps,Y,E,avgTC] = eventBoldContrast(tc,selectedConds,avgTC,unitFree,plotFlag)
% Calculate and plot the relative time series amplitudes for various
% experimental conditions
%
%   [relamps,Y,E,avgTC] = eventBoldContrast(tc,selectedConds,[avgTC],[unitFree],[plotFlag]);
%
% returns:
% relamps: the contrast values for every event
% Y: the average contrast value per condition
% E: standard error per condition
% avgTC: the mean time course of all the events in all conditions (r, baseline
% for contrast computation)
%
% The model is based on Ress & Heeger (Nature Neuroscience, 2000).  The
% relative amplitudes are  (after our recent edits) probably a misnomer.
% The values returned here are the percent contrast of the modulation.
%
% The model used here is this.  We assume that the basic time course of
% each condition has the same shape function, say r(t). They differ only by
% noise and a scalar, a*r(t) + N. We estimate the basic shape of the time
% course by averaging all of the conditions. This reduces the effect of the
% noise (N).  Then we compare the amplitude of the response by looking at
% the scalar (dot, inner) product of the measurements for one condition
% with respect to the estimated average.
%
% At present, the calculation of the relative amplitude is
%
%       a =  s . r / sqrt(r.r)
%
% where r is the average set to a mean of zero, and s is the time course
% for the selected condition, also set to a mean of zero.
%
% The logic for this is as follows.  The angle, theta, between s and r is
%
%     cos(theta) = s.r / |s||r|
%
% The length of the projection of s onto r is |x|/|s| = cos(theta).
% Putting these together,
%
%      |x| = |s|* (s.r) / |s||r|
%      |x| = (s.r)/|r|
%      |x| = (s.r)/sqrt(r.r)
%
% The |r| is the lenth of r and this is sqrt(dot(r,r,))
% Units, scaling: taking the sqrt in the denominator has the effect of keeping the
% original scale of the subject's signal. so, the units are %modulation now.
%
% unitFree flag: if this is set to 1, we do not take the square root in the denominator.
% this leaves us with a unit less measure, but one that is more comparable
% across subjects, which can be used for a group analysis, if overall signal
% variations are considered irrelevant.
%
% Limitations.
% If different conditions have different time courses, this too will be
% missed by this analysis.
%
% tc: a structure containing the time course information for an
% experiment.
%
% HISTORY
% 02/23/04 ras: wrote it. This is the first new functionality to
% timeCourseUI since incorporating it into mrLoadRet.
% 07/04 ras: now everything uses er_chopTSeries, so the relamps
% are already calculated. Nonetheless, I'm keeping the old calculations
% here, since with relamps in particular, having some conditions not
% selected may significantly change the mean time course, fundamentally
% changing the usefulness of the estimate. (e.g., if there's a baseline
% condition with a fundamentally diff't shape than other conds, disabling
% it, then recomputing makes sense.)
% 08/05 MBS, BW Comments. also changed the way relAmps are computed - avgTC
% is normalized to mean zero and unit size
%

if ieNotDefined('tc'), error('Time course structure required'); end
if ieNotDefined('selectedConds'), selectedConds = find(tc_selectedConds(tc)); end
if ieNotDefined('avgTC'),    avgTC = []; end
if ieNotDefined('unitFree'), unitFree = 0; end
if ieNotDefined('plotFlag'), plotFlag = 1; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up existing objects in figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if plotFlag
    otherAxes = findobj('Type','axes','Parent',gcf);
    delete(otherAxes);
    otherUiControls = findobj('Type','uicontrol','Parent',gcf);
    delete(otherUiControls);
    cla
end

% params
prestim = 1;

% tx contains the frame # in tc.wholeTc at which each trial starts
tx = round([(tc.trials.onsetSecs./tc.TR)+1 length(tc.wholeTc)]);
% intervals = diff(tx);

% figure out selected conditions, trials, labels, and colors

whichConds = tc.condNums(selectedConds);
nConds = length(whichConds);
trials = find(ismember(tc.trials.cond,whichConds));
nTrials = length(trials);

colors = tc.condColors(selectedConds);
labels = tc.condNames(selectedConds);

% initialize selTrials matrix to be NaNs of the size of the shortest trial
% minInt = min(intervals(trials)) + prestim;
% selTrials = NaN*ones(minInt,nTrials);

% get selected trials, avg time course of selected trials
selTrials = [];
for i = 1:nTrials
    rng = tx(trials(i))-prestim:tx(trials(i)+1)-1;
    rng = round(rng);
    rng = rng(rng>0 & rng<length(tc.wholeTc)); % clip to fit data
    selTrials(1:length(rng),i) = tc.wholeTc(rng)';
    conds(i) = tc.trials.cond(trials(i));
end

% normalize everything to have a mean zero
for i = 1:nTrials
    offset = mean(selTrials(:,i));
    selTrials(:,i) = selTrials(:,i) - offset;
end

% we calculate the avgTC after normalizing each trial, so its mean is set to zero too
if isempty(avgTC), avgTC = mean(selTrials,2);
else
    if length(avgTC) ~= size(selTrials,1)
        error('User supplied avgTC has length %.0f.  Should have length %.0f\n',...
            length(avgTC),size(selTrials,1));
    end
end

% take dot product of each trial's time course w/ avg
relamps = NaN*ones(1,nTrials);
for i = 1:nTrials
    relamps(i) = dot(avgTC,selTrials(:,i));
end

if unitFree % this will give Ress' measure, no units
    relamps = relamps / dot(avgTC,avgTC);
else
    % this will give a bold contrast measure
    relamps = relamps / sqrt(dot(avgTC,avgTC));
end

% get mean and stdev of each condition
for c = 1:length(whichConds)
    ind = (conds==whichConds(c));
	Y(c) = mean(relamps(ind));
	E(c) = std(relamps(ind)) ./ sqrt(sum(ind));
end

%%%%% plot color bars for each cond
% (gum to keep legend from being erased:)
if plotFlag

    % figure out nice axis range
    maxY = max([0 1.2*max(Y)]);
    minY = min([0 1.2*min(Y)]);
    AX = [0 nConds+1 minY maxY];

    set(gcf,'NextPlot','add');
    set(gca,'Visible','off');
    subplot('Position',get(gca,'Position'))
    mybar(Y,E,labels,labels,colors);
    ylabel('Relative amplitude');
    hold on

    % append calculated relamps to tc struct
    tc.relamps = relamps;
    if isfield(tc,'ui'),  set(gcf,'UserData',tc); end
end

return


text:
% AssignmentToIndex Convert assignment to index.
%
%   I = AssignmentToIndex(A, D) converts an assignment, A, over variables
%   with cardinality D to an index into the .val vector for a factor.
%   If A is a matrix then the function converts each row of A to an index.
%
%   See also IndexToAssignment.m and FactorTutorial.m

function I = AssignmentToIndex(A, D)

D = D(:)'; % ensure that D is a row vector
if (any(size(A) == 1)),
    I = cumprod([1, D(1:end - 1)]) * (A(:) - 1) + 1;
else
    I = sum(repmat(cumprod([1, D(1:end - 1)]), size(A, 1), 1) .* (A - 1), 2) + 1;
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% GETTEXT - This function prints a dialog box on screen and waits for
%             the user to enter a string. There is a cancel button which
%             returns a value of [].
% Usage:
%   >> out = gettext(label1,label2,...,label7);
%
% Author: Colin Humphries, CNL / Salk Institute, La Jolla, 1997

% Copyright (C) 1997 Colin Humphries, Salk Institute
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad

function out = gettext(label1,label2,label3,label4,label5,label6,label7);

global is_text_entered
is_text_entered = 0;

for i = 1:nargin
   eval(['leng(i) = length(label',num2str(i),');'])
end
figurelength = 15*1.0*max(leng);

if figurelength < 240
   figurelength = 290;
end

figureheight = 80+nargin*18;
% Set figure

figure('Position',[302 396 figurelength figureheight],'color',[.6 .6 .6],'NumberTitle','off')
f1 = gcf;

ax = axes('Units','Pixels','Position',[0 0 figurelength figureheight],'Visible','off','XLim',[0 figurelength],'YLim',[0 figureheight]);

for i = 1:nargin

   eval(['text(figurelength/2,',num2str(figureheight-(i-1)*18-20),',label',num2str(i),',''color'',''k'',''FontSize'',16,''HorizontalAlignment'',''center'')'])

end
% Set up uicontrols

TIMESTRING = ['global is_text_entered;is_text_entered = 1;clear is_text_entered'];

u = uicontrol('Style','Edit','Units','Pixels','Position',[15 20 figurelength-95 25],'HorizontalAlignment','left','Callback',TIMESTRING);

TIMESTRING = ['global is_text_entered;is_text_entered = 2;clear is_text_entered'];
v = uicontrol('Style','Pushbutton','Units','Pixels','Position',[figurelength-70 20 55 25],'String','Cancel','Callback',TIMESTRING);

while(is_text_entered == 0)
   drawnow
end

if is_text_entered == 1
   out = get(u,'string');
else
   out = [];
end

clear is_text_entered
delete(f1)

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function x = nonzeros(a)
%NONZEROS     Implements  nonzeros(a)  for sparse interval matrix
%
%   x = nonzeros(a)
%
%Functionality as in Matlab.
%

% written  08/09/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/25/07     S.M. Rump  huge indices for sparse matrices
% modified 02/18/08     S.M. Rump  no column vector output for ND arrays
% modified 10/18/08     S.M. Rump  huge arrays
% modified 05/31/09     S.M. Rump  multiple arrays
% modified 10/03/12     S.M. Rump  spones for ND-arrays
% modified 11/13/12     S.M. Rump  cure Matlab 7.0 bug
%

  if a.complex
    if isequal(a.rad,0)
      x = cintval(nonzeros(a.mid));
    else
      [m,n] = size(a.mid);          % n is product of remaining indices for multiple arrays
      if length(n)>1                % cure bug in Matlab 7.0 (R14)
        a.mid = reshape(a.mid,m,n);   % make sure dimension fits with sparse(...)
        a.rad = reshape(a.rad,m,n);   % a.rad cannot be 0
      end
      [I,J] = find(spones(a.mid)+spones(a.rad));    % correct index set
      xmid = complex( real(nonzeros(real(a.mid)+sparse(I,J,sqrt(-1),m,n))) , ...
                      real(nonzeros(imag(a.mid)+sparse(I,J,sqrt(-1),m,n))) );
      xrad = real( real(nonzeros(a.rad+sparse(I,J,sqrt(-1),m,n))) );
      x = intval(xmid,xrad,'midrad');
      if size(x.mid,1)==1
        x.mid = x.mid.';
        x.rad = x.rad';
      end
    end
  else
    [m,n] = size(a.inf);
    if prod(size(a.inf))<2^31       % not huge
      if issparse(a.inf)
        x = nonzeros(a.inf(:)+sqrt(-1)*a.sup(:));
      else
        x = nonzeros(complex(a.inf(:),a.sup(:)));
      end
      x = intval(real(x),imag(x),'infsup');
      return
      I = find(spones(a.inf)+spones(a.sup));
      x = intval(squeeze(a.inf(I)),squeeze(a.sup(I)),'infsup');
      sizexinf = size(x.inf);
      if ( length(sizexinf)<3) & ( size(x.inf,1)==1 )
        x.inf = x.inf';
        x.sup = x.sup';
      end
    else                            % huge array
      [I,J,asup] = find(a.sup);
      x = nonzeros(a.inf + sparse(I,J,complex(0,asup),m,n));
      x = intval(real(x),imag(x),'infsup');
    end
  end


text:
function Distance = CalDistance(PopObj,RefPoint)
% Calculate the distance between each solution to each adjusted reference
% point

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

    N  = size(PopObj,1);
    NR = size(RefPoint,1);

    %% Adjust the location of each reference point
    index=find(sum(PopObj.^2,2)==0);
    if ~isempty(index)
        PopObj(index,:)=PopObj(index,:)+1e-06*rand(size(PopObj(index,:)));
    end
    Cosine = 1 - pdist2(PopObj,RefPoint,'cosine');%The cosine of the individual and the reference point in the population, N*NR
    NormR  = sqrt(sum(RefPoint.^2,2)); %NR*1
    NormP  = sqrt(sum(PopObj.^2,2)); %N*1
    d1     = repmat(NormP,1,NR).*Cosine; %N*NR
    d2     = repmat(NormP,1,NR).*sqrt(1-Cosine.^2); %N*NR
    [~,nearest] = min(d2,[],1); %1*NR
    RefPoint    = RefPoint.*repmat(d1(N.*(0:NR-1)+nearest)'./NormR,1,size(RefPoint,2));

    %% Calculate the distance between each solution to each point
    Distance = pdist2(PopObj,RefPoint);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function r = cot(a)
%COT          Hessian (elementwise) cotangent
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = cot(a.x);
    f = -1 - sqr(r.x);
    r.dx = f * a.dx;
    r.hx = f * ( a.hx - reshape( (r.x*a.dx) * a.dx.' , size(a.hx) ) );

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = cot(full(a.x));
    if issparse(a.hx)               % input sparse

      ax = (-1) - sqr(full(r.x(:)));
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = -r.x(:);
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,r.dx);
      end
      [ia,ja,sa] = find(a.hx);      % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = cot(a.x);
      rx = r.x(:).';
      f = -1 - sqr(rx);
      f = f(ones(N*N,1),:);
      r.dx = a.dx .* f(1:N,:);
      adx = repmat(rx,N,1) .* a.dx;
      r.hx = f .* ( a.hx - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:) );

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end


text:
function f = simplify(f, tol)
%SIMPLIFY  Simplifies the smoothPart of a SINGFUN object.
%  F = SIMPLIFY(F) returns a SINGFUN which has a simplified smoothPart by
%  calling SIMPLIFY in SMOOTHFUN and has exponents less than 1 by absorbing
%  boundary roots with integer multiplicity into the smoothPart.
%
%  F = SIMPLIFY(F, TOL) does the same thing but uses TOL supplied by the user
%  and pass it to SIMPLIFY in SMOOTHFUN.

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

if ( nargin == 2 )
    f.smoothPart = simplify(f.smoothPart, tol);
elseif ( nargin == 1 )
    f.smoothPart = simplify(f.smoothPart);
end

% Replace the roots:
f = simplifyExponents(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [NUM,file,fileExt] = calculateNumber(imgPath)
    imgExt = {'*.bmp', '*.jpg', '*.png'};
    k=1;
    d1 = dir([imgPath char(imgExt(k))]);
    file = {d1(~[d1.isdir]).name};
    if isempty(file)
        k = k + 1;
        d1 = dir([imgPath char(imgExt(k))]);
        file = {d1(~[d1.isdir]).name};
        if isempty(file)
            k = k + 1;
            d1 = dir([imgPath char(imgExt(k))]);
            file = {d1(~[d1.isdir]).name};
            NUM = length(file);
        else
            NUM = length(file);
        end
    else
        NUM = length(file);
    end
    fileExt = char(imgExt(k));
end

repo: Saliency-Evaluation-Toolbox