sampleShell.cpp

#include <iostream>
#include <vector>
#include <tuple>
#include <random>
#include <algorithm>
#include <unordered_set>
#include <chrono>
#include <unordered_set>

using namespace std;

extern "C" {

vector<vector<int>> getNearestNeighbors(int dim) {
    if (dim == 2) {
        return {{-1, 0}, {0, 1}, {1, 0}, {0, -1}};
    } else if (dim == 3) {
        return {{-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {1, 0, 0}, {0, 0, -1}, {0, 0, 1}};
    }
    return {};
}

bool belongsToTheSpace(int newNode, std::unordered_set<int> spaceSet) {
    return spaceSet.find(newNode) != spaceSet.end();
}

bool isInGrid(const int node, const int labeledNodeSize) {
    if (node >= 0 && node < labeledNodeSize) {
        return true;
    }
    return false;
}

int to1DIndex(const vector<int>& Index, const int* shape, int dim) {
    int linearIndex = 0;
    int multiplier = 1;

    for (int i = dim - 1; i >= 0; --i) {
        linearIndex += Index[i] * multiplier;
        multiplier *= shape[i];
    }

    return linearIndex;
}

int convertCoord(const vector<int>& Index, const int* shape, int dim){
    if (dim == 2) {
        return Index[0] * shape[1] + Index[1];
    } 
    if (dim == 3) 
        return Index[0] * shape[1] * shape[2] + Index[1] * shape[2] + Index[2];
    return 0;
}

std::pair<float, int> neirestSumOfSamples(const int node, const int* sample, const int* labeledNodes, vector<int>& nearestNeighbors, const int sampleSize) {
    float NN_sum = 0.0;
    int count = 0;
    for (int& NN : nearestNeighbors) {
        int newNode = node + NN;
        if (isInGrid(newNode, sampleSize) && labeledNodes[newNode]) {
            NN_sum += sample[newNode];
            count++;
        }
    }
    return std::pair<float, int>(NN_sum/count, count);
}

float neighborInteraction(const int newNode, const float* p, const int* sample, const int* labeledNodes, vector<int>& nearestNeighbors, const int* shape, const int dim, const int sampleSize) {
    float NN_interaction = 0.0;
    int count = 0;
    for (int& NN : nearestNeighbors) {
        int sampledNode = newNode + NN;
        if (isInGrid(sampledNode, sampleSize) && labeledNodes[sampledNode] && p[newNode] ==  p[sampledNode]) { // When p[sampledNode] is equal to p[newNode] then the interaction is equal to the sum of the nearest neighbors samples
            if (p[newNode] >= p[sampledNode])
                NN_interaction += (p[newNode] - p[sampledNode])/(1 - p[sampledNode]) * (1 - sample[sampledNode]) + sample[sampledNode];
            else
                NN_interaction += p[newNode]/p[sampledNode] * sample[sampledNode];
            count++;
        }
    }
    return NN_interaction/count;
}

int multiplyElements(int* vector, int size){
    int result = 1;
    for (int i = 0; i < size; ++i) {
        result *= vector[i];
    }
    return result;
}

void printNoneZeroElements(int* sample, int sampleSize) {
    for (int i = 0; i < sampleSize; ++i) {
        if (sample[i] != 0) {
            cout << "Index " << i << ": " << sample[i] << '\t' << endl;
        }
    }
}

bool checkIfLabeled(int newNode, std::vector<int> labeledNodes) {

    return labeledNodes[newNode];
}

vector<int> getNearestNeighbors1D(vector<vector<int> > neirestNeighbors, int* sampleSize, int dim){
    vector<int> neirestNeighbors1D;
    for (const vector<int>& NN : neirestNeighbors) {
        neirestNeighbors1D.push_back(convertCoord(NN, sampleSize, dim));
    }
    return neirestNeighbors1D;
}

void sampleShell(int dim, int nOfElements, double q, int* sample, int* labeledNodes, float* p, int startIndex, int* shellSpace, int* sampleSize, int spaceSize) {
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<double> dis(0.0, 1.0);
    uniform_int_distribution<int> disInt(0, 1);

    vector<int> activeNodes;
    activeNodes.push_back(startIndex);
    int activeNodesCounter = 1;
    vector<vector<int> > neirestNeighbors = getNearestNeighbors(dim);
    vector<int> neirestNeighbors1D = getNearestNeighbors1D(neirestNeighbors, sampleSize, dim);
    const int sampleSizeInt = multiplyElements(sampleSize, dim);
    labeledNodes[startIndex] = 1;
    int maxActiveNodes = 0;
    for (int i = 0; i < nOfElements; ++i) {
        int Node = activeNodes[0];
        for (int& NN : neirestNeighbors1D) {
            int newNode = Node + NN;
            if (shellSpace[newNode] == 0 || labeledNodes[newNode] || newNode < 0 || newNode >= sampleSizeInt) {
                continue;
            }
            float p_samples = neighborInteraction(newNode, p, sample, labeledNodes, neirestNeighbors1D, sampleSize, dim, sampleSizeInt);
            sample[newNode] = dis(gen) < p_samples * q + (1-q) * p[newNode]; // p_samples * q + (1-q) * p
            activeNodes.push_back(newNode);
            activeNodesCounter++;
            labeledNodes[newNode] = 1;
        }
        if (maxActiveNodes < activeNodesCounter)
            maxActiveNodes = activeNodesCounter;

        if (activeNodesCounter!=0) {
        activeNodes.erase(activeNodes.begin());
            activeNodesCounter--;
        }
    }
}

void reSampleShell(int dim, int nOfElements, double q, int* sample, int* labeledNodes, float* p, int startIndex, int* shellSpace, int* sampleSize, int spaceSize) {
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<double> dis(0.0, 1.0);
    uniform_int_distribution<int> disInt(0, 1);

    vector<int> activeNodes;
    activeNodes.push_back(startIndex);
    int activeNodesCounter = 1;
    vector<vector<int> > neirestNeighbors = getNearestNeighbors(dim);
    vector<int> neirestNeighbors1D = getNearestNeighbors1D(neirestNeighbors, sampleSize, dim);
    const int sampleSizeInt = multiplyElements(sampleSize, dim);
    labeledNodes[startIndex] = 1;
    
    for (int i = 0; i < nOfElements; ++i) {
        int Node = activeNodes[0];
        for (int& NN : neirestNeighbors1D) {
            int newNode = Node + NN;
            if (shellSpace[newNode] == 0 || labeledNodes[newNode] || newNode < 0 || newNode >= sampleSizeInt) {
                continue;
            }
            float p_samples = neighborInteraction(newNode, p, sample, shellSpace, neirestNeighbors1D, sampleSize, dim, sampleSizeInt);

            sample[newNode] = dis(gen) < p_samples * q + (1-q) * p[newNode]; // p_samples * q + (1-q) * p

            activeNodes.push_back(newNode);
            activeNodesCounter++;
            labeledNodes[newNode] = 1;
        }
        if (activeNodesCounter!=0) {
        activeNodes.erase(activeNodes.begin());
            activeNodesCounter--;
        }
    }
}
}

float calculateProbState(float nearestNeighboorSum, int realization, float q, float p){
    return (1 - nearestNeighboorSum - realization * (1 - 2 * nearestNeighboorSum)) * q + (1 - q) * pow(p,realization) * pow(1 - p, 1 - realization);
}

bool switchState(const int* sample, int index, double q, float* p, vector<int>& nearestNeighbors, int sampleSize, const int* space){
    std::pair<float, int> sum = neirestSumOfSamples(index, sample, space, nearestNeighbors, sampleSize);
    float p_current = calculateProbState(sum.first, sample[index], q, p[index]);
    float p_new = calculateProbState(sum.first, 1 - sample[index], q, p[index]);
    for (int& NN : nearestNeighbors) {
        int sampledNode = index + NN;
        if (isInGrid(sampledNode, sampleSize) && space[sampledNode]) { 
            sum = neirestSumOfSamples(sampledNode, sample, space, nearestNeighbors, sampleSize);
            float sum_new = sum.first - sample[index] * 1.0/sum.second + (1 - sample[index]) * 1.0/sum.second;
            int realization = sample[sampledNode];
            p_current *= calculateProbState(sum.first, realization, q, p[sampledNode]);
            p_new *= calculateProbState(sum_new, realization, q, p[sampledNode]);
        }
    }
    return p_current/p_new < 1; 
}

class UniqueVector {
private:
    std::unordered_set<int> uniqueElements;
    std::vector<int> elements;

public:
    // Add an element to the collection
    bool add(int value) {
        if (uniqueElements.find(value) == uniqueElements.end()) {
            uniqueElements.insert(value);
            elements.push_back(value);
            return true; // Successfully added
        }
        return false; // Element already exists
    }

    // Get the element at a specific index
    int get(int index) const {
        if (index >= 0 && index < elements.size()) {
            return elements[index];
        }
        throw std::out_of_range("Index out of bounds");
    }

    // Get the number of elements
    int size() const {
        return elements.size();
    }
};

void addHeterogeneousVoxel(UniqueVector& indexesHeterogeneous, const int node, const int* sample, const int* labeledNodes, vector<int>& nearestNeighbors, const int sampleSize){

    std::pair<float, int> sum = neirestSumOfSamples(node, sample, labeledNodes, nearestNeighbors, sampleSize);
    if (sum.first != 0 && sum.first != 1) {
        indexesHeterogeneous.add(node);
    }
    else if (sum.first == 0 && sample[node] == 1) {
        indexesHeterogeneous.add(node);
    }
    else if(sum.first == 1 && sample[node] == 0) {
        indexesHeterogeneous.add(node);
    }
}

void reSampleShell_conv(int dim, int nOfElements, double q, int* sample, int* labeledNodes, float* p, int startIndex, int* shellSpace, int* sampleSize, int spaceSize, int numberOfIterations) {
    random_device rd;
    mt19937 gen(rd());
    uniform_real_distribution<double> dis(0.0, 1.0);
    uniform_int_distribution<int> disInt(0, 1);

    vector<int> activeNodes;
    activeNodes.push_back(startIndex);
    int activeNodesCounter = 1;
    vector<vector<int> > neirestNeighbors = getNearestNeighbors(dim);
    vector<int> neirestNeighbors1D = getNearestNeighbors1D(neirestNeighbors, sampleSize, dim);
    const int sampleSizeInt = multiplyElements(sampleSize, dim);
    labeledNodes[startIndex] = 1;
    UniqueVector indexesHeterogeneous;
    addHeterogeneousVoxel(indexesHeterogeneous, startIndex, sample, labeledNodes, neirestNeighbors1D, sampleSizeInt);
    for (int i = 0; i < nOfElements; ++i) {
        int Node = activeNodes[0];
        for (int& NN : neirestNeighbors1D) {
            int newNode = Node + NN;
            if (shellSpace[newNode] == 0 || labeledNodes[newNode] || newNode < 0 || newNode >= sampleSizeInt) {
                continue;
            }
            addHeterogeneousVoxel(indexesHeterogeneous, newNode, sample, labeledNodes, neirestNeighbors1D, sampleSizeInt);

            activeNodes.push_back(newNode);
            activeNodesCounter++;
            labeledNodes[newNode] = 1;
        }
        if (activeNodesCounter!=0) {
        activeNodes.erase(activeNodes.begin());
            activeNodesCounter--;
        }
    }
    printNoneZeroElements(sample,spaceSize);
    for (int i = 0; i < numberOfIterations; ++i) { // Try sorting the array and running through all the voxels.
        int index = indexesHeterogeneous.get(int(dis(gen) * indexesHeterogeneous.size()));
        bool switchStateBool = switchState(sample, index, q, p, neirestNeighbors1D, sampleSizeInt, shellSpace);
        if (switchStateBool == true)
            switchState(sample, index, q, p, neirestNeighbors1D, sampleSizeInt, shellSpace);
            sample[index] = 1 - sample[index];
            
        if (switchStateBool){
            for (int& NN : neirestNeighbors1D) {
                indexesHeterogeneous.add(index + NN);
            }
        }
    }
    std::cout<<"Aqii"<<std::endl;
    printNoneZeroElements(sample,spaceSize);
}




vector<int> createCube(int dim, int subDim) {
    vector<vector<vector<int>>> cube(dim, vector<vector<int>>(dim, vector<int>(dim, 0)));

    int start = (dim - subDim) / 2;
    int end = start + subDim;
    vector<int> space(dim * dim * dim, 0);
    int shape[3] = {dim, dim, dim};
    // Fill the 3x3x3 sub-cube with ones
    for (int i = start; i < end; ++i) {
        for (int j = start; j < end; ++j) {
            for (int k = start; k < end; ++k) {
                space[to1DIndex({i, j, k}, shape, 3)] = 1;
            }
        }
    }

    return space;
}

int main() {
    int dim = 3;
    // int nOfElements = 10;
    double q = 0.5;

    // Example 3D space
    int sampleSize[] = {11, 11, 11};
    int startIndex = to1DIndex({5,5,5}, sampleSize, dim);
    int sampleSizeInt = 11*11*11;

    int* sample = new int[sampleSizeInt];
    int* labeledNodes = new int[sampleSizeInt];
    float* p = new float[sampleSizeInt];
    std::vector<int> vec = createCube(11,11);
    int* shellSpace = vec.data();

    for (int i = 0; i < sampleSizeInt; ++i) {
        sample[i] = 0.0f;
        p[i] = 0.5f;
    }
    sample[startIndex] = 1;
    // sampleShell(dim, nOfElements, q, sample, labeledNodes, p, startIndex, shellSpace, sampleSize, vec.size());
    reSampleShell_conv(dim, sampleSizeInt, q, sample, labeledNodes, p, startIndex, shellSpace, sampleSize, vec.size(), 1000);
    // Output some results to verify
    // cout << "Sample values after running sampleShell:" << endl;
    // for (int i = 0; i < sampleSizeInt; ++i) {
    //     cout << "Index " << i << ": " << sample[i] << '\t' << endl;
    // }

    return 0;
}