interatomicDist-FW-SeqExample.cu

#include <stdio.h>
#include <cstdlib>

#define ATOMS_IN_MOLECULE 10

//These bounds breaks the physical realism of molecule sizes
#define UNKNOWN_UPPER_BOUND 1.0f
#define UNKNOWN_LOWER_BOUND 0.01f
#define UNKNOWN_ACTUAL_VALUE 0.0f

//Struct representing the weight of an edge in the graph
typedef struct __align__(16)
{
    float upper;
    float lower;
    float actual;
} weight;

float generateRandFloatInWeightRange();
int generateRandInt(int maxRageFromZero);
void generateGraph(weight adjMatrix[][ATOMS_IN_MOLECULE], int atomsRequested);
void printGraphMatrix(weight *adjMatrix, int atomsHeightWidth, bool fwPass);
void floydWarshallTheChemist(weight adjMatrix[][ATOMS_IN_MOLECULE], int atomsHeightWidth);

int main(){
    //Initializing the RNG seed:
    srand(time(NULL));

    //Create matrix of size Atoms x Atoms, for a molecule made of a Atoms
    weight moleculeGraph[ATOMS_IN_MOLECULE][ATOMS_IN_MOLECULE];
    generateGraph(moleculeGraph, ATOMS_IN_MOLECULE);

    printGraphMatrix(moleculeGraph[0], ATOMS_IN_MOLECULE, false);

    //Computer goes beep boop beep, and now we might have a better idea as to how the molecule looks
    floydWarshallTheChemist(moleculeGraph, ATOMS_IN_MOLECULE);
    printf("\n\nNewly Computed Graph: \n\n");
    printGraphMatrix(moleculeGraph[0], ATOMS_IN_MOLECULE, true);

    return 0;
}


float generateRandFloatInWeightRange(){
    int randNumPull = rand();
    while (randNumPull == 0){
        randNumPull = rand();
    }

    float randNum = ((float)(randNumPull % 101)) / 100;
    return randNum;
}

int generateRandInt(int maxRageFromZero){
    return (rand() % maxRageFromZero);
}

//Function to generate graph of molecule
void generateGraph(weight adjMatrix[][ATOMS_IN_MOLECULE], int atomsRequested){
    weight noEdgeWeight;
    noEdgeWeight.upper = 0.0f;
    noEdgeWeight.lower = 0.0f;
    noEdgeWeight.actual = 0.0f;

    weight unknownMessurementWeight;
    unknownMessurementWeight.upper = UNKNOWN_UPPER_BOUND;
    unknownMessurementWeight.lower = UNKNOWN_LOWER_BOUND;
    unknownMessurementWeight.actual = UNKNOWN_ACTUAL_VALUE;

    //Initializing the molecule/graph with unknown messurements/weights
    for(int row = 0; row < atomsRequested; row++){
        for (int column = 0; column < atomsRequested; column++)
        {
            //Since the graph is undirected, there should be no edges from a vertex to itself
            if (row == column)
            {
                //Performing 2D arr pointer arithmetic
                adjMatrix[row][column] = noEdgeWeight;
            }
            else
            {
                //Setting the rest of the edges to their unknow messurement boundaries
                adjMatrix[row][column] = unknownMessurementWeight;
            }
            
        }   
    }

    //For demo's sake lets say we know the actual distance between 20% of the atom pairs
    int atomPairsKnown = ((float)(atomsRequested * atomsRequested)) * 0.2f;
    printf("Atom Pairs Known: %d\n", atomPairsKnown);
    for (int i = 0; i < atomPairsKnown; i++)
    {
        int randRow = generateRandInt(atomsRequested);
        int randColumn = generateRandInt(atomsRequested);
        
        if (randRow != randColumn)
        {
            adjMatrix[randRow][randColumn].actual = generateRandFloatInWeightRange();
            adjMatrix[randRow][randColumn].upper = adjMatrix[randRow][randColumn].actual;
            adjMatrix[randRow][randColumn].lower = adjMatrix[randRow][randColumn].actual;

            //Enforcing Symmetry
            adjMatrix[randColumn][randRow].actual = adjMatrix[randRow][randColumn].actual;
            adjMatrix[randColumn][randRow].upper = adjMatrix[randRow][randColumn].actual;
            adjMatrix[randColumn][randRow].lower = adjMatrix[randRow][randColumn].actual;
        }
        
    }

}

void printGraphMatrix(weight *adjMatrix, int atomsHeightWidth, bool fwPass){
    int atomsInTotal = 0;
    int atomsPossiblyUnknown = 0;
    int atomsNotConnected = 0;

    for(int row = 0; row < atomsHeightWidth; row++){
        for (int column = 0; column < atomsHeightWidth; column++)
        {
            float currUpper = (*(adjMatrix + (row * atomsHeightWidth + column))).upper;
            float currLower = (*(adjMatrix + (row * atomsHeightWidth + column))).lower;
            float currActual = (*(adjMatrix + (row * atomsHeightWidth + column))).actual;

            float symCurrUpper = (*(adjMatrix + (column * atomsHeightWidth + row))).upper;
            float symCurrLower = (*(adjMatrix + (column * atomsHeightWidth + row))).lower;
            float symCurrActual = (*(adjMatrix + (column * atomsHeightWidth + row))).actual;

            printf("adjMatrix[%d][%d] => (Upper: %f, Lower: %f, Actual: %f)", row, column, currUpper, currLower, currActual);

            //Checks if we by accidently break physical bounds
            if (currLower > currUpper)
            {
                printf("[BREAKS PHYSICS]");
            }
            //Checks if symmetry in the undirected graph is upheld
            if ((currUpper != symCurrUpper) || (currLower != symCurrLower) || (currActual != symCurrActual))
            {
                printf("[BREAKS SYMMETRY]");
            }
            if (currUpper == 1.0f && currLower == 0.01f && currActual == 0.0f)
            {
                //Checks wheter the graph has has a pass in the bound-smoothing algorithm
                if (fwPass == false)
                {
                    printf("[POSSIBLY UNKNOWN]");
                    atomsPossiblyUnknown++;
                }
                else
                {
                    //Means that no direct or transitive data including this edge existed
                    printf("[UNABLE TO INFER FROM AVAILABLE DATA]");
                    atomsPossiblyUnknown++;
                }
                
                
            }
            if (currUpper == 0.0f && currLower == 0.0f && currActual == 0.0f)
            {
                atomsNotConnected++;
            }
            
            
            printf("\n");

            atomsInTotal++;
        }
        
    }
    printf("Atoms in Total: %d\n", atomsInTotal);
    printf("Atoms Possibly Unknown in Total: %d\n", atomsPossiblyUnknown);
    printf("Atom Pairs Not Connected: %d\n", atomsNotConnected);
}

void floydWarshallTheChemist(weight adjMatrix[][ATOMS_IN_MOLECULE], int atomsHeightWidth){
    for (int atom_k = 0; atom_k < atomsHeightWidth; atom_k++)
    {
        for (int atom_i = 0; atom_i < atomsHeightWidth; atom_i++)
        {
            for (int atom_j = 0; atom_j < atomsHeightWidth; atom_j++)
            {
                //Skipping the nodes attempting to compute distance bounds to themselves, or those who have true data
                if ((atom_i == atom_j) || (adjMatrix[atom_i][atom_j].actual != 0.0f))
                {
                    continue;
                }
                

                //Computing upper bound
                float currUpperDistanceFromItoJ = adjMatrix[atom_i][atom_j].upper;
                float upperDistanceFromItoK = adjMatrix[atom_i][atom_k].upper;
                float upperDistanceFromKtoJ = adjMatrix[atom_k][atom_j].upper;

                //New upper bound for the atom pair
                //Seeing wheter the current route is shorther or if the route through other atoms is shorter
                adjMatrix[atom_i][atom_j].upper = (currUpperDistanceFromItoJ <= (upperDistanceFromItoK + upperDistanceFromKtoJ)) ? currUpperDistanceFromItoJ : (upperDistanceFromItoK + upperDistanceFromKtoJ);

            }
            
        }
        
    }


    for (int atom_k = 0; atom_k < atomsHeightWidth; atom_k++)
    {
        for (int atom_i = 0; atom_i < atomsHeightWidth; atom_i++)
        {
            for (int atom_j = 0; atom_j < atomsHeightWidth; atom_j++)
            {
                //Skipping the nodes attempting to compute distance bounds to themselves, or those who have true data
                if ((atom_i == atom_j) || (adjMatrix[atom_i][atom_j].actual != 0.0f))
                {
                    continue;
                }

                //Computing lower bound
                float currLowerDistanceFromItoJ = adjMatrix[atom_i][atom_j].lower;

                //New lower bound for the atom pair
                //Seeing wheter the current route is longest or if the route through other atoms is longer
                adjMatrix[atom_i][atom_j].lower = max(currLowerDistanceFromItoJ, max((adjMatrix[atom_i][atom_k].lower - adjMatrix[atom_k][atom_j].upper), (adjMatrix[atom_k][atom_j].lower - adjMatrix[atom_i][atom_k].upper)));

                //Ensure that we dont break physics by surpassing the upper bounds
                if (adjMatrix[atom_i][atom_j].lower > adjMatrix[atom_i][atom_j].upper)
                {
                    adjMatrix[atom_i][atom_j].lower = adjMatrix[atom_i][atom_j].upper;
                }
                
            }
            
        }
        
    }
    
}