text:
function [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (dec)
% by Sundar Krishnan
% 2003, Edited in June, 2004
%
% Description :
% This function Fr_dec2bin.m will convert a POSITIVE Decimal system
% Fraction (dec) to Binary system Fraction Fr_bin.
% Matlab itself has bin2dec.m and dec2bin.m, but there seems to be
% no standard Matlab function when fractions are involved.
%
% This function Fr_bin2dec.m and it's companion / dual function Fr_dec2bin.m
% were developed mainly with a view to get  quick results
% while learning Arithmetic (Entropy) Coding in School.
% (Now, more comments have been added to better explain the programme.)
%
% The results of this function are limited in accuracy due to the
% "precision" used in the function num2str.m in addition to
% Floating Point limits and Rounding errors.
%
% Accumulation of errors due to these limits can be seen
% when Fr_bin2dec and Fr_dec2bin are tested back-to-back in pairs.
%
% After experiments, I observed that the best precision is 16.
% If all the digits of the input bin are used for a pure fraction,
% the results are likely to be more accurate since we have more margin
% wrt the limit of 16 digits.
%
% Given below under "Usage Eg" are the many cases
% that have been tested during the development of this program,
% together with the results obtained in each case.
%
% Pl do forward me any new case that breaks the code
% beyond the aforesaid limitations.
%
% Outputs str_Fr and Fr_dec are intermediate results.
%
% See also : [Fr_dec, str_Fr, Fr_bin] = Fr_bin2dec (bin)
%
% Additional Test Cases involving pairs of dual tests
% are given towards the end.
%
%                                   ********************
%
% Usage Eg : (The foll have been tried out.)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.6796875)         % 0.1010111
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (113.6796875)       % 1110001.1010111
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (113.68359374)
%   = 1110001.10101110111111111
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (1045.013671875)
%   = 10000010101.000000111
%
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.013671875)       % 0.000000111
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000131835937)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (10099300.131835937)
%   = 100110100001101001100100.0010000111000000
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (1.0450137e+018)
%
% Also try this ! and enjoy the result :
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (1.0450137e+100)
%
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (2987.120089)
%   % = 101110101011.0001111010111110
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (1167892987.120089)
%   % = 1000101100111001010000111111011.0001111010111110
%
%                                       &&&&&&&&&&&&
%
% Usage Eg : Check in pairs :
% Fr_dec = Fr_bin2dec (10000010100.0010000111)  % = 1.044125000000000e+003
% Fr_bin = Fr_dec2bin (1.044125000000000e+003)  % = 10000010100.001
% Fr_bin = Fr_dec2bin ( Fr_bin2dec (10000010100.0010000111) )
%   returns Fr_bin = 10000010100.001
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec (101110101011.00011111) )
%   returns Fr_bin = 101110101011.0001 (corr to 2987.0625)
%   instead of the expected (same) 101110101011.00011111
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ...
%          (1000101100111001010000111111011.0001111010111110) )
%   returns Fr_bin = 1000101100111001000000000000000.00000000000000000
%   (corr to 1167884288)
%   instead of the expected (same)
%   1000101100111001010000111111011.0001111010111110
%   which itself was obtained with Fr_dec2bin (1167892987.120089)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec (101110101011.0001111010111110) )
%   returns Fr_bin = 101110101011.0001 (corr to 2987.0625)
%   instead of the expected (same) 101110101011.0001111010111110
%   which itself was obtained with Fr_dec2bin (2987.120089)
%
%
%                                   ********************

% 1) Inits :
Fr_bin = 0 ;
exp_power  = 0 ;

%                                       &&&&&&&&&&&&

% 2) Use num2str to convert the input to string :
%
% After experiments, I observed that the best precision is 16.
% For eg, with precision >= 17,
% str_Fr = num2str ( .1010111, 17 )  =  0.10101110000000001
% str_Fr = num2str ( .1010111, 16 )  =  0.1010111
%
% num2str.m's output will also contain "0" prefix before the decimal dot "."
% which we remove later.

% 2-a) Check if the input is greater than 1.
% If yes, can we use higher precision ?
% NO, I have found problems with precision > 16 even when the input  > 1 !
% So, commenting out the foll code, and retaining precision = 16 only.
% str_Fr = num2str (dec) ;
% if str_Fr > 1
%     precision = 48 ;
% else
%     precision = 16 ;
% end

precision = 16 ;            % See the note above.
str_Fr = num2str (dec, precision) ;
% Some egs of dec = 1045.0137 , 1.0450137e+018 , 0.131835937 , 0.0000131835937

% NOTE : For long input dec strings, pl note that even with precision > 16,
% say, with precision = 48, the input itself is accurately read
% only for the first 16 digits ; or, if it is converted to an exp format,
% then the input is accurately read only till 15 decimals after the dot.
% For eg, if dec = 116789292349873465787.120089,
% the whole integer part is taken as :
% 116789292349873470000 = % 1.1678929234987347e+016
% So, this will by itself creep in errors !

% In general, it is observed errors will creep in
% if the whole integer part > 999999999999999

%                                       &&&&&&&&&&&&

% 3) Now, if str_Fr above is in exp format, as for eg,
% '2.987062500000000e+003', we would like to get it in the form = 2987.0625
%
% I have observed that if the input no < 0.0001 (ie, < 0.0001000...)
% num2str.m's output is in the exp form ie, with powers less than e-005.
% For eg, dec = 0.000100000001 gives str_Fr = 0.000100000001
% But dec = 0.0000999999999999 gives str_Fr = 9.9999999999900001e-005
%
% Also, with precision = 16, num2str.m's output for nos > 1, upto 1.0e+015,
% is WITHOUT the exp form of power. For eg,
% with dec = 999999999999999.9999999999999999
% str_Fr = num2str (dec, 16)    % gives = 1.0e+015 = 1 0000 0000 0000 000
%
% For nos > 1e+016, num2str.m's output is in exp form.

if ~isempty ( findstr ( str_Fr, 'e') )

    exp_power = 0 ;
    [str_Fr_Bef_Exp, exp] = strtok ( str_Fr, 'e' ) ;
        % Some egs = str_Fr_Bef_Exp = 1.119996810555458,   exp = e-005

    [exp_power, ign ] = strtok ( exp, 'e' ) ;
        % exp starts with 'e', hence see LHS

    exp_power = abs ( str2num (exp_power) ) ;

    % Remove the dot at the 2nd place : (as in 1.119996810555458)
    % However, there is no dot when it's a pure fraction,
    % and is an exact submultiple of 2 !
    if length (str_Fr_Bef_Exp) >= 2  ;
        str_Fr_Bef_Exp (2) = [] ;
    end

    if exp (2) == '-'           % < 1e-005
        for k = 1 : exp_power - 1
            str_Fr_Init_Zeros(k) = '0' ;
        end

        str_Fr = strcat (  '0.', str_Fr_Init_Zeros, str_Fr_Bef_Exp ) ;

    elseif exp (2) == '+'       % > 1.0e+015

        str_Fr = str_Fr_Bef_Exp ;

        % Normally, the foll "if" loop should not be necessary
        % since exp format does not occur for powers <= 1.0e+015. Still ...
        if length ( str_Fr ) > exp_power + 1

            str_Fr ( end + 1 ) = str_Fr (end) ;

            for j = length (str_Fr) - 1  :  -1 :  ...
                    length (str_Fr) - (exp_power + 1) + 1

                str_Fr ( j ) = str_Fr (j-1) ;
            end

            str_Fr (exp_power + 2) = '.' ;

        end

        % Foll logic when exp_power > 1.0e+015, like for eg,
        % str_Fr = '1.0450137e+018'
        % implies str_Fr_Bef_Exp = 10450137 (length = 8)
        % ie, str_Fr should become 10450137 0000 0000 000 (length = 19)
        % ie, padding with 0s at the end is reqd.
        if length ( str_Fr ) < exp_power
            str_Fr = strcat ( str_Fr,  ...
                repmat ( ['0'], 1, exp_power - (length ( str_Fr ) - 1) ) ) ;
        end

    end

end

%                                       &&&&&&&&&&&&

% 4) Separate the whole integer and fraction parts of str_Fr.
[bef_dec, Fr_dec] = strtok ( str_Fr, '.' ) ;

%                                       &&&&&&&&&&&&

% Now, we have bef_dec as the whole integer part, and
% the Fractional part starting "."

% 5) Convert first the whole integer part to binary
% by calling the std Matlab's fn dec2bin.m
bef_bin = dec2bin ( str2num (bef_dec) ) ;

%                                       &&&&&&&&&&&&

% 6) Now, finally, deal with the Fractional Part.

len_strFr = length (Fr_dec) ;
% eg of Fr_dec = '.123456789'  or  = '.000000001'  or  = '.12402343750000'

% The Fractional Part Fr_bin should start here with the dot :
% We will later concatenate bef_bin and Fr_bin
%
% Note : The part about the Fractional Part Fr_bin is not as starightforward
% as the Fractional part Fr_dec in the dual file Fr_bin2dec.m
% It is more complex due to the fact that we need to find the decreasing
% powers of 2 that will match with Fr_dec.

Fr_bin = '.' ;

Fr_dec_Current = str2num (Fr_dec) ;

for k = 1 : 16
    if Fr_dec_Current  >=  2^(-k)
        % Fr_bin = strcat ( Fr_bin,  repmat (['0'], 1, k - length(Fr_bin)), ...
        %     '1'  ) ; % Old round about code, but it seems it still works !
        Fr_bin = strcat ( Fr_bin,  '1'  ) ;

        Fr_dec_Current = Fr_dec_Current - 2^-(k) ;

        % Don't go beyond the pt where the current decremented balance
        % is 0 or negative. This will happen if input dec is <= 2^(-16) !
        if Fr_dec_Current <= 0    % Uncomment foll when you want to see details
            % fprintf ( '\n **********  Fr_dec_Current <= 0  ********** \n' ) ;
            % fprintf ( '\n *******  Pausing ... Prees any Key  ******* \n' ) ;
            % pause
            break ;
        end

    else
        % Fr_bin = strcat ( Fr_bin,  repmat (['0'], 1, k - length(Fr_bin)), ...
        %     '0'  ) ; % Old round about code, but it seems it still works !
        Fr_bin = strcat ( Fr_bin,  '0'  ) ;
    end

end % for k = 1 : 16

% k, Fr_bin, Fr_dec_Current     % Uncomment for testing

% Note that since precision is set to 16, the limit in our code is :
% 2^(-16) = 0.0000152587890625
% So, if a fraction is less than 2^(-16), we will have Fr_bin = "."
% at this point.

%                                       ++++++++++++

% 6-b) Also, check at the next level 2^(-k-1) ie, beyond the above k
% to add 1 at the end if Fr_dec_Current >= the half mark.!
% At the limit of k = 16 above, 2^(-17) = 0.00000762939453125

% However, we need to take caution if the no is lower than 2^(-16)
% in which case Fr_bin at this point, would be just '.0000000000000000'

if length(Fr_bin) == 17  &  all ( Fr_bin == '.0000000000000000' )
% Note for R13 : If short-circuiting double && were used (not in R12),
% the 2nd expr will NOT be evaluated if the 1st is false
% ie, if false AND X is always false, so X is not computed.

% However, it is observed that even with this single &,
% the 2nd expr is not computed if the 1st expr is false.

    if Fr_dec_Current >= 2^-(17)
        % Fr_bin = strcat ( Fr_bin,  repmat ( ['0'], 1, 16 ),  '1'  ) ; % Old
        Fr_bin = strcat ( Fr_bin,  '1'  ) ;

        Fr_dec_Current = Fr_dec_Current - 2^-(17) ;
        if Fr_dec_Current  >=  2^(-18)
            Fr_bin = strcat (  Fr_bin, '1' ) ;
        end

    else
        % Fr_bin = strcat ( Fr_bin,  repmat ( ['0'], 1, 17 ),  '1'  ) ;  % Old
        Fr_bin = strcat ( Fr_bin,  '0'  ) ;

        Fr_dec_Current = Fr_dec_Current - 2^-(18) ;
        if Fr_dec_Current  >=  2^(-19)
            Fr_bin = strcat (  Fr_bin, '1' ) ;
        end
    end

elseif Fr_dec_Current  >=  2^(-k-1)
    % At this point, normally, k should be 16
    % unless at some point above, Fr_dec_Current <= 0
    Fr_bin = strcat (  Fr_bin, '1' ) ;
    % fprintf ( '\n      ************  Last 1 added.  ************ \n' ) ;
end

%                                       &&&&&&&&&&&&

% 7) Concatenate the whole integer part and the fraction parts.
Fr_bin = strcat ( bef_bin, Fr_bin ) ;

% Fr_bin

% class_Fr_bin = class(Fr_bin)        % = char (Note)
% But note that the dual function :
% Fr_dec = Fr_bin2dec (bin) returns a double !

%                                   ********************

% 8) Some additional Test Cases :

% dec < 2^(-16) = 0.0000152587890625 (nearer to 2^-16 than 2^-17)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000131835937)
%   = 0.000000000000000011 (16 0s, 1, 1)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.000000000000000011 ) )
%   = 0.000000000000000011 (16 0s, 1, 1)
% Fr_bin2dec ( 0.000000000000000011 ) = 0.000011444091796875
% (= 2^-17 + 2^-18)         in place of 0.0000131835937

%                                       ++++++++++++

% dec < 2^(-16) = 0.0000152587890625 (nearer to 2^-17 than 2^-16)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000101835937)
%   = 0.00000000000000001  (16 0s, 1)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.00000000000000001 ) )
%   = 0.00000000000000001  (16 0s, 1)
% Fr_bin2dec ( 0.00000000000000001 ) = 0.00000762939453125
% (= 2^-17)                in place of 0.0000101835937

%                                       ++++++++++++

% Midway betn 2^-17 and 2^-18 = 0.0000057220458984375
% dec < 2^(-17) = 0.00000762939453125   (nearer to 2^-18 than 2^-17)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000056835937)
%   = 0.00000000000000000 (17 0s)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.00000000000000000 ) )
%   = 0.00000000000000000 (17 0s)
% Fr_bin2dec ( 0.00000000000000000 ) = 0.0
%                          in place of 0.0000056835937

%                                       ++++++++++++

% dec < 2^(-17) = 0.00000762939453125   (nearer to 2^-17 than 2^-18)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000070835937)
%   = 0.000000000000000001 (17 0s, 1)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.000000000000000001 ) )
%   = 0.000000000000000001 (17 0s, 1)
% Fr_bin2dec ( 0.000000000000000001 ) = 0.000003814697265625
% (= 2^-18)                 in place of 0.0000070835937

%                                       ++++++++++++

% dec < 2^(-17) = 0.00000762939453125   (nearer to 2^-18 than 2^-17)
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.0000039935937)
%   = 0.00000000000000000  (17 0s)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.00000000000000000 ) )
%   = 0.00000000000000000 (17 0s)
% Fr_bin2dec ( 0.00000000000000000 ) = 0.0
%                          in place of 0.0000039935937
%                          in place of anything < (2^-17 - 2^-19)
%                          ie, < Midway betn 2^-17 and 2^-18
%                          ie, < 0.0000057220458984375

%                                       ++++++++++++

% dec = 0.00001652587890625  very slightly >  2^(-16) = 0.0000152587890625
% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (0.00001652587890625)
%   = 0.0000000000000001   (15 0s, 1)
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 0.0000000000000001 ) )
%   = 0.0000000000000001   (15 0s, 1)
% Fr_bin2dec ( 0.0000000000000001 ) = 0.0000152587890625
%                         in place of 0.00001652587890625

%                                       ++++++++++++

% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (999 + 2^-11 + 2^-9)
%   = 1111100111.00000000101
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 1111100111.00000000101 ) )
%   = 1111100111.00000000000000000
%
% Fr_bin2dec ( 1111100111.00000000101 ) = 999
% [Fr_dec, str_Fr, Fr_bin] = Fr_bin2dec ( 1111100111.00000000101 )
% gives Fr_dec = 999 in place of 999.00244140625 ,
% str_Fr = 1111100111 and an empty Fr_bin
% because of the precision = 16 limit !
%
% However, Fr_bin2dec ( .00000000101 ) = 0.00244140625
% This shows that if the all the digits of the input bin are used
% for a pure fraction, the results are likely to be more accurate
% since we have more margin wrt the limit of 16 digits.

%                                       ++++++++++++

% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (879.0010365625)
%   = 1101101111.00000000010000111
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 1101101111.00000000010000111 ) )
%   = 1101101111.00000000000000000
% Fr_bin2dec ( 1101101111.00000000010000111 ) = 879
%                                   in place of 879.0010365625

%                                       ++++++++++++

% [Fr_bin, str_Fr, Fr_dec] = Fr_dec2bin (879.0012765625)
%   = 1101101111.00000000010100111
%
% Fr_bin = Fr_dec2bin ( Fr_bin2dec ( 1101101111.00000000010100111 ) )
%   = 1101101111.00000000000000000
% Fr_bin2dec ( 1101101111.00000000010100111 ) = 879
%                                   in place of 879.0012765625

%                                       ++++++++++++

%                                   ********************

% 9) Some useful values :
% (Pl note that the char length below in each line may cross 80 chars !
% But wrapping will not look nice nor easy to understand !)
%
% 2^-9    = 0.001953125
% 2^-10   = 0.0009765625
% 2^-11   = 0.00048828125
% 2^(-14) = 0.00006103515625
% 2^(-15) = 0.000030517578125
% 2^(-16) = 0.0000152587890625
% 2^(-17) = 0.00000762939453125
% 2^(-18) = 0.000003814697265625
% 2^(-19) = 0.0000019073486328125
%
% The pgm was tested with these values during development.
% These values can be spot-tested by testing the result of Fr_dec2bin ( dec). For eg :
% Fr_dec2bin ( 0.000285828865257397324183692319802554 ) ; = 0.00000000000100101
%
% Test Base = 2^(-15) :
% 2^(-15) + 2^(-16)    = 0.0000457763671875                         0.0000000000000011
%
% 2^(-15) + 2^(-16.9)  = 0.0000386945607188132718216934686662547    0.00000000000000101
% 2^(-15) + 2^(-17)    = 0.00003814697265625                        0.00000000000000101
% 2^(-15) + 2^(-17.1)  = 0.0000376360549281067460325725041667934    0.0000000000000010
%
% 2^(-15) + 2^(-18)    = 0.000034332275390625                       0.0000000000000010
% 2^(-15) + 2^(-18.01) = 0.0000343059253518563686429336937594913    0.0000000000000010
% 2^(-15)              = 0.000030517578125

% Test Base = 2^(-14) :
% 2^(-14) + 2^(-15)    = 0.000091552734375                          0.000000000000011
% 2^(-14) + 2^(-15.1)  = 0.0000895090634624269841302900166671735    0.00000000000001011
%
% 2^(-14) + 2^(-15.55) = 0.0000826143426875777442749281116365005    0.0000000000000101
% 2^(-14) + 2^(-16)    = 0.0000762939453125                         0.0000000000000101
%
% 2^(-14) + 2^(-16.55) = 0.0000714572163143493310459230799506385    0.00000000000001001
% 2^(-14) + 2^(-17)    = 0.00006866455078125                        0.00000000000001001
%
% 2^(-14) + 2^(-18)    = 0.000064849853515625                       0.0000000000000100
% 2^(-14) + 2^(-18.01) = 0.0000648235034768563686429336937594913    0.0000000000000100
% 2^(-14)              = 0.00006103515625                           0.00000000000001

% Test Base = 2^(-13) :

% 2^(-13) + 2^(-14)    = 0.00018310546875                           0.00000000000011
% 2^(-13) + 2^(-14.01) = 0.00018268386812970189828693910015186      0.00000000000010111
% 2^(-13) + 2^(-14.1)  = 0.000179018126924853968260580033334347     0.00000000000010111
% 2^(-13) + 2^(-14.45) = 0.000166750662107715619528003885783119     0.00000000000010101
%
% 2^(-13) + 2^(-14.75) = 0.000158362033538901399741134642656265     0.0000000000001010
% 2^(-13) + 2^(-15)    = 0.000152587890625                          0.000000000000101
%
% 2^(-13) + 2^(-15.75) = 0.000140216173019450699870567321328132     0.0000000000001001
% 2^(-13) + 2^(-16)    = 0.0001373291015625                         0.0000000000001001
% 2^(-13)              = 0.0001220703125                            0.0000000000001

% Test Base = 2^(-12) :
% 2^(-12) + 2^(-13)    = 0.0003662109375                            0.0000000000011
% 2^(-12) + 2^(-13.55) = 0.000327517105514794648367384639605108     0.0000000000010101
%
% 2^(-12) + 2^(-14)    = 0.00030517578125                           0.00000000000101
% 2^(-12) + 2^(-14.55) = 0.000285828865257397324183692319802554     0.00000000000100101
%
% 2^(-12) + 2^(-15)    = 0.000274658203125                          0.000000000001001
% 2^(-12)              = 0.000244140625                             0.000000000001

% Test Base = 2^(-1) :
% 2^(-1) + 2^(-16)     = 0.5000152587890625                         0.1000000000000001
% 2^(-1) + 2^(-16.1)   = 0.500014236953606213492065145008334        0.10000000000000001

% 2^(-1) + 2^(-16.9)   = 0.500008176982593813271821693468666        0.10000000000000001
% 2^(-1) + 2^(-17)     = 0.50000762939453125                        0.10000000000000001
% 2^(-1) + 2^(-17.1)   = 0.500007118476803106746032572504167        0.1000000000000000

% 2^(-1) + 2^(-17.99)  = 0.500003841230583407283053713600975        0.1000000000000000
% 2^(-1) + 2^(-18)     = 0.500003814697265625                       0.1000000000000000
% 2^(-1) + 2^(-18.01)  = 0.500003788347226856368642933693759        0.1000000000000000
% 2^(-1)               = 0.5                                        0.1

% Test Base = 2^(-16) :
% 2^(-16) + 2^(-16.99) = 0.0000229412502293145661074272019509372    0.00000000000000011
% 2^(-16) + 2^(-17)    = 0.00002288818359375                        0.00000000000000011
% 2^(-16) + 2^(-17.01) = 0.0000228354835162127372858673875189825    0.0000000000000001
%
% 2^(-16) + 2^(-17.99) = 0.0000191000196459072830537136009754686    0.0000000000000001
% 2^(-16) + 2^(-18)    = 0.000019073486328125                       0.0000000000000001
% 2^(-16) + 2^(-18.01) = 0.0000190471362893563686429336937594913    0.0000000000000001
%
% 2^(-16) + 2^(-19)    = 0.0000171661376953125                      0.0000000000000001
% 2^(-16)              = 0.0000152587890625                         0.0000000000000001

% Test Base = 2^(-17) :
% 2^(-17)              = 0.00000762939453125                        0.00000000000000001

% Test Base = 2^(-18) :
% 2^(-18)              = 0.000003814697265625                       0.0000000000000000

% Midway betn 2^-17 and 2^-18  =  0.0000057220458984375 is just above 0 ;
% 0.0000057220458984375  is the limit for this set of programmes.
% Anything < (2^-17 - 2^-19) ie, anything < 0.0000057220458984375 is 0.
%
%                                   ********************


text:
%function [path, filename] = splitpath(str)
function [path, filename] = splitpath(str)
  slashidx = strfind(str, '/');

  if isempty(slashidx)
    path = [];
    filename = str;
  else
    path = str(1:slashidx(end)-1);
    filename = str(slashidx(end)+1:end);
  end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% input: rvqs0: r s in s0, v s in s0, q s0 2 s, note s0 is a fixed local
% world frame, also called w-frame
% rqs02e: rs0 in e, q s0 2 e, a, acceleration by IMU, w, angular rate by
% gyro, dt, sample interval from the last measurement to the current
% measurement.
% gwomegaw, $g^w$, gravity in local world frame, 3-vector, and $w_{iw}^w$, 3-vector,
% if all of them set 0, this amount to integration in inertial frame

function rvqs0_new=strapdown_local_quat_bias(rvqs0, rqs02e, a, w, dt, gwomegaw)
rvqs0_new=zeros(size(rvqs0));
gs0=gwomegaw (1:3); % gravity in s0 frame
wie2s0=gwomegaw (4:6);
%Update attitude
%% method (1) second order integration
qe=rvec2quat_v000(wie2s0*dt);
vr_a=quatmult_v000(rvqs0(7:10),qe);
qb=rvec2quat_v000(-w*dt);
rvqs0_new(7:10)=quatmult_v000(qb,vr_a);
%% method (2) Runge-Kutta 4th order integration, empirically, this sometimes
% gives worse result than second order integration
% wie2s=quatrot_v000(rvqs0(7:10),wie2s0,0);
% omega=zeros(4,2);
% omega(1,2)=dt;
% omega(2:4,1)=w-wie2s;
% omega(2:4,2)=lastw-wie2s;
% qs2s0=rotationRK4( omega, [rvqs0(7); -rvqs0(8:10)]);
% rvqs0_new(7:10)=[qs2s0(1); -qs2s0(2:4)];

%% better velocity and position integration than first order rectanglar rule
%Update Velocity
vel_inc1=(quatrot_v000(rvqs0(7:10),a*dt,1)+quatrot_v000(rvqs0_new(7:10),a*dt,1))/2;
vel_inc2=(gs0+2*cross(rvqs0(4:6),wie2s0))*dt;
rvqs0_new(4:6)=rvqs0(4:6)+vel_inc1+vel_inc2;
%Update_pos
rvqs0_new(1:3)=rvqs0(1:3)+(rvqs0(4:6)+rvqs0_new(4:6))*dt/2;

repo: ekfmonoslam


text:
classdef SymmetricFourthOrder3DVoigtTensor <  AbstractTensor ...
                                              & FourthOrderDescriptor ...
                                              & VoigtRepresentation ...
                                              & Elasticity3dDescriptor

    methods (Access = public)

        function obj = SymmetricFourthOrder3DVoigtTensor()
        end

        function createRandomTensor(obj)
            obj.createRandomTensor@AbstractTensor();
            obj.makeSymmetrization();
        end
    end

    methods (Access = protected)

        function makeSymmetrization(obj)
            t = obj.getValue;
            ts = 0.5*(t + t');
            obj.setValue(ts);
        end

        function loadTensorSize(obj)
            obj.tensorSize = [6,6];
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
clear all, close all, clc
load testSys;  % previously saved system

q = 2;   % Number of inputs
p = 2;   % Number of outputs
n = 100; % State dimension
sysFull = drss(n,p,q); % Discrete random system
r = 10;  % Reduced model order

%% Plot Hankel singular values
hsvs = hsvd(sysFull); % Hankel singular values
figure
subplot(1,2,1)
semilogy(hsvs,'k','LineWidth',2)
hold on, grid on
semilogy(r,hsvs(r),'ro','LineWidth',2)
ylim([10^(-15) 100])
subplot(1,2,2)
plot(0:length(hsvs),[0; cumsum(hsvs)/sum(hsvs)],'k','LineWidth',2)
hold on, grid on
plot(r,sum(hsvs(1:r))/sum(hsvs),'ro','LineWidth',2)
set(gcf,'Position',[1 1 550 200])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/FIG_BT_HSVS');

%% Exact balanced truncation
sysBT = balred(sysFull,r);  % Balanced truncation

%% Compute BPOD
[yFull,t,xFull] = impulse(sysFull,0:1:(r*5)+1);
sysAdj = ss(sysFull.A',sysFull.C',sysFull.B',sysFull.D',-1);
[yAdj,t,xAdj] = impulse(sysAdj,0:1:(r*5)+1);
% Not the fastest way to compute, but illustrative
% Bot xAdj and xFull are size m x n x 2
HankelOC = [];  % Compute Hankel matrix H=OC
for i=2:size(xAdj,1) % start at 2 to avoid the D matrix
    Hrow = [];
    for j=2:size(xFull,1)
        Ystar = permute(squeeze(xAdj(i,:,:)),[2 1]);
        MarkovParameter = Ystar*squeeze(xFull(j,:,:));
        Hrow = [Hrow MarkovParameter];
    end
    HankelOC = [HankelOC; Hrow];
end
[U,Sig,V] = svd(HankelOC);
Xdata = [];
Ydata = [];
for i=2:size(xFull,1)  % start at 2 to avoid the D matrix
    Xdata = [Xdata squeeze(xFull(i,:,:))];
    Ydata = [Ydata squeeze(xAdj(i,:,:))];
end
Phi = Xdata*V*Sig^(-1/2);
Psi = Ydata*U*Sig^(-1/2);
Ar = Psi(:,1:r)'*sysFull.a*Phi(:,1:r);
Br = Psi(:,1:r)'*sysFull.b;
Cr = sysFull.c*Phi(:,1:r);
Dr = sysFull.d;
sysBPOD = ss(Ar,Br,Cr,Dr,-1);

%% Plot impulse responses for all methods
figure
impulse(sysFull,0:1:60), hold on;
impulse(sysBT,0:1:60)
impulse(sysBPOD,0:1:60)
legend('Full model, n=100','Balanced truncation, r=10','Balanced POD, r=10')

%% Plot impulse responses for all methods
figure
[y1,t1] = impulse(sysFull,0:1:200);
[y2,t2] = impulse(sysBT,0:1:100)
[y5,t5] = impulse(sysBPOD,0:1:100)
subplot(2,2,1)
stairs(y1(:,1,1),'LineWidth',2);
hold on
stairs(y2(:,1,1),'LineWidth',1.2);
stairs(y5(:,1,1),'LineWidth',1.);
ylabel('y_1')
title('u_1')
set(gca,'XLim',[0 60]);
grid on
subplot(2,2,2)
stairs(y1(:,1,2),'LineWidth',2);
hold on
stairs(y2(:,1,2),'LineWidth',1.2);
stairs(y5(:,1,2),'LineWidth',1.);
title('u_2')
set(gca,'XLim',[0 60]);
grid on
subplot(2,2,3)
stairs(y1(:,2,1),'LineWidth',2);
hold on
stairs(y2(:,2,1),'LineWidth',1.2);
stairs(y5(:,2,1),'LineWidth',1.);
xlabel('t')
ylabel('y_2')
set(gca,'XLim',[0 60]);
grid on
subplot(2,2,4)
stairs(y1(:,2,2),'LineWidth',2);
hold on
stairs(y2(:,2,2),'LineWidth',1.2);
stairs(y5(:,2,2),'LineWidth',1.);
xlabel('t')
set(gca,'XLim',[0 60]);
grid on
subplot(2,2,2)
legend('Full model, n=100',['Balanced truncation, r=',num2str(r)],['Balanced POD, r=',num2str(r)])
set(gcf,'Position',[100 100 550 350])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../figures/FIG_BT_IMPULSE');

author: dynamicslab


text:
% demo_geodet     December 16, 2012

% demonstrates the geocentric-to-geodetic functions

% Orbital Mechanics with Matlab

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

clear all;

global req flat

% conversion factor - degrees-to-radians

rtd = 180.0d0 / pi;

% Earth equatorial radius (kilometers)

req = 6378.1363;

% Earth flattening factor (non-dimensional)

flat = 1.0 / 298.257;

% Earth polar axis (kilometers)

rpolar = req * (1.0d0 - flat);

% eci position vector (kilometers)

reci(1) = -.586479273288D+04;
reci(2) = -.178173078828D+04;
reci(3) = -.215629990858D+04;

rmag = norm(reci);

dec = asin(reci(3) / rmag);

[alt1, lat1] = geodet1 (rmag, dec);

clc; home;

fprintf('geocentric declination    %14.8f  degrees \n\n', rtd * dec);

fprintf('geocentric radius         %14.8f  kilometers \n\n', rmag);

fprintf('\ngeodet1 function\n');
fprintf('================\n\n');

fprintf('geodetic latitude         %14.8f  degrees \n\n', rtd * lat1);

fprintf('geodetic altitude         %14.8f  kilometers \n\n', alt1);

[lat2, alt2] = geodet2(dec, rmag);

fprintf('\ngeodet2 function\n');
fprintf('================\n\n');

fprintf('geodetic latitude         %14.8f  degrees \n\n', rtd * lat2);

fprintf('geodetic altitude         %14.8f  kilometers \n\n', alt2);

[alt5, lat5] = geodet5 (req, rpolar, reci);

fprintf('\ngeodet5 function\n');
fprintf('================\n\n');

fprintf('geodetic latitude         %14.8f  degrees \n\n', rtd * lat5);

fprintf('geodetic altitude         %14.8f  kilometers \n\n', alt5);

[alt6, lat6] = geodet6 (req, rpolar, reci);

fprintf('\ngeodet6 function\n');
fprintf('================\n\n');

fprintf('geodetic latitude         %14.8f  degrees \n\n', rtd * lat6);

fprintf('geodetic altitude         %14.8f  kilometers \n\n', alt6);

author: Sable

repo: mcbench-benchmarks


text:
function [R2,out]=rsquared(x,y,kmax,method,h)

%RSQUARED calculates the R-squared value of the robust or classical PCR/PLS analysis. This function
% is used in rpcr.m and rsimpls.m.
%
% The Robust R-squared is described in:
%    Hubert, M., Verboven, S. (2003),
%    "A robust PCR method for high-dimensional regressors",
%    Journal of Chemometrics, 17, 438-452.
%
% The required input arguments
%    x    : the explanatory variables
%    y    : the response variables
%    kmax : the number of components to choose in the ROBPCA method
%  method : the method for which the R2 has to be computed. It can be 'RSIMPLS' or 'RPCR'.
%
%Optional input arguments
%    h    : the quantile used in RPCR/RSIMPLS
%
% I/O: R2=rsquared(x,y,10)
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S.Verboven on 31-10-2002
% Last updated on 18-02-2003

if nargin<4
    error('Missing input arguments')
end
[n,p]=size(x);
[n,q]=size(y);
cutoffWeights = sqrt(chi2inv(0.975,q));

if nargin<5
    h=floor(0.75*n);
end

count=1;
if q>1
    while count*q+q+(q*(q+1)/2) <= h
        count=count+1;
    end
else
    while count+2<h
        count=count+1;
    end
end
ktot=max(1,min(count-1,kmax));
if count-1<kmax
    disp(['Warning (rpcr): The maximal number of components is set to ', num2str(ktot),' to avoid overfitting.'])
end
weight=zeros(n,1);
bound=floor((n+ktot+q)/2);
h=max(h,bound);
switch method
case 'RPCR'
    outcv = cvRpcr(x,y,ktot,0,h);
case 'RSIMPLS'
    outcv = cvRsimpls(x,y,ktot,0,h);
end
R2 = outcv.R2;
W = outcv.outWeights.w_min;
rss = outcv.rss;

% place the two figures next to each other in the middle of the screen
bdwidth=5;
topbdwidth=30;
set(0,'Units','pixels');
scnsize=get(0,'ScreenSize');
pos1=[bdwidth, 1/3*scnsize(4)+bdwidth, scnsize(3)/2-2*bdwidth, scnsize(4)/2-(topbdwidth+bdwidth)];
pos2=[pos1(1)+scnsize(3)/2, pos1(2), pos1(3), pos1(4)];

figure('Position',pos1)
set(gcf,'Name', 'R-squared curve', 'NumberTitle', 'off');
plot(1:ktot,R2,'o-');
xlabel('Number of components');
ylabel('R2');
title(method)
ylabel('Robust R^2-value');
set(gca,'XTick',1:1:ktot)

figure('Position',pos2)
set(gcf,'Name', 'Square Root of Residual Sum of Squares curve', 'NumberTitle', 'off');
plot(1:ktot,sqrt(rss),'o-');
xlabel('Number of components');
title(method)
ylabel('Square Root of Robust RSS-value');
set(gca,'XTick',1:1:ktot)

k=input('How many components would you like to retain? ');
out.k=k;
out.weights=W(:,1:k);
out.rss = rss;

% closing the figures.
% close
% close

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function varargout = multiwaveplot(varargin)
% Plot stacked waves from a multichannel matrix
%
%   multiwaveplot(wave)
%   multiwaveplot(x,y,wave)
%   multiwaveplot(...,gain)
%   multiwaveplot(...,mode)
%   h = multiwaveplot(...)
%
%   Multiwaveplot draws a series of stacked waves (one on
%   top of the other) contained in the rows of an input 2-D
%   matrix. Each wave has a designated row on the plot; the
%   first row is plotted at the bottom of the plot.
%
%   multiwaveplot(wave) draws a series of waves
%   contained in the rows of wave.
%
%   multiwaveplot(x,y,wave) plot the waves against the data
%   in x and y, such that x specifies the common x data, and
%   y determines the vertical position of each row. Hence,
%   length(x)==size(wave,2) and length(y)==size(wave,1).
%
%   multiwaveplot(...,gain) scales the height of each wave
%   according to gain. With gain=1 (default), the height of
%   the tallest wave will be limited so as not to encroach
%   on the adjacent wave; all other waves are scaled
%   accordingly. This can be overridden by specifying gain >
%   1.
%
%   multiwaveplot(...,mode) plots the data with the
%   specified mode. The default mode depends on gain: for
%   gain <= 1, mode defaults to 'plot' and plots lines for
%   each row of wave; for gain > 1, mode defaults to 'fill'
%   and instead plots white patch objects for each row of
%   wave, covering the area under each wave, such that waves
%   with a lower row index mask those with a higher row
%   index. This behaviour can be overridden by specifying
%   MODE as a string (either 'plot' or 'fill').
%
%   h = multiwaveplot(...) returns a vector of handles to
%   patch (for mode = 'fill') or lineseries (for mode =
%   'plot') graphics objects, one handle per wave.
%
%   See also FILL, PATCH, PLOT, IMAGESC.

% !---
% ==========================================================
% Last changed:     $Date: 2013-06-24 17:45:41 +0100 (Mon, 24 Jun 2013) $
% Last committed:   $Revision: 249 $
% Last changed by:  $Author: ch0022 $
% ==========================================================
% !---

%% Derive inputs

input_spec = 'X and Y must be vectors, WAVE must be a matrix, GAIN must be a scalar and MODE must be a string';

% MODE
mode = find_inputs(@ischar,varargin,...
    ['Unknown parameter specified. ' input_spec]);

% GAIN
gain = find_inputs(@(x)(isscalar(x) & ~ischar(x)),varargin,...
    ['Unknown parameter specified. ' input_spec]);

if isempty(gain)
    gain = 1;
end

% Condtionally set MODE depending on GAIN
if isempty(mode)
    if gain > 1
        mode = 'fill';
    else
        mode = 'plot';
    end
end

% WAVE
wave = find_inputs(@(x)(isnumeric(x) & all(numel(x)>size(x)) & length(size(x))<3),varargin,...
    ['Both X and Y must be specified. ' input_spec]);

assert(~isempty(wave),'Please specify WAVE (which should be a matrix)')

[r,c] = size(wave);

% X and Y
[x,y] = find_inputs(@(x)(~isscalar(x) & isvector(x) & ~ischar(x)),varargin,...
    ['Both X and Y must be specified. ' input_spec]);

if isempty(x)
    x = 1:c;
    y = 1:r;
end

%% Plot

if min(wave(:))>=0
    adjust = 1; % for positive data (e.g. correlograms)
else
    adjust = 0.5; % for wave data varying on zero
end

for n = 1:r
    % just in case all values are zero
    if ~all(wave(n,:)==0)
        wave(n,:) = (adjust/max(abs(wave(:)))).*wave(n,:);
    end
end

scale = zeros(r,1); % this parameter allows for non-linear y-values, scaling each channel accordingly
h = zeros(r,1);

axis_min = y(1)-(adjust*(y(2)-y(1)));

hold on;

for n = r:-1:1
    try % calculate scaling factor
        scale(n) = y(n+1)-y(n);
    catch
        scale(n) = y(n)-y(n-1);
    end
    wave(n,:) = wave(n,:).*scale(n);
    switch mode % plot
        case 'plot'
            h(n) = plot(x,gain.*wave(n,:)+y(n),'k');
        case 'fill'
            xa=[x(1) x(1:c) x(c) x(1)];
            ya=[axis_min gain.*wave(n,:)+y(n) axis_min axis_min];
            h(n) = fill(xa,ya,'w');
        otherwise
            error('Unknown MODE specified: must be ''fill'' or ''plot''')
    end
end

hold off
axis([x(1) x(c) axis_min max(y(r)+(gain.*wave(r,:)))]); % set axis limits
set(gca,'layer','top') % move axis to top layer
box on

varargout(1:nargout) = {h};

% end of multiwaveplot()

% ----------------------------------------------------------
% Local functions:
% ----------------------------------------------------------

% ----------------------------------------------------------
% find_inputs: validate and provide inputs (if any)
% ----------------------------------------------------------
function varargout = find_inputs(fhandle,input,msg)

indices = cellfun(fhandle,input);

if any(indices) % inputs are specified
    if sum(indices)~=nargout
        error(msg) % return error message
    end
    varargout(1:nargout) = input(indices);
else % unspecified returns empty matrix
    varargout(1:nargout) = {[]};
end

% [EOF]


text:
function evalPCKbackproject(expidxs,bUseHeadSize,bSave,bHeadSizeFromRect)

fprintf('evalPCK()\n');

if (nargin < 2)
    bUseHeadSize = true;
end

if (nargin < 3)
    bSave = true;
end

if (nargin < 4)
    bHeadSizeFromRect = true;
end

if (bUseHeadSize)
    prefix = 'pck-head-';
    range = 0:0.05:0.5;
else
    prefix = 'pck-torso-';
    range = 0:0.02:0.2;
end

fontSize = 18;
figure(100); clf; hold on;
legendName = cell(0);
set(0,'DefaultAxesFontSize', fontSize)
set(0,'DefaultTextFontSize', fontSize)

[~,parts] = util_get_parts24();
nJoints = 14;
jidxsUpperBody = 7:12;

for expidx = expidxs
    fprintf('**********************************************************************************\n');
    % load experiment parameters
    p = rcnn_exp_params(expidx);
    % load ground truth
    [annolist, imgidxs_gt, rectidxs_gt, rect_ignore] = getAnnolist(expidx);
    annolist_gt = annolist(imgidxs_gt);

    fnameDist = [fileparts(p.evalTest) '/distAll'];

    mkdir_if_missing(p.latexDir);
    mkdir_if_missing(p.plotsDir);

    try
        assert(false); % no load
        load(fnameDist,'distAll','accAll','pidxs','keypointsAll');
        assert(exist('keypointsAll','var') && length(pidxs) == length(p.pidxs) && ...
               sum(pidxs ~= p.pidxs) == 0 && size(accAll,2) == nJoints+2);
    catch
        assert(length(p.pidxs) == nJoints);
        % load keypoints
        if (exist([p.evalTest '.mat'],'file')>0)
            fprintf('load %s\n',p.evalTest);
            load([p.evalTest '.mat']);
            assert(exist('keypointsAll','var')>0||exist('bboxAll','var')>0)
        else
            warning('file not found: %s\n',[p.evalTest '.mat']);
            keypointsAll = backprojectKeypointsImg(expidx,1,length(annolist_gt),false);
        end

        distAll = nan(sum(cellfun(@length,rectidxs_gt(imgidxs_gt))),nJoints);
        accAll = zeros(length(range),nJoints+2);

        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            % part is a joint
            assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
            jidx = parts(pidx+1).pos(1);

            % compute distance to the ground truth jidx
            distAll(:,i) = getNormGTJointDistBackproject(keypointsAll,annolist_gt,jidx,bUseHeadSize,rectidxs_gt(imgidxs_gt),bHeadSizeFromRect);
        end

        matchPCK = double(distAll <= range(end));

%         if (isfield(p,'singlePerson') && p.singlePerson)
%             issingle = getSinglePersonImages(annolist_gt);
%             distAll(issingle==0,:) = nan;
%         end
        matchPCK(isnan(distAll)) = nan;

        for i = 1:length(p.pidxs)
            dist = distAll(:,i);
            % remove the cases without the ground truth
            dist(isnan(dist)) = [];
            % compute accuracy for each threshold
            for k = 1:numel(range)
                accAll(k,i) = 100*mean(dist<=range(k));
            end
        end

        % compute avg PCKh upper body
        dist = reshape(distAll(:,jidxsUpperBody),size(distAll,1)*length(jidxsUpperBody),1);
        dist(isnan(dist)) = [];
        for k = 1:numel(range)
            accAll(k,end-1) = 100*mean(dist<=range(k));
        end

        % compute avg PCKh full body
        dist = reshape(distAll,size(distAll,1)*size(distAll,2),1);
        dist(isnan(dist)) = [];
        for k = 1:numel(range)
            accAll(k,end) = 100*mean(dist<=range(k));
        end

        pidxs = p.pidxs;
        try
            save(fnameDist,'distAll','accAll','pidxs','keypointsAll','matchPCK');
        catch
            warning('Unable to save %s\n',fnameDist);
        end
    end

    if (bSave)
        tableFilename = [p.latexDir '/' prefix 'expidx' num2str(expidx) '.tex'];
        [row, header] = genTable(accAll(end,:),p.name);
        fid = fopen(tableFilename,'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',row{1});fclose(fid);
        fid = fopen([p.latexDir '/' prefix 'header.tex'],'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',header);fclose(fid);
    end
    auc = area_under_curve(scale01(range),accAll(:,end));
    plot(range,accAll(:,end),'color',p.colorName,'LineStyle','-','LineWidth',3);
    legendName{end+1} = sprintf('%s, AUC: %1.1f%%', p.name, auc);

    fprintf('%s\n',legendName{end});
end

legend(legendName,'Location','NorthWest');
set(gca,'YLim',[0 100]);

xlabel('Normalized distance');
ylabel('Detection rate, %');

if (bSave)
    print(gcf, '-dpng', [p.plotsDir '/' prefix 'expidx' num2str(expidx) '.png']);
    printpdf([p.plotsDir '/' prefix 'expidx' num2str(expidx) '.pdf']);
end

end


text:
function M_hat = run_mc(params)
  M = params.M;
  Idx = params.Idx;
  Ms = sparse(M);
  [~,~,VALS] = find(Ms);
  Known = logical(Idx);
  [U,S,V,hist] = ScGrassMC(Known, VALS, 2,...
                    'tol', 1.e-6,...
                    'maxit',100,...
                    'grad_type','scaled',...
                    'beta_type','P-R',...
                    'sigma_type','approx',... %'tol_reschg', tol_reschg,...
                    'verbose', 1);
  M_hat = (U*S*V');
end

author: andrewssobral


text:
function loadSession(baseDir)
% Loads the mrSESSION and SCANS structures from mrSESSION.mat.
%
%  loadSession([baseDir])
%
% baseDir: directory containing the mrSESSION.mat file
%             defaults to current directory (pwd)
%
% djh, 2/17/2001
% 7/16/02 djh, update mrSESSION to version3.01
% - inhomoCorrect can now be 0, 1, or 2
%   previous default of 0 is now option 1 (divide by mean)
%   0 now means do nothing
% - eliminate mrSESSION.vAnatomyPath
%   replace with global vANATOMYPATH
%
% ras, 04/05 -- got rid of all references to versions earlier
% than 3.0, on the grounds that none of those sessions would
% work anyomre anyway, and times have changed...

mrGlobals

if ~exist('HOMEDIR', 'var') || isempty(HOMEDIR), HOMEDIR = pwd; end
if ~exist('baseDir', 'var') || isempty(baseDir), baseDir = HOMEDIR; end

mrSessPath = fullfile(baseDir, 'mrSESSION.mat');

% Load mrSESSION & dataTYPES
if exist(mrSessPath,'file'), load(mrSessPath)
else  myErrorDlg(sprintf('No mrSESSION.mat file in %s \n', baseDir))
end

% Old mrSESSIONs don't seem to have this field.  Can we do better than
% setting it empty?  Maybe we should write and call a mrSESSION update
% routine that queries the user? -- BW
if(~isfield(mrSESSION,'sessionCode')), mrSESSION.sessionCode = 'sCode'; end
if(~isfield(mrSESSION,'description')), mrSESSION.description = 'sDesc'; end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function output = callmexpress11(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
Q       = interfacedata.Q;
K       = interfacedata.K;
x0      = interfacedata.x0;
integer_variables = interfacedata.integer_variables;
binary_variables = interfacedata.binary_variables;
UB      = interfacedata.ub;
LB      = interfacedata.lb;

showprogress('Calling XPRESS',options.showprogress);

SENSE = 1;
C = full(c);
if isempty(F_struc)
    A = zeros(1,length(c));A(1)=1;
    B = 1e6;
else
    A =-(F_struc(:,2:end));
    B = full(F_struc(:,1));
end

% XPRESS complains about large bounds
if isempty(LB)
  %  LB = repmat(-2147483647,length(c),1);
else
    LB(LB<-2147483647) = -2147483647;
end
if isempty(UB)
  %  UB = repmat(2147483647,length(c),1);
else
    UB(UB>2147483647) = 2147483647;
end

CTYPE = repmat('L',K.l+K.f,1);  % Inequalities
CTYPE(1:K.f) = 'E';             % Equality constraints
VARTYPE = repmat('C',size(A,2),1);
VARTYPE(integer_variables)='I'; % Integer variables
VARTYPE(binary_variables) ='B'; % Binary variables

if nnz(Q)==0
    H = [];
else
    H = full(2*Q);
end

if options.verbose==0
    options.xpress.msglev = 0;
else
    options.xpress.msglev = 1;
end

if options.savedebug
    save mexpressdebug SENSE H C A B CTYPE LB UB VARTYPE options
end

% Call mex-interface
solvertime = tic;
[x,FMIN,STATUS,EXTRA] = xpress(H,C,A,B,LB,UB,CTYPE,VARTYPE,SENSE,options.xpress);
solvertime = toc(solvertime);
problem = 0;
if isstruct(EXTRA)
    D_struc = -EXTRA.lambda;
else
    D_struc = [];
end

% Check (error code depends on problem type!)
if isempty(union(binary_variables,integer_variables))
    switch STATUS
        case 1
            problem = 0;
        case 2
            problem = 1;
        case 5
            problem = 2;
        case 4
            problem = 3;
        case {0,3,6}
            problem = 11;
        otherwise
            problem = -1;
    end
else
    switch STATUS
        case 6
            problem = 0;
        case 5
            problem = 1;
        case 4
            problem = 2;
        case {4,3}
            problem = 3;
        case {0,1,2}
            problem = 11;
        otherwise
            problem = -1;
    end

end
infostr = yalmiperror(problem,'MEXPRESS');

% Save all data sent to solver?
if options.savesolverinput
	solverinput.H = H;
    solverinput.A = A;
	solverinput.C = C;
	solverinput.B = B;
	solverinput.CTYPE = CTYPE;
	solverinput.LB = LB;
	solverinput.UB = UB;
else
	solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
	solveroutput.x = x;
    solveroutput.FMIN = FMIN;
    solveroutput.STATUS = STATUS;
    solveroutput.EXTRA=EXTRA;
else
	solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [K] = ku1u1(x, xp, hyp, i)

logsigma = hyp(1);
logtheta = hyp(2);

n_x = size(x,1);
n_xp = size(xp,1);

x = repmat(x,1,n_xp);
xp = repmat(xp',n_x,1);

switch i

case 0

K=exp(1).^(logsigma+(-1/2).*exp(1).^((-1).*logtheta).*(x+(-1).*xp).^2);

case 1 % logsigma

K=exp(1).^(logsigma+(-1/2).*exp(1).^((-1).*logtheta).*(x+(-1).*xp).^2);

case 2 % logtheta

K=(1/2).*exp(1).^(logsigma+(-1).*logtheta+(-1/2).*exp(1).^((-1).*logtheta) ...
  .*(x+(-1).*xp).^2).*(x+(-1).*xp).^2;

otherwise

        K = zeros(n_x, n_xp);
end

if K == 0

    K = zeros(n_x, n_xp);

end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function FtI = warped_img_inv(I,u,v)

[M,N] = size(I);

% Stack the motion field into vector
vectors = zeros(length(xPos),2);
vectors(:,1) = u(sub2ind(size(u),yPos,xPos));
vectors(:,2) = v(sub2ind(size(u),yPos,xPos));

% Calculate frame2to1
[xPosv, yPosv] = meshgrid(1:1:n_size(2),1:1:n_size(1));
% vectors_full is the full size motion field
vectors_full = zeros(n_size(1),n_size(2),2);
vectors_full(:,:,1) = u;
vectors_full(:,:,2) = v;

xPosv = xPosv+vectors_full(:,:,1);
yPosv = yPosv+vectors_full(:,:,2);
xPosv = reshape(xPosv,n_size(1),n_size(2));
yPosv = reshape(yPosv,n_size(1),n_size(2));

xPosv(xPosv <= 1) = 1;
yPosv(yPosv <= 1) = 1;
xPosv(xPosv >= n_size(2)) = n_size(2);
yPosv(yPosv >= n_size(1)) = n_size(1);

FtI = interp2(I,xPosv,yPosv,'cubic');

author: thomas-koehler

repo: SupER


text:
function [ x, w ] = gen_hermite_dr_compute ( n, alpha )

%*****************************************************************************80
%
%% GEN_HERMITE_DR_COMPUTE computes a generalized Gauss-Hermite rule.
%
%  Discussion:
%
%    The integral to be approximated has the form:
%
%      Integral ( -oo < x < +oo ) |x|^ALPHA exp(-x^2) f(x) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Input, real ALPHA, the parameter.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  if ( n == 1 )
    arg = ( alpha + 1.0 ) / 2.0;
    x(1) = 0.0;
    w(1) = gamma ( arg );
    return
  end

  if ( mod ( n, 2 ) == 0 )
    n_laguerre = floor ( n / 2 );
    alpha_laguerre = ( alpha - 1.0 ) / 2.0;
  elseif ( mod ( n, 2 ) == 1 )
    n_laguerre = floor ( n - 1 ) / 2;
    alpha_laguerre = ( alpha + 1.0 ) / 2.0;
  end

  [ x_laguerre, w_laguerre ] = gen_laguerre_ss_compute ( n_laguerre,  alpha_laguerre );

  if ( mod ( n, 2 ) == 0 )

    x(1:n_laguerre) = - sqrt ( x_laguerre(n_laguerre:-1:1) );
    x(n_laguerre+1:n_laguerre+n_laguerre) = sqrt ( x_laguerre(1:n_laguerre) );

    w(1:n_laguerre) = 0.5 * w_laguerre(n_laguerre:-1:1);
    w(n_laguerre+1:n_laguerre+n_laguerre) = 0.5 * w_laguerre(1:n_laguerre);

  elseif ( mod ( n, 2 ) == 1 )

    x(1:n_laguerre) = - sqrt ( x_laguerre(n_laguerre:-1:1) );
    x(n_laguerre+1) = 0.0;
    x(n_laguerre+2:n_laguerre+n_laguerre+1) = sqrt ( x_laguerre(1:n_laguerre) );

    w(1:n_laguerre) = 0.5 * w_laguerre(n_laguerre:-1:1) ...
                             ./ x_laguerre(n_laguerre:-1:1);

    arg = ( alpha + 1.0 ) / 2.0;
    w(n_laguerre+1) = gamma ( arg ) ...
      - sum ( w_laguerre(1:n_laguerre) ./ x_laguerre(1:n_laguerre) );

    w(n_laguerre+2:n_laguerre+n_laguerre+1) = 0.5 * w_laguerre(1:n_laguerre) ...
                                                             ./ x_laguerre(1:n_laguerre);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = fliplr(f)
%FLIPLR   This has no effect on DELTAFUN object.

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

% This has no effect!

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function prob_test0275 ( )

%*****************************************************************************80
%
%% TEST0275 tests CARDIOID_CDF, CARDIOID_CDF_INV and CARDIOID_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  a = 0.0;
  b = 0.25;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0275\n' );
  fprintf ( 1, '  For the Cardioid PDF:\n' );
  fprintf ( 1, '  CARDIOID_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  CARDIOID_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  CARDIOID_PDF evaluates the PDF;\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %f\n', a );
  fprintf ( 1, '  PDF parameter B = %f\n', b );

  if ( ~cardioid_check ( a, b ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ x, seed ] = cardioid_sample ( a, b, seed );
    pdf = cardioid_pdf ( x, a, b );
    cdf = cardioid_cdf ( x, a, b );
    x2 = cardioid_cdf_inv ( cdf, a, b );

    fprintf ( 1,'  %12f  %12f  %12f  %12f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
%% Copyright (C) 2016, 2019 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym intersect (@var{A}, @var{B})
%% Return the common elements of two sets.
%%
%% Example:
%% @example
%% @group
%% A = finiteset(sym(1), 2, 3);
%% B = finiteset(sym(pi), 2);
%% intersect(A, B)
%%   @result{} ans = (sym) @{2@}
%% @end group
%% @end example
%%
%% The sets can also be intervals or a mixture of finite sets
%% and intervals:
%% @example
%% @group
%% C = interval(sym(2), 10);
%% intersect(A, C)
%%   @result{} ans = (sym) @{2, 3@}
%%
%% D = interval(0, sym(pi));
%% intersect(C, D)
%%   @result{} ans = (sym) [2, π]
%% @end group
%% @end example
%%
%% @seealso{@@sym/union, @@sym/setdiff, @@sym/setxor, @@sym/unique,
%%          @@sym/ismember, @@sym/finiteset, @@sym/interval}
%% @end defmethod

function r = intersect(a, b)

  if (nargin ~= 2)
    print_usage ();
  end

  cmd = {
         'a, b = _ins'
         'if isinstance(a, sp.Set) or isinstance(b, sp.Set):'
         '    return a & b,'
         ''
         'A = sp.FiniteSet(*(list(a) if isinstance(a, sp.MatrixBase) else [a]))'
         'B = sp.FiniteSet(*(list(b) if isinstance(b, sp.MatrixBase) else [b]))'
         'C = A & B'
         'return sp.Matrix([list(C)]),'
        };

    r = pycall_sympy__ (cmd, sym(a), sym(b));

end

%!test
%! A = sym([1 2 3]);
%! B = sym([1 2 4]);
%! C = intersect(A, B);
%! D = sym([1 2]);
%! assert (isequal (C, D))

%!test
%! % one nonsym
%! A = sym([1 2 3]);
%! B = [1 2 4];
%! C = intersect(A, B);
%! D = sym([1 2]);
%! assert (isequal (C, D))

%!test
%! % empty
%! A = sym([1 2 3]);
%! C = intersect(A, A);
%! assert (isequal (C, A))

%!test
%! % empty input
%! A = sym([1 2]);
%! C = intersect(A, []);
%! assert (isequal (C, sym([])))

%!test
%! % scalar
%! syms x
%! assert (isequal (intersect([x 1], x), x))
%! assert (isequal (intersect(x, x), x))

%!test
%! A = interval(sym(1), 3);
%! B = interval(sym(2), 5);
%! C = intersect(A, B);
%! assert( isequal( C, interval(sym(2), 3)))

repo: octsympy


text:
% LTFAT - Simple auditory processing
%
%  Peter L. Søndergaard, 2011 - 2018
%
%  Plots
%     SEMIAUDPLOT      - 2D plot on auditory scale.
%
%  Auditory scales
%     AUDTOFREQ        - Auditory unit to frequency conversion.
%     FREQTOAUD        - Frequency to auditory unit conversion.
%     AUDSPACE         - Auditory unit spaced vector
%     AUDSPACEBW       - Auditory unit spaced vector by equal bandwidth.
%     ERBTOFREQ        - Erb scale to frequency conversion.
%     FREQTOERB        - Frequency to erb scale conversion.
%     ERBSPACE         - Equidistant points on the erb scale.
%     ERBSPACEBW       - Equidistant points by equal bandwidth.
%     AUDFILTBW        - Bandwidth of audiory filters.
%
%  Range compression
%     RANGECOMPRESS    - Compress range of signal (mu-law etc).
%     RANGEEXPAND      - Expand range of signal.
%
%  Auditory filters
%     GAMMATONEFIR     - Gammatone FIR approximation.
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function B = BoundMirrorShrink(A)
% Shrink the matrix to remove the padded mirror boundaries
%
% for example
%
% A = [
%     5  4  5  6  12  6
%     2  1  2  3  11  3
%     5  4  5  6  12  6
%     8  7  8  9  13  9
%     5  4  5  6  12  6
%     ]
%
% B = BoundMirrorShrink(A) will yield
%
%     1  2  3  11
%     4  5  6  12
%     7  8  9  13

% Chenyang Xu and Jerry L. Prince, 9/9/1999
% http://iacl.ece.jhu.edu/projects/gvf

[m,n] = size(A);
yi = 2:m-1;
xi = 2:n-1;
B = A(yi,xi);


text:
function p = betaScores(r)
%BETASCORES Compute Beta scores for rank vector
%   p = BETASCORES(r)
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%   INPUTS:
%       r   vector of normalized rank values on interval [0,1]
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%   OUTPUTS:
%       p   a vector of p-values that corresponds to the sorted input
%           vector. The NaN-s are moved to the end.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%   See also CORRECTBETAPVALUES, THRESHOLDBETASCORE, AGGREGATERANKS.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%   Copyright (2013) Nejc Ilc <nejc.ilc@gmail.com>
%   Based on R package RobustRankAggreg written by Raivo Kolde.
%   Reference:
%     Kolde, R., Laur, S., Adler, P., & Vilo, J. (2012).
%     Robust rank aggregation for gene list integration and meta-analysis.
%     Bioinformatics, 28(4), 573-580
%
%   Revision: 1.0 Date: 2013/05/16
%--------------------------------------------------------------------------
    n = sum(~isnan(r));
    p = nan(1,length(r));
    % Sort the values.
    r = sort(r);
    % Get the order statistics and calculates p-values for each of the
    % order statistics. These are based on their expected distribution
    % under the null hypothesis of uniform distribution.
    p(1:n) = betacdf(r(1:n),1:n,n:-1:1);
end

author: BatzoglouLabSU

repo: SIMLR


text:
%%  Gaussian mask
%   k: the scaling factor
%   s: standard variance
%   M: mask

function M = gaussianMask(k,s)
M=[];
if s==0
    M=1;
    return;
end
R = ceil(3*s); % cutoff radius of the gaussian kernal
for i = -R:R,
    for j = -R:R,
        M(i+R+1,j+R+1) = k * exp(-(i*i+j*j)/2/s/s)/(2*pi*s*s);
    end
end


text:
function retPhases = atlasEstimatePhases(vw, areaCorners, scanNum, fieldName)
% Define retinotopy phases for an angle and wedge map
%
%   retPhases = atlasEstimatePhases(vw, areaCorners, scanNum, fieldName)
%
% Lets the user interactively define the data phases corresponding to the
% fovea, periphery, upper and lower vertical meridians for the atlas and
% data fits.
%
% scanNum is the scan number of the wedge and ring in this view (and for
% the view's selected data type). fieldName is a 1 x 2 cell array
% specifying the field containing the polar angle and eccentricity data,
% respectively. (For PRF analyses, these are not necessarily loaded into
% the phase slot.)
%
if notDefined('fieldName'),		fieldName = {'ph' 'ph'};			end

nAreas = length(areaCorners);
for ii=1:nAreas
    corners = areaCorners{ii};
    retPhases(ii,1) = atlasEstimateBoundaryPhases(vw, corners(1,:), ...
							corners(2,:), scanNum(1), fieldName{1});
    retPhases(ii,2) = atlasEstimateBoundaryPhases(vw, corners(3,:), ...
							corners(4,:), scanNum(1), fieldName{1});
    retPhases(ii,3) = atlasEstimateBoundaryPhases(vw, corners(1,:), ...
							corners(4,:), scanNum(2), fieldName{2});
    retPhases(ii,4) = atlasEstimateBoundaryPhases(vw, corners(2,:), ...
							corners(3,:), scanNum(2), fieldName{2});
end

if nAreas > 1, retPhases = meanPhase(retPhases); end

%% this part of the code seems to attempt to adjust the phase difference
%% between the stored data in the flat view, and real world units. A couple
%% comments / points of confusion here:
%% (1) It seems that things would be much simplified if the downstream code
%% to this function dealt exclusively with data already in real-world
%% units.
%% (2) for traveling-wave data, the functions polarAngle and eccentricity
%% map between phase and real-world units. Perhaps this should be used?
%% (3) for pRF data, the loaded data are already in real world units.
%% (4) what exactly do the retPhases do?
%% ras, 04/09.
prompt = {'Foveal', 'Peripheral', 'Lower Phase (angle map)'};
def = {num2str(retPhases(1)), num2str(retPhases(2)), num2str(retPhases(3))};
dlgTitle = 'Adjust retinal phase estimates';
lineNo = 1;
answer = inputdlg(prompt, dlgTitle, lineNo, def);

% adjust the retPhases based on the user response
angleShift = retPhases(3) - str2double(answer{3});
for ii=1:3, retPhases(ii) = str2double(answer{ii}); end

% We shift the fourth (UVM) by the same amount as the LVM
retPhases(4) = retPhases(4) - angleShift;

return;

%----------------------------------------------------
function meanPh = atlasEstimateBoundaryPhases(vw, p1, p2, scanNum, fieldName)
%
%   meanPh = atlasEstimateBoundaryPhases(vw, p1, p2, scanNum, fieldName)
%
% Author:  Wandell
% Purpose:
%     Estimate the average phase along a line between two points.  This
%     code is used to provide first estimates of the foveal, peripheral,
%     UVM and LVM phase in building atlases.
%
% Example:
%
%   scanNum = 1
%   retPhase(1)= atlasEstimateBoundaryPhases(vw,corners(1,:),corners(2,:),scanNum);
%
%
% ras 04/2009:  allows the field to be passed as a parameter -- not only
% confined to phase field, and doesn't require multiple different scans for
% polar angle / eccentricity. This is critical for use with pRF models.

curSlice = viewGet(vw, 'Current Slice');
[x, y] = findLinePoints([p1(1) p1(2)], [p2(1) p2(2)]);

newCoords = zeros(3,length(x));
newCoords(1,:) = y;
newCoords(2,:) = x;
newCoords(3,:) = curSlice*ones(1,length(x));

% Convert coords to canonical frame of reference
newCoords = curOri2CanOri(vw, newCoords);
ph = getCurDataROI(vw, fieldName, scanNum, newCoords);
if ~isequal( lower(fieldName), 'ph' )
	% put the data into the range [0 2*pi].
	% there may be a chance that it is not correct that the ph data span
	% this range -- for instance, if someone provides polar angle data in a
	% different field. I should deal with this contingency, but it might be
	% even better to have the downstream code deal with real-world-unit
	% data, and leave the conversion to these units up to the user before
	% entering this function (ras 04/09)...
	ph = normalize(ph, 0, 2*pi);
end
cxph = exp(sqrt(-1)*ph);

% We add pi to the output so that the variables run from [0,2pi] instead of
% from [-pi,pi].  This is consistent with mrLoadRet encoding of phase.
meanPh = angle(mean(cxph));

if meanPh < 0, meanPh = meanPh + 2*pi; end
if meanPh > 2*pi, meanPh = meanPh - 2*pi; end

return;


text:
function DEM = pad(DEM,varargin)

%PAD add or remove a border of pixels around a GRIDobj
%
% Syntax
%
%     DEMp = pad(DEM);
%     DEMp = pad(DEM,px,val)
%
% Description
%
%     This function adds or removes a border of constant values along all
%     edges of a GRIDobj. By default, pad adds a one-pixel wide border with
%     zeros. To control the value and amount of padding, use the arguments
%     val and px, respectively. px can be negative. In this case, the DEM
%     is cropped by the number of pixels at each grid border.
%
% Input arguments
%
%     DEM     GRIDobj
%     px      amount of padding on each of the four edges of the grid.
%             For px=1 the resulting size of DEMp is DEM.size+2.
%     val     pad value (default=0). This value does not have any effects
%             if px<0.
%
% Output arguments
%
%     DEMp    padded or cropped GRIDobj
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     DEMp20 = pad(DEM,20);
%     DEMm100 = pad(DEM,-100);
%     subplot(1,2,1)
%     imageschs(DEMp20)
%     subplot(1,2,2)
%     imageschs(DEMm100)
%
% See also: GRIDobj/crop, padarray
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 18. August, 2017

% check input arguments
narginchk(1,3)
p = inputParser;
p.FunctionName = 'GRIDobj/pad';

defaultval = 0;
defaultpx  = 1;
addRequired(p,'DEM')
addOptional(p,'px',defaultpx, @(x) isnumeric(x) && isscalar(x));
addOptional(p,'val',defaultval,@(x) isnumeric(x) && isscalar(x));
parse(p,DEM,varargin{:});

px  = sign(p.Results.px)*ceil(abs(p.Results.px));
val = p.Results.val;

% new size of output grid
newsize = DEM.size + px*2;

% check if size is negative and issue an error
if any(newsize < 0)
    error('the number of pixels removed from the grid edges is too large')
end

% create new array
if val == 0
    Znew = zeros(newsize,class(DEM.Z));
elseif isnan(val)
    Znew = nan(newsize,class(DEM.Z),class(DEM.Z));
else
    Znew = zeros(newsize,class(DEM.Z))+val;
end

% transfer values to new array
if px>0;
    Znew(px+1:end-px,px+1:end-px) = DEM.Z;
elseif px<0;
    abspx = abs(px);
    Znew = DEM.Z(abspx+1:end-abspx,abspx+1:end-abspx);
end

% adjust referencing
DEM.Z           = Znew;
DEM.refmat(3,:) = DEM.refmat(3,:)-px*[DEM.refmat(2,1) DEM.refmat(1,2)];
DEM.size        = size(DEM.Z);

if ~isempty(DEM.georef);
DEM.georef.SpatialRef = refmatToMapRasterReference(DEM.refmat,size(Znew));
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
% SCAT_ENERGY Calculate scattering energy
%
% Usage
%    energy = scat_energy(S, U)
%
% Input
%    S (cell): The scattering transform.
%    U (cell): The wavelet modulus coefficients (optional).
%
% Output
%    energy (numeric): The energy of the scattering transform (the sum of the
%       squares of the coefficients). If both S and U are given, the energy of
%       both are computed and summed.

function energy = scat_energy(S, U)
	if nargin < 2
		U = [];
	end

	if ~isempty(U)
		energy = scat_energy(S) + scat_energy(U);
	else
		if iscell(S)
			energy = 0;
			for m = 1:numel(S)
				energy = energy + scat_energy(S{m});
			end
		else
			energy = 0;
			for p = 1:numel(S.signal)
				sig = S.signal{p};
				% TODO: fix so that multiple signals work!!
				%sz = size(sig);
				%sig = reshape(sig,[prod(sz(1:end-1)) sz(end)]);
				energy = energy + sum(abs(sig(:)).^2,1);
			end
		end
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [x]=bias_nodes1d(x,f)

% function [x]=bias_nodes1d(x,f)
% ------------------------------------------------------------------------
% This function biases the spacing for the entries in x using the factor f
% such that the output array goes from x(1) to x(end) but biased using x^f
% (i.e. the spacing decreasing depending on the magnitude of x).
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/09/25
%------------------------------------------------------------------------

min_x=min(x,[],2)*ones(1,size(x,2));
x=x-min_x;
max_x=max(x,[],2)*ones(1,size(x,2));
x=max_x.*((x.^f)./(max_x.^f));
x=x+min_x;

%%
% _*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 varargout = finish_processing_step(this, module, varargin)
% finishes current processing step by deleting SPM-related temporary data
% on disk and writing back processed data from disk into data matrix
% results of processing step
%
%   Y = MrImage()
%   Y.finish_processing_step(module, varargin)
%
% This is a method of class MrImage.
%
% IN
%
% OUT
%
% EXAMPLE
%   Y.finish_processing_step('smooth', fileOutputSpm);
%
%   See also MrImage

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-03
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC 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/>.

itemsSave = this.parameters.save.keepCreatedFiles;

switch itemsSave
    case 0
        itemsSave = 'none';
    case 1
        itemsSave = 'all';
end

doSave = ~strcmpi(itemsSave, 'none');
doSaveRaw = ismember(itemsSave, {'all'});
doSaveNifti = ismember(itemsSave, {'nii', 'all', 'processed'});
doSaveObject = ismember(itemsSave, {'object', 'all', 'processed'});
doSaveProcessed = ismember(itemsSave, {'processed'});

pathSave        = this.parameters.save.path;
fileRaw         = this.get_filename('prefix', 'raw');
pathRaw         = fileparts(fileRaw);
fileProcessed   = this.get_filename;
nameImage       = this.name;

% for file prefixes
isSuffix = false;
isMixedCase = true;
doKeepCreatedFiles = (this.parameters.save.keepCreatedFiles == 1) & ...
    ~strcmpi(this.parameters.save.keepCreatedFiles, 'none');
hasMatlabbatch = ismember(module, this.get_all_matlabbatch_methods());
varargout = {};

% nifti leftovers exist and resulting files have to be renamed
if hasMatlabbatch
    fileMatlabbatch = fullfile(this.parameters.save.path, ...
        'matlabbatch.mat');
    filesCreated = {
        fileRaw
        fileProcessed
        fileMatlabbatch
        };
    switch module
        case 'apply_transformation_field'
            prefixOutput = 'w';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);

        case 'coregister_to'
            % has matlabbatch, but does not create unnecessary files...,
            % since matlabbatch not executed...
            prefixOutput = '';
            fileStationaryImage = varargin{1};
            tapas_uniqc_delete_with_hdr(fileStationaryImage);
            fileOutputSpm = varargin{2}; % SPM output of reslice
            filesCreated = [
                filesCreated
                {fileStationaryImage}
                ];

        case 'apply_realign'

            prefixOutput = 'r';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);
            filesCreated = {
                fileRaw
                fileProcessed
                };

        case 'realign'

            fileWeightingRealign = varargin{1};

            prefixOutput = 'r';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);

            fileRealignmentParameters = regexprep(...
                tapas_uniqc_prefix_files(fileRaw, 'rp_'), '\.nii', '\.txt');
            fileRealignMean = tapas_uniqc_prefix_files(fileRaw, 'mean');
            filesCreated = [
                filesCreated
                {fileRealignmentParameters; fileRealignMean;}
                fileWeightingRealign
                ];

            varargout{1} = load(fileRealignmentParameters, '-ascii');

        case 'reslice'

            prefixOutput = 'r';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);

            fnTargetGeometry = varargin{1};

            % dummy image of target geometry always deleted
            tapas_uniqc_delete_with_hdr(fnTargetGeometry);

        case 'segment'
            % parse inputs
            tissueTypes = varargin{1};
            imageOutputSpace = varargin{2};
            deformationFieldDirection = varargin{3};
            splitSuffix = varargin{4};

            % the bias-corrected anatomical is the output
            prefixOutput = 'm';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);
            tpmOutputSpm = MrImage(tapas_uniqc_prefix_files(fileOutputSpm, '*', 1), ...
                'updateProperties', 'save');
            delete(tapas_uniqc_prefix_files(fileOutputSpm, '*', 1));
            tpmOutputSpm.save();

            % get current and new tissue probability map file names
            allTissueTypes = {'GM', 'WM', 'CSF', 'bone', 'fat', 'air'};
            indTissueTypes = find(ismember(lower(allTissueTypes), ...
                lower(tissueTypes)));
            nTissues = numel(indTissueTypes);

            filesTpm = cell(nTissues,1);
            filesTpmProcessed = cell(nTissues,1);

            for iTissue = 1:nTissues
                indTissue = indTissueTypes(iTissue);
                filesTpm{iTissue} = tapas_uniqc_prefix_files(fileRaw, ...
                    sprintf('c%d%', indTissue));
                % add wildcard for multi-channel segmentation
                filesTpm{iTissue} = tapas_uniqc_prefix_files(filesTpm{iTissue}, splitSuffix, 1);
                filesTpmProcessed{iTissue} = tapas_uniqc_prefix_files( ...
                    fullfile(this.parameters.save.path, sprintf('%s.nii', ...
                    lower(allTissueTypes{indTissue}))), ...
                    'tissueProbabilityMap', isSuffix, isMixedCase);
            end

            if ismember(imageOutputSpace, {'mni', 'standard', 'template', 'warped'})
                filesTpm = tapas_uniqc_prefix_files(filesTpm, 'w');
                filesTpmProcessed = tapas_uniqc_prefix_files(filesTpmProcessed, ...
                    'warped', isSuffix, isMixedCase);
            end

            % deformation field file names, if saved
            filesDeformationField = {};
            filesDeformationFieldProcessed = {};

            hasForwardField = ismember(deformationFieldDirection, {'forward', 'both', 'all'});
            hasBackwardField = ismember(deformationFieldDirection, {'backward', 'both', 'all'});
            if hasForwardField
                filesDeformationField{end+1,1} = tapas_uniqc_prefix_files(fileRaw, 'y_');
                filesDeformationField{end,1} = tapas_uniqc_prefix_files(filesDeformationField{end,1}, splitSuffix, 1);
                filesDeformationFieldProcessed{end+1,1} = ...
                    fullfile(pathSave, 'forwardDeformationField.nii');
            end
            if hasBackwardField
                filesDeformationField{end+1, 1} = tapas_uniqc_prefix_files(fileRaw, 'iy_');
                filesDeformationField{end,1} = tapas_uniqc_prefix_files(filesDeformationField{end,1}, splitSuffix, 1);
                filesDeformationFieldProcessed{end+1, 1} = ...
                    fullfile(pathSave, 'backwardDeformationField.nii');
            end

            % bias field names
            fileBiasField = cellstr(tapas_uniqc_prefix_files(fileRaw, ...
                'BiasField_'));
            fileBiasFieldProcessed = cellstr(fullfile(pathSave, ...
                'biasField.nii'));

            % join bias field and bias field corrected outputs
            % add wildcard for multi-channel segmentation
            biasFieldPath = dir(tapas_uniqc_prefix_files(fileBiasField{1}, '*', 1));
            if ~isempty(biasFieldPath)
                tmpBiasField = MrImage(tapas_uniqc_prefix_files(fileBiasField{1}, '*', 1), ...
                    'updateProperties', 'save');
                delete(tapas_uniqc_prefix_files(fileBiasField{1}, '*', 1));
                tmpBiasField.save();
            else
                this.save('fileName', fullfile(fileBiasField{1}));
            end

            % move all image files to their final names
            filesMoveSource = [
                filesTpm
                filesDeformationField
                fileBiasField
                ];
            filesMoveTarget = [
                filesTpmProcessed
                filesDeformationFieldProcessed
                fileBiasFieldProcessed
                ];
            tapas_uniqc_move_with_hdr(filesMoveSource, filesMoveTarget);

            % now load all output variables

            % load tissue probability maps
            doLoadTissueProbabilityMaps = nargin >= 1;
            if doLoadTissueProbabilityMaps
                nTissues = numel(filesTpmProcessed);
                for iTissue = 1:nTissues
                    varargout{1}{iTissue,1} = MrImage(filesTpmProcessed{iTissue}, ...
                        'updateProperties', 'all');
                end
            end

            % load deformation fields, if wanted
            doLoadDeformationFields = nargout >= 2;
            if doLoadDeformationFields
                varargout{2} = {};
                if hasForwardField
                    varargout{2}{1,1} = MrImage(...
                        filesDeformationFieldProcessed{1}, ...
                        'updateProperties', 'all');
                    % 5th dimension is the deformation dimensions
                    allDimLabels = varargout{2}{1,1}.dimInfo.dimLabels;
                    varargout{2}{1,1}.dimInfo.dimLabels{5} = ['d_', ...
                        allDimLabels{1}, allDimLabels{2}, allDimLabels{3}];
                end

                if hasBackwardField
                    varargout{2}{end+1,1} = MrImage(...
                        filesDeformationFieldProcessed{end}, ...
                        'updateProperties', 'all');
                    % 5th dimension is the deformation dimensions
                    allDimLabels = varargout{2}{end,1}.dimInfo.dimLabels;
                    varargout{2}{end,1}.dimInfo.dimLabels{5} = ['d_', ...
                        allDimLabels{1}, allDimLabels{2}, allDimLabels{3}];
                end

            end

            doLoadBiasField = nargout >= 3;
            if doLoadBiasField
                varargout{3}{1} = MrImage(fileBiasFieldProcessed, ...
                    'updateProperties', 'all');
            else
                varargout{3}{1} = {};
            end

            % other file with normalization information for old
            % segmentation
            fileSeg8 = regexprep(fileRaw, '\.nii$', '_seg8\.mat');

            % files to be deleted, if specified
            if doSaveProcessed
                filesCreated = [
                    filesCreated
                    fileSeg8
                    ];
            else
                filesCreated = [
                    filesCreated
                    filesMoveTarget
                    fileSeg8
                    ];
            end

        case 'smooth'
            prefixOutput = 's';
            fileOutputSpm = tapas_uniqc_prefix_files(fileRaw, prefixOutput);
    end

    % copy dimInfo to SPM-output file, if it exists
    % coregister has already written new file incl. dimInfo
    % segment does not change the image
    if ~any(strcmp(module, {'coregister_to'}))
        fileDimInfoRaw = this.get_filename('prefix', 'dimInfoRaw');
        if exist(fileDimInfoRaw, 'file')
            copyfile(fileDimInfoRaw, tapas_uniqc_prefix_files(fileDimInfoRaw, prefixOutput))
        end
    end

    tapas_uniqc_move_with_hdr(fileOutputSpm, fileProcessed);

    % load back data into matrix
    this.load(fileProcessed);

    % check whether dimInfo is available as well
    [pathstr,name] = fileparts(fileProcessed);
    dimInfoFileName = dir(fullfile(pathstr, [name, '_dimInfo.mat']));
    if ~isempty(dimInfoFileName) && ~isempty(dimInfoFileName.name)
        loadDimInfo = load(fullfile(dimInfoFileName.folder, dimInfoFileName.name));
        % add only non-geometry related properties, as they are already set
        % from the loaded file (and change e.g. for reslice)
        newDimInfo = MrDimInfo(loadDimInfo.objectAsStruct);
        nonGeomDims = find(~ismember(newDimInfo.dimLabels, {'x', 'y', 'z'}));
        for nUpdateDims = nonGeomDims
            dimLabel = newDimInfo.dimLabels{nUpdateDims};
            % dim exists but needs to be set to correct value
            if ismember(dimLabel, this.dimInfo.dimLabels)
                this.dimInfo.set_dims(dimLabel, ...
                    'units', newDimInfo.(dimLabel).units, ...
                    'samplingPoints', newDimInfo.(dimLabel).samplingPoints, ...
                    'samplingWidths', newDimInfo.(dimLabel).samplingWidths);
                % dim does not exist and needs to be added
            else
                this.dimInfo.add_dims(dimLabel, ...
                    'units', newDimInfo.(dimLabel).units{:}, ...
                    'samplingPoints', newDimInfo.(dimLabel).samplingPoints{:}, ...
                    'samplingWidths', newDimInfo.(dimLabel).samplingWidths);
            end
        end
        % also add dimInfo to the bias field
        if any(strcmp(module, {'segment'})) && ~isempty(biasFieldPath)
            varargout{3}{1}.dimInfo = newDimInfo.copyobj();
        end
    end

    % remove NaNs
    this.data(isnan(this.data)) = 0;

    % delete all unwanted files
    if ~doSaveRaw
        tapas_uniqc_delete_with_hdr(filesCreated);

        % delete raw sub-folder of current processing step
        if exist(pathRaw, 'dir')
            delete(fullfile(pathRaw, '*'));
            rmdir(pathRaw);
        end

        [stat, mess, id] = rmdir(this.parameters.save.path);
    end
else % no matlabbatch created

    % NOTE: this saving is important, e.g. for compute_masks to allow for
    % multiple processing steps within a bigger step to be saved
    if doSave
        this.save();
    end
end

% to preserve name after reloading, because not in nifti-file
this.name = nameImage;

end

author: translationalneuromodeling

repo: tapas


text:
% MLF2PDF prints a PDF figure for inclusion in LaTeX documents, with
%    LaTeX fonts. It creates a complete pdf-file with annotation of the
%    figure (titles, labels and texts) and graphics (lines,arrows,markers,
%    ...). This function uses the matlabfrag function and makes simply a
%    concatenation of the text part and the graphical part in a single PDF
%    file. It requires:
%
%                  MATLAB        |      OTHER
%               -----------------|----------------------------
%                 matlabfrag.m   |      pdflatex
%
%
%    USAGE:
%    ------------------------
%
%    mlf2pdf creates the PDF file of the current graphical figure (CGF)
%    named LaTeXfile.pdf
%
%    mlf2pdf(H) creates the PDF file from the graphical figure with
%    handle H, named LaTeXfile.pdf
%
%    mlf2pdf(H,FILENAME) creates the PDF file from the graphical figure
%    with handle H, named FILENAME.PDF (FILENAME is a character array of
%    the filename, with or without the extension .PDF)
%
%    mlf2pdf(H,FILENAME,PACKAGES) to adding
%    packages such as fonts, mathfont,...
%    ex: PACKAGES = 'amssymb, times'
%

function mlf2pdf(h,FileName,Packages)

if nargin < 2
    FileName='LaTeXfile';
end

if nargin < 1
    h=gcf;
end

TempName = strcat('TEMP',num2str(round(rand*10000))); %Generate random file name

if ~exist('matlabfrag.m')
    disp('MatLabFrag M-file does not exist or is not in the MATLAB''s search path.');
    disp('This file can be downloaded at: http://www.mathworks.com/matlabcentral/fileexchange');
    disp('                                                         Try again...');
    return;
end

 matlabfrag(TempName,'handle',h);  %call matlabfrag to export figure to .eps and .tex file.

%-------------------------------------------------------
% Temporary LaTeX file
%-------------------------------------------------------
fid = fopen(strcat(TempName,'2.tex'),'w');

fprintf(fid,'\\documentclass[11pt, oneside]{article}\n');
fprintf(fid,'\\usepackage{graphicx}\n');
fprintf(fid,'\\usepackage{amsmath}\n');
fprintf(fid,'\\usepackage[T1]{fontenc}\n');
fprintf(fid,'\\usepackage[latin1]{inputenc}\n');
fprintf(fid,'\\usepackage{ae}\n');
fprintf(fid,'\\usepackage{psfrag}\n');
fprintf(fid,'\\usepackage{color}\n');
fprintf(fid,'\\usepackage{pstool}\n');

if nargin==3
    fprintf(fid,'\\usepackage{%s}\n',Packages); % Suplementary packages
end
fprintf(fid,'\\pagestyle{empty}\n');
fprintf(fid,' \n');
fprintf(fid,'\\begin{document}\n');
fprintf(fid,'    \\begin{figure}\n');
fprintf(fid,'        \\centering\n');
fprintf(fid,'        \\psfragfig{%s}\n',TempName);
fprintf(fid,'    \\end{figure}\n');
fprintf(fid,' \n');
fprintf(fid,'\\end{document}\n');
fclose(fid);

%-------------------------------------------------------
% LaTeX Command
%-------------------------------------------------------
Str=sprintf('pdflatex -shell-escape --src -interaction=nonstopmode %s2.tex',TempName);
disp(sprintf('\n[LaTeX Command] %s',Str));
[hdos,wdos]=system(Str);

if hdos ~=0
   if isunix==0
     dos(sprintf('del %s*',TempName));
   else
     unix(sprintf('rm %s*',TempName));
   end
   error('Error %d -- LATEX:\n%s',hdos ,wdos);
   return;
end

%-------------------------------------------------------
% Rename TempFile to FileName and delete FileName before (if it exists)
%-------------------------------------------------------
if isunix==0
   dos(sprintf('del %s.pdf',FileName));
   dos(sprintf('ren %s.pdf %s.pdf',TempName,FileName));
else
   unix(sprintf('rm %s.pdf',FileName));
   unix(sprintf('mv %s.pdf %s.pdf',TempName,FileName));
end
%-------------------------------------------------------
% Success
%-------------------------------------------------------
  disp(sprintf('... OK!\nPDF file [%s.pdf] has been created in the current directory\n',FileName));
 %-------------------------------------------------------
% Delete all the temporary files
%-------------------------------------------------------
if isunix==0
  dos(sprintf('del %s*',TempName));
else
  unix(sprintf('rm %s*',TempName));
end
return;


text:
function predict_label = CNN_SVM_prediction(imdb,OP)
%
% predict_label:
%   if OP.class = 1 then: predict_label = 1:Wake, 2:REM sleep, 3:NREM light sleep, 4:NREM deep sleep
%   if OP.class = 2 then: predict_label = 1:NREM sleep (light + deep), 2:REM sleep, 3: wake
%   if OP.class = 3 then: predict_label = 1:NREM deep sleep, 2:NREM light sleep, 3:Wake + REM sleep
%   if OP.class = 4 then: predict_label = 1:NREM sleep (light + deep), 2: Wake + REM sleep

imdb.meta.features(isinf(imdb.meta.features))=NaN;
sqi=nanmean(imdb.meta.sqi,1);
sqi(find(isnan(sqi)))=nanmean(sqi);

% load model
load(OP.model);

net=stat.net;
net.layers{1,end}.type='softmax';
mean_features=stat.mean_features;
std_features=stat.std_features;

mean_sqi=stat.mean_sqi;
std_sqi=stat.std_sqi;

model=stat.model;

% standardization
for i=8:13
    for j=1:size(imdb.meta.features,1)
        imdb.meta.features(j,i)=log(imdb.meta.features(j,i));
    end
end
imdb.meta.features(find(isinf(imdb.meta.features)))=NaN;

% subtract data_mean
for k=1:size(imdb.images.data,4)
    imdb.images.data(:,:,1,k)=imdb.images.data(:,:,1,k)-stat.data_mean;
end

for i=1:size(imdb.meta.features,2)
    imdb.meta.features(find(isnan(imdb.meta.features(:,i))),i)=mean_features(i);
    imdb.meta.features(:,i)=bsxfun(@minus, imdb.meta.features(:,i),mean_features(i));
    imdb.meta.features(:,i)=bsxfun(@rdivide, imdb.meta.features(:,i),std_features(i));
end
sqi=bsxfun(@minus,sqi,mean_sqi);
sqi=bsxfun(@rdivide,sqi,std_sqi);

% CNN prediction
if size(imdb.images.data,4)>1000
    split=round(size(imdb.images.data,4)/1000);
else
    split=1;
end

retall=[];
for kk=1:split
    each=ceil(size(imdb.images.data,4)/split);
    res=vl_simplenn(net,imdb.images.data(:,:,:,(kk-1)*each+1:min(kk*each,size(imdb.images.data,4))));

    ret=[];
    for i=1:OP.nclass
        for j=1:size(res(end).x,4)
            ret(i,j)=res(end).x(1,1,i,j);
        end
    end
    retall=[retall,ret];
end
% probability of CNN output
ret=retall;

xtest=[];
for i=1:size(ret,1)
    xtest(:,end+1)=ret(i,:);
end
for i=1:13
    if i~=6 % remove DFA feature
        xtest(:,end+1)=imdb.meta.features(:,i);
    end
end
xtest(:,end+1)=sqi;

% SVM prediction
[predict_label, accuracy_test, ytest_out] = svmpredict(ones(size(xtest,1),1), xtest, model,'-b 1');
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [ x, y, z, w ] = ld1454 ( )

%*****************************************************************************80
%
%% LD1454 computes the 1454 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(1454,1);
  y = zeros(1454,1);
  z = zeros(1454,1);
  w = zeros(1454,1);
  a = 0.0;
  b = 0.0;
  v = 0.7777160743261247E-04;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.7557646413004701E-03;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.3229290663413854E-01;
  v = 0.2841633806090617E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.8036733271462222E-01;
  v = 0.4374419127053555E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1354289960531653;
  v = 0.5417174740872172E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1938963861114426;
  v = 0.6148000891358593E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2537343715011275;
  v = 0.6664394485800705E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3135251434752570;
  v = 0.7025039356923220E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3721558339375338;
  v = 0.7268511789249627E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4286809575195696;
  v = 0.7422637534208629E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4822510128282994;
  v = 0.7509545035841214E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5320679333566263;
  v = 0.7548535057718401E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6172998195394274;
  v = 0.7554088969774001E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6510679849127481;
  v = 0.7553147174442808E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6777315251687360;
  v = 0.7564767653292297E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6963109410648741;
  v = 0.7587991808518730E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.7058935009831749;
  v = 0.7608261832033027E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.9955546194091857;
  v = 0.4021680447874916E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.9734115901794209;
  v = 0.5804871793945964E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.9275693732388626;
  v = 0.6792151955945159E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.8568022422795103;
  v = 0.7336741211286294E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.7623495553719372;
  v = 0.7581866300989608E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.5707522908892223;
  b = 0.4387028039889501;
  v = 0.7538257859800743E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5196463388403083;
  b = 0.3858908414762617;
  v = 0.7483517247053123E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4646337531215351;
  b = 0.3301937372343854;
  v = 0.7371763661112059E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4063901697557691;
  b = 0.2725423573563777;
  v = 0.7183448895756934E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3456329466643087;
  b = 0.2139510237495250;
  v = 0.6895815529822191E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2831395121050332;
  b = 0.1555922309786647;
  v = 0.6480105801792886E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2197682022925330;
  b = 0.9892878979686097E-01;
  v = 0.5897558896594636E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1564696098650355;
  b = 0.4598642910675510E-01;
  v = 0.5095708849247346E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6027356673721295;
  b = 0.3376625140173426;
  v = 0.7536906428909755E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5496032320255096;
  b = 0.2822301309727988;
  v = 0.7472505965575118E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4921707755234567;
  b = 0.2248632342592540;
  v = 0.7343017132279698E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4309422998598483;
  b = 0.1666224723456479;
  v = 0.7130871582177445E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3664108182313672;
  b = 0.1086964901822169;
  v = 0.6817022032112776E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2990189057758436;
  b = 0.5251989784120085E-01;
  v = 0.6380941145604121E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6268724013144998;
  b = 0.2297523657550023;
  v = 0.7550381377920310E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5707324144834607;
  b = 0.1723080607093800;
  v = 0.7478646640144802E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5096360901960365;
  b = 0.1140238465390513;
  v = 0.7335918720601220E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.4438729938312456;
  b = 0.5611522095882537E-01;
  v = 0.7110120527658118E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.6419978471082389;
  b = 0.1164174423140873;
  v = 0.7571363978689501E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.5817218061802611;
  b = 0.5797589531445219E-01;
  v = 0.7489908329079234E-03;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end


text:
%%%%%%%%%%%%%%%%%%%% RECOMPUTES THE REPROJECTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%

check_active_images;

% Reproject the patterns on the images, and compute the pixel errors:

ex = []; % Global error vector
x = []; % Detected corners on the image plane
y = []; % Reprojected points

if ~exist('alpha_c'),
   alpha_c = 0;
end;

for kk = 1:n_ima,

   eval(['omckk = omc_' num2str(kk) ';']);
   eval(['Tckk = Tc_' num2str(kk) ';']);

   if active_images(kk) & (~isnan(omckk(1,1))),

      %Rkk = rodrigues(omckk);

      eval(['y_' num2str(kk) '  = project_points_fisheye(X_' num2str(kk) ',omckk,Tckk,fc,cc,kc,alpha_c);']);

      eval(['ex_' num2str(kk) ' = x_' num2str(kk) ' - y_' num2str(kk) ';']);

      eval(['x_kk = x_' num2str(kk) ';']);

      eval(['ex = [ex ex_' num2str(kk) '];']);
      eval(['x = [x x_' num2str(kk) '];']);
      eval(['y = [y y_' num2str(kk) '];']);

   else

      %	eval(['y_' num2str(kk) '  = NaN*ones(2,1);']);

      % If inactivated image, the error does not make sense:
      eval(['ex_' num2str(kk) ' = NaN*ones(2,1);']);

   end;

end;

err_std = std(ex')';

repo: ekfmonoslam


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, space_v, vel, space_p, press] = ...
                          solve_stokes_2d (problem_data, method_data)

warning ('geopdes:obsolete', 'Function SOLVE_STOKES_2D is obsolete. Using SOLVE_STOKES instead')

[geometry, msh, space_v, vel, space_p, press] = solve_stokes (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [inform, x] = SteepDescent(fun, x, sdparams)
%  Implements steepest descent using simple Wolfe conditions via StepSize.m.
%
%  Implementation Parameters:
%    * STOP when
%      - Objective function gradient has norm less than gtol.
%      - maxit steps have been taken.
%    * ftol = 1.0e-20
%    * gtol = 1.0e-4
%    * xtol = 1.0e-20
%    * maxit = 1000
%    * p_k = -grad(f(x_k))
%    * alfa_k = max(10*xtol, f'(x_{k-1})*p_{k-1}*alfa_{k-1}/(f'(x_k)*p_k))
%
%  Input:
%    fun      - a pointer to a function
%    x        - the following structure:
%               * x.p - the starting point values
%               * x.f - the function value of x.p
%               * x.g - the gradient value of x.p
%    sdparams - the following structure, as an example:
%         sdparams = struct('maxit',1000,'toler',1.0e-4);
%
%  Output:
%    inform - structure containing two fields:
%      * inform.status - 1 if gradient tolerance was
%                        achieved;
%                        0 if not.
%      * inform.iter   - the number of steps taken
%    x      - the solution structure, with point, function,
%             and gradient evaluations at the solution point.

%  Number of function and gradient evaluations.
global numf numg
numf = 0;
numg = 0;

%  Populate local caching of nparams parameters, and params parameters.
toler = sdparams.toler;  % Set gradient tolerance.
maxit = sdparams.maxit;  % Set maximum number of allowed iterations.
xtol = 1.0e-20;  % Set point tolerance.
ftol = 1.0e-20;  % Set function tolerance.
gtol = toler;  % Set gradient tolerance.

%  Initialize parameter structure for StepSize function call.
params = struct('ftol', ftol, 'gtol', gtol, 'xtol', xtol);

% Below was tailored to geodesic functions.
% If calling geodesic function...
% if isfield(sdparams, 'geoparams') && isstruct(sdparams.geoparams)
%     params.geoparams = sdparams.geoparams;
%     x.f = feval(fun, x.p, 1, params.geoparams);
%     x.g = feval(fun, x.p, 2, params.geoparams);
% end
% params.geoparams = -1;

%  Compute function and gradient at starting point.
x.f = feval(fun, x.p, 1);
x.g = feval(fun, x.p, 2);

alfa = 1;  % Initial alpha value.
iter = 0;  % Number of iterations.
while iter < maxit && norm(x.g) >= toler
    pg = x.g;  % Store previous gradient value.
    d = -x.g;  % Steepest descent direction.

    %  Get step size satisfying simple Wolfe conditions.
    [alfa, x] = StepSize(fun, x, d, alfa, params);
    alfa = max(10*xtol, pg'*d*alfa / (x.g'*(-x.g)));  % Calculate new alpha.
    iter = iter + 1;  % Increment step counter.
end
inform.iter = iter;  % Number of iterations = iter at this point.

%  If norm of gradient is less than tolerance, the method succeeded.
if norm(x.g) < toler
    inform.status = 1;  % Update status to success indicator, 1.
else
    inform.status = 0;  % Update status to failure indicator, 0.
end
return;  % Return inform and final point x
end

author: clarkzinzow

repo: Nonlinear-Optimization-Algorithms