Кто-то в качестве своего первого серьезного проекта пишет калькулятор, другие веб-сайт и в принципе идея ограничивается только скилами и вашей собственной фантазией (или фантазией преподавателя, если у вас такой есть). Калькулятор мне писать не хотелось т.к. слишком уж это тривиально, а для веба на тот момент я был слишком зелен и ничего не понимал в протоколах, html, фреимворках и т.д.

Но месяцы за изучением Python были пройдены, сотни задач на Stepik, Leetcode и подобных ресурсах решены, десятки часов Хирьянова и других святил на Youtube с упоением просмотрены. И руки невероятно чесались написать что-то «своё, масштабное и нужное».

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

И тут через «сарафанное радио» мы находим друг друга!

Если кратко, то у одной из Самарских лабораторий есть программа на С#, которая была написана лет 5-10 назад силами одного из учёных (не профессиональным программистом) и на момент её написания всё было круто, но время шло и количество заказов, как и нагрузка увеличивалась. И наступило время, когда данная программа стала узким местом в общем рабочем процессе.

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

Звучит просто, не так ли?)

По крайней мере мне так показалось.

Вроде бы с xlsx я ещё не работал, но ведь работал с другими форматами, вряд ли с ним будет сложнее. Но в середине марта пришло т.з. и я понял, что работа с xlsx, будет пожалуй самым простым, что предстоит выполнить.

Вот парочка к пунктов для понимания (в моём случае непонимания):

3. Найти петли гистерезиса, при этом следует учитывать, что их число может быть произвольным. Способ поиска – смена знака скорости (т.е. перемещение начинается в другую сторону).

4. Для каждой петли гистерезиса определить следующие параметры:

a. Жёсткость – определяется как отношение разницы силы к разнице перемещения в начале и конце петли.

b. Коэффициент рассеяния – определяется как умноженное на 4 отношение площади петли к площади под срединной линией. Площадь петли – разность между площадью под кривой нагрузки и разгрузки. Площадь под срединной линией – полусумма площадей под двумя этими кривыми. Площадь под кривой можно определить, например, методом Симпсона. Следует учесть, что петля может быть незамкнутой.

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

Было чувство, что человек заказавший мне написание и отправивший ТЗ, был уверен, что я хорошо разбираюсь в физике. Да, я разбираюсь в физике. Но в физике уровня скорее 7-го класса, в котором она у меня собсна началась. Я смогу дать ответ почему небо голубое или откуда ветер берётся. Но интегральные же вычисления, которые требовались для подсчёта были от меня невероятно далеки.

Опять же для понимания, входные данные представляли из себя следующее, примерно до 6000 контрольных точек:

Как всё это выглядело на графике по шкале удлинения (пройденного расстояния):

Сами же вычисления должны начинаться с конца плоского участка («полки»).

И график, где по одной шкале пройденное расстояние, а по другой оказываемое воздействие (такой красивый он получается, если строить его с конца «полки»):

И вот эти самые петли (Гистерезиса), что получаются в процессе трёх минутного тестирования изделия и нужно обсчитывать.

В общем нужно с чего то начинать и я начал.

# Импорт библиотеки openpyxl для чтения и записи файлов формата .xlsx.
import openpyxl

def read_and_write_xlsx(link: str, num_sheet: int) -> (list[float], list[float]):
    """
    Функция для четения файла формата xlsx и запись его данных по отдельным именнованым контейнерам типа list.
    :param link: Ссылка на файл.
    :param num_sheet: Номер листа для чтения.
    :return: Два именнованных контейнера с данными об удлинении и нагрузки испытуемого образца.
    """
    elongation: list = []  # Удлинение (см)
    burden: list = []  # Нагрузка (N)
    book = openpyxl.open(link, read_only=True)  # Объект принимающий ссылку с целью чтения.
    sheet = book.worksheets[num_sheet]  # Объект для уточнения листа для чтения.
    # Цикл, на запись данных по контейнерам.
    for e, b in sheet.iter_rows(min_row=2, min_col=3, max_col=4):
        elongation.append(abs(e.value))
        burden.append(abs(b.value))
    return elongation, burden

Контейнерами было решено выбрать list (хотя ни каких изменений в исходные данные вносить не планируется), т.к. объём не настолько велик, чтобы можно было почувствовать разницу с тем же tuple.

Далее нужно было выработать некий алгоритм для нахождения упомянутой ранее полки.

Небольшая ремарка. В проекте используется лишь одна библиотека — openpyxl, служащая для считывания и записи данных. И не потому что я хотел сделать всё своими силами, принципиально не пользуясь чужими по умному сконструированными велосипедами. Я бы с удовольствием импортировал тот же numpy и itertools (да хоть через * прости госпади), дабы упростить себе задачу, но я ни в них ни в других прекрасных библиотеках и модулях, не нашёл ничего, чем бы можно было воспользоваться «из коробки». Единственное исключение возможно numpy.trapz, которым я безуспешно пытался посчитать площадь петель Гистерезиса и по итогу навелосипедил своего.

Для поиска полки же был написан следующий алгоритм:

def search_shelf(burden: list[float]) -> tuple[int, int]:
    """
    Функция для нахождения начала и конца полки.
    :param burden: Данные по нагрузке.
    :return: Индексы начала и конца полки.
    """
    start_shelf = 0  # Индекс элемента начала полки.
    end_shelf = 0  # Индекс элемента конца полки.
    counter = 10  # Счётчик для отхода от индекса начала полки.
    new_start = None  # Индекс самого первого экстремума.
    for i in range(1000, len(burden)):
        # Стартует с 1000-го значения, т.к. до этого мы точно ничего не ждём.
        if all(burden[i] > j for j in burden[i - 10:i]) and all(burden[i] > j for j in burden[i + 1 : i + 10]):
            new_start = i
            break

    for i in range(new_start, len(burden)):  # Цикл для поиска начала и конца полки.
        if not start_shelf:  # Пока не задано значение начала полки.
            if len({el // 100 for el in burden[i:i + 10]}) == 2:
            # Если длина можества из 10-ти следущих значений равна 2.
                start_shelf = i  # Значит, спустя 10 значений будет начало полки.

        else:  # Если мы знаем значение начала полки.
            if counter:  # Пропускаем несколько циклов.
                counter -= 1
                continue
            if len({el // 100 for el in burden[i:i + 10]}) > 3:
            # Если длина множества из 10-ти следующих значений больше 3.
                end_shelf = i  # Значит, найден конец полки.
                break  # Остановка цикла после нахождения конца полки.

    return start_shelf, end_shelf

Вроде пока всё выполнимо, если немного посидеть и подумать.

Дальше шёл следующий этап по поиску (корректное название подсказал отец) экстремумов. Если простым языком, то экстремум — это крайняя точка, когда график идёт вверх или вниз.

def search_extremes(elng: list[float], end_shelf: int) -> list[int]:
    """
    Функция для поиска экстемумов.
    :param elng: Список со значениями удлинения.
    :param end_shelf: Индекс конца полки.
    :return: Скисок индексов экстремумов.
    """
    inds_extremes = []  # Список для точек экстремума.
    counter = 0  # Счётчик для отхода от экстремума, на случей если значение экстремума на графике повторяется.

    for i, el in enumerate(elng[end_shelf:-10], start=end_shelf):  # Цикл для нахождения экстремумов.
        if counter:
            counter -= 1
        # Проверка условия, что экстремум найден, даже если рядом точка с идентичными данными.
        elif ((el <= elng[i - 1] and el < elng[i - 2] and el < elng[i - 3] and el < elng[i - 4] and
               el <= elng[i + 1] and el < elng[i + 2] and el < elng[i + 3] and el < elng[i + 4]) or
              (el >= elng[i - 1] and el > elng[i - 2] and el > elng[i - 3] and el > elng[i - 4] and
               el >= elng[i + 1] and el > elng[i + 2] and el > elng[i + 3] and el > elng[i + 4])):
            inds_extremes.append(i)
            counter = 10

    inds_extremes.append(len(elng))  # Добавление к экстремумам индекса последнего значения.

    return inds_extremes

Внимательный читатель спросит, «а для чего такое всратое запутанное условие?» И будет прав!

И тут начинаются первые звоночки, касательно входных данных.

Скажу так, это не первая вариация данной функции и до определённого момента она была достаточно компактная и лаконичная. Но потом обнаружилось, что на одном из 14ти экстремумов одного из 30 листов с данными (то бишь примерно 1 из 400 раз), значение дублируется до 5-го знака после запятой и экстремум превращается в миркополку так сказать.

Но это ещё пол беды.

Возможно помните строчку из ТЗ: «Способ поиска – смена знака скорости». Так вот — забудьте! Вторая же часть такого условия, в том, что по ходу пути к очередному экстремуму, линия может резко поменять вектор буквально на одно значение, а после вновь вернуться на исходное направление. Так проявляется дрожание датчика.

И дабы всё это учесть, функция была реализована именно таким вот образом.

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

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

Но всё хорошо в меру.

Идём дальше.

А дальше была ещё одна не особо примечательная функция, после которой я и пытался применять numpy.trapz, но безуспешно.

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

def search_area_under_line(start_ind: int, end_ind: int, elongation: list[float], burden: list[float]) -> int:
    """
    Функция считает площадь под кривой.
    :param start_ind: Индекс начала кривой.
    :param end_ind: Индекс конца кривой.
    :param elongation: Список удлинения.
    :param burden: Список нагрузки.
    :return: Площадь под заданой кривой.
    """
    need_burden = [i - min(burden[start_ind:end_ind]) for i in
                   burden[start_ind:end_ind]]  # Скорректированные занчения нагрузки.
    need_elongation = elongation[start_ind:end_ind]  # Копия значений удлинения.
    area = 0  # Площадь под кривой.
    for i in range(len(need_elongation) - 1):  # Цикл для подсчёта площадь под кривой методом интегрирования.
        area += (need_burden[i] + need_burden[i + 1]) / 2 * (need_elongation[i + 1] - need_elongation[i])

    return area


def dissipation_coefficient_and_rigidity(total_list: list, elongation: list[float], burden: list[float]) -> (
        float, int):
    """
    Функция для расчёта коэффициента рассеивания и жёсткости по петлям гистерезиса, начиная с амплитуды 0,2.
    :param total_list: список контрольных точек (экстремумов и антиподов).
    :param elongation: Список удлинения.
    :param burden: Список нагрузки.
    :return: Коэффициент рассеивания и жёсткость.
    """
    list_dissipation_coefficient = []  # Коэффициент рассеивания.
    list_rigidity = []  # Жёсткость.
    for i in range(0, len(total_list), 3):  # Цикл для вычисления коэффициента рассеивания и жётскости.
        var_1 = abs(search_area_under_line(total_list[i], total_list[i + 1], elongation, burden))  # Площадь разгружения.
        var_2 = search_area_under_line(total_list[i + 1], total_list[i + 2], elongation, burden)  # Площадь нагружения.
        loop_area = var_2 - var_1  # Площадь петли.
        sum_area = var_1 + var_2  # Сумма площадей под линией разгружения и нагружения.
        dissipation_coefficient = 4 * loop_area / (sum_area * 0.5)  # Вычисление коэффициента рассеивания.

        rigidity = ((burden[total_list[i]] - burden[total_list[i + 1]]) /       # Вычисление жёсткости.
                    (elongation[total_list[i]] - elongation[total_list[i + 1]]))
        list_dissipation_coefficient.append(round(dissipation_coefficient, 4))
        list_rigidity.append(int(rigidity))

    return list_dissipation_coefficient, list_rigidity

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

Финал

И вот, спустя 14 дней разработки и календарный месяц, к середине апреля программа была готова. С изначальных двух часов обработки 64-ёх листов с данными на С#, которые нужно было последовательно отправлять на расчёт, время сократилось до 30 секунд и одного клика, что не может не радовать. Оставалось только написание текстового интерфейса взаимодействия с пользователем (что бы даже школьник смог разобраться, но ничего случайно при этом не уронить) и дополнение документации.

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

Очень внимательный читатель спросит: «а как можно считать площадь незамкнутой фигуры?» Кусочек ТЗ: «Следует учесть, что петля может быть незамкнутой.»

А вот и я не знаю, а возможно и никто не знает, но вот что точно известно, что если приблизить, любая петля, на любом графике будет выглядеть примерно вот так:

Но моя программа на такие мелочи внимания не обращала и продолжала считать всё с некоторой долей погрешности.

Красным отмечено место, где петля должна замыкаться, как бы перекрывая саму себя и идя внахлёст.

Узнав данное обстоятельство, я написал алгоритм для программного замыкания петель, но учёные умы сказали, что не православному это, менять исходные данные. Отчего пришлось от данного решения проблемы отказаться.

И тут мы плавно возвращаемся к теме конфликта интересов и зон ответственности.

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

И всё таки, насколько должны быть точны входные данные? Ответа на данный вопрос у меня до сих пор нет и ооооочень внимательный читатель может заметить, что в оглавлении данной статьи, слово «писал», а не «написал». Т.к. технически проект есть и теперь любой желающий может воспользоваться им GitHub и произвести подсчёт прекрасных и завораживающих петель Гистерезиса, но на «вооружение» его так и не приняли, отчего несколько печально.

Если вы до читали до этого момента, то хочу выразить вам свою благодарность, надеюсь, было не очень скучно и местами даже интересно.) На момент публикации я с радость принимаю заказы, т.к. всё ещё верю в силу «сарафанного радио»!

Всем добра!)