text:
function [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
           Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
           str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
           Diophantine_1 ( X_Real, X_Imag,  Y_Real, Y_Imag,  Ga2b2, Gc2d2,  ...
                           str_Det_fprt_Init )
%
% Subject  : A programme to generate more examples of Ramanujam's
%            Diophantine Equation Numbers
% Release  : R13
% Author   : Sundar Krishnan
% email id : sundark100@yahoo.com
% Date     : Read the article MATHEMATICAL MINIATURE 9.pdf in Aug 2003 ;
%            Skeleton Trial Programmes in Sep 2003 ;
%            Revisit and Fine Tuning in Sep 2004
%
% Introduction :
% --------------
% This programme is about finding many examples of Ramanujam's
% Diophantine Equation.
% For mathematical formulas used in making this programme, I have referred
% to the article MATHEMATICAL MINIATURE 9.pdf
% by John Butcher, butcher@math.auckland.ac.nz.
%
% To circumvent the problem of editing the Greek symbols, I have replaced
% them with suitable alphabets.
%
% The basic equation we are trying to solve is to get the numbers :
% x & y and u & v and N that satisfy the Diophantine equation :
% x^3 + y^3 = u^3 + v^3 = N
% where the four integers x, y, u, v have no common factor.
%
% For eg, the "lowest" Diophantine Number N of Ramanujam is :
% 1729 = 9^3 + 10^3 = 1^3 + 12^3          % see Usage Eg 1 below
%
% Other examples given in Butcher's article are :
%  9^3 + 15^3    =  2^3 + 16^3 =  4104    % see Usage Eg 3 below
% 33^3 + 15^3    =  2^3 + 34^3 = 39312    % see Usage Eg 4 below
%  6^3 + (-3)^3  =  5^3 + 4^3  =   189    % see Usage Eg 5 below
% (ie, 3^3 + 4^3 + 5^3 = 6^3)
%
% All these examples are covered in the various Usage Examples given below.
%
%                                    &&&&&&&&&&&&&&
%
% Derivation of Formulae, Notations followed and Step by Step Procedure :
% -----------------------------------------------------------------------
%
% If N is even, let us have :
% a2 = (x + y) / 2 ;    b2 = (x - y) / 2        % a', b' in the article
% c2 = (u + v) / 2 ;    d2 = (u - v) / 2        % c', d' in the article
%
% If N is odd, let us have :
% a2 = (x + y) ;        b2 = (x - y)
% c2 = (u + v) ;        d2 = (u - v)
%
% gcd (x, y,  u, v) = 1 => gcd (a2, b2,  c2, d2) = 1  % Theoretical Concept 1
%
% Let gcd (a2, b2) = Ga2b2                      % Myu  in the article
% and gcd (c2, d2) = Gc2d2                      % Veta in the article
%
% Note that gcd ( Ga2b2, Gc2d2 ) = 1                    % FULCRUM 1
%
% a = a2 / ( Ga2b2 * Gc2d2^3 ) ;    b = b2 / Ga2b2
% c = c2 / ( Gc2d2 * Ga2b2^3 ) ;    d = d2 / Gc2d2
%
%                                       ++++++++
%
% X_Mag  =  X_Real +  j * (sqrt(3) * X_Imag) ;        X_Sq  =  X_Mag^2 ;
%
% Y_Mag  =  Y_Real +  j * (sqrt(3) * Y_Imag) ;        Y_Sq  =  Y_Mag^2 ;
%
% Z_Mag  =  Z_Real +  j * (sqrt(3) * Z_Imag) ;        Z_Sq  =  Z_Mag^2 ;
%
% My interpretation : Initially assume the ratio "r" as 1 ;
% later, adjust r to get b and d as whole numbers.      % Logical Concept 1
% a = r * Y_Sq ;        c = r * X_Sq ;
%
% The "tricky" part here is to first assume r = 1
% and solve for Z_Real and Z_Imag using these 2 eqs.
% Gc2d2^3 * ( Y_Real^2 + 3 * Y_Imag^2 ) * r  =  ...
%     X_Real * Z_Real - 3 * X_Imag * Z_Imag
%
% Ga2b2^3 * ( X_Real^2 + 3 * X_Imag^2 ) * r  =  ...
%    Y_Real * Z_Real - 3 * Y_Imag * Z_Imag
%
% ie, in matrix form, solve for Z_Real and Z_Imag :
% [ X_Real,  - 3 * X_Imag ;    Y_Real,  - 3 * Y_Imag ] * [ Z_Real ;  Z_Imag ]
% = [ Gc2d2^3 * ( Y_Real^2 + 3 * Y_Imag^2 ) * r ;
%     Ga2b2^3 * ( X_Real^2 + 3 * X_Imag^2 ) * r ]
%
% Obtain b and d from :
% Gc2d2^6 * a^2  +  3 * b^2  =  X_Sq * Z_Sq
%
% Ga2b2^6 * c^2  +  3 * d^2  =  Y_Sq * Z_Sq
%
% Then, choose r in a way as to convert b and d to whole numbers.
% Also, correct a and c to reflect the change in r.
% Note that the "appropriate" roots b and d may be negative in some cases.
%
% Then, back-calculate to get a2, b2, c2, d2 :
%
% If a2 + b2 is odd, N is even, else N is odd.
% use appropriate formulae for obtaining x, y, u, v.
%
% If x^3 + y^3 and u^3 + v^3 don't match, try with negative roots -b, -d etc.
% But this probably may never happen - see below.
%
%                                    &&&&&&&&&&&&&&
%
% I wish to thank Prof John Butcher for his article which triggered
% and enabled me to write this Matlab code for "generating"
% Ramanujam's Diophantine Equation Numbers.
% I hope that this programme will be useful to many the world over.
% Time permitting, I may be improving on this programme to make it suitable
% for generating a series of Diophantine Numbers automatically.
% But the base groundwork is already laid now, and we need
% only to build further.
%
% Q1a : I am curious to know why Prof John Butcher has said that
% gcd (a, b) and gcd (c, d) should be 1.
% In the case of his own example : 9^3 + 15^3  =  2^3 + 16^3 =  4104,
% intermediate calculations give a = 12, b = 3, c = 9, d = 7.
% See Usage Eg Case 4 below. Obviously, gcd (a, b) = 3 (not 1).
% But, gcd (x, y,  u, v) = 1 (var name used in my programme is G_Gxy_Guv.)
%
% Also, in his second example : 33^3 + 15^3  =  2^3 + 34^3 = 39312,
% calculations give a = 3, b = 9, c = 9, d = 8. See Usage Eg Case 5 below.
% Obviously, gcd (a, b) = 3 (not 1). But, gcd (x, y,  u, v) = 1.
%
% Similarly, in his 3rd example too, gcd (a, b) = 3 (not 1) !
%
% In many of my own examples with randomly chosen values, I have mixed results:
% In Usage Eg 2, gcd (c, d) = 9 (not 1), but, gcd (x, y,  u, v) = 1.
% In Usage Eg 3, gcd (a, b) = 3 (not 1), but, gcd (x, y,  u, v) = 3 (not 1) !
% Egs 7 and 8 are similar to Eg 2 ; Egs 9, 10, and 11 are similar to Eg 3.
%
% Q1b) Therefore, I would like to know the condition or constraint,
% which when fulfilled, will ensure that we will certainly obtain
% x, y,  u, v such that gcd (x, y,  u, v) = 1
%
% Q2 : I would also like to know how the article's Theorem 1 is used
% in deriving the formulae for Diophantine Numbers.
% It says : ... "although beneath the surface" ...
%
%                                    &&&&&&&&&&&&&&
%
% % Usage Egs :
%
% % Usage Egs Case 6, Case 7 and Case 8 can be used as Regression Test Cases
% % after any modifications are done to this programme.
%
% % Case 1) : Ramanujam's Diophantine Equation : Lowest No 1729 :
% % 1729 = 9^3 + 10^3 = 1^3 + 12^3
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 1, 2,  4, 1,  1, 1,   'Y' )
%
% % 10, 9,  12, 1,  1729, 1729, 1729,   1,  19, 1, 13, 11,   19, 1, 13, 11,
% % 1, -3,   1, 1, 1,   1, 1,   N_ODD, Pos_Pos,   0, 1.729e-7
% %               ---  ---
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 2) : With my own randomly chosen values : (G_Gxy_Guv = 1)
% % 321236^3  +  (-78236)^3  =  418185^3  +  (-343305)^3  =  32670295158384000
% % But gcd (x, y,  u, v) = 1 ;    gcd (c, d) = 9 (NOT 1)
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 1, -2,  0, -3,  4, 5,   'Y' )
%
% % Note that this is a good example for 2 reasons :
% % a) We had to toggle our first choice of N_ODD to N_EVN.
% % b) error_Limit had to be used because the computed Calc_Error is 8,
% % but checking with the Calculator gives a difference of 0.
% % But I corrected this later by taking factored cube.
%
% % Note :
% Gcd          % = gcd ( c, d )     = 9 (not 1)      % See Q1a, b above.
% G_Gxy_Guv    % = gcd ( Gxy, Guv ) = 1              % But, = 3 in Case 3
%
% % 321236, -78236,  418185, -343305,
% % 32670295158384000, 32670295158384000,
% % 3.267029515838401e16 (32670295158384000),   9,
% % (RHS error of +8 later corrected to 0)
% % 243, 49934, 117, 76149,   121500, 199736, 37440, 380745,
% % 2820.333333333334, 92.44444444444444,   4, 45, 1,   1, 9,
% %                                               ---     ***
% % N_EVN, Pos_Pos,   8 (0), 3000    % Calc_Error later corrected to 0
% %                  *******
% %
% %                                    &&&&&&&&&&&&&&
%
% % Case 3) : With my own randomly chosen values : (G_Gxy_Guv = 3)
% % 115821^3  +  (-72045Z)^3  =  267420^3  +  (-261804)^3  =  1179732995041536
% % But gcd (x, y,  u, v) = 3 (NOT 1) ;    gcd (a, b) = 3 (NOT 1)
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 1, 2,  3, 4,  3, 4,   'Y' )
%
% % Note :
% Gab          % = gcd ( a, b )     = 3 (not 1)      % See Q1a, b above.
% G_Gxy_Guv    % = gcd ( Gxy, Guv ) = 3              % But, = 1 in Case 1, 4, 5, 6
%
% % 115821, -72045,  267420, -261804,
% % 1179732995041536, 1179732995041536, 1179732995041536,   2,
% % 114, 31311, 26, 66153,   21888, 93933, 2808, 264612,
% % -6945, -1765.5,   9, 12, 3,   3, 1,   N_EVN, Pos_Pos,   0, 3000
% %                         ***  ***
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 4) : Butcher's 2nd example :
% %  9^3 + 15^3  =  2^3 + 16^3 =  4104
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 0, 1,  1, 1,  1, 1,   'Y' )
%
% % Note :
% Gab          % = gcd ( a, b )     = 3 (not 1)      % See Q1a, b above.
% % But here, G_Gxy_Guv    % = gcd ( Gxy, Guv ) = 1  % But, = 3 in Case 3
%
% % 15, 9,  16, 2,  4104, 4104, 4104,  3,   12, 3, 9, 7,   12, 3, 9, 7,
% % -1, -1.33333333333333,   3, 2, 1,   3, 1,   N_EVN, Pos_Pos,   0, 4.104e-7
% %                               ---  ***
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 5) : Butcher's 3rd example :
% % 33^3 + 15^3  =  2^3 + 34^3 = 39312
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 0, 1,  1, 0,  1, 2,   'Y' )
%
% % Note :
% Gab          % = gcd ( a, b )     = 3 (not 1)      % See Q1a, b above.
% % But here, G_Gxy_Guv    % = gcd ( Gxy, Guv ) = 1  % But, = 3 in Case 3
%
% % 33, 15,  34, 2,  39312, 39312, 39312,  3,   3, 9, 9, 8,   24, 9, 18, 16,
% % 3, -2.66666666666667,   3, 2, 1,   3, 1,   N_EVN, Pos_Pos,   0, 3.9312e-6
% %                              ---  ***
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 6) : Butcher's 4th example :
% % 6^3 + (-3)^3  =  5^3 + 4^3  =   189        % ie, 3^3 + 4^3 + 5^3 = 6^3
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 0, 1,  -1, 0,  1, 1,   'Y' )
%
% % Note :
% Gab          % = gcd ( a, b )     = 3 (not 1)      % See Q1a, b above.
% % But here, G_Gxy_Guv    % = gcd ( Gxy, Guv ) = 1  % But, = 3 in Case 3
%
% % 6, -3,  5, 4,  189, 189, 189,  3,   3, 9, 9, 1,  3, 9, 9, 1,
% % -3, -0.33333333333333,   3, 1, 1,   3, 1,   N_ODD, Pos_Pos,   0, 1.89e-8
% %                               ---  ***
%
% %                                    &&&&&&&&&&&&&&
%
% % Many more Usage Egs are given at the end.
%
%                                 ********************

% 1) Inputs' Defaults and Checks :

if nargin < 4
    error ('Diophantine_1.m needs atleast 4 input args.' ) ;
end

if nargin < 5
    Ga2b2 = 1 ;    % Default Value of gcd (a2, b2)
    Gc2d2 = 1 ;    % Default Value of gcd (c2, d2)

    % Default is that we do NOT want extra Print-outs.
    str_Det_fprt_Init = 'N' ;
end

if nargin < 6
    Gc2d2 = 1 ;    % Default Value of gcd (c2, d2)

    str_Det_fprt_Init = 'N' ;
end

if nargin < 7
    str_Det_fprt_Init = 'N' ;
end

G_Ga2b2_Gc2d2 = gcd ( Ga2b2, Gc2d2 ) ;        % FULCRUM 1

if G_Ga2b2_Gc2d2 ~= 1
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    fprintf ( '\n**** Diophantine_1.m : ' ) ;
    fprintf ( 'gcd ( Ga2b2, Gc2d2 ) = %d **** \n', G_Ga2b2_Gc2d2 ) ;
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    % Note : Actually, we should exit here with an error stmt.
end

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

% 2) Inits :
format long

% Constant repeatedly used in the programme :
root_3 = sqrt (3) ;

% Initially assume the ratio "r" as 1        % Logical Concept 1
r = 1 ;

% If there is no convergence, let us show it as N = -1
N = -1 ;

LHS = -1 ;
RHS = -1 ;

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

% 3) Calculate X_Sq and Y_Sq :

X_Mag  =  abs ( X_Real  +  j * X_Imag * root_3 ) ;
X_Sq   =  X_Mag^2 ;

Y_Mag  =  abs ( Y_Real  +  j * Y_Imag * root_3 ) ;
Y_Sq   =  Y_Mag^2 ;

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

% 4) Calculate a and c :

a = r * Y_Sq ;
c = r * X_Sq ;

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

% 5a) Solve for s (Z_Real) and t (Z_Imag) :

% NOTE that if you choose a case like :
% ... = Diophantine_1 ( -3, -4,  -3, -4,  11, 7 )
% A_st will be singular !
% This will also lead to a problem below in : r = lcm ( Den_b, Den_d ) ;

A_st = [ X_Real,  -3 * X_Imag ;
         Y_Real,  -3 * Y_Imag ] ;

B_st = [ Gc2d2^3 * ( Y_Real^2 + 3 * Y_Imag^2 ) * r ;
         Ga2b2^3 * ( X_Real^2 + 3 * X_Imag^2 ) * r ] ;

% The 'mldivide' operator may only produce a single output.
% So, we do this in 3 steps.
st_Soln = A_st \ B_st ;

Z_Real = st_Soln(1) ;
Z_Imag = st_Soln(2) ;

%                                       ++++++++

% 5b) Compute Z_Mag and Z_Sq :

Z_Mag  =  abs ( Z_Real  +  j * Z_Imag * root_3 ) ;
Z_Sq   =  Z_Mag^2 ;

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

% 6) The important logical step in this programme :

% 6a) Obtain b and d :

b_pos = sqrt ( ( X_Sq * Z_Sq  -  Gc2d2^6 * a^2 ) / 3 ) ;

d_pos = sqrt ( ( Y_Sq * Z_Sq  -  Ga2b2^6 * c^2 ) / 3 ) ;

% Because of sqrt, b_pos and d_pos may not come out exactly as
% whole integers in the calculations due to FP errors.
% We correct this below.

%                                       ++++++++

% 5b) Choose r in a way as to convert b and d to whole numbers.

[ Num_b, Den_b ]  =  rat (b_pos) ;
[ Num_d, Den_d ]  =  rat (d_pos) ;

r = lcm ( Den_b, Den_d ) ;        % Logical Concept 1

% Try with positive roots first.
% Note that the "appropriate" roots b and d may be negative in some cases.
b_pos = r * b_pos ;
d_pos = r * d_pos ;

%                                       ++++++++

% 5c) Correct a and c to reflect the change in r :
a = r * Y_Sq ;
c = r * X_Sq ;

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

% 6) Because of sqrt, b_pos and d_pos may not come out exactly as
% whole integers in the calculations due to FP errors.

if abs ( a - round(a) ) < 1e-6 * abs (a)
    a = round(a) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** a is NOT a whole integer. ' ) ;
        fprintf ( 'a = %d **** \n' ) ;
    end
end

if abs ( c - round(c) ) < 1e-6 * abs (c)
    c = round(c) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** c is NOT a whole integer. ' ) ;
        fprintf ( 'c = %d **** \n' ) ;
    end
end

%                                       ++++++++

if abs ( b_pos - round(b_pos) ) < 1e-6 * abs (b_pos)
    b_pos = round(b_pos) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** b_pos is NOT a whole integer. ' ) ;
        fprintf ( 'b_pos = %d **** \n' ) ;
    end
end

if abs ( d_pos - round(d_pos) ) < 1e-6 * abs (d_pos)
    d_pos = round(d_pos) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** d_pos is NOT a whole integer. ' ) ;
        fprintf ( 'd_pos = %d **** \n' ) ;
    end
end

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

% 7a) As per Butcher's article, gcd(a, b) and gcd(c, d) should be = 1.

Gab = gcd (a, b_pos) ;    % should be the same as gcd (a, b) calculated below
Gcd = gcd (c, d_pos) ;    % should be the same as gcd (c, d) calculated below

if Gab ~= 1 | Gcd ~= 1
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    fprintf ( '\n** Diophantine_1.m : gcd ( a, b_pos ) = %d ** \n', Gab ) ;
    fprintf ( '\n** Diophantine_1.m : gcd ( c, d_pos ) = %d ** \n', Gcd ) ;
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    % Note : Actually, we should exit here with an error stmt.
end

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

% 8) Back calculate a2, b2, c2, d2 :
a2 = a * Ga2b2 * Gc2d2^3 ;
b2_pos = b_pos * Ga2b2 ;

c2 = c * Gc2d2 * Ga2b2^3 ;
d2_pos = d_pos * Gc2d2 ;

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

% 9) Because of sqrt in b_pos and d_pos, let us check
% a2, c2, b2_pos, d2_pos also for wholeness.

if abs ( a2 - round(a2) ) < 1e-6 * abs (a2)
    a2 = round(a2) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** a2 is NOT a whole integer. ' ) ;
        fprintf ( 'a2 = %d **** \n' ) ;
    end
end

if abs ( c2 - round(c2) ) < 1e-6 * abs (c2)
    c2 = round(c2) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** c2 is NOT a whole integer. ' ) ;
        fprintf ( 'c2 = %d **** \n' ) ;
    end
end

%                                       ++++++++

if abs ( b2_pos - round(b2_pos) ) < 1e-6 * abs (b2_pos)
    b2_pos = round(b2_pos) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** b2_pos is NOT a whole integer. ' ) ;
        fprintf ( 'b2_pos = %d **** \n' ) ;
    end
end

if abs ( d2_pos - round(d2_pos) ) < 1e-6 * abs (d2_pos)
    d2_pos = round(d2_pos) ;
else
    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n    **** d2_pos is NOT a whole integer. ' ) ;
        fprintf ( 'd2_pos = %d **** \n' ) ;
    end
end

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

% 10) As a first choice, if a2 + b2 is even, N is odd, else N is even.
% N_ODD or N_EVN may be reqd to be decided again in xyuv_from_a2b2c2d2
% For eg, it is necessary to toggle the first choice in the case of Usage Eg 6.
%
% Note that a and b are positive because X_Sq and Y_Sq are pos,
% and r = 1 or lcm is also pos.
% So, a2 and b2 must also be pos.
if mod ( (a2 + b2_pos), 2 ) < 1e-6 * min ( abs (a2), abs (b2_pos) )
% mod is preferable to rem here : rem (-5, 2) = -1 ;  mod (-5, 2) = -1 + 2 = +1
% Note that mod ( 19.999999, 2 ) = 1.999999
    str_N_odd_even = 'N_ODD' ;
else
    str_N_odd_even = 'N_EVN' ;

end

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

% 11a) Calculate x & y and u & v using b2_pos and d2_pos :
b2 = b2_pos ;
d2 = d2_pos ;

[ x, y,  u, v ] = xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even ) ;

% If the above choice of N_ODD or N_EVN was not correct, we will have
% a factor of 0.5 in either or both pairs : x, y or u, v
% So, we need to change the choice.
% To account for cascading FP errors, I have taken 0.25 instead of 0.5
if abs ( x - round(x) ) > 0.25 |  abs ( u - round(u) ) > 0.25 | ...
   abs ( y - round(y) ) > 0.25 |  abs ( v - round(v) ) > 0.25

    if str_N_odd_even == 'N_ODD'
        str_N_odd_even = 'N_EVN' ;
    elseif str_N_odd_even == 'N_EVN'
        str_N_odd_even = 'N_ODD' ;
    end

    if str_Det_fprt_Init == 'Y'
        fprintf ( '\n**** Choice of str_N_odd_even toggled. **** \n' ) ;
        % Usage Eg 6 needs this toggling.
        str_N_odd_even
    end

    % Note that this final choice of N_EVN or N_ODD will stay valid even for
    % the other match tests between x^3 + y^3 and u^3 + v^3 below.

    % Don't forget to call xyuv_from_a2b2c2d2 again.
    [ x, y,  u, v ] = xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even ) ;

end

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

% 11b) Check with pos roots b2_pos and d2_pos if x^3 + y^3 and u^3 + v^3 match:

% LHS = x^3 + y^3
% RHS = u^3 + v^3

% As values of x and y increases, errors creep in the computation
% of x^3 + y^3 due to cascading FP errors.
% Instead, we factorise and try.

LHS = ( x + y ) * ( x^2 + y^2 - x*y ) ;    % instead of : LHS = x^3 + y^3
RHS = ( u + v ) * ( u^2 + v^2 - u*v ) ;    % instead of : RHS = u^3 + v^3

Calc_Error  =  abs ( LHS - RHS ) ;

% If LHS or RHS becomes very large like say, of the order of 1e20,
% a multiplying factor 1e-4 is not enough.
% For eg, in the case of Usage Eg 6, Calc_Error computed above is 8,
% but checking with the Calculator gives a difference of 0.
error_Limit = compute_error_Limit ( LHS, RHS ) ;

if str_Det_fprt_Init == 'Y'
    fprintf ( '\nCalc_Error with pos roots b_pos and d_pos = %d \n', ...
              Calc_Error ) ;
    error_Limit
end

% if x^3 + y^3 == u^3 + v^3        % ideal if no FP errors creep in !
if Calc_Error < error_Limit        % if b_pos, d_pos

    b = b_pos ;
    d = d_pos ;

    str_b_d = 'Pos_Pos' ;

else
    % I discovered later during development that the solution will "converge"
    % in all (almost or really all ?) cases with pos roots b_pos and d_pos ;
    % it (probably) really makes no difference to the end result at all.
    % But since I had already written more than 50 % of the foll code, and
    % also tested it partly, I have let it remain - ofcourse after testing it
    % sufficiently - no harm in retaining it.
    % One secret ? : It was actually an error in my hand calculations which
    % initially prompted me to write the foll code for different signs
    % for b and d !

    % 11c) Check if x^3 + y^3 and u^3 + v^3 match with b2_neg and d2_pos :
    b2_neg = -b2_pos ;    b2 = b2_neg ;
    d2 = d2_pos ;

    [ x, y,  u, v ] = ...
        xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even ) ;

    LHS = ( x + y ) * ( x^2 + y^2 - x*y ) ;    % instead of : LHS = x^3 + y^3
    RHS = ( u + v ) * ( u^2 + v^2 - u*v ) ;    % instead of : RHS = u^3 + v^3

    Calc_Error  =  abs ( LHS - RHS ) ;

    error_Limit = compute_error_Limit ( LHS, RHS ) ;

    if str_Det_fprt_Init == 'Y'
        fprintf ( '\nCalc_Error with b_neg and d_pos = %d \n', ...
                  Calc_Error ) ;
        error_Limit
    end

    if Calc_Error < error_Limit       % if b_neg, d_pos

        b = -b_pos ;    % b_neg = -b_pos ;
        d = d_pos ;

        str_b_d = 'Neg_Pos' ;

    else

        % 11d) Check if x^3 + y^3 and u^3 + v^3 match with b2_pos and d2_neg :
        b2 = b2_pos ;
        d2_neg = -d2_pos ;    d2 = d2_neg ;

        [ x, y,  u, v ] = ...
            xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even ) ;

        LHS = ( x + y ) * ( x^2 + y^2 - x*y ) ;
        RHS = ( u + v ) * ( u^2 + v^2 - u*v ) ;

        Calc_Error  =  abs ( LHS - RHS ) ;

        error_Limit = compute_error_Limit ( LHS, RHS ) ;

        if str_Det_fprt_Init == 'Y'
            fprintf ( '\nCalc_Error with b_pos and d_neg = %d \n', ...
                      Calc_Error ) ;
            error_Limit
        end

        if Calc_Error < error_Limit  % if b_pos, d_neg

            b = b_pos ;
            d = -d_pos ;    % d_neg = -d_pos ;

            str_b_d = 'Pos_Neg' ;

        else

            % 11e) Check if x^3 + y^3 and u^3 + v^3 match
            % with b2_neg and d2_neg :
            b2_neg = -b2_pos ;    b2 = b2_neg ;
            d2_neg = -d2_pos ;    d2 = d2_neg ;

            b = -b_pos ;    % b_neg = -b_pos ;
            d = -d_pos ;    % d_neg = -d_pos ;

            str_b_d = 'Neg_Neg' ;

            [ x, y,  u, v ] = ...
                xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even ) ;

            LHS = ( x + y ) * ( x^2 + y^2 - x*y ) ;
            RHS = ( u + v ) * ( u^2 + v^2 - u*v ) ;

            Calc_Error  =  abs ( LHS - RHS ) ;

            error_Limit = compute_error_Limit ( LHS, RHS ) ;

            if str_Det_fprt_Init == 'Y'
                fprintf ( '\nCalc_Error with b_neg and d_neg = %d \n', ...
                          Calc_Error ) ;
                error_Limit
            end

            if Calc_Error > error_Limit    % if b_neg, d_neg

                fprintf ( '\n\n**** NO MATCH found. INVESTIGATE. **** \n' ) ;

                fprintf ( '\n**** Values calculated so far :   **** \n' ) ;
                x, y,  u, v,  N, LHS, RHS,  r
                a, b, c, d,   a2, b2, c2, d2
                str_N_odd_even, str_b_d,   Calc_Error, error_Limit
                Gab = gcd (a, b)
                Gcd = gcd (c, d)

                Gxy = gcd (x, y)
                Guv = gcd (u, v)

                fprintf ( '\n************************************ \n\n' ) ;

                error ( 'Can''t find matching x & y and u & v ! INVESTIGATE.');

            end    % if b_neg, d_neg

        end    % if b_pos, d_neg

    end    % if b_neg, d_pos

end    % if b_pos, d_pos

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

% 12) Compute N :

N = x^3 + y^3 ;

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

% 13) Compute gcd of x, y  &  u, v  and  a, b  &  c, d :
% a, c,  b, d
% x, y,  u, v

Gxy = gcd (x, y) ;    % If these are not integers, we could get this error :
Guv = gcd (u, v) ;    % Error using ==> gcd : Requires integer input arguments.

Gab = gcd (a, b) ;    % should be the same as gcd (a, b_pos) calculated above
Gcd = gcd (c, d) ;    % should be the same as gcd (c, d_pos) calculated above

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

% 14) Check if gcd of (x, y,  u, v ) = 1
% gcd (x, y,  u, v) = 1 => gcd (a2, b2,  c2, d2) = 1  % Theoretical Concept 1

G_Gxy_Guv = gcd ( Gxy, Guv ) ;

if G_Gxy_Guv ~= 1
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    fprintf ( '\n**** Diophantine_1.m : gcd ( Gxy, Guv ) = %d **** \n', ...
              G_Gxy_Guv ) ;
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    % Note : Actually, we should exit here with an error stmt.
end

% 7b) Similarly, we should exit if Gab ~= 1 | Gcd ~= 1   % Yes, it is NOT 14b !
% Note that this condition (eqvt) has already been addressed above as :
% gcd (a, b_pos) and gcd (c, d_pos).
if Gab ~= 1 | Gcd ~= 1
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    fprintf ( '\n** Diophantine_1.m : gcd ( a, b ) = %d ** \n', Gab ) ;
    fprintf ( '\n** Diophantine_1.m : gcd ( c, d ) = %d ** \n', Gcd ) ;
    fprintf ( '\n********************* NOTE ********************* \n' ) ;
    % Note : Actually, we should exit here with an error stmt.
end

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

fprintf ( '\n\n************  Pgm output begins here.  ************ \n ' ) ;

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

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

% 16)
function [ error_Limit ] = compute_error_Limit ( LHS, RHS )

min_L_R = min ( abs (LHS), abs (RHS) ) ;

if min_L_R < 1e6
    error_Limit = 1e-10 * min_L_R ;
    % max = 1e-4 for 1e6 - 1

elseif min_L_R < 1e10
    error_Limit = max ( 1e-8 * min_L_R, 5e-4 ) ;
    % 5e-4 for 1e6 ; 100 for 1e10 - 1

elseif min_L_R < 1e15
    error_Limit = max ( 1e-12 * min_L_R, 300 ) ;
    % 300 for 1e10 ; 1000 for 1e15 - 1

elseif min_L_R < 1e20
    error_Limit = max ( 1e-16 * min_L_R, 3000 ) ;
    % 3000 for 1e15 ; 10000 for 1e20 - 1

elseif min_L_R < 1e25
    error_Limit = max ( 1e-19 * min_L_R, 30000 ) ;
    % 30000 for 1e20 ; 100000 for 1e25 - 1

else
    error_Limit = 2e6 ;

end

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

% 17)
function [ x, y,  u, v ] = ....
    xyuv_from_a2b2c2d2 ( a2, b2,  c2, d2,  str_N_odd_even )

% Note : As a first choice, if a2 + b2 is even, str_N_odd_even is supplied
% as N_ODD to this subfunction, and vice-versa.
% N_ODD or N_EVN may be reqd to be decided again based on the results.

if str_N_odd_even == 'N_ODD'

    x = ( a2 + b2 ) / 2 ;
    y = ( a2 - b2 ) / 2 ;

    u = ( c2 + d2 ) / 2 ;
    v = ( c2 - d2 ) / 2 ;

else

    x = ( a2 + b2 ) ;
    y = ( a2 - b2 ) ;

    u = ( c2 + d2 ) ;
    v = ( c2 - d2 ) ;

end

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

% Some more Usage Egs are given below.
%
% % Case 7) : With my own randomly chosen values : (G_Gxy_Guv = 1)
% % 116244^3  +  (-115872)^3  =  37249^3  +  (-33217)^3  =  15031929519936
%
% clear, clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( -3, 5,  2, -3,  2, 1, 'Y' )
%
% Gab        % = gcd ( a, b )     = 3 (not 1)     % See Q1a, b above.
% G_Gxy_Guv  % = gcd ( Gxy, Guv ) = 1
%
% r % = 3
%
% Z_Real % = 3453
% Z_Imag % = -692.6666666666666
%
% % N_EVN, Pos_Pos
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 8) : With my own randomly chosen values : (G_Gxy_Guv = 1)
% % (-107766)^3  +  (-634932)^3  =  (-2013055)^3  +  1991671^3
% % = -257217167508536664 (Calculator)   { -2.572171675085367e+017 (Matlab) }
% % ie, 2013055^3 = 107766^3  +  634932^3  +  1991671^3
%
% % Since we have the successful special case of 3 negatives and 1 positive,
% % this is an important Regression Test Case.
%
% clear, clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( -3, 1,  2, -3,  3, -11, 'Y' )
%
% % -107766, -634932, -2013055, 1991671,
% % -257217167508536664, -257217167508536664, -257217167508536664,
% % 3,   93, 87861, 36, 182033,   -371349, 263583, -10692, -2002363,
% % 17637, -3883.333333333333,   18, 11, 1,   3, 1, N_EVN, Pos_Pos, 3000
%
% Gab        % = gcd ( a, b )     = 3 (not 1)     % See Q1a, b above.
% G_Gxy_Guv  % = gcd ( Gxy, Guv ) = 1
%
% r % = 3
%
% Z_Real % = 17637
% Z_Imag % = -3883.333333333333
%
% % str_N_odd_even toggled from N_ODD to N_EVN, Pos_Pos
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 9) : With my own randomly chosen values : (G_Gxy_Guv = 343)
% % 327908^3  +  (-72716)^3  =  360493^3  +  (-228781)^3
% % =  3.487337281103962e+016
% % But gcd (x, y,  u, v) = 343 (NOT 1) ;    gcd (c, d) = 49 (NOT 1)
%
% clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( 1, -4,  2, -3,  4, 7,   'Y' )
%
% Gcd        % = gcd ( c, d )     = 49 (not 1)     % See Q1a, b above.
%
% G_Gxy_Guv  % = gcd ( Gxy, Guv ) = 343  % = 3 in Case 2 ; = 1 in Case 1, 4, 5, 6
%
% r % = 3
%
% Z_Real % = 3871
% Z_Imag % = 1208.666666666667
%
% % N_EVN, Pos_Pos
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 10) : With my own randomly chosen values : (G_Gxy_Guv = 7)
% % 50932^3  +  (-21448)^3  =  102543^3  +  (-98511)^3  =  1.222546648901760e+014
%
% clear, clc
% [ x, y,  u, v,  N, LHS, RHS,  r,   a, b, c, d,   a2, b2, c2, d2,   ...
%   Z_Real, Z_Imag,   Gxy, Guv, G_Gxy_Guv,   Gab, Gcd,               ...
%   str_N_odd_even, str_b_d,   Calc_Error, error_Limit ] =           ...
%       Diophantine_1 ( -3, -5,  9, 8,  2, 3 )
%
% % Note :
% Gab        % = gcd ( a, b )     = 7 (not 1)     % See Q1a, b above.
% Gcd        % = gcd ( c, d )     = 7 (not 1)     % See Q1a, b above.
%
% G_Gxy_Guv  % = gcd ( Gxy, Guv ) = 7  % = 3 in Case 2 ; = 1 in Case 1, 4, 5, 6
%
% r % = 1
%
% Z_Real % = 2968
% Z_Imag % = 1085
%
% % N_EVN, Pos_Pos
%
% %                                    &&&&&&&&&&&&&&
%
% % Case 11) :
% Diophantine_1 (  1, -2,  -3, -4,  3,  8,   'Y' )
% % 1483821^3 + 267219^3 = 4444392^3 + (-4388232)^3 = 3.286046473469061e+018
% % G_Gxy_Guv % = 3 ;
% % Gab % = 1
%
% %                                    &&&&&&&&&&&&&&
%
% % Cases of type : Can't find matching x & y and u & v ! INVESTIGATE.
% Diophantine_1 (  1, -2,   5, -3,  4,  7,   'Y' )
% Diophantine_1 ( -3, -4,   4,  3,  11, 7 )
%
% %                                    &&&&&&&&&&&&&&
%
% % A_st is singular for this case.
% Diophantine_1 ( -3, -4,  -3, -4,  11, 7 )
%
%                                 ********************
%                                 ********************


text:
function [trQYtY, YQ] = evaluate_objective(Y, problem_data, use_Cholesky)
% function [trQYtY, YQ] = evaluate_objective(Y, problem_data, use_Cholesky)
%
% This function computes and returns the value of the objective function
% tr(Q Y^T Y).  Optionally, it returns the product YQ as the second
% argument

% Copyright (C) 2016 by David M. Rosen

if nargin < 3
    use_Cholesky = true;
end

Yt = Y';
YQ = Qproduct(Yt, problem_data, use_Cholesky)';

trQYtY = trace(YQ * Yt);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [points, pos, faceInds] = intersectEdgeMesh3d(edge, varargin)
%INTERSECTEDGEMESH3D Intersection points of a 3D edge with a mesh.
%
%   INTERS = intersectEdgeMesh3d(EDGE, VERTICES, FACES)
%   Compute the intersection points between a 3D edge and a 3D mesh defined
%   by vertices and faces.
%
%   [INTERS, POS, INDS] = intersectEdgeMesh3d(EDGE, VERTICES, FACES)
%   Also returns the position of each intersection point on the input edge,
%   and the index of the intersected faces.
%   For edges, the values of POS are expected to be comprised between 0 and
%   1.
%
%   Example
%     [V, F] = createCube;
%     edge = [-1 0.5 0.5  +3 0.5 0.5];
%     pts = intersectEdgeMesh3d(edge, V, F)
%     pts =
%         1.0000    0.5000    0.5000
%              0    0.5000    0.5000
%
%   See also
%     meshes3d, interesectLineMesh3d, triangulateFaces
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2021-02-24, using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2021-2022 INRA - Cepia Software Platform

% perform computation on supporting line
line = edgeToLine3d(edge);
[points, pos, faceInds] = intersectLineMesh3d(line, varargin{:});

% identifies intersection points within parameterization bounds
inds = pos >= 0 & pos <= 1;

% select relevant results
points = points(inds, :);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [X, K, cluster_sizes, true_labels] = prepare_data(noise)
    if nargin < 1
        noise = false;
    end
% Ambient space dimension
M = 4000;
% Number of subspaces
K = 4;
% common dimension for each subspace
D = 10;
% dimensions of each subspace
Ds = D * ones(1, K);
% Noise level
sigma = 0.8;
fprintf('Generating synthetic data.\n');
bases = spx.data.synthetic.subspaces.random_subspaces(M, K, Ds);
% Number of points on each subspace
Sk = 16 * D;
cluster_sizes = Sk * ones(1, K);
% total number of points
S = sum(cluster_sizes);
% Let us generate uniformly distributed points in each subspace
points_result = spx.data.synthetic.subspaces.uniform_points_on_subspaces(bases, cluster_sizes);
X0 = points_result.X;
start_indices = points_result.start_indices;
end_indices = points_result.end_indices;
if noise
    Noise = sigma * spx.data.synthetic.uniform(M, S);
    % Add noise to signal
    X = X0 + Noise;
else
    X = X0;
end
% Normalize noisy signals.
X = spx.norm.normalize_l2(X);
% true labels for each of the data samples
true_labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes);
fprintf('Size of the data matrix: %d x %d\n', M, S);
end

author: indigits

repo: sparse-plex


text:
function [fitresult, gof,stats] = evaluate_spatial_scale(data_obj,parcel_obj,varargin)
% Evaluate information coding at multiple spatial scales, using a prespecified parcellation
% and the predict function (and associated optional inputs). This function evaluates the
% spatial scale of information coding by constructing predictive models using random draws
% of features from different parcels to test whether information is contained 1) at the
% whole-brain scale, 2) within individual parcels, or 3) across all parcels (may be similar
% to whole-brain if parcels cover most cortical and subcortical areas)
% Usage:
% ::
%
%    [fitresult, gof,stats] = evaluate_spatial_scale(obj, parcel_obj, varargin);
%
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2018 Philip Kragel
%
%     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/>.
% ..
%
% :Inputs:
%
%   **data_obj:**
%   **parcel_obj:**
%
% :Optional inputs:
%
%   **colors:**
%       Cell array of RGB values that indicate colors for each parcel
%
%   **labels:**
%       Cell array of strings that indicate names for each parcel
%
%   **num_voxels:**
%      Vector of integers that indicate how many voxels to randomly sample
%      these should be smaller than the total number of voxels in the
%      smallest parcel
%
%
%   **do_plot:**
%      Automatically make plot of sampling curves
%
%   **render_brain**
%      Make a brain rendering of parcels on plot
%
%
% :Outputs:
%
%   **fitresult:**
%        cell array of fit objects for each parcel
%
%   **gof:**
%        cell array of with performance metrics for each parcel
%
%   **stats:**
%        structure with normal approximation Z-stats and P-values for
%        pairwise comparisons of effect size at the second largest spatial
%        scale
%
% :Examples:
% ::
%
%  % Example 1: Load bmrk3 data and predict pain intensity sampling
%  within/across different resting-state networks
% % -------------------------------------------------------------------------
%   canlab_help_set_up_pain_prediction_walkthrough;
%   [mask_obj, networknames, imagenames] = load_image_set('bucknerlab');
%   mycolors={[120 18 134]/255 [70 130 180]/255  [0 118 14]/255  [196 58 250]/255  [220 248 164]/255  [230 148 34]/255  [205 62 78]/255 };
%   [fitresult, gof] = evaluate_spatial_scale(image_obj,mask_obj,'cv_lassopcr', 'nfolds',[rem(subject_id,2)+1],'verbose',0,'labels',networknames,'colors',mycolors);
%
%
%  % Example 1: Load bmrk3 data and predict pain intensity sampling
%  within/across different brainstem/thalamus regions
% % -------------------------------------------------------------------------
%   canlab_help_set_up_pain_prediction_walkthrough;
%   [mask_obj, networknames, imagenames] = load_atlas('brainstem');
%   [fitresult, gof] = evaluate_spatial_scale(image_obj,mask_obj,'cv_lassopcr', 'nfolds',[rem(subject_id,2)+1],'verbose',0,'labels',networknames,'colors',mycolors);
%
%
% :References:
%   Kragel, P. A., Koban, L., Barrett, L. F., & Wager, T. D. (2018). Representation, pattern information, and brain signatures: from neurons to neuroimaging. Neuron, 99(2), 257-273.
%   Kragel, P. A., Reddan, M., LaBar, K. S., & Wager, T. D. (2018). Emotion Schemas are Embedded in the Human Visual System. bioRxiv, 470237.
%
% :See also:
%
% predict, load_atlas

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
%
%   Phil: created, Dec 2018
%
%   TODOs: - change outcome metric to optionally incude accuracy/AUC/etc.
%          - change outcome to optionally reflect % of noise ceiling
%          - account for regions that have small number of voxels
% ..

%if we have a 1d atlas object convert to matrix for
parcel_obj=resample_space(parcel_obj,data_obj,'nearest');

if size(parcel_obj.dat,2)==1
    parcel_obj.dat=condf2indic(parcel_obj.dat);
    parcel_obj.dat=parcel_obj.dat(:,2:end);

end

% defaults
num_iterations = 1000;
num_parcels = size(parcel_obj.dat,2);
num_vox = [ 50 150 250 500 750 1000 1500 5000 10000];

do_plot=true;
render_brain=false;

mycolors=scn_standard_colors(num_parcels);

atlas_labels=cell(num_parcels,1);
for i=1:length(atlas_labels)
    atlas_labels{i}='';
end

% optional inputs
if any(strcmp(varargin,'colors'))
    mycolors=varargin{find(strcmp(varargin,'colors'))+1};
end

if any(strcmp(varargin,'labels'))
    atlas_labels=varargin{find(strcmp(varargin,'labels'))+1};
end

if any(strcmp(varargin,'num_voxels'))
    num_vox=varargin{find(strcmp(varargin,'num_voxels'))+1};
end

if any(strcmp(varargin,'iterations'))
    num_iterations=varargin{find(strcmp(varargin,'iterations'))+1};
end

%% adjust number of voxels to sample based on size of parcels

for n=1:num_parcels %for each parcel
    mask=parcel_obj;
    mask.dat=mask.dat(:,n);
    nf(n)=length(find(mask.dat>0));
end

num_vox(num_vox>(min(nf)))=[];
%% test with increasing number of features randomly sampled from each parcel
num_feats_within = nan(num_parcels,length(num_vox)+1);
pred_outcome_r=nan(num_iterations,num_parcels,length(num_vox)+1);
for n=1:num_parcels %for each parcel
    mask=parcel_obj;
    mask.dat=mask.dat(:,n);
    masked_data=apply_mask(data_obj,mask);

    num_feats_within(n,:)=[num_vox size(masked_data.dat,1)];

    for it=1:num_iterations
        temp_inds=randperm(size(masked_data.dat,1));
        cs=1;
        for nf=num_feats_within(n,:)
            try
                subsampled_masked_data=masked_data;

                if nf==size(subsampled_masked_data.dat,1) %if subsampling all data, do nothing
                else
                    subsampled_masked_data.dat=subsampled_masked_data.dat(temp_inds(1:nf),:); %otherwise grab first 1:nf reordered features
                end

                if ~any(strcmp(varargin,'nfolds')) %if cv not specified do 10-fold
                    [~, stats]=predict(subsampled_masked_data,varargin{:},'nfolds',10);
                else
                    [~, stats]=predict(subsampled_masked_data,varargin{:});
                end

                try
                    pred_outcome_r(it,n,cs)=stats.pred_outcome_r;

                catch
                    pred_outcome_r(it,n,cs)=1-stats.cverr;

                end
                cs=cs+1;
            catch
            end

        end
    end
end

%% random samples from entire brain

num_feats_full=[num_vox size(data_obj.dat,1)];
pred_outcome_r_full=nan(num_iterations,length(num_vox)+1);

for it=1:num_iterations
    temp_inds=randperm(size(data_obj.dat,1));
    cs=1;
    for nf=num_feats_full
        try
            subsampled_masked_data=data_obj;

            if nf==size(subsampled_masked_data.dat,1) %if subsampling all data, do nothing
            else
                subsampled_masked_data.dat=subsampled_masked_data.dat(temp_inds(1:nf),:); %otherwise grab first 1:nf reordered features
            end

            if ~any(strcmp(varargin,'nfolds')) %if cv not specified do 10-fold
                [~, stats]=predict(subsampled_masked_data,varargin{:},'nfolds',10);
            else
                [~, stats]=predict(subsampled_masked_data,varargin{:});
            end

            try
                pred_outcome_r_full(it,cs)=stats.pred_outcome_r;

            catch
                pred_outcome_r_full(it,cs)=1-stats.cverr;

            end

            cs=cs+1;
        catch
        end

    end
end

%% random samples across all parcels

mask=parcel_obj;

mask.dat=double(any(mask.dat')');

masked_data=apply_mask(data_obj,mask);
num_feats_all_parcels=[num_vox size(masked_data.dat,1)];
pred_outcome_r_all_parcels=nan(num_iterations,length(nf));

for it=1:num_iterations
    temp_inds=randperm(size(masked_data.dat,1));
    cs=1;
    for nf=num_feats_all_parcels
        try
            subsampled_masked_data=masked_data;

            if nf==size(subsampled_masked_data.dat,1) %if subsampling all data, do nothing
            else
                subsampled_masked_data.dat=subsampled_masked_data.dat(temp_inds(1:nf),:); %otherwise grab first 1:nf reordered features
            end

            if ~any(strcmp(varargin,'nfolds')) %if cv not specified do 10-fold
                [~, stats]=predict(subsampled_masked_data,varargin{:},'nfolds',10);
            else
                [~, stats]=predict(subsampled_masked_data,varargin{:});
            end

            try
                pred_outcome_r_all_parcels(it,cs)=stats.pred_outcome_r;

            catch
                pred_outcome_r_all_parcels(it,cs)=1-stats.cverr;

            end

            cs=cs+1;
        catch
        end

    end
end

%% do curve fitting and make plots

if do_plot

    [fitresult{1}, gof{1},x_mean_wb,y_mean_wb] = createFit(num_feats_full(:)', mean(pred_outcome_r_full),'color',[.5 .5 .5],'linewidth',2);
    [fitresult{2}, gof{2},x_mean_parcels,y_mean_parcels] = createFit(num_feats_all_parcels(:)', mean(pred_outcome_r_all_parcels),'color',[.1 .1 .1],'linewidth',2);

    for r=1:num_parcels

        mean_r_within=squeeze(mean(pred_outcome_r(:,r,1:size(pred_outcome_r_all_parcels,2))));
        [fitresult{2+r}, gof{2+r},x_mean_parcel{r},y_mean_parcel{r}] = createFit(num_feats_within(r,:)', mean_r_within(:),'color',mycolors{r},'linewidth',2);

    end

    for it=1:num_iterations
        [~, ~,x_wb(it,:),y_wb(it,:)] = createFit(num_feats_full(:)', pred_outcome_r_full(it,:),'color',[.5 .5 .5],'linewidth',2);
        [~, ~,x_all_parcels(it,:),y_all_parcels(it,:)] = createFit(num_feats_all_parcels(:)', pred_outcome_r_all_parcels(it,:),'color',[.1 .1 .1],'linewidth',2);

        for r=1:num_parcels

            r_within=squeeze((pred_outcome_r(it,r,1:size(pred_outcome_r_all_parcels,2))));
            [~, ~,x{it,r},y{it,r}] = createFit(num_feats_within(r,:)', r_within(:),'color',mycolors{r},'linewidth',2);

        end
    end

    figure;
    hold on;
    boundedline(mean(x_wb),mean(y_wb), std(y_wb),'alpha', 'cmap',[.5 .5 .5]);
    boundedline(mean(x_all_parcels), mean(y_all_parcels), std(y_all_parcels),'alpha', 'cmap',[.1 .1 .1]);

    for r=1:num_parcels
        clear to_plot_x to_plot_y
        for i=1:size(x,1)
            to_plot_x(i,:)=x{i,r};
            to_plot_y(i,:)=y{i,r};
        end
        boundedline(mean(to_plot_x), mean(to_plot_y), std(to_plot_y),'alpha', 'cmap',mycolors{r});

    end
    h=findobj(gca,'type','line');
    set(h,'linewidth',2);
    set(gca,'XScale','log')
    grid on
    legend(h(1:end),{atlas_labels{:} 'All Parcels' 'Whole-brain'});

    xlabel 'Number of voxels'
    ylabel 'Model performance'

    % add brain render to plot

    axis tight;
    if render_brain

        clear o2
        o2=fmridisplay();
        o2 = surface(o2, 'axes', [0.5 0 .5 .5], 'direction', 'surface left', 'orientation', 'lateral','trans');
        o2 = surface(o2, 'axes', [0.5 .5 .4 .4], 'direction', 'surface right', 'orientation', 'lateral');

        pname = 'R.pial_MSMAll_2_d41_WRN_DeDrift.32k.mat'; % from Glasser_et_al_2016_HCP

        try
            load(which(pname))
        catch
            error('Unable to find ''R.pial_MSMAll_2_d41_WRN_DeDrift.32k.mat'' please check your path');
        end

        set(o2.surface{1}.object_handle(1),'FaceVertexCData',cdata)
        set(o2.surface{2}.object_handle(1),'FaceVertexCData',cdata)
        colormap gray

        %add mask regions
        for i=1:size(parcel_obj.dat,2)
            tv=parcel_obj;
            tv.dat=tv.dat(:,i);

            o2=addblobs(o2,region(tv),'splitcolor',{mycolors{i} mycolors{i} mycolors{i} mycolors{i}});

        end

        %clear brainstem
        delete(o2.surface{1}.object_handle(2))
        delete(o2.surface{2}.object_handle(2))

        %position axes
        ha=findobj(gcf,'type','axes');
        newpos=[[.315 .1 .3 .3];[.61 .1 .3 .3]];
        for i=1:2
            set(ha(i),'Position',newpos(i,:))
        end
        set(ha([1]),'view',[-90 0]);
        set(ha([2]),'view',[90 0]);

    end
end

%%  stats across regions
random_subsampled_outcomes=cat(2,permute(cat(3,pred_outcome_r_full,pred_outcome_r_all_parcels),[1,3,2]),pred_outcome_r);

for i=1:size(random_subsampled_outcomes,2)
    for j=1:size(random_subsampled_outcomes,2)

        if i<j

            b_mean = squeeze(nanmean(random_subsampled_outcomes(:,i,size(random_subsampled_outcomes,3)-1)- random_subsampled_outcomes(:,j,size(random_subsampled_outcomes,3)-1)));
            b_ste = squeeze(nanstd(random_subsampled_outcomes(:,i,size(random_subsampled_outcomes,3)-1) - random_subsampled_outcomes(:,j,size(random_subsampled_outcomes,3)-1)));
            b_ste(b_ste == 0) = Inf;
            b_Z(i,j) = b_mean ./ b_ste;
            b_P(i,j) = 2*normcdf(-1*abs(b_Z(i,j)),0,1);

        end
    end
end

stats.pairwise_comparisons.Z = b_Z;
stats.pairwise_comparisons.P = b_P;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Locate place fields in a 1D firing map.
%
% Locates the place fields and calculates the field sizes and spacing
% between fields.
%
%  USAGE
%   [fieldsCurve, fields, spacing] = analyses.placefield1D(map, <options>)
%   map           Firing rate map either structure obtained using analyses.map
%                 or a matrix that represents firing map. Note that if you want to filter place fields based on number
%                 of spikes, then you should use map structure.
%   <options>     optional list of property-value pairs (see table below)
%
%   ======================================================================================
%    Properties     Values
%   --------------------------------------------------------------------------------------
%    'thresholdType' Type of the firing rate threshold. '%' specifies that 'threshold'
%                    value is in range [0; 1] and rates lower than 'threshold * curveAmplitude'
%                    are not considered as place fields. curveAmplitude is the range of a firing
%                    curve, i.e. max(curve) - min(curve).
%                    'rate' specifies that 'threshold' value is an absolute rate. Values lower than
%                    'threshold' are not considered as place fields. Default it '%'.
%
%    'threshold'    Threshold for a firing rate. How the value is treated depends on options
%                   'thresholdType'. Default = 0.2 and 'thresholdType' = '%'.
%
%    'minRows'      Minimum number of row bins in a place field. Fields with
%                   fewer bins are not considered as place fields. Remember to
%                   adjust this value when you change the bin width. Default = 3.
%                   If minRows == 0, then minimum number of bins is not checked.
%
%    'binWidth'     Bin length in centimetres. Used to calculate size of fields in cm.
%                   Default is 1.
%
%    'minSpikes'    Minimum number of spikes in a place field. Default is 0, which means
%                   that minimum number of spikes is not checked. To use this value 'map'
%                   must be a result of function analyses.map.
%
%    'minDistance'  Minimum number of bins between adjacent fields. Distance is calculated between
%                   end of one field and start of another field. If distance is smaller or
%                   equals to 'minDistance', then two fields are merged together. Note that more
%                   than two fields can be merged if they are all adjacent. Default is 0, which
%                   means that fields should have a common border.
%
%    'debug'        'on' plot turning curve, place fields and their centres of mass. 'off' do not
%                   plot anything. Default is 'off'.
%
%    'pos'          Position samples. Used to calculate posInd. If not provided,
%                   then posInd will be an empty matrix.
%   ======================================================================================
%
%   fieldsCurve     Vector of the same size as underline firing rate curve. The elements of fieldsCurve
%                   are integer values greater or equal to 0. 0 elements correspond to background (not
%                   a place field). Values of 1 constitute first field. Values of 2 constitute second
%                   field; and so on.
%   fields          Structure with information about each field. Structure fields are:
%       col         vector of columns that constitute this field.
%       row         vector of rows (artificial). Same size as col, all ones;
%       peak        field peak;
%       x, y        field centre of mass point. See NOTES for details about centre of mass calculation;
%       meanRate    mean firing rate;
%       width       field width in cm;
%       size        size of the field. Do not rely on it! It's from old code and could be misleading.
%                   Calculated as width * <mean firing rate of a field>.
%       posInd      Indices of position samples that correspond to this field. In case position sample
%                   matrix is of size Nx5 ([t x y x1 y1]), posInd corresponds to the left most position
%                   columns ([x y]). If pos argument is not provided, then posInd will be an empty matrix.
%   spacing         Matrix of spacing between fields. Entry (i, j) equals to distance between fields i and j.
%
%  NOTE
%   Since centre of mass for a curve is somewhat undefined, it is calculated for plate bounded by a firing curve
%   and a function y = 0 for all x. See http://tutorial.math.lamar.edu/Classes/CalcII/CenterOfMass.aspx
%
function [fieldsCurve, fields, spacing] = placefield1D(mapS, varargin)
    if nargin < 1 || mod(length(varargin), 2) ~= 0
        error('BNT:numArgs', 'Incorrect number of parameters (type ''help analyses.placefield1D'' for details).');
    end

    % Default values
    threshold = 0.1;
    minRows = 2;
    binWidth = 1; % [cm];
    minSpikes = 0;
    minFieldDistance = 0; % [bins]
    debug = false;
    isThrshPercentage = true;
    pos = [];

    % Parse parameter list
    i = 1;
    while i < length(varargin)
        if ~ischar(varargin{i}),
            error(['Parameter ' num2str(i+2) ' is not a property (type ''help analyses.placefield1D'' for details).']);
        end

        switch(lower(varargin{i})),
            case 'threshold'
                threshold = varargin{i+1};
                i = i + 2;

            case 'minrows'
                minRows = varargin{i+1};
                if ~helpers.isdscalar(minRows, '>=0')
                    error('Incorrect value for property ''minRows'' (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            case 'binwidth'
                binWidth = varargin{i+1};
                if ~helpers.isdscalar(binWidth, '>0')
                    error('Incorrect value for property ''binWidth'' (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            case 'minspikes'
                minSpikes = varargin{i+1};
                if ~helpers.isdscalar(minSpikes, '>=0')
                    error('Incorrect value for property ''minSpikes'' (type ''help analyses.placefield1D'' for details).');
                end
                if minSpikes > 0 && ~isstruct(mapS)
                    error('You should only use minSpikes with a structure-like firing curve (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            case 'mindistance'
                minFieldDistance = varargin{i+1};
                if ~helpers.isdscalar(minFieldDistance, '>=0')
                    error('Incorrect value for property ''minDistance'' (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            case 'debug'
                debug = strcmpi(varargin{i+1}, 'on');
                i = i + 2;

            case 'thresholdtype'
                isThrshPercentage = strcmpi(varargin{i+1}, '%');
                if ~helpers.isstring(varargin{i+1}, '%', 'rate'),
                    error('Incorrect value for property ''thresholdType'' (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            case 'pos'
                pos = varargin{i+1};
                if ~ismatrix(pos) || size(pos, 2) < 2
                    error('Incorrect value for property ''pos'' (type ''help analyses.placefield1D'' for details).');
                end
                i = i + 2;

            otherwise,
                error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help analyses.placefield1D'' for details).']);
        end
    end
    if isThrshPercentage
        if ~helpers.isdscalar(threshold, '>=0', '<=1')
            error('Incorrect value for property ''threshold'' (type ''help analyses.placefield1D'' for details).');
        end
    else
        if ~helpers.isdscalar(threshold, '>=0')
            error('Incorrect value for property ''threshold'' (type ''help analyses.placefield1D'' for details).');
        end
    end

    if isstruct(mapS)
%         mapAxis = mapS.x;
        if size(mapS.z, 1) > 1 && size(mapS.z, 2) > 1
            error('Incorrect value for property ''map''. Firing curve is 2D, but should be a vector (type ''help analyses.placefield1D'' for details).');
        end

        curve = mapS.z;
    else
        if size(mapS, 1) > 1 && size(mapS, 2) > 1
            error('Incorrect value for property ''map''. It should be a vector (type ''help analyses.placefield1D'' for details).');
        end
%         mapAxis = 1:length(mapS);

        curve = mapS;
    end
    curveLen = length(curve);

    fieldsCurve = zeros(1, curveLen);
    fields = struct('row', {}, 'col', {}, ...
        'size', {}, 'peak', {}, ...
        'x', {}, 'y', {}, ...
        'meanRate', {}, 'width', {}, ...
        'posInd', {} ...
        );
    spacing = [];

    if curveLen == 0
        return
    end

    % extend curve, so that findpeaks function finds all peaks
    map = zeros(1, curveLen + 2);
    map(2:2+curveLen-1) = curve;

    [~, locs] = findpeaks(map);
    dMap = diff(map);
    dMap(end+1) = dMap(end);

    if isThrshPercentage
        threshold = min(curve) + ((max(curve) - min(curve)) * threshold);
    end

    startField = zeros(1, 100);
    stopField = zeros(1, 100);
    fieldWidth = zeros(1, 100);
    fieldSize = zeros(1, 100);
    curField = 1;

    for i = 1:length(locs)
        if map(locs(i)) < threshold
            continue;
        end
        curStopField = locs(i);
        while curStopField < length(dMap) && dMap(curStopField) < 0 && map(curStopField) >= threshold
            curStopField = curStopField + 1;
        end
        % substract one because of the loop
        curStopField = curStopField -1;

        curStartField = locs(i) - 1;
        while curStartField > 0 && dMap(curStartField) > 0 && map(curStartField) >= threshold
            curStartField = curStartField - 1;
        end
        if curStartField == 0
            curStartField = 1;
        end
        curStopField = curStopField - 1; % substract one cause of map -> curve translation

        if curStopField > length(curve)
            curStopField = length(curve);
        end

        startField(curField) = curStartField;
        stopField(curField) = curStopField;

        fieldWidth(curField) = (curStopField - curStartField + 1) * binWidth;
        fieldSize(curField) = fieldWidth(curField) * nanmean(curve(curStartField:curStopField));
        curField = curField + 1;
    end

    startField(curField:end) = [];
    stopField(curField:end) = [];
    fieldWidth(curField:end) = [];
    fieldSize(curField:end) = [];

    nFields = length(fieldWidth);
    removeField = false(1, nFields);

    fieldsCurve = zeros(1, curveLen);
    if nFields == 0
        fields = struct();
        fields(1) = [];
        spacing = [];
        return;
    end

    % merge neighbouring fields
    curField = 1;
    for i = 1:nFields
        if i + 1 > nFields
            break;
        end

        fieldsDist = abs(startField(i+1) - stopField(i));
        if fieldsDist <= minFieldDistance
            stopField(curField) = stopField(i+1);
            removeField(i+1) = true;
        else
            curField = i + 1;
        end
    end
    startField(removeField) = [];
    stopField(removeField) = [];
    fieldWidth(removeField) = [];
    fieldSize(removeField) = [];

    if ~isempty(find(removeField, 1))
        % adjust fields width and size after merge
        nFields = length(startField);
        for i = 1:nFields
            fieldWidth(i) = (stopField(i) - startField(i) + 1) * binWidth;
            fieldSize(i) = fieldWidth(i) * nanmean(curve(startField(i):stopField(i)));
        end
    end

%     fields(nFields) = struct('row', {}, 'col', {}, ...
%         'size', {}, 'peak', {}, ...
%         'x', {}, 'y', {}, ...
%         'meanRate', {}, 'width', {}, ...
%         'posInd', {} ...
%         );
%     fields(nFields) = struct('row', [], 'peak', 0, 'x', -1, 'y', -1, 'meanRate', -1, 'posInd', []);

    curField = 1;
    removeField = false(1, nFields);
    for i = 1:nFields
        fieldData = curve(startField(i):stopField(i));
        peak = max(fieldData);

        % check field peak rate
        if peak < threshold
            removeField(i) = true;
            continue;
        end

        % Number of spikes in the field
        if minSpikes > 0
            numSpikes = sum(mapS.Nspikes(startField(i):stopField(i)));
            if numSpikes < minSpikes
                % Marks field as having to few spikes
                removeField(i) = true;
                continue;
            end
        end

        % check number of bins
        if length(startField(i):stopField(i)) <= minRows
            removeField(i) = true;
            continue;
        end

        fields(curField).col = startField(i):stopField(i);
        fields(curField).row = ones(1, length(fields(curField).col));
        fields(curField).peak = peak;
        fields(curField).meanRate = nanmean(fieldData);
        fields(curField).width = fieldWidth(i);
        fields(curField).size = fieldSize(i);
        fieldsCurve(startField(i):stopField(i)) = curField;

        if ~isempty(pos)
            if isstruct(mapS)
                xSpace = mapS.x;
            else
                nBins = length(curve);
                limitsX = [nanmin(pos(:, bntConstants.PosX)) nanmax(pos(:, bntConstants.PosX))];
                xSpace = linspace(limitsX(1), limitsX(2), nBins);
            end
            xMin = xSpace(startField(i));
            xMax = xSpace(stopField(i));
            posIndX = pos(:, bntConstants.PosX) >= xMin & pos(:, bntConstants.PosX) <= xMax;
            fields(curField).posInd = find(posIndX);
        end

        curField = curField + 1;
    end

    fields(curField:end) = [];
    nFields = length(fields);

    if debug && nFields > 0
        figure, plot(curve);
        v = axis();
        ymin = v(3);
        ymax = v(4);
        color = [1 0 0];
        fieldSize = zeros(1, nFields);
        for i = 1:nFields
            startField = fields(i).col(1);
            stopField = fields(i).col(end);

            p = patch([startField startField stopField stopField], [ymax ymin ymin ymax], color);
            set(p, 'FaceAlpha', 0.5);
            color(1) = color(1) - 0.1;
            color(2) = color(2) + 0.2;
%             color(3) = color(3) + 0.2;
            if color(1) < 0
                color(1) = 0;
            end
            if color(2) > 1
                color(2) = 1;
            end

            fieldSize(i) = fields(i).size;
        end
        xlabel('Bins'), ylabel('Rate');
        titleStr = sprintf('Field(s) with size(s) %s', sprintf('%f, ', fieldSize));
        titleStr(end-1:end) = [];
        title(titleStr);
    end

    % Calculate the spacing using centre of mass (COM). Since centre of mass for a curve is somewhat
    % undefined I'll calculate COM for palte bounded by curve and a function y = 0 for all x.
    % http://tutorial.math.lamar.edu/Classes/CalcII/CenterOfMass.aspx
    distMat = zeros(nFields, 2);

    % store points for debug purposes
    pointsX = zeros(1, nFields);
    pointsY = zeros(1, nFields);

    for i = 1:nFields
        x = fields(i).col(:);
        y = curve(x);
        y = y - min(y); % make it from 0. This gives better estimation

        y1 = x .* y;
        y2 = (y.^2)/2;
        A = trapz(x, y);

        Mx = trapz(x, y1)/A;
        My = min(curve(x)) + trapz(x, y2)/A;

        pointsX(i) = Mx;
        pointsY(i) = My;

        fields(i).x = Mx;
        fields(i).y = My;

        distMat(i, :) = [Mx My];

%             % Raymond's code
%             posMass = sum(mapAxis(x).*curve(x));
%             mass = sum(curve(x));
%             com(i) = posMass/mass;
    end

    if debug
        hold on;
        plot(pointsX, pointsY, 'o');
    end

    if nFields <= 1
        spacing = [];
    else
        D = pdist(distMat);
        spacing = squareform(D);
    end
end

author: kavli-ntnu


text:
function results = vl_test_kmeans(varargin)
% VL_TEST_KMEANS

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

vl_test_init ;

function s = setup()
randn('state',0) ;
s.X = randn(128, 100) ;

function test_basic(s)
[centers, assignments, en] = vl_kmeans(s.X, 10, 'NumRepetitions', 10) ;
[centers_, assignments_, en_] = simpleKMeans(s.X, 10) ;
assert(en_ <= 1.1 * en, 'vl_kmeans did not optimize enough') ;

function test_algorithms(s)
distances = {'l1', 'l2'} ;
dataTypes = {'single','double'} ;

for dataType = dataTypes
  for distance = distances
    distance = char(distance) ;
    conversion = str2func(char(dataType)) ;
    X = conversion(s.X) ;
    vl_twister('state',0) ;
    [centers, assignments, en] = vl_kmeans(X, 10, ...
                                           'NumRepetitions', 1, ...
                                           'MaxNumIterations', 10, ...
                                           'Algorithm', 'Lloyd', ...
                                           'Distance', distance) ;
    vl_twister('state',0) ;
    [centers_, assignments_, en_] = vl_kmeans(X, 10, ...
                                              'NumRepetitions', 1, ...
                                              'MaxNumIterations', 10, ...
                                              'Algorithm', 'Elkan', ...
                                              'Distance', distance) ;
    vl_assert_almost_equal(centers, centers_, 1e-5) ;
    vl_assert_almost_equal(assignments, assignments_, 1e-5) ;
    vl_assert_almost_equal(en, en_, 1e-5) ;
  end
end

function test_patterns(s)
distances = {'l1', 'l2'} ;
dataTypes = {'single','double'} ;
for dataType = dataTypes
  for distance = distances
    distance = char(distance) ;
    conversion = str2func(char(dataType)) ;
    data = [1 1 0 0 ;
            1 0 1 0] ;
    data = conversion(data) ;
    [centers, assignments, en] = vl_kmeans(data, 4, ...
                                           'NumRepetitions', 100, ...
                                           'Distance', distance) ;
    assert(isempty(setdiff(data', centers', 'rows'))) ;
  end
end

function [centers, assignments, en] = simpleKMeans(X, numCenters)
[dimension, numData] = size(X) ;
centers = randn(dimension, numCenters) ;

for iter = 1:10
  [dists, assignments] = min(vl_alldist(centers, X)) ;
  en = sum(dists) ;
  centers = [zeros(dimension, numCenters) ; ones(1, numCenters)] ;
  centers = vl_binsum(centers, ...
                      [X ; ones(1,numData)], ...
                      repmat(assignments, dimension+1, 1), 2) ;
  centers = centers(1:end-1, :) ./ repmat(centers(end,:), dimension, 1) ;
end

author: Cloud-CV

repo: object-proposals


text:
%% Create an artificial set of QRS detections, based on multilead QRS detections created with wavedet algorithm.
% This is an auxiliar function for the ECGtask_QRS_detections_post_process.
% *wavedetMix* use the output of wavedet algorithm to perform a  multilead
% composition. The result of this algorithm is a new set of detections
% based on the concatenation of the ""best"" detections found for each
% 20-seconds window in a recording. So, this algorithm generates *new* QRS
% detection series, as described in [add reference].
%
%
%   all_detections = wavedetMix(struct_in, ECG_header, start_end_this_segment )
%
% Arguments:
%
%      + struct_in: the structure which results from loading the result of
%                   the QRS detection task used for invoking wavedet:
%
%                   cached_filenames = ECG_w.GetCahchedFileName('QRS_detection');
%                   struct_in = load(cached_filenames{1});
%
%      +ECG_header: [struct] OPTIONAL.
%
%             Description of the ECG typically available in the
%             ECG_header. Structure with fields:
%
%               -freq: Sampling rate in Hz. (1)
%
%               -nsig: Number of ECG leads. (size(ECG,2))
%
%               -nsamp: Number of ECG samples. (size(ECG,1))
%
%      + start_end_this_segment: an array with the first and last sample
%                                indexes.
%
% Output:
%
%      + all_detections: struct with the artificial detections named with
%                        prefix "wavedetMix_ECGmix".
%
% Example
%
%     ECG_w.ECGtaskHandle = 'QRS_detections_post_process';
%     % Mix QRS detections strategy function
%     ECG_w.ECGtaskHandle.post_proc_func = 'wavedet_QRS_detection_mix';
%     ECG_w.ECGtaskHandle.payload = load(cached_filenames{1});
%     ECG_w.ECGtaskHandle.CalculatePerformance = true;
%     ECG_w.Run;
%
% See also ECGtask_QRS_detections_post_process, mixartif, best_m_lead, qrs_detection_and_correction
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Birthdate: 01/01/2014
% Last update: 14/07/2015
% Copyright 2008-2015
%
function all_detections = wavedetMix(struct_in, ECG_header, start_end_this_segment )

    all_detections = [];

    % attemp to build a better detection from single-lead detections.

    AnnNames = struct_in.series_quality.AnnNames(:,1);
    detector_name = cellfun( @(a)( strtok(a, '_')), AnnNames, 'UniformOutput', false);

    aux_idx = find(strcmpi(detector_name, 'wavedet'));

    if( isempty(aux_idx) )
        return
    end

    all_annotations = {};
    for ii = rowvec(aux_idx)
        all_annotations = [all_annotations; {struct_in.(struct_in.series_quality.AnnNames{ii,1}).time}];
    end

    [ ratios, estimated_labs ] = CalcRRserieRatio(all_annotations, ECG_header, start_end_this_segment);

    [~, best_detections_idx] = sort(ratios, 'descend');

    % generate artificial annotations combining K best annotations
    aux_idx = best_detections_idx(1:min(10, length(best_detections_idx)) );
    artificial_annotations = combine_anns(all_annotations(aux_idx), estimated_labs(aux_idx), ECG_header );

    for ii = 1:length(artificial_annotations)
        aux_str = ['wavedetMix_ECGmix' num2str(ii)];
        all_detections.(aux_str) = artificial_annotations(ii);
    end

author: marianux

repo: ecg-kit

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


text:
function varargout = FA2(varargin)
% FA2 M-file for FA2.fig
%      FA2, by itself, creates a new FA2 or raises the existing
%      singleton*.
%
%      H = FA2 returns the handle to a new FA2 or the handle to
%      the existing singleton*.
%
%      FA2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FA2.M with the given input arguments.
%
%      FA2('Property','Value',...) creates a new FA2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before FA2_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to FA2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help FA2

% Last Modified by GUIDE v2.5 07-Feb-2008 06:51:04

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @FA2_OpeningFcn, ...
                   'gui_OutputFcn',  @FA2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before FA2 is made visible.
function FA2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to FA2 (see VARARGIN)

% Choose default command line output for FA2
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes FA2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = FA2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)

set(handles.panel2,'visible','on');
set(handles.panel3,'visible','on');
set(handles.panel4,'visible','off');

% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)

set(handles.panel2,'visible','on');
set(handles.panel4,'visible','on');
set(handles.panel3,'visible','off');
set(handles.text23,'visible','on');

% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
set(handles.panel4,'visible','off');
global x ;
x=linspace(-10,10,101);
a=zeros(1,4);
v=str2num(get(handles.edit1,'string'));
a=[v];
if isnumeric(a)
     set(handles.text13,'visible','on');
end
if length(a) < 4 || length(a) > 4
     set(handles.text12,'visible','on');
end

if a(1,1)==a(1,2) && a(1,2)==a(1,3) && a(1,3)==a(1,4)
     set(handles.text11,'visible','on');
end

if a(1,1)==a(1,2) || a(1,1)==a(1,3) || a(1,1)==a(1,4)
    set(handles.text18,'visible','on');
end

if a(1,1) > a(1,2)
    set(handles.text15,'visible','on');
end
if a(1,3) > a(1,4)
    set(handles.text14,'visible','on');
end

if a(1,1) < -10 || a(1,1) > 10 || a(1,2) < -10 || a(1,2) > 10 || a(1,3) < -10 || a(1,3) > 10 || a(1,4) < -10 || a(1,4) > 10

  set(handles.text10,'visible','on');
end
global A;
A=trapmf(x,a);

% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
global x;
x=linspace(-10,10,101);
a=zeros(1,3);
v=str2num(get(handles.edit2,'string'));
a=[v];
if isnumeric(a)
     set(handles.text13,'visible','on');
end
if length(a) < 3 || length(a) > 3
     set(handles.text17,'visible','on');
end

if a(1,1)==a(1,2) && a(1,2)==a(1,3)
     set(handles.text11,'visible','on');
end

if a(1,1)==a(1,2) || a(1,1)==a(1,3)
    set(handles.text18,'visible','on');
end

if a(1,1) > a(1,2)
    set(handles.text15,'visible','on');
end
if a(1,2) > a(1,3)
    set(handles.text24,'visible','on');
end

global A;
A=trimf(x,a);

% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: Sable

repo: mcbench-benchmarks


text:
function M = mtkronprod(T,U,n,transpose)
%MTKRONPROD Compute a matricized tensor Kronecker product.
%   mtkronprod(T,U,n) computes the product
%
%      tens2mat(T,n)*conj(kron(U([end:-1:n+1 n-1:-1:1])))
%
%   without permuting the tensor T. Note that for improved performance, it
%   is advisable for the largest two dimensions of T to be the first and
%   last modes of T.
%
%   mtkronprod(T,U,n,'T') and mtkronprod(T,U,n,'H') transpose or complex
%   conjugate transpose the matrices U{n} before computing the matricized
%   tensor Kronecker product.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 4, transpose = 0; end
if ischar(transpose)
    if strcmp(transpose,'T'), transpose = 1;
    elseif strcmp(transpose,'H'), transpose = 1i; end
end
U = U(:).';
if transpose == 0, size_tens = cellfun('size',U,1);
else size_tens = cellfun('size',U,2); end
data = T;
if isstruct(T)
    % Incomplete or sparse tensor.
    T = data.matrix;
    alloc = @sparse;
else
    % Full tensor.
    alloc = @zeros;
end
N = length(size_tens);
if transpose == 0, size_core = cellfun('size',U,2);
else size_core = cellfun('size',U,1); end
if n > 0, size_core(n) = size_tens(n); end
R = prod(size_core([1:n-1 n+1:N]));

% Determine if large-scale version of the algorithm should be executed.
ratio = size_core./size_tens;
perm = zeros(1,N);
l = 1; r = N;
for i = 1:N
    if r == n || (l < n && ratio(l) < ratio(r)), perm(i) = l; l = l+1;
    else perm(i) = r; r = r-1;
    end
end
if n == 0, mem = 8*prod(size_tens);
else mem = 8*prod(size_tens([1:perm(1)-1 perm(1)+1:N]))*size_core(perm(1));
end
largescale = mem > 2e9 || (isstruct(data) && isempty(data.matrix));
if largescale && ~isstruct(data)
    error('mtkronprod:mem',['Intermediate result is too large. Try ' ...
        'to format data with fmt first so that the large-scale ' ...
        'version of the algorithm can be applied.']);
end

if largescale

    % The large scale implementation fires if the required memory would
    % otherwise be larger than 2GB and assumes the dataset has been
    % formatted by fmt.
    idx = [1:n-1 n+1:N];
    jdx = cell(1,N);
    switch transpose
        case 0,  U = cellfun(@(u)conj(u),U,'UniformOutput',false);
        case 1,  U = cellfun(@(u)u',U,'UniformOutput',false);
        case 1i, U = cellfun(@(u)u.',U,'UniformOutput',false);
    end
    if n == 0, M = zeros(1,R);
    else M = zeros(size_tens(n),R); end
    for r = 1:R
        s = data.val;
        [jdx{:}] = ind2sub(size_core([1:n-1 n+1:N]),r);
        for j = 1:length(idx)
            s = s.*U{idx(j)}(data.sub{idx(j)},jdx{j});
        end
        if n == 0, M(r) = sum(s);
        else M(:,r) = accumarray(double(data.sub{n}),s,[size_tens(n) 1]);
        end
    end

else

    % Apply structure-exploiting matricized tensor Kronecker product.
    M = T;
    cpl = cumprod([1 size_core(perm(perm < n))]); l = 1;
    cpr = cumprod([1 size_core(perm(perm > n))]); r = 1;
    for i = 1:length(perm)

        mode = perm(i);
        if mode < n

            tmp = reshape(M,cpl(l)*size_tens(mode),[]);
            M = alloc(cpl(l),size(tmp,2)*size_core(mode));
            for j = 1:cpl(l)
                idx = j:cpl(l):size(tmp,1);
                switch transpose
                    case 0,  tmp2 = U{mode}'*tmp(idx,:);
                    case 1,  tmp2 = conj(U{mode})*tmp(idx,:);
                    case 1i, tmp2 = U{mode}*tmp(idx,:);
                end
                M(j,:) = tmp2(:);
            end
            l = l+1;

        elseif mode > n

            tmp = reshape(M,[],cpr(r)*size_tens(mode));
            M = alloc(size(tmp,1)*size_core(mode),cpr(r));
            for j = 1:cpr(r)
                idx = (1:size_tens(mode))+(j-1)*size_tens(mode);
                switch transpose
                    case 0,  tmp2 = tmp(:,idx)*conj(U{mode});
                    case 1,  tmp2 = tmp(:,idx)*U{mode}';
                    case 1i, tmp2 = tmp(:,idx)*U{mode}.';
                end
                M(:,j) = tmp2(:);
            end
            r = r+1;

        end

    end

    % Permute and reshape output.
    if n > 1
        if issparse(M)
            idx = find(M);
            [i,j,k] = ind2sub([cpl(end) size_tens(n) cpr(end)],idx);
            ik = sub2ind([cpl(end) cpr(end)],i,k);
            M = sparse(j,ik,full(M(idx)),size_tens(n),cpl(end)*cpr(end));
        else
            M = reshape(M,[cpl(end) size_tens(n) cpr(end)]);
            M = reshape(permute(M,[2 1 3]),size_tens(n),[]);
        end
    else
        if n == 0, M = reshape(M,1,[]);
        else M = reshape(M,size_tens(n),[]); end
    end

end

end

author: andrewssobral

repo: lrslibrary


text:
function [rate]=recognition(Training_matrix,Vectors,Y)

% recognition process
sucess=0;
error=0;

psi=mean(Training_matrix,2);

%--------------recognition from s01-s40-------------------------------
for i=1:40
   for j=4:10
     test=double(imread(['D:\face database\orl_faces(olivetti)\s' num2str(i) '\' num2str(j) '.bmp' ]));
     test=extract_fea(test);
     test=reshape(test,size(test,1)*size(test,2),1);
     %Projecting the test image onto the eigenface space...
     newomega=Vectors'*(test-psi);
     %Calculating the distance of the new face from each sample..NN
     %classifier
     distances = dist(Y',newomega);
     %Choosing the nearest sample : the smallest distance...
     indice=find(distances==min(distances));
     k=indice(1);  % k is the label of the class,correspond the column of the trainingmatrix.
     if mod(k,3)==0
         if (k/3==i)
             sucess=sucess+1;
         else
             error=error+1;
         end
     else
         if (floor(k/3)+1)==i
             sucess=sucess+1;
         else
           error=error+1;
         end
     end
end
end

rate=sucess/(sucess+error);

author: sfvsfv

repo: Mathematical-modeling


text:
function tec_write ( tec_file_name, dim_num, node_num, element_num, ...
  element_order, node_data_num, node_coord, element_node, node_data )

%*****************************************************************************80
%
%% TEC_WRITE writes finite element data to a TEC file.
%
%  Discussion:
%
%    This program writes the node, element and data files that define
%    a finite element geometry and data based on that geometry:
%    * a set of nodes,
%    * a set of elements based on those nodes,
%    * a set of data values associated with each node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string TEC_FILE_NAME, the name of the TEC file.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer NODE_DATA_NUM, the number of data items per node.
%
%    Input, real NODE_COORD(DIM_NUM,NODE_NUM), the coordinates of nodes.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    the global index of local node I in element J.
%
%    Input, real NODE_DATA(NODE_DATA_NUM,NODE_NUM), the data values associated
%    with each node.
%
  tec_file_unit = fopen ( tec_file_name, 'wt' );
%
%  Write the title.
%
  fprintf ( tec_file_unit, 'TITLE = "%s"\n', tec_file_name );
%
%  Write the variable names.
%
  fprintf ( tec_file_unit, 'VARIABLES = ' );

  name = 'X';
  for dim = 1 : dim_num
    if ( dim == 1 )
      fprintf ( tec_file_unit, '"%s"', name );
    else
      fprintf ( tec_file_unit, ', "%s"', name );
    end
    name = s_inc ( name );
  end

  name = 'data_001';
  for dim = 1 : node_data_num
    fprintf ( tec_file_unit, ', "%s"', name );
    name = file_name_inc ( name );
  end

  fprintf ( tec_file_unit, '\n' );
%
%  Write the ZONE record.
%
  if ( dim_num == 2 & element_order == 3 )
    zonetype = 'FETRIANGLE';
  elseif ( dim_num == 2 & element_order == 4 )
    zonetype = 'FEQUADRILATERAL';
  elseif ( dim_num == 3 & element_order == 4 )
    zonetype = 'FETETRAHEDRON';
  elseif ( dim_num == 3 & element_order == 8 )
    zonetype = 'FEBRICK';
  else
    zonetype = 'FEUNKNOWN';
  end

  fprintf ( tec_file_unit, 'ZONE N = %d, E = %d, ', node_num, element_num );
  fprintf ( tec_file_unit, 'DATAPACKING = POINT, ZONETYPE = %s\n', ...
    zonetype );
%
%  Write the node coordinates and node data.
%
  for node = 1 : node_num
    for dim = 1 : dim_num
      fprintf ( tec_file_unit, '  %10f', node_coord(dim,node) );
    end
    for data = 1 : node_data_num
      fprintf ( tec_file_unit, '  %10f', node_data(data,node) );
    end
    fprintf ( tec_file_unit, '\n' );
  end
%
%  Write the element-node connectivity.
%
  for element = 1 : element_num
    for order = 1 : element_order
      fprintf ( tec_file_unit, '  %6d', element_node(order,element) );
    end
    fprintf ( tec_file_unit, '\n' );
  end

  fclose ( tec_file_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEC_WRITE wrote all data to "%s".\n', ...
    tec_file_name );

  return
end

author: johannesgerer

repo: jburkardt-m


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

curr_dir = pwd;
cd(model_dir);
model_obj_rec.net = caffe_load_model(model_obj_rec.net_def_file, model_obj_rec.net_weights_file, model_phase);

if use_detection_svms
    assert(exist(model_obj_rec.svm_weights_file, 'file')>0);
    weights = read_svm_detection_weights(model_obj_rec.svm_weights_file);
    model_obj_rec.net = set_net_layer_weights(model_obj_rec.net, model_obj_rec.svm_layer_name, weights);
end

cd(curr_dir);
end

author: gidariss

repo: mrcnn-object-detection


text:
%GRIDROTATE Rotate grid.
%
%   [ GRID ] = GRIDROTATE( GRID, TH, AX ) Applies rotation angle TH
%   (radians) to grid points around the origin (in 3D around axis AX:
%   1 = x-axis (default), 2 = y-axis, 3 = z-axis). The input argument
%   GRID can either be a grid struct or simply an array of
%   coordinates. Normals will also be recalculated in the case of a
%   grid struct with boundary information.
%
%   See also GRIDEXTRUDE, GRIDMERGE, GRIDREVOLVE, GRIDSCALE

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
classdef RWMOP11 < PROBLEM
% <multi> <real> <constrained>
% Water resource management problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 5;
            obj.D        = 3;
            obj.lower    = [0.01 0.01 0.01];
            obj.upper    = [0.45 0.1 0.1];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1);
            x2 = x(:,2);
            x3 = x(:,3);
            % Objectives
            f(:,1) = 106780.37 .* (x2 + x3) + 61704.67 ;
            f(:,2) = 3000 .* x1 ;
            f(:,3) = 305700 .* 2289 .* x2 ./ power(0.06.*2289, 0.65) ;
            f(:,4) = 250 .* 2289 .* exp(-39.75.*x2+9.9.*x3+2.74) ;
            f(:,5) = 25 .* (1.39 ./(x1.*x2) + 4940.*x3 -80) ;
            % Constraints
            g(:,1) = 1 - (0.00139./(x1.*x2)+4.94.*x3-0.08);
            g(:,2) = 1 - (0.000306./(x1.*x2)+1.082.*x3-0.0986);
            g(:,3) = 50000 - (12.307./(x1.*x2) + 49408.24.*x3+4051.02);
            g(:,4) = 16000 - (2.098./(x1.*x2)+8046.33.*x3-696.71);
            g(:,5) = 10000 - (2.138./(x1.*x2)+7883.39.*x3-705.04);
            g(:,6) = 2000 - (0.417.*x1.*x2 + 1721.26.*x3-136.54);
            g(:,7) = 550 - (0.164./(x1.*x2)+631.13.*x3-54.48);
            g      = -g;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
         %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [7.3450511e+04   1.3500000e+03   2.8534690e+06   6.6200320e+06   2.5000000e+04];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function writeG2oPose(file_id, robot_idx, frame_idx, T_W_C)

frame_id = gtsamFrameId(robot_idx, frame_idx);

x = T_W_C(1, 4); y = T_W_C(2, 4); z = T_W_C(3, 4);
R = T_W_C(1:3, 1:3);
assert(sum(sum(imag(R))) == 0);
% Seems to be necessary to avoid complex quaternions
% (to avoid trace(R) < -1)
R = R / (det(R) + 1e-5);
q = rot2quat(R, 1e-4);
%quat2rot(q) - R
%assert(all(all(quat2rot(q) - R < 1e-5)));
assert(norm(q)>1e-3);
q = q/norm(q);
assert(norm(q)>1e-3);
qw = q(1); qx = q(2); qy = q(3); qz = q(4);
if (sum(imag(q)) ~= 0)
    q
    assert(false)
end

fprintf(file_id, ...
    'VERTEX_SE3:QUAT %d %f %f %f %f %f %f %f\n', ...
    frame_id, x, y, z, qx, qy, qz, qw);

end

author: uzh-rpg