codebase / stimlib

Commit ea2056c0efdd9b0b6b0fa8dba8a2559b76677595

Authored by cherub
1 parent 5687cf1b

added flow.h

Showing 1 changed file with 408 additions and 0 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/flow.h 0 → 100644
  Show/Hide comments   View file @ea2056c
  1 +#pragma once
  2 +#include <fstream> // Required for ofstream, etc.
  3 +#include <iomanip> // Required for setw
  4 +#include <iostream> // Required for cout, cin, etc.
  5 +#include <tuple> // Required for returning multiple values from a function
  6 +//#include <vector>
  7 +//#include <string>
  8 +//#include <stdlib.h> // Required to remove ambiguous error for cout, cin, etc.
  9 +
  10 +using namespace std;
  11 +
  12 +
  13 +class flow
  14 +{
  15 +public:
  16 + void backupToTxt(unsigned int nL, float **D, char filename[]);
  17 + tuple<int, int> copySrcDesRadLen(char filename[]);
  18 + void copyToArray(int *src, int *dest, float *radii, float *len);
  19 + int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes);
  20 + void inversion(float **a, int n, float **b);
  21 + float determinant(float **a, int n);
  22 +
  23 +protected:
  24 + int minor(float **src, float **dest, int row, int col, int order);
  25 +};
  26 +
  27 +/* Function to find the dangle nodes in a network */
  28 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  29 +// Modified by Cherub P. Harder on 8/12/2015
  30 +int flow::getDangleNodes(int datarow, int numNodes, int *column1, int *column2, int *dangleNodes)
  31 +{
  32 + int count = datarow, diff1 = 0, diff2 = 0, numPress = 0, st = 0;
  33 +
  34 + // Find matching nodes within column2
  35 + for( int i = 0; i < count; i++ )
  36 + {
  37 + for( int y = i+1; y < datarow; y++ )
  38 + {
  39 + if( column2[i] == column2[y] ) // Is there a match?
  40 + {
  41 + st = column2[i]; // Save the matching node
  42 + cout << endl << column2[i] << " = " << column2[y] << endl; // Display the matching nodes
  43 + memmove(column2+i, column2+i+1, (datarow-(i+1)) * sizeof(column2[0])); // Move up the rows
  44 + // taking the places of the rows before them starting
  45 + // with where the matching node is located
  46 + column2[datarow-1] = st; // Place the matching node at the very end of the array--
  47 + // this is for comparison purpose so that the other match-
  48 + // ing node will be moved as well and then omitted later.
  49 + diff1++; // Count the matching node
  50 +
  51 + // Display the updated array (with the matching node moved to the bottommost row)
  52 + cout << "Updated array:" << endl;
  53 + for( int k = 0; k < datarow; k++ )
  54 + cout << column2[k] << endl;
  55 +
  56 + // Decrement the counters
  57 + // NOTE: The counters need to be decremented because the rows had been moved up, so the same
  58 + // locations need to be read again because they contain different values now after the move.
  59 + i--; // Decrement i to read the node that took over the place
  60 + // of the matching node. Otherwise, it will be skipped.
  61 + y--; // Decrement y to read the following node for comparison
  62 + count--; // The maximum count need to be decremented so that the
  63 + // matching nodes that had been moved will not be read again.
  64 + // However, the maximum count (datarow) for finding a match
  65 + // will not be decremented because the remaining matching
  66 + // node that has not been moved yet needs to be moved and
  67 + // the only way to do that is to match it with its duplicate.
  68 + }
  69 + }
  70 + }
  71 +
  72 + // Store the nodes that have no duplicates
  73 + // NOTE: This will only save the nodes that have not been moved to the bottom.
  74 + cout << "\ndangleNodes array:" << endl;
  75 + for( int j = 0; j < datarow-diff1; j++ )
  76 + {
  77 + dangleNodes[numPress] = column2[j];
  78 + cout << dangleNodes[j] << endl; // DELETE!!!
  79 + numPress++; // Count the non-duplicated node
  80 + }
  81 +
  82 + // Find if the non-duplicated nodes have a match from column1
  83 + count = datarow-diff1; // Reinitialize the counter
  84 +
  85 + for( int i = 0; i < count; i++ )
  86 + {
  87 + for( int j = 0; j < datarow; j++ )
  88 + {
  89 + if( dangleNodes[i] == column1[j] ) // Is there a match?
  90 + {
  91 + st = column1[j]; // Save the matching node
  92 + cout << endl << dangleNodes[i] << " = " << column1[j] << endl; // Display the matching nodes
  93 + memmove(dangleNodes+i, dangleNodes+i+1, (datarow-diff1-(i+1)) * sizeof(dangleNodes[0]));
  94 + dangleNodes[count-1] = st; // Move the matching node to the bottom of the array
  95 + diff2++; // Count the matching node
  96 +
  97 + // Display the updated array
  98 + cout << "Updated dangleNodes array:" << endl;
  99 + for( int k = 0; k < count-1; k++ )
  100 + {
  101 + cout << dangleNodes[k] << endl;
  102 + }
  103 +
  104 + // Decrement the counters
  105 + i--;
  106 + j--;
  107 + count--;
  108 + numPress--; // Decrement to the exact number of dangle nodes
  109 + }
  110 + }
  111 + }
  112 +
  113 + return numPress; // Number of dangle nodes
  114 +}
  115 +
  116 +
  117 +// Function to make a backup copy of the contents of a matrix to a .txt file
  118 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  119 +void flow::backupToTxt(unsigned int nL, float **D, char filename[])
  120 +{
  121 + ofstream output_file(filename);
  122 +
  123 + for( unsigned int i = 0; i < nL; i++ )
  124 + {
  125 + for( int j = 0; j < 4; j++ )
  126 + {
  127 + if( j < 3 )
  128 + output_file << D[i][j] << "\t";
  129 +
  130 + else
  131 + output_file << D[i][j];
  132 + }
  133 +
  134 + output_file << "\n";
  135 + }
  136 +
  137 + output_file.close( );
  138 +}
  139 +
  140 +
  141 +// Function to make separate copies of the source nodes, destination nodes, radii, and lengths
  142 +// Created by Cherub P. Harder (8/10/2015), U of Houston
  143 +tuple<int, int> flow::copySrcDesRadLen(char filename[])
  144 +{
  145 + int cnt = 0, numElements = 0, numNodes = 0;
  146 + float number = 0.0;
  147 + ofstream srcData("srcCol.txt"); // A .txt file to store the source nodes
  148 + ofstream desData("destCol.txt"); // A .txt file to store the destination nodes
  149 + ofstream radiiData("radii.txt"); // A .txt file to store the radii
  150 + ofstream lenData("lengths.txt"); // A .txt file to store the lengths
  151 + FILE *fp = fopen(filename, "r"); // Create a variable of type FILE* and open the file using
  152 + // the fopen function and assign the file to the variable
  153 + // Check if the file exists
  154 + if(fp == NULL) // Alternative: if(!fp)
  155 + {
  156 + printf("Error! File does not exist.\n");
  157 + getchar( ); // Pause
  158 + exit(-1); // NOTE: Must include stdlib.h.
  159 + }
  160 +
  161 + // Store data to their respective .txt files
  162 + while(fscanf(fp, "%f", &number) == 1)
  163 + {
  164 + cnt++; // Increment counter
  165 +
  166 + // Store to srcCol.txt
  167 + if(cnt == 1)
  168 + srcData << number << endl;
  169 +
  170 + // Store to destCol.txt
  171 + if(cnt == 2)
  172 + desData << number << endl;
  173 +
  174 + // Save the current number of nodes
  175 + if(cnt < 3)
  176 + {
  177 + if(number > numNodes)
  178 + numNodes = (int)number;
  179 + }
  180 +
  181 + // Store to radii.txt
  182 + if(cnt == 3)
  183 + radiiData << number << endl;
  184 +
  185 + // Store to lengths.txt
  186 + if(cnt == 4)
  187 + {
  188 + lenData << number << endl;
  189 +
  190 + numElements++; // Count the elements
  191 + cnt = 0; // Reset counter
  192 + }
  193 + }
  194 +
  195 + srcData.close( );
  196 + desData.close( );
  197 + radiiData.close( );
  198 + lenData.close( );
  199 +
  200 + return make_tuple(numNodes, numElements); // Return two values
  201 +}
  202 +
  203 +
  204 +// Function to copy data for .txt files to their respective arrays
  205 +// Created by Cherub P. Harder (8/11/2015), U of Houston
  206 +void flow::copyToArray(int *src, int *dest, float *radii, float *len)
  207 +{
  208 + int v = 0;
  209 + float tmp = 0, R = 0, L = 0;
  210 +
  211 + // Store source node values to the array src
  212 + ifstream readSrc("srcCol.txt");
  213 +
  214 + while( readSrc >> tmp )
  215 + {
  216 + src[v] = (int)tmp;
  217 + v++;
  218 + }
  219 +
  220 + readSrc.close( );
  221 +
  222 + // Store destination node values to the array dest
  223 + v = 0; // Reset counter
  224 + ifstream readDest("destCol.txt");
  225 +
  226 + while( readDest >> tmp )
  227 + {
  228 + dest[v] = (int)tmp;
  229 + v++;
  230 + }
  231 +
  232 + readDest.close( );
  233 +
  234 + // Store radius values to the array radii
  235 + v = 0; // Reset counter
  236 + ifstream readRad("radii.txt");
  237 +
  238 + while( readRad >> tmp )
  239 + {
  240 + radii[v] = tmp;
  241 + v++;
  242 + }
  243 +
  244 + readRad.close( );
  245 +
  246 + // Store length values to the array len
  247 + v = 0; // Reset counter
  248 + ifstream readLen("lengths.txt");
  249 +
  250 + while( readLen >> tmp )
  251 + {
  252 + len[v] = tmp;
  253 + v++;
  254 + }
  255 +
  256 + readLen.close( );
  257 +}
  258 +
  259 +
  260 +// Function to find the inverse of a square matrix
  261 +void flow::inversion(float **a, int n, float **b)
  262 +{
  263 + // Get the determinant of A
  264 + float det = (float)(1.0/determinant(a, n));
  265 +
  266 + // Memory allocation
  267 + float *tmp = new float[(n-1) * (n-1)];
  268 + float **m = new float * [n-1];
  269 + for( int i = 0; i < n-1; i++ )
  270 + m[i] = tmp + ( i * (n-1) );
  271 +
  272 + for( int j = 0; j < n; j++)
  273 + {
  274 + for( int i = 0; i < n; i++ )
  275 + {
  276 + // Get the cofactor (matrix) of a(j,i)
  277 + minor(a, m, j, i, n);
  278 +// coFactor(a, 3, b);
  279 + b[i][j] = det * determinant( m, n-1 );
  280 + if( (i+j)%2 == 1 )
  281 + b[i][j] = -b[i][j];
  282 + }
  283 + }
  284 +
  285 + // Release memory
  286 + // Delete [] minor[0];
  287 + delete [] tmp;
  288 + delete [] m;
  289 +}
  290 +
  291 +
  292 +// Function to find the determinant of a matrix using recursion
  293 +// Contribution by Edward Popko
  294 +// Modified by Cherub P. Harder (7/15/2015), U of Houston
  295 +// Arguments: a(double **) - pointer to a pointer of an arbitrary square matrix
  296 +// n(int) - dimension of the square matrix
  297 +float flow::determinant(float **a, int n)
  298 +{
  299 + int i, j, j1, j2; // General loop and matrix subscripts
  300 + float det = 0; // Initialize determinant
  301 + float **m = NULL; // Pointer to pointer to implement 2D square array
  302 +
  303 + // Display contents of matrix C (DELETE!!!)
  304 +/* std::cout << "\nThe updated matrix C:\n";
  305 + for( int j = 0; j < n; ++j )
  306 + {
  307 + std::cerr << "\t";
  308 +
  309 + for( int k = 0; k < n; ++k )
  310 + std::cerr << left << setw(15) << a[j][k];
  311 +
  312 + std::cerr << endl;
  313 + }
  314 +
  315 + getchar(); // DELETE!!!*/
  316 +
  317 + if(n < 1) { } // Error condition - should never get here
  318 +
  319 + else if (n == 1) // Should never get here
  320 + {
  321 + det = a[0][0];
  322 + }
  323 +
  324 + else if(n == 2) // Basic 2x2 sub-matrix determinate definition
  325 + { // When n == 2, this ends the recursion series
  326 + det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
  327 + }
  328 + // Recursion continues, solve next sub-matrix
  329 + else // Solve the next minor by building a sub-matrix
  330 + {
  331 + det = 0; // Initialize determinant of sub-matrix
  332 +
  333 + for (j1 = 0; j1 < n; j1++) // For each column in sub-matrix get space for the
  334 + { // pointer list
  335 + m = (float **) malloc((n-1) * sizeof(float *));
  336 +
  337 + for (i = 0; i < n-1; i++)
  338 + m[i] = (float *) malloc((n-1)* sizeof(float));
  339 + // i[0][1][2][3] first malloc
  340 + // m -> + + + + space for 4 pointers
  341 + // | | | | j second malloc
  342 + // | | | +-> _ _ _ [0] pointers to
  343 + // | | +----> _ _ _ [1] and memory for
  344 + // | +-------> _ a _ [2] 4 doubles
  345 + // +----------> _ _ _ [3]
  346 + //
  347 + // a[1][2]
  348 + // Build sub-matrix with minor elements excluded
  349 +
  350 + for (i = 1; i < n; i++)
  351 + {
  352 + j2 = 0 ; // Start at first sum-matrix column position
  353 + // Loop to copy source matrix less one column
  354 + for (j = 0; j < n; j++)
  355 + {
  356 + if (j == j1) continue; // Do not copy the minor column element
  357 +
  358 + m[i-1][j2] = a[i][j]; // Copy source element into new sub-matrix
  359 + // i-1 because new sub-matrix is one row
  360 + // (and column) smaller with excluded minors
  361 + j2++; // Move to next sub-matrix column position
  362 + }
  363 + }
  364 +
  365 + det += (float)pow(-1.0, 1.0 + j1 + 1.0) * a[0][j1] * determinant(m, n-1);
  366 + // Sum x raised to y power
  367 + // recursively get determinant of next
  368 + // sub-matrix which is now one
  369 + // row & column smaller
  370 +
  371 + for (i = 0; i < n-1; i++) free(m[i]); // Free the storage allocated to
  372 + // this minor's set of pointers
  373 + free(m); // Free the storage for the original
  374 + // pointer to pointer
  375 + }
  376 + }
  377 +
  378 + return(det);
  379 +}
  380 +
  381 +
  382 +// Function to calculate the cofactor of element (row, col)
  383 +int flow::minor(float **src, float **dest, int row, int col, int order)
  384 +{
  385 + // Indicate which col and row is being copied to dest
  386 + int colCount=0,rowCount=0;
  387 +
  388 + for(int i = 0; i < order; i++ )
  389 + {
  390 + if( i != row )
  391 + {
  392 + colCount = 0;
  393 + for(int j = 0; j < order; j++ )
  394 + {
  395 + // When j is not the element
  396 + if( j != col )
  397 + {
  398 + dest[rowCount][colCount] = src[i][j];
  399 + colCount++;
  400 + }
  401 + }
  402 +
  403 + rowCount++;
  404 + }
  405 + }
  406 +
  407 + return 1;
  408 +}
0 409 \ No newline at end of file
... ...