function conix = cpatt_script(n_nod, ConsPattern, v_walk, dscheme, bsymm, model) 
% returns the "coloring map" of constraints Jacobian using the 
% sparsity pattern and utilizing the TOMLAB function findpatt.m
% Note1: this script was written to allow for the use of pre-computed 
%       coloring maps to decrease the computation time required to
%       compute the constraints Jacobian by finite differences avoiding
%       the occurance of a memory error.
% Note2: the logic is based on the folowing observations:
%       1) Matlab crashes for a number of nodes greater than about 70,
%       reporting "insufficient memory" when a coloring map in Prob.ConIx
%       is specified. This occurs because of an expensive computation 
%       performed in the TOMLAB FDJac.m function (in line 369);
%       2) Matlab crashes for a number of nodes grater than about 157,
%       reporting "insufficient memory" when running the TOMLAB-function
%       findpatt.m, which computes a coloring map for the constraints
%       Jacobian
% Note3: Proposed solution
%       1) Since the expensive computation in line 369 of FDJac.m is not 
%       required in our optimizations, this line is simply commented out in
%       a local version of FDJac.m. This also substantially decreases
%       computation time.
%       2) Prob.ConsDiff is set to 1 (compute Jacobian by FD) instead of 11
%       (compute Jacobian with finite differences computing coloring map
%       internally by calling findpatt.m). This allows for the computation
%       of the coloring map here, either by using findpatt.m when possible
%       (for n_nod < limit) or by constructing a coloring map based on a
%       coloring map obtained for a reduced number of nodes (for n_nod > 
%       limit).
% Note4: for a number of nodes superior to n_limit2 this function is 
%        intended to be used only when bilateral symmetry is
%        activated. It will have to be changed to allow for the use with
%        bilateral symmetry for a number of nodes greater than n_limit2
% Note5: for a number of nodes greater than n_limit2 the computation is
%        specific for the legs2d model
%      
% Author: Marko Ackermann
% Last modification: MA (01.29.2009)

% For Euler discretization scheme and legs2d model
limit1_eu = 70; % nodes upper limit for computation of line 168 in FDJac.m when coloring map is provided (70 in Marko's machine, 79 in Ton's machine)
limit2_eu = 161; % nodes upper limit for computation of coloring map using findpatt.m 

% For trapezoidal discretization scheme and legs2d model
limit1_tr = 70; % nodes upper limit for computation of line 168 in FDJac.m when coloring map is provided (70 in Marko's machine, 79 in Ton's machine)
limit2_tr = 161; % % nodes upper limit for conputation of coloring map using findpatt.m 

% For Euler discretization scheme and legs3d model
limit2_eu = 35; % nodes upper limit for computation of coloring map using findpatt.m 

% For trapezoidal discretization scheme and legs3d model
limit2_tr = 35; % nodes upper limit for computation of coloring map using findpatt.m 


if model == 0 % legs2d
    if dscheme == 0 % Euler discretization scheme
        if n_nod <= limit2_eu
            % Coloring map is computed using findpatt.m and the sparsity 
            % pattern and stored in Prob.ConIx.
            % Note: line 369 of FDJac.m, which performs expensive computation, 
            % should be commented out to decese computation time for 
            % n_nod < limit1_eu, and must be commented out for 
            % n_nod > limit1_eu
            display('Ensure that Prob.NumDiff=1!')
            if n_nod < limit1_eu
                display('Line 369 in FDJac.m should be commented out to decrease computation time')
            else
                display('Ensure line 369 in FDJac.m is commented out!')
            end
            display('Coloring map being computed by findpatt.m ..')   
            conix = findpatt(ConsPattern);

        else % n_nod > limit2_eu
            % Coloring map cannot be computed by findpatt.m because of memory
            % limitations.           
            display('Ensure that Prob.NumDiff=1!')
            display('Ensure that line 369 in FDJac.m is commented out!')             
            
            if bsymm == 0 % symmetry deactivated
                % Use a modificaton of the routine provided by Tomlab, which avoids
                % crashing due to memory limitations
                display('Using a modification of findpatt.m, findpatt_mod.m, which avoids crashing due to memory limitations.');
                display('Coloring map being computed by findpatt_mod.m ..')
                conix = findpatt_mod(ConsPattern);
                
            else % symmetry activated
                % Coloring map is constructed based on coloring map
                % obtained for 25 nodes. This can be done for a number of nodes 
                % 6*p+1, where p is an integer.
                % Note: here line 369 of FDJac.m, which performs 
                % expensive computation, should be c

                % Check whether number of nodes is n_nod=6*p+1, with p integer.
                if mod(n_nod-1,6) ~= 0
                    n_nod_sugg = n_nod - mod(n_nod-1,6);
                    error(['(n_nod-1) is not a multiple of 6! For n_nod > ' int2str(limit2_eu) ' use n_nod=p*6+1, with p integer. Use e.g. n_nod = ' int2str(n_nod_sugg)])  
                end

                % Computing coloring map
                display('Coloring map being assembled ..');

                % Initial part (first 6*66 elements) 
                BC =  [1     2     3     4     5     6     7     8     9    10    11    12 ...
                      13    14    15    16    17    18    19    20    21    22    23    24 ...
                      25    26    27    28    29    30    31    32    33    34     1     1 ...
                      1     1     1     1     1     1     1     1     1     1     1     1 ...
                      1     1     2     2     2     2     2     2     2     2     2     2 ...
                      2     2     2     2     2     2     2     1     4     3     6     5 ...
                      8     7    19     9    20    21    22    23    24    25    26    27 ...
                      28    29    30    31    32    33    34    35    36    37    38    39 ...
                      40    41    42    43     2     2     2     2     2     2     2     2 ...
                      2     2     2     2     2     2     2     2     3     3     3     3 ...
                      3     3     3     3     3     3     3     3     3     3     3     3 ...
                      1     2     3     4     5     6     7     8    10    11    12    13 ...
                      14    15    16    17    18    28    19    30    29    32    31    34 ...
                      33    36    35    38    37    40    39    42    41    44     1     1 ...
                      1     1     1     1     1     1     1     1     1     1     1     1 ...
                      1     1     3     3     3     3     3     3     3     3     3     3 ...
                      3     3     3     3     3     3     2     1     4     3     6     5 ...
                      8     7     9    10    19    20    21    22    23    24    25    26 ...
                      27    29    30    31    32    33    34    35    36    37    38    39 ...
                      40    41    42    43     2     2     2     2     2     2     2     2 ...
                      2     2     2     2     2     2     2     2     3     3     3     3 ...
                      3     3     3     3     3     3     3     3     3     3     3     3 ...
                      1     2     3     4     5     6     7     8    11     9    12    13 ...
                      14    15    16    17    18    27    28    30    29    32    31    34 ...
                      33    36    35    38    37    40    39    42    41    44     1     1 ...
                      1     1     1     1     1     1     1     1     1     1     1     1 ...
                      1     1     3     3     3     3     3     3     3     3     3     3 ...
                      3     3     3     3     3     3     2     1     4     3     6     5 ...
                      8     7    10    11    19    20    21    22    23    24    25    26 ...
                      29    28    30    31    32    33    34    35    36    37    38    39 ...
                      40    41    42    43     2     2     2     2     2     2     2     2 ...
                      2     2     2     2     2     2     2     2     3     3     3     3 ...
                      3     3     3     3     3     3     3     3     3     3     3     3];

                % Repeating patterns (6*66 elements repeated nr times)
                RC = [1     2     3     4     5     6     7     8     9    10    12    13 ...
                    14    15    16    17    18    27    28    29    31    30    33    32 ...
                    35    34    37    36    39    38    41    40    43    42     1     1 ...
                    1     1     1     1     1     1     1     1     1     1     1     1 ...
                    1     1     3     3     3     3     3     3     3     3     3     3 ...
                    3     3     3     3     3     3     2     1     4     3     6     5 ...
                    8     7    11     9    19    20    21    22    23    24    25    26 ...
                    29    28    30    31    32    33    34    35    36    37    38    39 ...
                    40    41    42    43     2     2     2     2     2     2     2     2 ...
                    2     2     2     2     2     2     2     2     3     3     3     3 ...
                    3     3     3     3     3     3     3     3     3     3     3     3 ...
                    1     2     3     4     5     6     7     8    10    11    12    13 ...
                    14    15    16    17    18    27    28    29    31    30    33    32 ...
                    35    34    37    36    39    38    41    40    43    42     1     1 ...
                    1     1     1     1     1     1     1     1     1     1     1     1 ...
                    1     1     3     3     3     3     3     3     3     3     3     3 ...
                    3     3     3     3     3     3     2     1     4     3     6     5 ...
                    8     7     9    10    19    20    21    22    23    24    25    26 ...
                    29    28    30    31    32    33    34    35    36    37    38    39 ...
                    40    41    42    43     2     2     2     2     2     2     2     2 ...
                    2     2     2     2     2     2     2     2     3     3     3     3 ...
                    3     3     3     3     3     3     3     3     3     3     3     3 ...
                    1     2     3     4     5     6     7     8    11     9    12    13 ...
                    14    15    16    17    18    27    28    29    31    30    33    32 ...
                    35    34    37    36    39    38    41    40    43    42     1     1 ...
                    1     1     1     1     1     1     1     1     1     1     1     1 ...
                    1     1     3     3     3     3     3     3     3     3     3     3 ...
                    3     3     3     3     3     3     2     1     4     3     6     5 ...
                    8     7    10    11    19    20    21    22    23    24    25    26 ...
                    29    28    30    31    32    33    34    35    36    37    38    39 ...
                    40    41    42    43     2     2     2     2     2     2     2     2 ...
                    2     2     2     2     2     2     2     2     3     3     3     3 ...
                    3     3     3     3     3     3     3     3     3     3     3     3];

                % final part (6*66+1 elements, if containing only T, or 6*66+2 elements
                % if containing T and v_walk)
                if isempty(v_walk) % locomotion speed NOT prescribed
                    FC = [1     2     3     4     5     6     7     8     9    10    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7    11     9    19    20    21    22    23    24    25    26 ...
                        29    28    30    31    32    33    34    35    36    37    38    39 ...
                        40    41    42    43     2     2     2     2     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     3     3     3     3 ...
                        3     3     3     3     3     3     3     3     3     3     3     3 ...
                        1     2     3     4     5     6     7     8    10    11    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7     9    10    19    20    21    22    23    24    25    26 ...
                        29    28    30    31    32    33    34    35    36    37    38    39 ...
                        40    41    42    43     2     2     2     2     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     3     3     3     3 ...
                        3     3     3     3     3     3     3     3     3     3     3     3 ...
                        1     2     3     4     5     6     7     8    11     9    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7    10    35    36    37    38    39    40    41    42    43 ...
                        11    19    20    21    22    23    24    25    26    28    29    30 ...
                        31    32    33    34     3     3     3     3     3     3     3     3 ...
                        4     4     4     4     3     3     3     3     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     2     2     2     2 ...
                        3 45];

                else % locomotion speed prescribed
                    FC = [1     2     3     4     5     6     7     8     9    10    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7    11     9    19    20    21    22    23    24    25    26 ...
                        29    28    30    31    32    33    34    35    36    37    38    39 ...
                        40    41    42    43     2     2     2     2     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     3     3     3     3 ...
                        3     3     3     3     3     3     3     3     3     3     3     3 ...
                        1     2     3     4     5     6     7     8    10    11    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7     9    10    19    20    21    22    23    24    25    26 ...
                        29    28    30    31    32    33    34    35    36    37    38    39 ...
                        40    41    42    43     2     2     2     2     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     3     3     3     3 ...
                        3     3     3     3     3     3     3     3     3     3     3     3 ...
                        1     2     3     4     5     6     7     8    11     9    12    13 ...
                        14    15    16    17    18    27    28    29    31    30    33    32 ...
                        35    34    37    36    39    38    41    40    43    42     1     1 ...
                        1     1     1     1     1     1     1     1     1     1     1     1 ...
                        1     1     3     3     3     3     3     3     3     3     3     3 ...
                        3     3     3     3     3     3     2     1     4     3     6     5 ...
                        8     7    10    35    36    37    38    39    40    41    42    43 ...
                        11    19    20    21    22    23    24    25    26    28    29    30 ...
                        31    32    33    34     3     3     3     3     3     3     3     3 ...
                        4     4     4     4     3     3     3     3     2     2     2     2 ...
                        2     2     2     2     2     2     2     2     2     2     2     2 ...
                        45];
                end

                % Assembling coloring map
                nr = ((n_nod-1)-12)/6;

                Cpatt = BC;
                for i = 1:nr
                    Cpatt = [Cpatt RC];
                end
                Cpatt = [Cpatt FC];

                conix=Cpatt;
            end
        end
        
    elseif dscheme == 1 % trapezoidal discretization scheme
        % Note: odd number of nodes seems to result in more efficient
        % coloring maps, n_nod=2p+1, where p is integer
       if n_nod <= limit2_tr
            % Coloring map is computed using findpatt.m and the sparsity 
            % pattern and stored in Prob.ConIx.
            % Note: line 369 of FDJac.m, which performs expensive computation, 
            % should be commented out to decrease computation time for 
            % n_nod < limit1_tr, and must be commented out for 
            % n_nod > limit1_tr
            display('Ensure that Prob.NumDiff=1!')
            if n_nod < limit1_tr
                display('Line 369 in FDJac.m should be commented out to decrease computation time')
            else
                display('Ensure line 369 in FDJac.m is commented out!')
            end
            display('Coloring map being computed by findpatt.m ..')   
            conix = findpatt(ConsPattern);

        else % n_nod > limit2_tr
            % Coloring map cannot be computed by findpatt.m because of memory
            % limitations.           
            display('Ensure that Prob.NumDiff=1!')
            display('Ensure that line 369 in FDJac.m is commented out!')             
            
            if bsymm == 0 % symmetry deactivated
                % Use a modificaton of the routine provided by Tomlab, which avoids
                % crashing due to memory limitations
                display('Using a modification of findpatt.m, findpatt_mod.m, which avoids crashing due to memory limitations.');
                display('Coloring map being computed by findpatt_mod.m ..')
                conix = findpatt_mod(ConsPattern);
                
            else         
                % Coloring map is constructed based on coloring map
                % obtained for 25 nodes. This can be done for an odd number of
                % nodes, i.e. n_nod=2*p+1, where p is an integer.
                % Note: here line 369 of FDJac.m, which performs expensive 
                % computation, must be commented out.

                % Check whether number of nodes is n_nod=2*p+1, with p integer.
                if mod(n_nod-1,2) ~= 0
                    error('Use an odd number of nodes!')  
                end

                % Computing coloring map
                display('Coloring map being assembled ..');

                % Repeating patterns (2*66 elements repeated nr times)
                RC = [1     2     3     4     5     6     7     8     9    10    11    12 ...
                      13    14    15    16    17    18    19    20    21    22    23    24 ...
                      25    26    27    28    29    30    31    32    33    34     1     1 ...
                      1     1     1     1     1     1     1     1     1     1     1     1 ...
                      1     1     2      2     2     2     2     2     2     2     2     2 ...
                      2     2     2     2     2     2    35    36    37    38    39    40 ...
                      41    42    43    44    45    46    47    48    49    50    51    52 ...
                      53    54    55    56    57    58    59    60    61    62    63    64 ...
                      65    66    67    68     3     3     3     3     3     3     3     3 ...
                      3     3     3     3     3     3     3     3     4     4     4     4 ...
                      4     4     4     4     4     4     4     4     4     4     4     4];

                % final part (2*66+1 elements, if containing only T, or 2*66+2 elements
                % if containing T and v_walk)
                if isempty(v_walk) % locomotion speed NOT prescribed
                    FC = [1     2     3     4     5     6     7     8     9    10    11    12 ...
                          13    14    15    16    17    18    19    20    21    22    23    24 ...
                          25    26    27    28    29    30    31    32    33    34     1     1 ...
                          1     1     1     1     1     1     1     1     1     1     1     1 ...
                          1     1     2     2     2     2     2     2     2     2     2     2 ...
                          2     2     2     2     2     2    35    36    37    38    39    40 ...
                          41    42    43    44    45    46    47    48    49    50    51    52 ...
                          53    54    55    56    57    58    59    60    61    62    63    64 ...
                          65    66    67    68     3     3     3     3     3     3     3     3 ...
                          3     3     3     3     3     3     3     3     4     4     4     4 ...
                          4     4     4     4     4     4     4     4     4     4     4     4 ...
                          2    69];

                else % locomotion speed prescribed
                    FC = [1     2     3     4     5     6     7     8     9    10    11    12 ...
                          13    14    15    16    17    18    19    20    21    22    23    24 ...
                          25    26    27    28    29    30    31    32    33    34     1     1 ...
                          1     1     1     1     1     1     1     1     1     1     1     1 ...
                          1     1     2     2     2     2     2     2     2     2     2     2 ...
                          2     2     2     2     2     2    35    36    37    38    39    40 ...
                          41    42    43    44    45    46    47    48    49    50    51    52 ...
                          53    54    55    56    57    58    59    60    61    62    63    64 ...
                          65    66    67    68     3     3     3     3     3     3     3     3 ...
                          3     3     3     3     3     3     3     3     4     4     4     4 ...
                          4     4     4     4     4     4     4     4     4     4     4     4 ...
                          69];
                end

                % Assembling coloring map
                nr = ((n_nod-1)-4)/2;

                Cpatt = RC;
                for i = 1:nr
                    Cpatt = [Cpatt RC];
                end
                Cpatt = [Cpatt FC];

                conix=Cpatt;
           end
       end
    else % do nothing
    end
    
elseif model == 1 % legs3d
    display('Ensure that Prob.NumDiff=1!')
    display('Ensure line 369 in FDJac.m is commented out!')
    display('Coloring map being computed by findpatt.m ..')
    if (n_nod <= limit2_eu) && (n_nod <= limit2_tr)
        % Use the routine provided by Tomlab
        conix = findpatt(ConsPattern);
    else
        % Use a modificaton of the routine provided by Tomlab, which avoids
        % crashing due to memory problems
        display('Using a modification of findpatt.m, findpatt_mod.m, which avoids crashing due to memory limitations.');
        conix = findpatt_mod(ConsPattern);
    end
else % other models
end