Расчёт токов в цепи постоянного тока

/*
Дана цепь постоянного тока с сопротивлениями и источниками (напряжения, ЭДС и тока). Найти токи и напряжения во
всех элементах. Программа составляет систему уравнений по 1 и 2 законам Кирхгофа и решает её с проверкой по
балансу мощностей.
Последовательность действий:
п.1. Пронумеровать все узлы схемы (нумерация с 1).
п.2. Напечатать в файл in.txt все элементы схемы в виде:
2 3 r5 3
2 1 u2 4
7 2 e13 4
2 4 j5 4
где первое и второе число это номера узлов к которым подключен элемент, r5, u2, e13, j5 - это сам элемент схемы
(сопротивление или источник напряжения, ЭДС, тока), и последнее число это величина в омах, вольтах или амперах. Везде пробелы.
Для напряжения, эдс и тока узлы должны быть по порядку: первый - от которого напряжение/эдс/ток направлен, второй - к которому
он направлен(направлен внутри источника а не во внешней цепи); для сопротивления направление безразлично. Пустые ветви можно
вводить как ветви с сопротивлением 0. Если ветвь содержит более одного элемента то между каждыми соседними элементами ветви
надо вставить узел и присвоить узлу номер. Допустим если в одной ветви два сопротивления то между ними рисуем узел и
присваиваем ему номер следующий за максимальным номером на схеме. Таким образом каждый элемент будет находиться между двумя
узлами и этот элемент можно будет ввести в файл.
п.3. Поместить файл in.txt в папку с программой и запустить программу.
В результате для каждого элемента ветви будет вычислен его ток и напряжение.
Найденное напряжение или ток в элементе считается направленным от меньшего узла элемента к большему.
Пример файла содержащего описание схемы:
 
1 2 r4  1
3 1 r   2
2 1 u12 -4
2 3 r2  3
 
Расшифровка файла:
К первому и второму узлу подключено сопротивление r4=1 ом.
К третьему и первому узлу подключено сопротивление r=2 ом.
К второму и первому узлу подключено напряжение u12=-4 вольт (направлено от узла 2 к узлу 1).
К второму и третьему узлу подключено сопротивление r2=3 ом.
В итоге файл с описанием схемы получается как в задачнике Бычкова только без номеров элементов.
*/
#include <iostream>
#include <algorithm>
#include <utility>
#include <fstream>
#include <vector>
#include <string>
#include <regex>
#include <ranges>
#include <set>
#include <clocale>
#include <random>
#include <iomanip>
#include <iterator>
#include <queue>
#include <optional>
#include <cassert>
#include <numeric>
#include <cctype>
#include <sstream>
#include <cmath>
#include <functional>
#include <map>
#include <ctime>
#include <format>
#include <unordered_set>
#include <climits>
#include <cfloat>
namespace rng = std::ranges;
 
double const EPS = 1e-10; // порог при котором число типа double считается нулём
using md = std::vector<std::vector<double>>;
 
class Elem;
class Smd;
std::map<int, std::vector<Elem>> mNodes, mNodesDfs;
std::set<std::set<Elem>> setContTot;
std::unordered_set<int> visNod;
std::vector<Elem> visEl;
std::vector<Elem> vEls;
int nodeStart = -1;
bool foundDfs = false;
 
// Вывод матрицы в поток 
auto& operator<<(std::ostream& ost, md const& mtx)
{
    ost << "\n------------------------ " << mtx.size() << " x " << mtx[0].size() << " ------------------------------\n\n";
    for (auto const& v : mtx)
    {
        for (auto val : v)
            ost << std::setw(8) << std::left << std::setprecision(4) << val;
        ost << "\n\n";
    }
    ost << "----------------------------------------------------------------------\n";
    return ost;
}
 
bool zero(double rhs)
{
    return std::fabs(rhs) < EPS;
}
bool notZero(double rhs)
{
    return std::fabs(rhs) > EPS;
}
 
// Прямой ход - приведение матрицы mtx к верхнетреугольному виду.
void triangularMatrix(md& mtx)
{
    for (size_t j = 0; j < std::min(mtx.size() - 1, mtx[0].size() - 2); ++j)
    {
        size_t maxrowindex = j;
        for (size_t i = j + 1; i < mtx.size(); ++i)
        {
            if (std::fabs(mtx[i][j]) > std::fabs(mtx[maxrowindex][j]))
                maxrowindex = i;
        }
        std::swap(mtx[maxrowindex], mtx[j]);
        if (notZero(mtx[j][j]))
        {
            auto mtxjj = mtx[j][j];
            rng::for_each(mtx[j], [mtxjj](auto& val) { val /= mtxjj; }); // Делим j-ю строку на её коэф-т с гл. диагонали
            for (size_t i = j + 1; i < mtx.size(); ++i)
            {
                auto mtxij = mtx[i][j];
                rng::transform(mtx[i], mtx[j], mtx[i].begin(), [mtxij](auto a, auto b) { return a - b * mtxij; }); // Вычитаем из i-й строки j-ю умноженную на j-й коэф-т нижней i-й строки
            }
        }
    }
    // обнуляем все коэф-ты меньшие EPS
    for (auto& v : mtx)
        rng::transform(v, v.begin(), [](auto val) { return std::abs(val) < EPS ? 0. : val; });
}
 
// элемент схемы
struct Elem
{
    Elem(int node1, int node2, std::string name, double val) : _node1(node1), _node2(node2), _name(name), _val(val)
    {
        --_node1;
        --_node2;
        char type = std::tolower(name.at(0));
        if (type == 'j') type = 'i';
        if (type != 'r' && _node1 > _node2)
        {
            std::swap(_node1, _node2);
            _val = -_val;
        }
        if (type == 'e')
        {
            type = 'u';
            _val = -_val; // превращаем в источник напряжения
        }
        _type = type;
    }
    auto operator<=>(Elem const& rhs) const
    {
        return _id <=> rhs._id;
    }
    bool operator==(Elem const& rhs) const
    {
        return _id <=> rhs._id == 0;
    }
 
    int _node1{ -1 }, _node2{ -1 };
    int _id{ INT_MAX };
    char _type = CHAR_MAX;
    std::string _name;
    double _val{ DBL_MAX };
};
 
// слагаемое уравнения
class Summand
{
public:
    Summand() = default;
    Summand(int varid, double coef, std::string const& elName) : _varid(varid), _coef(coef), _elName(elName)
    {
    }
    friend std::ostream& operator<<(std::ostream& ost, Summand const& smd)
    {
        std::string str;
        if (smd._coef != 1)
            str += std::format("{}{}{}", '(', smd._coef, ')');
        if (smd._varid != -1)
            str += std::format("{}{}{}{}", (smd.type() == 'i' ? 'u' : 'i'), '(', smd._elName, ')');
        ost << std::setw(11) << std::left << str;
        return ost;
    }
    char type() const
    {
        char c = 'f'; // free member
        if (auto itEl = rng::find(vEls, _varid, &::Elem::_id); itEl != vEls.end())
            c = itEl->_type;
        return c;
    }
    double _coef{ 1 };
    std::string _elName;// равен name элемента схемы
    int _varid{ -1 }; // равен ид элемента схемы. Для своб. слагаемых равен -1
};
 
// уравнение
class Equation
{
public:
    // добавляет слагаемое в уравнение с приведением подобных и удалением нулевых
    void Add(Summand const& smd)
    {
        if (zero(smd._coef))
            return;
        if (auto it = rng::find_if(_v, [&smd](auto const& val) { return smd._varid == val._varid; });
            it != _v.end())
        {
            it->_coef += smd._coef;
            if (zero(it->_coef))
                _v.erase(it);
        }
        else
            _v.push_back(smd);
    }
    Summand& operator[](int rhs) { return _v[rhs]; }
    Summand const& operator[](int rhs) const { return _v[rhs]; }
 
    friend std::ostream& operator<<(std::ostream& ost, Equation const& eq)
    {
        for (auto& smd : eq._v)
            ost << '+' << smd << ' ';
        ost << "= 0";
        return ost;
    }
    std::vector<Summand> _v;
};
 
// обход в глубину для поиска контуров
void dfsGetCircle(int nodeCur)
{
    visNod.insert(nodeCur);
    for (auto& elCur : mNodesDfs[nodeCur])
    {
        int nodeNext = elCur._node1 == nodeCur ? elCur._node2 : elCur._node1;
        if (!visNod.contains(nodeNext))
        {
            visEl.push_back(elCur);
            dfsGetCircle(nodeNext);
            visEl.pop_back();
        }
        else if (nodeNext == nodeStart && rng::find(visEl, elCur) == visEl.end())
        {
            std::set<Elem> setCont;
            for (auto& el : visEl)
            {
                setCont.insert(el);
            }
            setCont.insert(elCur);
            setContTot.insert(setCont);
        }
    }
    visNod.extract(nodeCur);
}
 
// обход в ширину для поиска дерева графа
std::queue<int> qBfs;
void bfsGetTree()
{
    while (!qBfs.empty())
    {
        int nodeCur = qBfs.front();
        qBfs.pop();
        for (auto& el : mNodes[nodeCur])
        {
            int nodeRem = el._node1 == nodeCur ? el._node2 : el._node1;
            if (!visNod.contains(nodeRem))
            {
                qBfs.push(nodeRem);
                visNod.insert(nodeRem);
                visEl.push_back(el);
            }
        }
    }
}
 
int main()
{
    std::setlocale(LC_CTYPE, ".1251");
    std::ifstream ifs{ "in.txt" };
    if (!ifs.is_open())
        throw 1;
 
    std::unordered_set<std::string> usNames;
    //double n = 4.3; // gggggggggggggggggg
    std::string str;
    std::stringstream ss;
    int node1{}, node2{};
    double val{};
    std::string name;
    for (int i = 0; std::getline(ifs, str); ++i)
    {
        str = std::regex_replace(str, std::regex{ "^\\s" }, "");
        if (str.empty() && mNodes.empty())
        {
            i = -1; continue;
        }
        else if (str.empty()) break;
        ss.clear();
        ss.str(str);
        node1 = -1; node2 = -1; name.clear();  val = DBL_MAX;
        ss >> node1 >> node2 >> name >> val;
        if (node1 == -1 || node2 == -1 || name.empty() || val == DBL_MAX)
        {
            std::cout << "Неверные данные в файле\n"; return 0;
        }
        //if (!usNames.insert(name).second)
        //{ std::cout << "Имена элементов схемы не должны совпадать\n"; return 0; }
        Elem el{ node1, node2, name, val };
        el._id = i;
        //el._val *= n;
        vEls.push_back(el);
        mNodes[el._node1].push_back(el);
        mNodes[el._node2].push_back(el);
    }
    ifs.close();
 
    // поиск дерева
    qBfs.push(0);
    visNod.insert(0);
    bfsGetTree();
    std::vector<Elem> vTree{ visEl };
 
    // поиск контуров
    std::vector<Elem> vAbsent;
    for (auto& el : vEls)
    {
        if (rng::find(vTree, el) == vTree.end())
        {
            mNodesDfs = mNodes;
            for (auto& [n, v] : mNodesDfs)
            {
                for (int i = 0; i < v.size(); ++i)
                {
                    if (rng::find(vTree, v[i]) == vTree.end())
                    {
                        v.erase(v.begin() + i);
                        --i;
                    }
                }
            }
            visEl.clear();
            visNod.clear();
            vAbsent.push_back(el);
            mNodesDfs[el._node1].push_back(el);
            mNodesDfs[el._node2].push_back(el);
            nodeStart = el._node1;
            dfsGetCircle(nodeStart);
            setContTot;
            vAbsent.pop_back();
        }
    }
 
    std::cout << "Независимые контуры\n";
    for (auto& st : setContTot)
    {
        for (auto& el : st)
        {
            std::cout << el._name << " ";
        }
        std::cout << "\n";
    }
    std::cout << "\n\n";
 
    std::set<std::set<Elem>> vContIndep;
    vContIndep = setContTot;
 
 
    // ур-я по 1 з. Кирхгофа
    std::vector<Equation> vSysz1;
    for (auto& [node, v] : mNodes)
    {
        Equation vEq;
        for (auto& el : v)
        {
            int remoteNode = el._node1 == node ? el._node2 : el._node1;
            Summand smd;
            if (el._type != 'i')
            {
                smd._varid = el._id;
                smd._elName = el._name;
            }
            if (el._type == 'i')
                smd._coef = el._val;
            if (node < remoteNode)
                smd._coef *= -1;
            vEq.Add(smd);
        }
        vSysz1.push_back(vEq);
    }
    //vSysz1.erase(vSysz1.begin());
 
    // ур-я по 2 з. Кирхгофа
    std::vector<Equation> vSysz2;
    for (auto& st : vContIndep)
    {
        Equation vEq;
        Elem elPrev = *st.begin();
        int nodePrev = elPrev._node1;
        Summand smd;
        if (elPrev._type != 'u')
        {
            smd._varid = elPrev._id;
            smd._elName = elPrev._name;
        }
        if (elPrev._type != 'i')
            smd._coef = elPrev._val;
        vEq.Add(smd);
        bool prevPosDir = true;
        for (int i = 0; i < st.size() - 1; ++i)
        {
            auto itCur = rng::find_if(st, [&](auto& el)
                { return (el._node1 == nodePrev || el._node2 == nodePrev) && el != elPrev; });
            Elem elCur = *itCur;
            Summand smd;
            if (elCur._type != 'u')
            {
                smd._varid = elCur._id;
                smd._elName = elCur._name;
            }
            if (elCur._type != 'i')
                smd._coef = elCur._val;
            int nodeConn = elCur._node1 == nodePrev ? elCur._node1 : elCur._node2;
            int maxNodePrev = std::max(elPrev._node1, elPrev._node2);
            int maxNodeCur = std::max(elCur._node1, elCur._node2);
            if ((nodeConn == maxNodePrev && nodeConn == maxNodeCur) || (nodeConn != maxNodePrev && nodeConn != maxNodeCur)) // вразнобой
            {
                if (prevPosDir)
                {
                    prevPosDir = false;
                    smd._coef *= -1;
                }
                else prevPosDir = true;
            }
            else
            {
                if (!prevPosDir)
                    smd._coef *= -1;
            }
            vEq.Add(smd);
            elPrev = elCur;
            nodePrev = elCur._node1 == nodePrev ? elCur._node2 : elCur._node1;
        }
        vSysz2.push_back(vEq);
    }
    vSysz1;
    std::cout << "Система уравнений\n";
    for (auto& eq : vSysz1)
        std::cout << eq << "\n";
    std::cout << "\n";
    for (auto& eq : vSysz2)
        std::cout << eq << "\n";
    std::cout << "\n";
 
    std::vector<Equation> vSys{ vSysz1 };
    vSys.append_range(vSysz2);
    std::set<int> setVars{ -1 }; //ид всех переменных входящих в систему (включая своб. член)
    for (int i = 0; i < vSys.size(); ++i)
        for (auto& smd : vSys[i]._v)
            setVars.insert(smd._varid);
 
    // приведение подобных слагаемых
    for (auto& eq : vSys)
    {
        rng::sort(eq._v, std::greater{}, &Summand::_varid); // чтобы своб. член попал назад
        for (int id : setVars)
        {
            auto er = rng::equal_range(eq._v, id, std::greater{}, &Summand::_varid);
            if (er.size() > 1)
            {
                for (auto ib = er.begin() + 1; ib != er.end(); ++ib)
                {
                    er.begin()->_coef += ib->_coef;
                }
                eq._v.erase(er.begin() + 1, er.end());
            }
            else if (!er.size())
            {
                // вставка нулевых слагаемых не нарушая сортировку вектора
                auto it = rng::find(vEls, id, &Elem::_id);
                eq._v.emplace(er.begin(), id, 0., it == vEls.end() ? "" : it->_name);
            }
        }
    }
    for (auto& eq : vSys)
    {
        rng::sort(eq._v.begin(), eq._v.end() - 1, {}, &Summand::_varid);
        eq._v.back()._elName.clear();
        if (notZero(eq._v.back()._coef))
            eq._v.back()._coef *= -1;
    }
 
    std::vector mtx{ vSys.size(), std::vector<double>(vSys[0]._v.size()) }; // расширенная матрица коэффициентов
    for (int i = 0; i < mtx.size(); ++i)
        for (int j = 0; j < mtx[i].size(); ++j)
            mtx[i][j] = vSys[i][j]._coef;
 
    std::cout << "\n\nИсходная матрица:\n" << mtx;
    triangularMatrix(mtx);
    std::cout << "\n\nТреугольная матрица:\n";
    std::cout << mtx << "\n\n";
 
    // обратный ход
    std::vector<std::optional<double>> vRoots(mtx[0].size() - 1);
    bool found = false;
    do
    {
        found = false;
        for (auto& v : mtx)
        {
            if (rng::count_if(v.begin(), v.begin() + vRoots.size(), [&](auto const& val)
                { return notZero(val) && !vRoots[&val - v.data()]; }) == 1)
            {
                found = true;
                double sum = 0;
                for (int i = 0; i < v.size() - 1; ++i)
                {
                    if (vRoots[&v[i] - v.data()])
                        sum += v[i] * *vRoots[&v[i] - v.data()];
                }
                auto it = rng::find_if(v, [&](auto const& val) { return notZero(val) && !vRoots[&val - v.data()]; });
                vRoots[&*it - v.data()] = (v.back() - sum) / *it;
            }
        }
    } while (found);
    for (auto& opt : vRoots)
        if (opt && zero(*opt))
            *opt = 0.;
 
    if (rng::all_of(vRoots, [](auto& opt) { return bool(opt); }))
    {
        std::cout << "\nПроверка уравнений подстановкой:\n\n";
        for (auto& eq : vSys)
        {
            double sum{};
            for (int i = 0; i < vRoots.size(); ++i)
            {
                sum += eq._v[i]._coef * *vRoots[i];
                if (zero(sum)) sum = 0.;
            }
            std::cout << sum << "==" << eq._v.back()._coef << "\n";
            if (notZero(sum - eq._v.back()._coef))
                throw;
        }
    }
 
    std::cout << "\nТоки и напряж направлены от меньшего узла к большему(см схему):\n\n";
    for (int i = 0; i < vRoots.size(); ++i)
        for (auto& v : vSys)
        {
            if (vRoots[i])
            {
                char varName = 'i';
                if (v[i].type() == 'i')
                    varName = 'u';
                std::cout << std::format("{:<9}", std::format("{}({}): ", varName, v[i]._elName)) << *vRoots[i] << '\n';
                break;
            }
        }
    std::cout << "\n";
 
    int cntSol = 1; // кол-во решений системы(2==бесконеч.)
    int rankExt = rng::count_if(mtx, [](auto const& v) // ранг расширенной м-цы
        { return rng::any_of(v, [](auto val) { return val != 0; }); });
    int rankNotExt = rng::count_if(mtx, [](auto const& v) // ранг нерасш. м-цы
        { return rng::any_of(v.begin(), v.end() - 1, [](auto val) { return val != 0; }); });
    if (rankExt != rankNotExt) cntSol = 0;
    else if (mtx[0].size() - 1 > rankExt) cntSol = 2; // ранги равны но неполны
 
    if (cntSol == 2)
        std::cout << "\nCистема имеет МНОЖЕСТВО РЕШЕНИЙ\n\n";
    else if (cntSol == 0)
        std::cout << "\nCистема НЕСОВМЕСТНА(не имеет решений)\n\n";
 
 
    std::cout << "\nНапряжения на сопротивлениях:\n\n";
    vEls;
    for (int i = 0; i < vRoots.size(); ++i)
    {
        if (!vRoots[i]) continue;
        for (auto& v : vSys)
        {
            if (!v[i]._elName.empty())
            {
                auto it = rng::find(vEls, v[i]._elName, &Elem::_name);
                if (it != vEls.end() && it->_type == 'r')
                {
                    std::cout << std::format("{:<23}", std::format("{:<4}{}{:.4}{}{:.4}{}", it->_name, ": ", *vRoots[i], " * ", it->_val, "  = "))
                        << *vRoots[i] * it->_val << "\n";
                }
                break;
            }
        }
    }
 
    std::cout << "\nБаланс мощностей:\n";
    double pGen{}, pCons{};
    for (int i = 0; i < vRoots.size(); ++i)
    {
        if (!vRoots[i])
        {
            std::cout << "Не все корни найдены\n"; break;
        }
        if (vSys[0][i]._elName.empty())
            throw;
        auto itEl = rng::find(vEls, vSys[0][i]._elName, &Elem::_name);
        switch (itEl->_type)
        {
        case 'r':
        {
            auto power = *vRoots[i] * *vRoots[i] * itEl->_val;
            pCons += power;
        }
        break;
        default: // j, u
        {
            auto power = *vRoots[i] * itEl->_val;
            if (power > 0) pCons += power;
            else pGen += -power;
        }
        break;
        }
    }
 
    std::cout << "(ген: " << pGen << ") == (" << "потр: " << pCons << ")\n\n";
    if (notZero(pGen - pCons))
        throw;
    int y = 0;
}