codebase / stimlib

Commit 28b882ce020c50db81326f56c344c2aa0f4d1b84

Authored by Jiaming Guo
1 parent 800ff264

new flow header

Showing 2 changed files with 645 additions and 372 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/flow.h
  Show/Hide comments   View file @28b882c
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;
  1 +#ifndef STIM_FLOW_H
  2 +#define STIM_FLOW_H
  3 +
  4 +#include <vector>
  5 +#include <algorithm>
  6 +
  7 +//STIM include
  8 +#include <stim/math/vec3.h>
  9 +#include <stim/parser/arguments.h>
  10 +#include <stim/biomodels/network.h>
  11 +
  12 +#ifdef __CUDACC__
  13 +#include <cublas_v2.h>
  14 +#include <stim/cuda/cudatools/error.h>
  15 +#endif
  16 +
  17 +namespace stim {
  18 + template <typename A, typename B, typename C>
  19 + struct triple {
  20 + A first;
  21 + B second;
  22 + C third;
  23 + };
  24 +
  25 + template <typename T>
  26 + struct bridge {
  27 + std::vector<unsigned> v; // vertices' indices
  28 + std::vector<typename stim::vec3<T> > V; // vertices' coordinates
  29 + T l; // length
  30 + T r; // radii
  31 + T deltaP; // pressure drop
  32 + T Q; // volume flow rate
  33 + };
  34 +
  35 + template <typename T>
  36 + class flow {
  37 +
  38 + private:
  39 +
  40 + // calculate the cofactor of elemen[row][col]
  41 + void get_minor(T** src, T** dest, int row, int col, int order) {
  42 +
  43 + // index of element to be copied
  44 + int rowCount = 0;
  45 + int colCount = 0;
  46 +
  47 + for (int i = 0; i < order; i++) {
  48 + if (i != row) {
  49 + colCount = 0;
  50 + for (int j = 0; j < order; j++) {
  51 + // when j is not the element
  52 + if (j != col) {
  53 + dest[rowCount][colCount] = src[i][j];
  54 + colCount++;
  55 + }
  56 + }
  57 + rowCount++;
99 58 }
100   -*/
101   - // Decrement the counters
102   - i--;
103   - j--;
104   - count--;
105   - numPress--; // Decrement to the exact number of dangle nodes
106 59 }
107 60 }
108   - }
109 61  
110   - return numPress; // Number of dangle nodes
111   -}
  62 + // calculate the det()
  63 + T determinant(T** mat, int order) {
112 64  
  65 + // degenate case when n = 1
  66 + if (order == 1)
  67 + return mat[0][0];
113 68  
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   - }
  69 + T det = 0.0; // determinant value
130 70  
131   - output_file << "\n";
132   - }
  71 + // allocate the cofactor matrix
  72 + T** minor = (T**)malloc((order - 1) * sizeof(T*));
  73 + for (int i = 0; i < order - 1; i++)
  74 + minor[i] = (T*)malloc((order - 1) * sizeof(T));
133 75  
134   - output_file.close( );
135   -}
136 76  
  77 + for (int i = 0; i < order; i++) {
137 78  
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   - }
  79 + // get minor of element(0, i)
  80 + get_minor(mat, minor, 0, i, order);
177 81  
178   - // Store to radii.txt
179   - if(cnt == 3)
180   - radiiData << number << endl;
  82 + // recursion
  83 + det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1);
  84 + }
181 85  
182   - // Store to lengths.txt
183   - if(cnt == 4)
184   - {
185   - lenData << number << endl;
  86 + // release memory
  87 + for (int i = 0; i < order - 1; i++)
  88 + free(minor[i]);
  89 + free(minor);
186 90  
187   - numElements++; // Count the elements
188   - cnt = 0; // Reset counter
  91 + return det;
189 92 }
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 93  
219   - // Store destination node values to the array dest
220   - v = 0; // Reset counter
221   - ifstream readDest("destCol.txt");
  94 + public:
  95 + std::vector<typename stim::triple<unsigned, unsigned, T> > V; // velocity
  96 + std::vector<typename stim::triple<unsigned, unsigned, T> > Q; // volume flow rate
  97 + std::vector<typename stim::triple<unsigned, unsigned, T> > deltaP; // pressure drop
  98 + std::vector<T> pressure; // pressure
  99 +
  100 + flow() {} // default constructor
  101 +
  102 + // calculate the flow rate of 3D model(circle cross section)
  103 + void calculate_flow_rate(unsigned e, T r) {
  104 + stim::triple<unsigned, unsigned, T> tmp_Q;
  105 + tmp_Q.first = V[e].first; // copy the vertices information
  106 + tmp_Q.second = V[e].second;
  107 + tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s
  108 + Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  109 + }
222 110  
223   - while( readDest >> tmp )
224   - {
225   - dest[v] = (int)tmp;
226   - v++;
227   - }
  111 + // calculate the flow rate of 2D model(rectangular cross section)
  112 + void calculate_flow_rate(unsigned e, T r, T h) {
  113 + stim::triple<unsigned, unsigned, T> tmp_Q;
  114 + tmp_Q.first = V[e].first; // copy the vertices information
  115 + tmp_Q.second = V[e].second;
  116 + tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s
  117 + Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  118 + }
228 119  
229   - readDest.close( );
  120 + // calculate the pressure drop of 3D model(circle cross section)
  121 + void calculate_deltaP(unsigned e, T u, T l, T r) {
  122 + stim::triple<unsigned, unsigned, T> tmp_deltaP;
  123 + tmp_deltaP.first = V[e].first; // copy the vertices information
  124 + tmp_deltaP.second = V[e].second;
  125 + tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa
  126 + deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  127 + }
230 128  
231   - // Store radius values to the array radii
232   - v = 0; // Reset counter
233   - ifstream readRad("radii.txt");
  129 + // calculate the pressure drop of 2D model(rectangular cross section)
  130 + void calculate_deltaP(unsigned e, T u, T l, T r, T h) {
  131 + stim::triple<unsigned, unsigned, T> tmp_deltaP;
  132 + tmp_deltaP.first = V[e].first; // copy the vertices information
  133 + tmp_deltaP.second = V[e].second;
  134 + tmp_deltaP.third = (12 * u * l * Q[e].third) / (stim::PI * h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa
  135 + deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  136 + }
234 137  
235   - while( readRad >> tmp )
236   - {
237   - radii[v] = tmp;
238   - v++;
239   - }
  138 + // better way to do this???
  139 + // find the maximum and minimum pressure positions
  140 + void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) {
  141 + std::vector<T> P(n_v, FLT_MAX);
  142 + // set one to reference
  143 + P[Q[0].first] = 0.0;
  144 + unsigned first = 0;
  145 + unsigned second = 0;
  146 + // calculate all the relative pressure in brute force manner
  147 + for (unsigned e = 0; e < n_e; e++) {
  148 + // assuming the obj file stores in a straight order, in other words, like swc file
  149 + first = Q[e].first;
  150 + second = Q[e].second;
  151 + if (P[first] != FLT_MAX) // if pressure at start vertex is known
  152 + P[second] = P[first] - deltaP[e].third;
  153 + else if (P[second] != FLT_MAX) // if pressure at end vertex is known
  154 + P[first] = P[second] + deltaP[e].third;
  155 + }
240 156  
241   - readRad.close( );
  157 + // find the maximum and minimum pressure position
  158 + auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number
  159 + auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number
242 160  
243   - // Store length values to the array len
244   - v = 0; // Reset counter
245   - ifstream readLen("lengths.txt");
  161 + max = std::distance(P.begin(), m1);
  162 + min = std::distance(P.begin(), m2);
246 163  
247   - while( readLen >> tmp )
248   - {
249   - len[v] = tmp;
250   - v++;
251   - }
  164 + T tmp_m = *m2;
  165 + // Now set the lowest pressure port to reference pressure(0.0 Pa)
  166 + for (unsigned i = 0; i < n_v; i++)
  167 + P[i] -= tmp_m;
252 168  
253   - readLen.close( );
254   -}
  169 + for (unsigned i = 0; i < n_v; i++)
  170 + pressure.push_back(P[i]);
  171 + }
255 172  
  173 + void inversion(T** A, int order, T** C) {
  174 +
  175 +#ifdef __CUDACC__
  176 +
  177 + // convert from double pointer to single pointer, make it flat
  178 + T* Aflat = (T*)malloc(order * order * sizeof(T));
  179 + for (unsigned i = 0; i < order; i++)
  180 + for (unsigned j = 0; j < order; j++)
  181 + Aflat[i * order + j] = A[i][j];
  182 + T* Cflat = (T*)malloc(order * order * sizeof(T));
  183 +
  184 + // create device pointer
  185 + T* d_Aflat; // flat original matrix
  186 + T* d_Cflat; // flat inverse matrix
  187 + T** d_A; // put the flat original matrix into another array of pointer
  188 + T** d_C;
  189 + int *d_P;
  190 + int *d_INFO;
  191 +
  192 + // allocate memory on device
  193 + HANDLE_ERROR(cudaMalloc((void**)&d_Aflat, order * order * sizeof(T)));
  194 + HANDLE_ERROR(cudaMalloc((void**)&d_Cflat, order * order * sizeof(T)));
  195 + HANDLE_ERROR(cudaMalloc((void**)&d_A, sizeof(T*)));
  196 + HANDLE_ERROR(cudaMalloc((void**)&d_C, sizeof(T*)));
  197 + HANDLE_ERROR(cudaMalloc((void**)&d_P, order * 1 * sizeof(int)));
  198 + HANDLE_ERROR(cudaMalloc((void**)&d_INFO, 1 * sizeof(int)));
  199 +
  200 + // copy matrix from host to device
  201 + HANDLE_ERROR(cudaMemcpy(d_Aflat, Aflat, order * order * sizeof(T), cudaMemcpyHostToDevice));
  202 +
  203 + // copy matrix from device to device
  204 + HANDLE_ERROR(cudaMemcpy(d_A, &d_Aflat, sizeof(T*), cudaMemcpyHostToDevice));
  205 + HANDLE_ERROR(cudaMemcpy(d_C, &d_Cflat, sizeof(T*), cudaMemcpyHostToDevice));
256 206  
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   -}
  207 + // calculate the inverse of matrix based on cuBLAS
  208 + cublasHandle_t handle;
  209 + CUBLAS_HANDLE_ERROR(cublasCreate_v2(&handle)); // create cuBLAS handle object
287 210  
  211 + CUBLAS_HANDLE_ERROR(cublasSgetrfBatched(handle, order, d_A, order, d_P, d_INFO, 1));
288 212  
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++)
  213 + int INFO = 0;
  214 + HANDLE_ERROR(cudaMemcpy(&INFO, d_INFO, sizeof(int), cudaMemcpyDeviceToHost));
  215 + if (INFO == order)
348 216 {
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
  217 + std::cout << "Factorization Failed : Matrix is singular." << std::endl;
  218 + cudaDeviceReset();
  219 + exit(1);
  220 + }
  221 +
  222 + CUBLAS_HANDLE_ERROR(cublasSgetriBatched(handle, order, (const T **)d_A, order, d_P, d_C, order, d_INFO, 1));
  223 +
  224 + CUBLAS_HANDLE_ERROR(cublasDestroy_v2(handle));
  225 +
  226 + // copy inverse matrix from device to device
  227 + HANDLE_ERROR(cudaMemcpy(&d_Cflat, d_C, sizeof(T*), cudaMemcpyDeviceToHost));
  228 +
  229 + // copy inverse matrix from device to host
  230 + HANDLE_ERROR(cudaMemcpy(Cflat, d_Cflat, order * order * sizeof(T), cudaMemcpyDeviceToHost));
  231 +
  232 + for(unsigned i = 0; i < order; i++)
  233 + memcpy(C[i], &Cflat[i*order], order * sizeof(T*));
  234 +
  235 + // clear up
  236 + free(Aflat);
  237 + free(Cflat);
  238 + HANDLE_ERROR(cudaFree(d_Aflat));
  239 + HANDLE_ERROR(cudaFree(d_Cflat));
  240 + HANDLE_ERROR(cudaFree(d_A));
  241 + HANDLE_ERROR(cudaFree(d_C));
  242 + HANDLE_ERROR(cudaFree(d_P));
  243 + HANDLE_ERROR(cudaFree(d_INFO));
  244 +
  245 +#else
  246 + // get the determinant of a
  247 + double det = 1.0 / determinant(A, order);
  248 +
  249 + // memory allocation
  250 + T* tmp = (T*)malloc((order - 1)*(order - 1) * sizeof(T));
  251 + T** minor = (T**)malloc((order - 1) * sizeof(T*));
  252 + for (int i = 0; i < order - 1; i++)
  253 + minor[i] = tmp + (i * (order - 1));
  254 +
  255 + for (int j = 0; j < order; j++) {
  256 + for (int i = 0; i < order; i++) {
  257 + // get the co-factor (matrix) of A(j,i)
  258 + get_minor(A, minor, j, i, order);
  259 + C[i][j] = det * determinant(minor, order - 1);
  260 + if ((i + j) % 2 == 1)
  261 + C[i][j] = -C[i][j];
359 262 }
360 263 }
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
  264 +
  265 + // release memory
  266 + free(tmp);
  267 + free(minor);
  268 +#endif
372 269 }
373   - }
374   -
375   - return(det);
  270 + };
376 271 }
377 272  
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   -}
406 273 \ No newline at end of file
  274 +#endif
407 275 \ No newline at end of file
... ...
stim/biomodels/flow_dep.h 0 → 100644
  Show/Hide comments   View file @28b882c
  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
... ...