count_number_of_modes = function( data, h ){
    #
    # Counts the number of density humps found when we estimate the
    # density of data with a KDE and a spacing of h
    #
    de = density( data, bw=h, adjust=1, kernel="gaussian" )

    # Look for local maximum (what I'm defining as a peak/hump/mode):
    #
    mask = diff( diff( de$y )>=0 ) < 0

    nm = sum(mask)
    return(nm)
}