Сжатие нескольких вложенных циклов `for`

Question

Подобно этому и многим другим вопросам, у меня есть много вложенных циклов (до 16) одной и той же структуры.

Проблема : у меня 4-буквы алфавита и хочу получить все возможные слова длиной 16. Мне нужно отфильтровать эти слова .Это последовательности ДНК (следовательно, 4 буквы: ATGC), правила фильтрации довольно просты:

• нет подстрок XXXX (т.е. не может иметь одну и ту же букву подряд более 3 раз, ATGCAT GGGG CTA - это «плохо»)
• конкретное содержание GC, то есть количество Gs + количество Cs должно находиться в определенном диапазоне (40-50%).ATATATATATATA и GCGCGCGCGCGC - плохие слова

itertools.product будет работать для этого, но структура данных здесь будет гигантской (4 ^ 16 = 4 * 10 ^ 9 слов)

Более важноЕсли я использую product, то мне все равно придется пройти через каждый элемент, чтобы отфильтровать его.Таким образом, у меня будет 4 миллиарда шагов 2

Мое текущее решение является вложенным for loop

alphabet = ['a','t','g','c']
for p1 in alphabet:
    for p2 in alphabet:
       for p3 in alphabet:
       ...skip...
          for p16 in alphabet:
             word = p1+p2+p3+...+p16
             if word_is_good(word):
                 good_words.append(word)
                 counter+=1

Есть ли хороший шаблон для программирования без 16 вложенных циклов?Есть ли способ эффективно его распараллелить (на многоядерных или нескольких узлах EC2)? Также с этим шаблоном я могу подключить word_is_good? проверить в середине цикла: плохо запускается слово

...skip...
for p3 in alphabet:
   word_3 = p1+p2+p3
   if not word_is_good(word_3):
     break
   for p4 in alphabet:
     ...skip...

Answer 1

from itertools import product, islice
from time import time

length = 16

def generate(start, alphabet):
    """
    A recursive generator function which works like itertools.product
    but restricts the alphabet as it goes based on the letters accumulated so far.
    """

    if len(start) == length:
        yield start
        return

    gcs = start.count('g') + start.count('c')
    if gcs >= length * 0.5:
        alphabet = 'at'

    # consider the maximum number of Gs and Cs we can have in the end
    # if we add one more A/T now
    elif length - len(start) - 1 + gcs < length * 0.4:
        alphabet = 'gc'

    for c in alphabet:
        if start.endswith(c * 3):
            continue

        for string in generate(start + c, alphabet):
            yield string

def brute_force():
    """ Straightforward method for comparison """
    lower = length * 0.4
    upper = length * 0.5
    for s in product('atgc', repeat=length):
        if lower <= s.count('g') + s.count('c') <= upper:
            s = ''.join(s)
            if not ('aaaa' in s or
                    'tttt' in s or
                    'cccc' in s or
                    'gggg' in s):
                yield s

def main():
    funcs = [
        lambda: generate('', 'atgc'),
        brute_force
    ]

    # Testing performance
    for func in funcs:

        # This needs to be big to get an accurate measure,
        # otherwise `brute_force` seems slower than it really is.
        # This is probably because of how `itertools.product`
        # is implemented.
        count = 100000000
        start = time()
        for _ in islice(func(), count):
            pass
        print(time() - start)

    # Testing correctness
    global length
    length = 12
    for x, y in zip(*[func() for func in funcs]):
        assert x == y, (x, y)

main()

На моей машине generate был немного быстрее, чем brute_force, примерно 390 секунд против 425. Это было почти так же быстро, как я мог их сделать.Я думаю, что все это займет около 2 часов.Конечно, на самом деле их обработка займет гораздо больше времени.Проблема в том, что ваши ограничения не сильно уменьшают полный набор.

Вот пример того, как использовать это параллельно в 16 процессах:

from multiprocessing.pool import Pool

alpha = 'atgc'

def generate_worker(start):
    start = ''.join(start)
    for s in generate(start, alpha):
        print(s)

Pool(16).map(generate_worker, product(alpha, repeat=2))

Answer 2

Поскольку у вас есть алфавит длиной 4 (или любой "степень целого числа 2" ), идея использования целочисленного идентификатора и побитовых операций приходит на ум вместо проверки последовательныхсимволы в строках.Мы можем присвоить целое значение каждому из символов в alphabet, для простоты давайте используем индекс, соответствующий каждой букве.

Пример:

65463543 ₁₀ = 3321232103313 ₄ = 'aaaddcbcdcbaddbd'

Следующая функция преобразует целое число из 10 в слово, используя alphabet.

def id_to_word(word_id, word_len):
    word = ''
    while word_id:
        rem = word_id & 0x3  # 2 bits pet letter
        word = ALPHABET[rem] + word
        word_id >>= 2  # Bit shift to the next letter
    return '{2:{0}>{1}}'.format(ALPHABET[0], word_len, word)

Теперь дляфункция для проверки того, является ли слово «хорошим» на основе целочисленного идентификатора.Следующий метод имеет формат, аналогичный id_to_word, за исключением того, что счетчик используется для отслеживания последовательных символов.Функция вернет False, если будет превышено максимальное количество идентичных последовательных символов, в противном случае она вернет True.

def check_word(word_id, max_consecutive):
    consecutive = 0
    previous = None
    while word_id:
        rem = word_id & 0x3
        if rem != previous:
                consecutive = 0
        consecutive += 1
        if consecutive == max_consecutive + 1:
            return False
        word_id >>= 2
        previous = rem
    return True

Мы фактически рассматриваем каждое слово как целое число с основанием 4. Еслидлина алфавита не была значением "степень 2" , тогда по модулю % alpha_len и целочисленное деление // alpha_len можно было бы использовать вместо & log2(alpha_len) и >> log2(alpha_len) соответственно, хотя для этого потребовалось бы многобольше.

Наконец , найти все хорошие слова для данного word_len.Преимущество использования диапазона целочисленных значений заключается в том, что вы можете уменьшить количество for-loop с в вашем коде с word_len до 2, хотя внешний цикл очень велик.Это может позволить более дружественную многопроцессорную работу по поиску хорошего слова.Я также добавил в быстрый расчет, чтобы определить наименьшие и наибольшие идентификаторы, соответствующие хорошим словам, что помогает значительно сузить поиск хороших слов

ALPHABET = ('a', 'b', 'c', 'd')

def find_good_words(word_len):
    max_consecutive = 3
    alpha_len = len(ALPHABET)

    # Determine the words corresponding to the smallest and largest ids
    smallest_word = ''  # aaabaaabaaabaaab
    largest_word = ''   # dddcdddcdddcdddc
    for i in range(word_len):
        if (i + 1) % (max_consecutive + 1):
            smallest_word = ALPHABET[0] + smallest_word
            largest_word = ALPHABET[-1] + largest_word
        else:
            smallest_word = ALPHABET[1] + smallest_word
            largest_word = ALPHABET[-2] + largest_word

    # Determine the integer ids of said words
    trans_table = str.maketrans({c: str(i) for i, c in enumerate(ALPHABET)})

    smallest_id = int(smallest_word.translate(trans_table), alpha_len)  # 1077952576
    largest_id = int(largest_word.translate(trans_table), alpha_len)    # 3217014720

    # Find and store the id's of "good" words
    counter = 0
    goodies = []

    for i in range(smallest_id, largest_id + 1):
        if check_word(i, max_consecutive):
            goodies.append(i)
            counter += 1

В этом цикле я специально сохранил идентификатор слова какВ отличие от самого слова, если вы собираетесь использовать слова для дальнейшей обработки.Однако, если вы используете только слова, измените вторую на последнюю строку, чтобы прочитать goodies.append(id_to_word(i, word_len)).

ПРИМЕЧАНИЕ: Я получаю MemoryError при попытке сохранить все хорошие идентификаторы для word_len >= 14.Я предлагаю записать эти идентификаторы / слова в какой-нибудь файл!