codebase / stimlib

Commit 47b35dbebda8e9756b8b66b40cfc2f12371d0476

Authored by David Mayerich
2 parents faef7718 dae97842

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

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