Несколько советов о том, как написать лучший код для науки о данных, используя python и numpy. Мы используем индивидуальный случай, чтобы показать хорошие шаблоны.

Введение

Большая часть любых научных усилий заключается в написании кода. Будь то типичное моделирование машинного обучения, анализ или участие в проекте данных, значительная часть времени уходит на прототипирование новых функций. Ожидается, что, будучи исследовательской, некоторые части программы будут заменены или изменены, зачастую сверх того, что было запланировано изначально.

В отличие от «потребительского программного обеспечения», изменения часто вызваны не требованиями заказчика, а скорее результатами процесса на этом пути. По этой причине чрезвычайно важно спроектировать его таким образом, чтобы не потребовалась «полная реконструкция», если экспериментальные данные указывают на другой путь.

Написание научного кода связано с двумя (дополнительными) специфическими проблемами: первая связана с математическими ошибками. Ошибки в вычислениях часто трудно отследить, особенно когда код семантически корректен. Багов не обнаружено. Никаких исключений не возникает. Все выглядит хорошо, но (числовой) результат неверен. В частности, при реализации вероятностных моделей результаты иногда могут выглядеть хорошо, в зависимости от некоторых начальных условий или случайного фактора.

Второй исходит из факта, описанного ранее. Всегда будут экспериментальные части, и они постоянно меняются. Таким образом, ключевым моментом является проектирование каждого компонента, чтобы большая часть работы могла остаться и служить надежной основой для следующего этапа разработки.

В этой статье мы сосредоточимся на шаблонах, которые могут сделать код более надежным, понятным и в целом более простым в обращении. Вы увидите, как простые улучшения могут привести к большему количеству повторно используемого кода и меньшему количеству ошибок.

Игрушечная задача

Наша задача, ради демонстрации, состоит в том, чтобы вычислить математическое ожидание результатов случайного процесса. Математически это сводится к следующей формуле:

где - функция массы вероятности (PMC), или для непрерывной переменной:

где – функция плотности вероятности (PDF).

Как вам может быть интересно, проблема здесь в том, что вы хотите, чтобы он работал с любым распределением: непрерывным или дискретным. Или, если это неправдоподобно, вы хотите, по крайней мере, распознать природу проблемы, чтобы вы могли плавно адаптировать код.

Плохой код — отправная точка

Начнем с не очень хорошего кода. Допустим, вы хотите бросить шестигранный кубик. Поскольку каждый исход равновероятен, все сводится к вычислению среднего значения по выборочным результатам.

import random
import statistics


def die(sides = 6):
    return random.randint(1, sides)

def expected_value(n_samples):
    samples = [die() for _ in range(n_samples)]
    return statistics.mean(samples)

Что не так с этим кодом? Несколько вещей…

Во-первых, функция die возвращает по одной выборке за раз. Его нужно вызывать N раза, чтобы получить N выборки, что очень медленно.

Во-вторых, функция expected_value сильно зависит от функции die, производящей выборки. Это очевидно, если вы решите использовать другую матрицу, скажем, 12-гранную. В этом дизайне вам нужно «открыть» expected_value, чтобы принять дополнительный параметр sides, только чтобы передать его die, чтобы расширить его для более общего случая. Хотя это сработало бы, но делает интерфейс expected_value нелогичным, но решение по-прежнему основано на использовании die в качестве источника образцов, что затрудняет рассмотрение других дистрибутивов.

Средства защиты…

Рассмотрим варианты, которые у вас есть:

Идея №1

Вы можете использовать внешнюю переменную для хранения образцов:

def expected_value(samples):
    return statistics.mean(samples)

samples = [die(12) for _ in range(n_samples)]

ev = expected_value(samples)

Это довольно очевидно, но вы просто перемещаете проблему в другое место… Теперь samples становится новой сущностью, хранилищем данных (пусть даже очень больших), и оно довольно анонимно. Функция expected_value ожидает его получения, но его подготовка лежит на вас.

Идея №2

Другой способ — сохранить die внутри, передав его expected_value как объект.

from functools import partial

twelve_sided_die = partial(die, 12)

def expected_values(die, n_samples):
    samples = [die() for _ in range(n_samples)]
    return statistics.mean(samples)


ev = expected_values(twelve_sided_die, 10000)

Идея использует подготовленную «версию» die и заставляет expected_value использовать ее в качестве источника сэмплов. Однако возникает новая проблема: expected_value совместим только с die. Он не может вычислить результат с помощью любого другого «генератора образцов», или, по крайней мере, не гарантируется, что он сделает это правильно.

Идея №3

Третья идея заключается в том, чтобы распознать проблему на более абстрактном уровне и разработать более совершенные интерфейсы.

На абстрактном уровне у нас есть два понятия:

  • Существует распределение вероятностей, из которого мы можем сделать выборку. (Это может быть кубик, монета, нормальное распределение — не важно).
  • Существует математическая операция, которая использует и преобразует данные. (Например, вычисление среднего значения, дисперсии и т. д.).

Давайте уделим больше внимания тому, как мы можем построить правильные абстракции и как контролировать их взаимную совместимость.

Распространение (данные)

Математически распределение вероятностей может быть функцией — непрерывной или дискретной, конечной или бесконечной, из которой мы можем извлекать выборки. «Рецепт» для этой функции может сильно различаться в зависимости от задачи. Мы можем использовать «существующую» формулу, такую ​​как распределение Гаусса или Пуассона, но это также может быть «нестандартное» создание, полученное, например, из гистограмма.

Имея это в виду, давайте установим следующую абстракцию:

from abc import ABC, abstractmethod

class Distribution(ABC):
    
    @abstractmethod
    def sample(self):
        ...

Выполнение

Из-за @abstractmethod наш дистрибутив требует реализации метода sample для любого дочернего класса, производного от этой абстракции. Для нашего штампа это может быть:

import numpy as np


class Die(Distribution):
    def __init__(self, sides = 6):
        self.sides = sides

    def sample(self, n_samples):
        return np.random.randint(1, self.sides + 1, size=n_samples)

Здесь доставка образцов происходит путем вызова метода sample, специфичного для подбрасывания честной кости: Die(12).sample(10000). Более того, благодаря numpy мы можем очень быстро создавать большие объемы данных, заменив понимание списка на np.ndarray.

На самом деле, все можно улучшить еще больше. В настоящее время вызов Die() возвращает что-то вроде этого <__main__.Die at 0x7f43f4448400>, что не является информативным. Также Die() == Die() оценивается как False, поскольку для python это два разных экземпляра объекта одного и того же класса. Чтобы исправить это, нам нужно реализовать еще два метода:

def __repr__(self):
        return f"{self.__class__.__name__}(sides={self.sides})"

    def __eq__(self, other):
        if isinstance(other, self.__class__):
            return self.sides == other.sides
        return False

Метод __repr__ делает визуализацию объекта приятной, а __eq__ вернет True только в том случае, если игральные кости «равносторонние».

Классы данных

Внедрение четырех методов каждый раз может быть утомительным. Кроме того, текущая реализация Die не предотвращает изменение объекта, даже случайное, путем присвоения атрибута существующему объекту, такому как этот die.sides = 20.

Следовательно, мы можем переопределить класс Die, используя класс dataclasses Python.

from dataclasses import dataclass


@dataclass(frozen=True)
class Die(Distribution):
    sides: int = 6

    def sample(self, n):
        return np.random.randint(1, self.sides + 1, size=n)

Поведение этого примера такое же, как и в предыдущем. Кроме того, установка frozen=True, присвоение нового значения die.sides вызовет исключение. Если нам нужен новый кубик, мы должны создать новый объект.

Теперь наша функция expected_value, скорее всего, возьмет die в качестве объекта распределения и выполнит математические действия, вызвав свой метод sample.

def expected_value(distribution, n):
    return distribution.sample(n=n).mean()

Типы

Приведенный выше пример аккуратен. Мы точно знаем, что делает expected_value, и это легко проверить. Однако n-гранная кость — не единственное распределение, для которого мы можем захотеть вычислить ожидаемое значение. Например, при подбрасывании монеты результаты не являются числовыми (если только мы не установим соглашение и не будем его придерживаться). Естественно, имеет смысл дать несколько советов о том, какие интерфейсы можно использовать вместе и как.

Для языка динамически типизированного, такого как python, вы не обязаны придерживаться типов переменных. Однако при использовании различных IDE и инструментов, таких как mypy, ввод текста может помочь выявить потенциальные точки отказа и сделать код более прозрачным.

Давайте переопределим наши классы, используя typing, и создадим еще два конкретных дистрибутива.

from typing import Generic, Sequence, TypeVar


D = TypeVar("D")


class Distribution(ABC, Generic[D]):
    
    @abstractmethod
    def sample(self, n: int) -> Sequence[D]:
        ...


@dataclass(frozen=True)
class Die(Distribution[int]):
    sides: int = 6

    def sample(self, n: int) -> Sequence[int]:
        return np.random.randint(1, self.sides + 1, size=n)


@dataclass(frozen=True):
class Coin(Distribution[str]):
    outcomes: tuple = ("H", "T")
    fairness: float = 0.5

    def sample(self, n: int) -> Sequence[str]:
        p = (self.fairness, 1.0 - self.fairness)
        return np.random.choice(self.outcomes, size=n, p=p)


@dataclass(frozen=True):
class Gaussian(Distribution[float]):
    mu: float = 0.0
    sigma: float = 1.0

    def sample(self, n: int) -> Sequence[float]:
        np.random.normal(loc=self.mu, scale=self.sigma, size=n)

Несколько вещей, которые происходят здесь. Благодаря D = TypeVar("D") теперь мы можем определить новый тип переменной, с помощью которого мы параметризуем тип каждого распределения. Вы можете заметить, что Distribution наследуется не только от абстрактного базового класса, но и от Generic[D], что также превращает его в новый (параметризованный) тип. Теперь это становится своего рода идентичностью, составляющей новый тип данных.

Ожидается, что каждая версия sample будет возвращать последовательность определенного типа, которая имеет смысл в контексте каждого отдельного дистрибутива. Таким образом, у нас есть единый интерфейс, который также параметризуется. Мы можем использовать это, чтобы обеспечить правильное поведение expected_value:

def expected_value(
    distribution: Distribution[float | int], n: int
) -> float:
    return distribution.sample(n=n).mean()

При прохождении, например. die = Die() или gaussian = Gaussian() до expected_value будут работать (поскольку int и float являются числовыми), передача coin = Coin() будет помечена, например. mypy, заявив, что

> error: Argument 1 to "expected_value" has incompatible type "Coin";
> expected "Distribution[Union[float, int]]"

Это может дать нам раннее предупреждение еще до того, как мы запустим код.

Повышение чистоты

Как видите, проектирование интерфейсов с использованием typing помогает формализовать намерения и выявить ошибки на ранней стадии. Вы даже можете вывести его на новый уровень, используя dtype numpy. Таким образом, вы не только убедитесь, что различные элементы подходят друг к другу, но и будете более осознанно относиться к объему памяти, занимаемому вашими данными.

Например:

import numpy.typing as npt


class Distribution(ABC, Generic[D]):
    
    @abstractmethod
    def sample(self, n: int) -> np.NDArray[np.generic]:
        ...


class Die(Distrinution[int]):
    sides: int = 6

    def sample(self, n: int) -> npt.NDArray[np.uint8]:
        return np.random.randint(
            1, self.sides + 1, size=n, dtype=np.uint8
        )

Таким образом, вы даже будете проинформированы, если метод die.sample возвращает числа, отличные от строго беззнаковых 8-битных целых чисел. Вопрос в том, хотите ли вы идти так глубоко? Есть над чем подумать.

Расчеты

Вернемся к разработке расчетной части. На данный момент у нас есть expected_value, готовый к работе с числовыми распределениями. Естественно, мы можем вычислить ожидаемое значение для Die и Gaussian, но не для Coin. Не с текущим дизайном.

Чтобы исправить это, есть два варианта:

  • Мы можем создать дистрибутив proxy, сопоставив, например, ("H", "T") -> (0, 1) или
  • Мы включаем сопоставление в expected_value, давая возможный «адаптер».

Первый подход создает искусственное тело, идея которого опирается на условность. Это не мешает никому определять другой прокси с ("H", "T") -> (1, 0), что затрудняет обнаружение ошибки.

Вместо этого мы можем изменить expected_value и дать ему возможность использовать пользовательские адаптеры.

def expected_value(
    d: Distribution[D],
    f: Callable[[D], Any] = lambda x: x,
    n: int = 1000
) -> float:
    return np.mean(np.apply_along_axis(f, axis=0, arr=d.sample(n)))

Второй аргумент expected_value — это вызываемая функция (функция), которую мы можем дополнительно использовать для перевода результатов, например. Coin() раздача. Однако по умолчанию он оставляет результаты нетронутыми.

die = Die(12)
expected_values(die, n=10000)

gaussian = Gaussian(mu=4.0, sigma=2.0)
expected_value(gaussian, n=100000)

# but
coin = Coin(fairness=0.75)
expected_value(coin, f=lambda x: np.where(x == "H", 1.0, 0.0))

Здесь мы не только избежали создания прокси-дистрибутива, но и сумели избежать привязки expected_value к какому-либо конкретному способу преобразования данных. Функция expected_value делает только то, что обещает: вычисляет ожидаемое значение. Если требуется какая-либо адаптация или преобразование, это обеспечивается извне. Обратите внимание, что здесь у нас также есть вариант: мы можем определить именованную функцию (например, coin_conversion) на случай, если мы планируем использовать ее повторно, или придерживаться lambda, когда отдельное определение не добавляет ценности.

Состав и итерация

Абстрагирование от математических вычислений оказывается полезным, особенно при разработке итерационных алгоритмов. Часто, помимо основных вычислений, нам приходится следить за некоторыми побочными результатами, такими как сходимость, ранняя остановка (максимальное количество итераций), метрики и так далее.

Возьмем, к примеру, константу. Математически мы можем получить его значение, приняв следующий предел:

Чем выше, тем точнее приближение.

Как мы можем решить ее адекватно?

Не лучший способ…

Начнем с подхода бедняка с петлей.

def approximate_e(
    initial_guess: float,
    max_iter: int = 10,
    epsilon: float = 0.01
) -> float:
    e_value = initial_guess
    for n in range(1, max_iter + 1):
        new_e_value = (1.0 + 1.0 / n) ** n
        if abs(new_e_value - e_value) < epsilon:
            return new_e_value
        e_value = new_e_value

    return new_e_value

Опять же, что не так с этим подходом?

Прежде всего, функция делает три вещи вместо одной. Строка 8. является абсолютной сутью расчета. Тем не менее, с условиями ранней остановки и сходимости у нас остается много накладных расходов кода, которые тесно связаны с фактическими вычислениями. Хотя эти два условия выглядят как нечто более общее, если мы решим заменить объект аппроксимации (например, на квадратный корень), нам придется скопировать и вставить этот дополнительный код и убедиться, что он не нарушит новый алгоритм. .

Во-вторых, наши единственные варианты параметризации двух условий — либо жестко закодировать значения для max_iter и epsilon, либо позволить пользователям предоставлять их в качестве аргументов. Это портит интерфейс и усложняет тестирование.

Наконец, алгоритм генерирует данные «с нетерпением». Вместо того, чтобы сосредоточиться на математике и предоставлять значения «по запросу», он выдает вам данные. Для больших объемов данных это может вызвать проблемы с памятью.

Абстракция над вычислениями

Теперь давайте решим сразу все три проблемы, разделив ответственность между разными частями. У нас есть три вещи:

  • что-то, что мы ожидаем вернуть нам правильные числа (фактические расчеты),
  • что-то, чтобы остановить процесс, если пройдет достаточное количество итераций (ранняя остановка),
  • что-то, чтобы остановить процесс, если значения больше не улучшаются заметно (конвергенция).
from typing import Iterator
import itertools


def approximate_e() -> Iterator[float]:
    n = 1
    while True:
        yield (1.0 + 1.0 / n) ** n
        n += 1

def early_stop(
    values: Iterator[float], max_iter: int = 100
) -> Iterator[float]:
    return itertools.islice(values, max_iter)

def convergence(
    values: Iterator[float], epsilon: float = 0.01
) -> Iterator[float]:
    for a, b in itertools.pairwise(values):
        yield a

        if (a - b) < epsilon:
            break

В этом дизайне используются итераторы, реализующие «ленивую» загрузку. Элементы данных возвращаются только один за другим при запросе (отсюда и ключевое слово yield). Благодаря этому нам (почти) не нужно беспокоиться о памяти.

Кроме того, каждая из трех функций может существовать отдельно. У них есть свой специфический интерфейс, и их можно тестировать отдельно.

Наконец (и самое главное), окончательный результат можно получить, связав их вместе!

values = approximate_e()
values = early_stop(values, max_iter=50)
values = convergence(values, epsilon=0.001)

for value in values:
    print("e becomes:", value)

Старый питон

Для Python старше 3.10 itertools.pairwise можно записать:

def pairwise(values: Iterator[D]) -> Iterator[Tuple[D, D]]:
    a = next(values, None)
    if a is None:
        return
    for b in values:
        yield a, b
        a = b

Выводы

Научное программирование сопряжено с дополнительными трудностями. Это захватывающая дисциплина, но количество ловушек, кажется, растет, по крайней мере, квадратично с увеличением сложности задачи.

В этой статье мы обсудили два основных компонента, которые снова и снова встречаются в каждой задаче научного программирования: данные и вычисления. Надеемся, что, правильно абстрагируя их, вы сможете не только повысить эффективность своего кода и уменьшить количество ошибок, но и сделать программирование намного удобнее. Для вас и вашей команды!

Первоначально опубликовано на https://zerowithdot.com.