From ea2056c0efdd9b0b6b0fa8dba8a2559b76677595 Mon Sep 17 00:00:00 2001 From: cherub Date: Thu, 13 Aug 2015 15:35:22 -0500 Subject: [PATCH] added flow.h --- stim/biomodels/flow.h | 408 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 408 insertions(+), 0 deletions(-) create mode 100644 stim/biomodels/flow.h diff --git a/stim/biomodels/flow.h b/stim/biomodels/flow.h new file mode 100644 index 0000000..8f080fd --- /dev/null +++ b/stim/biomodels/flow.h @@ -0,0 +1,408 @@ +#pragma once +#include // Required for ofstream, etc. +#include // Required for setw +#include // Required for cout, cin, etc. +#include // Required for returning multiple values from a function +//#include +//#include +//#include // Required to remove ambiguous error for cout, cin, etc. + +using namespace std; + + +class flow +{ +public: + void backupToTxt(unsigned int nL, float **D, char filename[]); + tuple copySrcDesRadLen(char filename[]); + void copyToArray(int *src, int *dest, float *radii, float *len); + int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes); + void inversion(float **a, int n, float **b); + float determinant(float **a, int n); + +protected: + int minor(float **src, float **dest, int row, int col, int order); +}; + +/* Function to find the dangle nodes in a network */ +// Created by Cherub P. Harder (8/10/2015), U of Houston +// Modified by Cherub P. Harder on 8/12/2015 +int flow::getDangleNodes(int datarow, int numNodes, int *column1, int *column2, int *dangleNodes) +{ + int count = datarow, diff1 = 0, diff2 = 0, numPress = 0, st = 0; + + // Find matching nodes within column2 + for( int i = 0; i < count; i++ ) + { + for( int y = i+1; y < datarow; y++ ) + { + if( column2[i] == column2[y] ) // Is there a match? + { + st = column2[i]; // Save the matching node + cout << endl << column2[i] << " = " << column2[y] << endl; // Display the matching nodes + memmove(column2+i, column2+i+1, (datarow-(i+1)) * sizeof(column2[0])); // Move up the rows + // taking the places of the rows before them starting + // with where the matching node is located + column2[datarow-1] = st; // Place the matching node at the very end of the array-- + // this is for comparison purpose so that the other match- + // ing node will be moved as well and then omitted later. + diff1++; // Count the matching node + + // Display the updated array (with the matching node moved to the bottommost row) + cout << "Updated array:" << endl; + for( int k = 0; k < datarow; k++ ) + cout << column2[k] << endl; + + // Decrement the counters + // NOTE: The counters need to be decremented because the rows had been moved up, so the same + // locations need to be read again because they contain different values now after the move. + i--; // Decrement i to read the node that took over the place + // of the matching node. Otherwise, it will be skipped. + y--; // Decrement y to read the following node for comparison + count--; // The maximum count need to be decremented so that the + // matching nodes that had been moved will not be read again. + // However, the maximum count (datarow) for finding a match + // will not be decremented because the remaining matching + // node that has not been moved yet needs to be moved and + // the only way to do that is to match it with its duplicate. + } + } + } + + // Store the nodes that have no duplicates + // NOTE: This will only save the nodes that have not been moved to the bottom. + cout << "\ndangleNodes array:" << endl; + for( int j = 0; j < datarow-diff1; j++ ) + { + dangleNodes[numPress] = column2[j]; + cout << dangleNodes[j] << endl; // DELETE!!! + numPress++; // Count the non-duplicated node + } + + // Find if the non-duplicated nodes have a match from column1 + count = datarow-diff1; // Reinitialize the counter + + for( int i = 0; i < count; i++ ) + { + for( int j = 0; j < datarow; j++ ) + { + if( dangleNodes[i] == column1[j] ) // Is there a match? + { + st = column1[j]; // Save the matching node + cout << endl << dangleNodes[i] << " = " << column1[j] << endl; // Display the matching nodes + memmove(dangleNodes+i, dangleNodes+i+1, (datarow-diff1-(i+1)) * sizeof(dangleNodes[0])); + dangleNodes[count-1] = st; // Move the matching node to the bottom of the array + diff2++; // Count the matching node + + // Display the updated array + cout << "Updated dangleNodes array:" << endl; + for( int k = 0; k < count-1; k++ ) + { + cout << dangleNodes[k] << endl; + } + + // Decrement the counters + i--; + j--; + count--; + numPress--; // Decrement to the exact number of dangle nodes + } + } + } + + return numPress; // Number of dangle nodes +} + + +// Function to make a backup copy of the contents of a matrix to a .txt file +// Created by Cherub P. Harder (8/10/2015), U of Houston +void flow::backupToTxt(unsigned int nL, float **D, char filename[]) +{ + ofstream output_file(filename); + + for( unsigned int i = 0; i < nL; i++ ) + { + for( int j = 0; j < 4; j++ ) + { + if( j < 3 ) + output_file << D[i][j] << "\t"; + + else + output_file << D[i][j]; + } + + output_file << "\n"; + } + + output_file.close( ); +} + + +// Function to make separate copies of the source nodes, destination nodes, radii, and lengths +// Created by Cherub P. Harder (8/10/2015), U of Houston +tuple flow::copySrcDesRadLen(char filename[]) +{ + int cnt = 0, numElements = 0, numNodes = 0; + float number = 0.0; + ofstream srcData("srcCol.txt"); // A .txt file to store the source nodes + ofstream desData("destCol.txt"); // A .txt file to store the destination nodes + ofstream radiiData("radii.txt"); // A .txt file to store the radii + ofstream lenData("lengths.txt"); // A .txt file to store the lengths + FILE *fp = fopen(filename, "r"); // Create a variable of type FILE* and open the file using + // the fopen function and assign the file to the variable + // Check if the file exists + if(fp == NULL) // Alternative: if(!fp) + { + printf("Error! File does not exist.\n"); + getchar( ); // Pause + exit(-1); // NOTE: Must include stdlib.h. + } + + // Store data to their respective .txt files + while(fscanf(fp, "%f", &number) == 1) + { + cnt++; // Increment counter + + // Store to srcCol.txt + if(cnt == 1) + srcData << number << endl; + + // Store to destCol.txt + if(cnt == 2) + desData << number << endl; + + // Save the current number of nodes + if(cnt < 3) + { + if(number > numNodes) + numNodes = (int)number; + } + + // Store to radii.txt + if(cnt == 3) + radiiData << number << endl; + + // Store to lengths.txt + if(cnt == 4) + { + lenData << number << endl; + + numElements++; // Count the elements + cnt = 0; // Reset counter + } + } + + srcData.close( ); + desData.close( ); + radiiData.close( ); + lenData.close( ); + + return make_tuple(numNodes, numElements); // Return two values +} + + +// Function to copy data for .txt files to their respective arrays +// Created by Cherub P. Harder (8/11/2015), U of Houston +void flow::copyToArray(int *src, int *dest, float *radii, float *len) +{ + int v = 0; + float tmp = 0, R = 0, L = 0; + + // Store source node values to the array src + ifstream readSrc("srcCol.txt"); + + while( readSrc >> tmp ) + { + src[v] = (int)tmp; + v++; + } + + readSrc.close( ); + + // Store destination node values to the array dest + v = 0; // Reset counter + ifstream readDest("destCol.txt"); + + while( readDest >> tmp ) + { + dest[v] = (int)tmp; + v++; + } + + readDest.close( ); + + // Store radius values to the array radii + v = 0; // Reset counter + ifstream readRad("radii.txt"); + + while( readRad >> tmp ) + { + radii[v] = tmp; + v++; + } + + readRad.close( ); + + // Store length values to the array len + v = 0; // Reset counter + ifstream readLen("lengths.txt"); + + while( readLen >> tmp ) + { + len[v] = tmp; + v++; + } + + readLen.close( ); +} + + +// Function to find the inverse of a square matrix +void flow::inversion(float **a, int n, float **b) +{ + // Get the determinant of A + float det = (float)(1.0/determinant(a, n)); + + // Memory allocation + float *tmp = new float[(n-1) * (n-1)]; + float **m = new float * [n-1]; + for( int i = 0; i < n-1; i++ ) + m[i] = tmp + ( i * (n-1) ); + + for( int j = 0; j < n; j++) + { + for( int i = 0; i < n; i++ ) + { + // Get the cofactor (matrix) of a(j,i) + minor(a, m, j, i, n); +// coFactor(a, 3, b); + b[i][j] = det * determinant( m, n-1 ); + if( (i+j)%2 == 1 ) + b[i][j] = -b[i][j]; + } + } + + // Release memory + // Delete [] minor[0]; + delete [] tmp; + delete [] m; +} + + +// Function to find the determinant of a matrix using recursion +// Contribution by Edward Popko +// Modified by Cherub P. Harder (7/15/2015), U of Houston +// Arguments: a(double **) - pointer to a pointer of an arbitrary square matrix +// n(int) - dimension of the square matrix +float flow::determinant(float **a, int n) +{ + int i, j, j1, j2; // General loop and matrix subscripts + float det = 0; // Initialize determinant + float **m = NULL; // Pointer to pointer to implement 2D square array + + // Display contents of matrix C (DELETE!!!) +/* std::cout << "\nThe updated matrix C:\n"; + for( int j = 0; j < n; ++j ) + { + std::cerr << "\t"; + + for( int k = 0; k < n; ++k ) + std::cerr << left << setw(15) << a[j][k]; + + std::cerr << endl; + } + + getchar(); // DELETE!!!*/ + + if(n < 1) { } // Error condition - should never get here + + else if (n == 1) // Should never get here + { + det = a[0][0]; + } + + else if(n == 2) // Basic 2x2 sub-matrix determinate definition + { // When n == 2, this ends the recursion series + det = a[0][0] * a[1][1] - a[1][0] * a[0][1]; + } + // Recursion continues, solve next sub-matrix + else // Solve the next minor by building a sub-matrix + { + det = 0; // Initialize determinant of sub-matrix + + for (j1 = 0; j1 < n; j1++) // For each column in sub-matrix get space for the + { // pointer list + m = (float **) malloc((n-1) * sizeof(float *)); + + for (i = 0; i < n-1; i++) + m[i] = (float *) malloc((n-1)* sizeof(float)); + // i[0][1][2][3] first malloc + // m -> + + + + space for 4 pointers + // | | | | j second malloc + // | | | +-> _ _ _ [0] pointers to + // | | +----> _ _ _ [1] and memory for + // | +-------> _ a _ [2] 4 doubles + // +----------> _ _ _ [3] + // + // a[1][2] + // Build sub-matrix with minor elements excluded + + for (i = 1; i < n; i++) + { + j2 = 0 ; // Start at first sum-matrix column position + // Loop to copy source matrix less one column + for (j = 0; j < n; j++) + { + if (j == j1) continue; // Do not copy the minor column element + + m[i-1][j2] = a[i][j]; // Copy source element into new sub-matrix + // i-1 because new sub-matrix is one row + // (and column) smaller with excluded minors + j2++; // Move to next sub-matrix column position + } + } + + det += (float)pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * determinant(m, n-1); + // Sum x raised to y power + // recursively get determinant of next + // sub-matrix which is now one + // row & column smaller + + for (i = 0; i < n-1; i++) free(m[i]); // Free the storage allocated to + // this minor's set of pointers + free(m); // Free the storage for the original + // pointer to pointer + } + } + + return(det); +} + + +// Function to calculate the cofactor of element (row, col) +int flow::minor(float **src, float **dest, int row, int col, int order) +{ + // Indicate which col and row is being copied to dest + int colCount=0,rowCount=0; + + for(int i = 0; i < order; i++ ) + { + if( i != row ) + { + colCount = 0; + for(int j = 0; j < order; j++ ) + { + // When j is not the element + if( j != col ) + { + dest[rowCount][colCount] = src[i][j]; + colCount++; + } + } + + rowCount++; + } + } + + return 1; +} \ No newline at end of file -- libgit2 0.21.4