function extrenum_data = find_min_max_locations( sampling_t, fdobj, degree, error_level )
    %
    % find_min_max_locations returns a structure of maxima and minima in fdobj
    % in the degree'th and degree+1'th derivative functions. Also 
    % includes the end points (boundaries) of the vector.
    %
    % This method works by looking for changing signs in the velocity 
    % (degree+1) and acceleration (degree+2) of the curve described in the
    % fdobj object:
    % Function maximum: degree+1th derivative = 0 and degree+2th derivative < 0
    % Function minimum: degree+1th derivative = 0 and degree+2th derivative > 0
    % Derivative maximum: degree+2th derivative = 0 and degree+3th derivative < 0
    % Derivative minimum: degree+2th derivative = 0 and degree+3th derivative > 0
    %
    % Note that cases have been found (e.g. subjectId = 11 of the growth
    % data) where both an extrenum in velocity and an extrenum in
    % acceleration occur between the same two x values. For this reason,
    % we look for both velocity and acceleration extrenum between each pair
    % of sampling values.
    %
    % Input parameters
    % ----------------
    % sampling_t    t values where where the function (and derivative) 
    %               signs are sampled.
    %               This sampling must be dense enough not to miss any 
    %               extrenum event. It is recommended to use the x values 
    %               of the raw data. 
    % fdobj         Functional object that fits the raw data
    % degree        In which derivative are we looking for maxima and minima. To
    %               find max in the function, use degree = 0, and to find spurts (that is,
    %               max in velocity), use degree = 1
    % error_level   indicates the precision that is required in the
    %               estimation of the zero (i.e., the xvalue selected will be between:
    %               [- error_level, + error_level]
    %
    % Output value
    % ------------
    % extrenum_data  a structure of information about the extrenum found. 
    %                Fields are: 
    %                * location, where the extrenum event is located
    %                * type: type of extrenum (see global_constants)
    %                (One entry per extrenum).

    % loading global parameters
    global_constants;

    tolerance = 0.00;   % sign changes will have to be at least 
                        % tolerance % of the total y range  
                        % for instance, if y varies from -4 to +6
                        % and tolerance = 0.02 (2%), then when a sign
                        % change is detected, the difference between the two
                        % points will have to be at least 0.02 x 6 - (-4) = 0.2


    optimize = true;   % for debug purposes only. Leave to true
    current_extrenum_index = 1;

    % include first point data (as boundary value)
    first_point.location = sampling_t(1);
    first_point.type     = BOUNDARY;

    extrenum_data(current_extrenum_index) = first_point;
    current_extrenum_index = current_extrenum_index + 1;

    sampling_t_size = size(sampling_t);
    last_point_index      = sampling_t_size(1);

    % initialize deriv and accel values
    prev_deriv  = eval_fd(sampling_t(1), fdobj, degree+1);
    prev_acc    = eval_fd(sampling_t(1), fdobj, degree+2);
    prev_xvalue = sampling_t(1);

    % computing ranges
    all_deriv   = eval_fd(sampling_t, fdobj, degree+1);
    all_acc     = eval_fd(sampling_t, fdobj, degree+2);
    range_deriv = max(all_deriv) - min(all_deriv);
    range_acc   = max(all_acc)   - min(all_acc);

    % look for min and max in non-boundary points
    for i = 2:last_point_index-1
        current_xvalue = sampling_t(i);

        current_deriv  = eval_fd(current_xvalue, fdobj, degree+1);
        current_acc    = eval_fd(current_xvalue, fdobj, degree+2);
        current_jerk   = eval_fd(current_xvalue, fdobj, degree+3);

        % creating empty structures
        new_deriv_event = struct('location', {}, 'type', {});
        new_accel_event = struct('location', {}, 'type', {});

        % Has velocity (1st derivative) changed sign? --> max or min of
        % function
        if ((sign(current_deriv) ~= sign(prev_deriv)) && ... % sign has changed
            ((abs(current_deriv - prev_deriv) / range_deriv) > tolerance)) % the difference was large enough (i.e., not due to noise)

            % a zero crossing occurs for an x in the interval current_xvalue - prev_xvalue
            % Does it correspond to a maxima or a minima?
            if (current_acc < 0)
                % a maxima corresponds to a negative value of the derivative
                type = MAX_FCT;
            else
                type = MIN_FCT;
            end

            % now use Newton method to find a better approximation of this

            % limit the search to the interval between current value and
            % previous value because we know the sign has changed between
            % those two values, so the zero crossing has to be within that
            % range
            interval_width = current_xvalue - prev_xvalue; 

            % Begin with the current x value as the estimated location 
            if (optimize == true)
                estimated_location = newton_optimize( current_xvalue, interval_width, ...
                    sampling_t(1), sampling_t(last_point_index), ... % function is defined over this interval
                    false, fdobj, degree+1, error_level );
                % store the x value estimate
                new_deriv_event(1).location = estimated_location;
                new_deriv_event(1).type     = type;
            else
                % no optimization -- just use sampled value
                new_deriv_event(1).location = current_xvalue;
                new_deriv_event(1).type     = type;
            end
        end

        % Has acceleration (2nd derivative) changed sign? --> max or min of
        % velocity
        if ((sign(current_acc) ~= sign(prev_acc)) && ... % sign has changed
            ((abs(current_acc - prev_acc) / range_acc) > tolerance)) % the difference was large enough (i.e., not due to noise)

            % same thing as above, but in the acceleration (degree + 2)
            if (current_jerk < 0)
                type = MAX_DERIV;
            else
                type = MIN_DERIV;
            end

            interval_width = current_xvalue - prev_xvalue;
            if (optimize == true)
                estimated_location = newton_optimize( current_xvalue, interval_width, ...
                    sampling_t(1), sampling_t(last_point_index), ... % function is defined over this interval
                    false, fdobj, degree+2, error_level );
                % store the x value estimate
                new_accel_event(1).location = estimated_location;
                new_accel_event(1).type     = type;
            else
                % no optimization -- just use sampled value
                new_accel_event(1).location = current_xvalue;
                new_accel_event(1).type     = type;
            end
        end

        % insert deriv and accel events in the right order
        % first, check if both events occured in the current interval
        if (size(new_deriv_event,1) > 0 && size(new_accel_event,1) > 0)
            % yes, we have both events, check which one came first

            if (new_deriv_event.location < new_accel_event.location)
                % velocity event happens first
                extrenum_data(current_extrenum_index) = new_deriv_event;
                current_extrenum_index = current_extrenum_index + 1;

                extrenum_data(current_extrenum_index) = new_accel_event;
                current_extrenum_index = current_extrenum_index + 1;
            else
                % acceleration event happens first
                extrenum_data(current_extrenum_index) = new_accel_event;
                current_extrenum_index = current_extrenum_index + 1;

                extrenum_data(current_extrenum_index) = new_deriv_event;
                current_extrenum_index = current_extrenum_index + 1;
            end
        else
            % either an acceleration or a velocity event (extrenum) 
            % occured, but not both. Insert that event in the
            % extrenum_data structure
            if (size(new_deriv_event) > 0)
                extrenum_data(current_extrenum_index) = new_deriv_event;
                current_extrenum_index = current_extrenum_index + 1;
            end
            if (size(new_accel_event) > 0)
                extrenum_data(current_extrenum_index) = new_accel_event;
                current_extrenum_index = current_extrenum_index + 1;
            end
        end

        % to move to next point, current estimates become the previous ones
        prev_deriv  = current_deriv;
        prev_acc    = current_acc;
        prev_xvalue = current_xvalue;
    end

    % append info of the boundary
    last_point.location = sampling_t(last_point_index);
    last_point.type = BOUNDARY;
    extrenum_data(current_extrenum_index) = last_point;
    current_extrenum_index = current_extrenum_index + 1;

    % Now perform the extrenum alternation test
    extrenum_data = eliminate_inflection_points(extrenum_data);
return;