stim/biomodels/flow.h

#pragma once
#include <fstream>									// Required for ofstream, etc.
#include <iomanip>									// Required for setw
#include <iostream>									// Required for cout, cin, etc.
#include <tuple>									// Required for returning multiple values from a function
#include <vector>
//#include <string>
//#include <stdlib.h>								// Required to remove ambiguous error for cout, cin, etc.
using namespace std;
class flow
{
public:
	void backupToTxt(unsigned int nL, double **D, char filename[]);
	tuple<int, int> copySrcDesRadLen(char filename[]);
	void copyToArray(int *src, int *dest, double *radii, double *len);
	int getDangleNodes(int datarow, int numNodes, int *row, int *column, int *dangleNodes);
	void inversion(double **a, int n, double **b);
protected:
	float determinant(double **a, int n);
	int minor(double **src, double **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, double **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<int, int> 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, double *radii, double *len)
{
	int v = 0;
	double 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(double **a, int n, double **b)
{
	// Get 1 over the determinant of A
	double det = (double)(1.0/determinant(a, n));
	cerr << "\n1/det(C) = " << det << endl;				// DELETE!!!
    // Memory allocation
    double *tmp = new double[(n-1) * (n-1)];
    double **m = new double * [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);
            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(double **a, int n)
{
	int i, j, j1, j2;								// General loop and matrix subscripts
    double det = 0;									// Initialize determinant
    double **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 = (double **) malloc((n-1) * sizeof(double *));
            for (i = 0; i < n-1; i++)
				m[i] = (double *) malloc((n-1)* sizeof(double));
                       //     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 += (double)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(double **src, double **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;
}