I have a little program that implements an ecological network analysis
algorithm. (It looks at the importance of indirect energy or matter
flows between species.)
01 //dea.cpp -- The fundamental DEA algorithm
02
03 #include
04 #include
05 #include //for test
06 #include "storage_adaptors.hpp" //for test
07
08 using boost::numeric::ublas::matrix;
09
10 void dea(const matrix<double>[], int, int, matrix<double>[]);
11 matrix<double> indflow(const matrix<double> & Ndea, const
matrix<double> & Gdea);
12
13 int main() {
14 using boost::numeric::ublas::make_matrix_from_pointer;
15 matrix<double> Garr[4];
16 matrix<double> Narr[2];
17 double G1[3][3] = {0, 0, 0, 0.33333, 0, 0, 0.66667, 1, 0};
18 double G2[3][3] = {0, 0, 0, 0.6, 0, 0, 0.4, 1, 0};
19 double G3[3][3] = {0, 0, 0, 0.55, 0, 0, 0.45, 1, 0};
20 double G4[3][3] = {0, 0, 0, 0.5, 0, 0, 0.5, 1, 0};
21 Garr[0] = make_matrix_from_pointer(G1);
22 Garr[1] = make_matrix_from_pointer(G2);
23 Garr[2] = make_matrix_from_pointer(G3);
24 Garr[3] = make_matrix_from_pointer(G4);
25 dea(Garr, 2, 4, Narr);
26 std::cout << Narr[0] << std::endl;
27 matrix<double> ratio = indflow(Narr[0], Garr[0]);
28 std::cout << ratio << std::endl;
29 return 0;
30 }
31
32 //The fundamental DEA algorithm
33 //Takes array of G matrices, window size, array length and empty
array for N_DEA matrices as input
34 void dea(const matrix<double> Garr[], int win, int arrLen,
matrix<double> Narr[]) {
35 using namespace boost::numeric::ublas;
36 int matSize = Garr[0].size1();
37 for (int i = 0; i < arrLen - win; i++) {
38 matrix<double> Gprod = identity_matrix<double>(matSize);
39 Narr[i] = identity_matrix<double>(matSize);
40 //the core loop
41 for (int j = 0; j < win; j++) {
42 Gprod = prod(Gprod, Garr[i+j]);
43 Narr[i] += Gprod;
44 }
45 }
46 }
47
48 //computes indirect flow fraction
49 matrix<double> indflow(const matrix<double> & Ndea, const
matrix<double> & Gdea) {
50 using namespace boost::numeric::ublas;
51 int matSize = Ndea.size1();
52 matrix<double> iflow = Ndea - Gdea -
identity_matrix<double>(matSize); //Ind = N-G-I
53 matrix<double> frac = zero_matrix<double>(matSize, matSize);
//inialize indirect flow fraction matrix
54 //perform division for nonzero, off-diagonal elements
55 for (int i=0; i
Check out my blog, http://perceivingwholes.blogspot.comPerceiving Wholes
"The whole person must have both the humility to nurture the
Earth and the pride to go to Mars." --Wyn Wachhorst, The Dream
of Spaceflight
