#include "Clusterizator.h"
#include "Utils.h"

bio::HiCMatrix & bio::Clusterizator::prepareHiCMatrix(HiCMatrix & matrix)
{
    for (auto it = matrix.getData().begin(); it != matrix.getData().end();)
    {
        bool val = false;
        for (size_t j = 0; j < it->size(); ++j)
        {
            if ((*it)[j] > 0.0)
            {
                val = true;
                break;
            }
        }

        if (!val)
            it = matrix.getData().erase(it);
        else
            ++it;
    }

    for (size_t j = 0; j < matrix[0].size();)
    {
        bool val = false;
        for (int i = 0; i < matrix.getSize(); ++i)
        {
            if (matrix[i][j] > 0.0)
            {
                val = true;
                break;
            }
        }


        if (!val)
        {
            for (int i = 0; i < matrix.getSize(); ++i)
                matrix[i].erase(matrix[i].begin() + j);
        }
        else
            ++j;
    }

    return matrix;
}

std::vector<double> bio::Clusterizator::getAvgFlow(const HiCMatrix & matrix)
{
    double *buf[2];
    buf[0] = new double[matrix.getSize()];
    buf[1] = new double[matrix.getSize()];
    int cur = 0;

    for (size_t line = 0; line < matrix.getSize() - 1; ++line)
    {
        for (size_t pos = 0; pos <= line; ++pos)
        {
            // (p + q + 2 * a) / 2 = (p + q) / 2 + a
            const double& a = matrix[pos][matrix[pos].size() - 1 - line + pos];
            if (pos == 0 && pos == line)
                buf[cur][pos] = a;
            else if (pos == 0)
                buf[cur][pos] = (buf[1 - cur][pos] * (line)+a) / (line + 1.0);
            else if (pos == line)
                buf[cur][pos] = (buf[1 - cur][pos - 1] * (line)+a) / (line + 1.0);
            else
            {
                double leftAvg = (buf[1 - cur][pos - 1] * (line)+a) / (line + 1.0);
                double rightAvg = (buf[1 - cur][pos] * (line)+a) / (line + 1.0);
                buf[cur][pos] = (leftAvg + combinationRatio(line - 1, pos) * rightAvg) / (1.0 + combinationRatio(line - 1, pos));
            }
        }

        cur = 1 - cur;
    }

    std::vector<double> res(matrix.getSize() - 1);
    for (int i = 0; i < res.size(); ++i)
        res[i] = buf[1 - cur][i];

    delete[] buf[0];
    delete[] buf[1];

    return res;
}

size_t bio::Clusterizator::findMinValuePos(const double * func, size_t leftBin, size_t rightBin)
{
    size_t left = leftBin + 1;
    size_t right = rightBin - 1;
    // Ñäâèãàåì ëåâóþ ãðàíèöó, ïîêà íå ïåðåâàëèì çà ëîêàëüíûé ìàêñèìóì ñëåâà
    while (func[left] < func[left - 1] && left != right)
        ++left;

    if (left == right)
        return (leftBin + rightBin) / 2;

    // Ñäâèãàåì ïðàâóþ ãðàíèöó, ïîêà íå ïåðåâàëèì çà ëîêàëüíûé ìàêñèìóì ñïðàâà
    while (func[right] < func[right + 1] && right != left)
        --right;

    if (left == right)
        return (leftBin + rightBin) / 2;

    double min = func[left];
    size_t minPos = left;

    for (int i = left + 1; i <= right; ++i)
        if (func[i] < min)
        {
            min = func[i];
            minPos = i;
        }

    return minPos;
}

void bio::Clusterizator::addChilds(const double * func, ClusterNode * parent)
{
    if (parent->getCluster().leftBin == parent->getCluster().rightBin) return;

    size_t minPos;
    if (parent->getCluster().rightBin - parent->getCluster().leftBin == 1)
        minPos = (parent->getCluster().leftBin + parent->getCluster().rightBin) / 2;
    else
        minPos = findMinValuePos(func, parent->getCluster().leftBin, parent->getCluster().rightBin);
    parent->setLeft(new ClusterNode(
        Cluster(parent->getCluster().leftBin, minPos, 0.0),
        &parent->getRoot(),
        parent,
        nullptr, nullptr));

    parent->setRight(new ClusterNode(Cluster(minPos + 1, parent->getCluster().rightBin, 0.0), &parent->getRoot(), parent, nullptr, nullptr));

    addChilds(func, &parent->getLeft());
    addChilds(func, &parent->getRight());
}

bio::ClusterNode * bio::Clusterizator::buildHierarchy(const std::vector<double>& func)
{
    ClusterNode *root = new ClusterNode(Cluster(0, func.size() - 1, 0.0), nullptr, nullptr, nullptr, nullptr);
    root->setRoot(root);

    addChilds(func.data(), root);


    /*std::vector<ClusterNode*> nodes(func.size());

    size_t curCluster = 1;
    for (;;)
    {
        for(int i = 0; )
    }*/
    return root;
}