Zhenghua Gong / Structures and MB

Download as
Browse Code ยป

Commit 7f982ea62aaac9de6a72455eb7e43ad484908594

Authored by Zhenghua Gong
1 parent c643e466ce
Exists in master

Updated MyOrder.cpp. So it can calculate factorial of 100.

Showing 1 changed file with 115 additions and 86 deletions   Show diff stats
1 #include <string> 1 #include <string>
2 #include <vector> 2 #include <vector>
3 #include <numeric> 3 #include <numeric>
4 #include <fstream> 4 #include <fstream>
5 #include <iterator> 5 #include <iterator>
6 #include <iostream> 6 #include <iostream>
7 #include <utility> 7 #include <utility>
8 #include <iomanip> 8 #include <iomanip>
9 #include <ctime> 9 #include <ctime>
10 #include <boost/math/special_functions.hpp> 10 #include <boost/math/special_functions.hpp>
11 #include <boost/lexical_cast.hpp> 11 #include <boost/lexical_cast.hpp>
12 #include <boost/algorithm/string/replace.hpp> 12 #include <boost/algorithm/string/replace.hpp>
13 #include <boost/multiprecision/cpp_int.hpp>
13 14
14 #include <stdio.h> 15 #include <stdio.h>
15 #include <math.h> 16 #include <math.h>
16 #include <algorithm> // std::find 17 #include <algorithm> // std::find
17 #include <map> 18 #include <map>
18 #include <iterator> 19 #include <iterator>
19 using namespace std; 20 #include <bits/stdc++.h>
21 #include <tuple>
20 22
23 using namespace std;
21 24
25 //Function for calculating the binomial coefficient
26 boost::multiprecision::cpp_int BinomialCoefficient(unsigned int n, unsigned int k) {
27 if (k == 0) { return 1; }
28 else { return (n * BinomialCoefficient(n - 1, k - 1)) / k; }
29 }
22 struct compare { 30 struct compare {
23 bool operator()(const std::string& first, const std::string& second) { 31 bool operator()(const std::string& first, const std::string& second) {
24 if(first.size() == second.size()) 32 if(first.size() == second.size())
25 return first < second; 33 return first < second;
26 else 34 else
27 return first.size() < second.size(); 35 return first.size() < second.size();
28 } 36 }
29 }; 37 };
30 38
31 string int_to_str(int num) 39 string int_to_str(int num)
32 { 40 {
33 stringstream ss; 41 stringstream ss;
34 42
35 ss << num; 43 ss << num;
36 44
37 return ss.str(); 45 return ss.str();
38 }; 46 };
39 47
40 48
41 int str_to_int(string st) 49 int str_to_int(string st)
42 { 50 {
43 int result; 51 int result;
44 52
45 stringstream(st) >> result; 53 stringstream(st) >> result;
46 54
47 return result; 55 return result;
48 }; 56 };
49 57
50 bool compareI(const pair<string, double>&i, const pair<string, double>&j) 58 bool compareI(const pair<string, double>&i, const pair<string, double>&j)
51 { 59 {
52 return i.second < j.second; 60 return i.second < j.second;
53 } 61 }
54 62
55 bool compareD(const pair<string, double>&i, const pair<string, double>&j) 63 bool compareD(const pair<string, double>&i, const pair<string, double>&j)
56 { 64 {
57 return i.second > j.second; 65 return i.second > j.second;
58 } 66 }
59 67
60 bool compareIn(const pair<double, string>&i, const pair<double, string>&j) 68 bool compareIn(const pair<double, string>&i, const pair<double, string>&j)
61 { 69 {
62 return i.first < j.first; 70 return i.first < j.first;
63 } 71 }
64 72
65 bool compareDe(const pair<double, string>&i, const pair<double, string>&j) 73 bool compareDe(const pair<double, string>&i, const pair<double, string>&j)
66 { 74 {
67 return i.first > j.first; 75 return i.first > j.first;
68 } 76 }
69 77
70 double logAB (double x, double y) 78 double logAB (double x, double y)
71 { 79 {
72 double result; 80 double result;
73 double maxVal = max(x,y); 81 double maxVal = max(x,y);
74 82
75 if(maxVal == x) 83 if(maxVal == x)
76 { 84 {
77 result = maxVal + log(1+exp(y-maxVal)); 85 result = maxVal + log(1+exp(y-maxVal));
78 } 86 }
79 else 87 else
80 { 88 {
81 result = maxVal + log(1+exp(x-maxVal)); 89 result = maxVal + log(1+exp(x-maxVal));
82 } 90 }
83 return result; 91 return result;
84 } 92 }
85 93
86 double persentageXofY (double newS, double oldS) 94 double persentageXofY (double newS, double oldS)
87 { 95 {
88 double result; 96 double result;
89 result = exp(oldS-newS)*100; 97 result = exp(oldS-newS)*100;
90 return result; 98 return result;
91 } 99 }
92 100
93 void findAndReplaceAll(std::string & data, std::string toSearch, std::string replaceStr) 101 void findAndReplaceAll(std::string & data, std::string toSearch, std::string replaceStr)
94 { 102 {
95 // Get the first occurrence 103 // Get the first occurrence
96 size_t pos = data.find(toSearch); 104 size_t pos = data.find(toSearch);
97 105
98 // Repeat till end is reached 106 // Repeat till end is reached
99 while( pos != std::string::npos) 107 while( pos != std::string::npos)
100 { 108 {
101 // Replace this occurrence of Sub String 109 // Replace this occurrence of Sub String
102 data.replace(pos, toSearch.size(), replaceStr); 110 data.replace(pos, toSearch.size(), replaceStr);
103 // Get the next occurrence from the current position 111 // Get the next occurrence from the current position
104 pos =data.find(toSearch, pos + toSearch.size()); 112 pos =data.find(toSearch, pos + toSearch.size());
105 } 113 }
106 } 114 }
107 115
116 string findMB(string structures, string variable){
117 string MB;
118
119
120 return MB;
121 }
122
108 int main() { 123 int main() {
109 /* 124 /*
110 //TEMPORARY INPUT FILE 125 //TEMPORARY INPUT FILE
111 vector< vector<int> > DAT; 126 vector< vector<int> > DAT;
112 vector <int> temp, temp2, temp3, temp4, temp5, temp6, temp7, temp8, temp9, temp10, temp11, temp12, temp13, temp14, temp15, temp16; 127 vector <int> temp, temp2, temp3, temp4, temp5, temp6, temp7, temp8, temp9, temp10, temp11, temp12, temp13, temp14, temp15, temp16;
113 //Filling up some of the test variables needed for the data set 128 //Filling up some of the test variables needed for the data set
114 temp.push_back(0); 129 temp.push_back(0);
115 temp.push_back(0); 130 temp.push_back(0);
116 temp.push_back(1); 131 temp.push_back(1);
117 //set 2 132 //set 2
118 temp2.push_back(0); 133 temp2.push_back(0);
119 temp2.push_back(1); 134 temp2.push_back(1);
120 temp2.push_back(0); 135 temp2.push_back(0);
121 //set3 136 //set3
122 temp3.push_back(0); 137 temp3.push_back(0);
123 temp3.push_back(1); 138 temp3.push_back(1);
124 temp3.push_back(1); 139 temp3.push_back(1);
125 //set4 140 //set4
126 temp4.push_back(1); 141 temp4.push_back(1);
127 temp4.push_back(0); 142 temp4.push_back(0);
128 temp4.push_back(0); 143 temp4.push_back(0);
129 //set5 144 //set5
130 temp5.push_back(1); 145 temp5.push_back(1);
131 temp5.push_back(1); 146 temp5.push_back(1);
132 temp5.push_back(1); 147 temp5.push_back(1);
133 //set6 148 //set6
134 temp6.push_back(1); 149 temp6.push_back(1);
135 temp6.push_back(1); 150 temp6.push_back(1);
136 temp6.push_back(0); 151 temp6.push_back(0);
137 //set7 152 //set7
138 temp7.push_back(1); 153 temp7.push_back(1);
139 temp7.push_back(0); 154 temp7.push_back(0);
140 temp7.push_back(1); 155 temp7.push_back(1);
141 //set8 156 //set8
142 temp8.push_back(1); 157 temp8.push_back(1);
143 temp8.push_back(0); 158 temp8.push_back(0);
144 temp8.push_back(1); 159 temp8.push_back(1);
145 //set9 160 //set9
146 temp9.push_back(1); 161 temp9.push_back(1);
147 temp9.push_back(0); 162 temp9.push_back(0);
148 temp9.push_back(1); 163 temp9.push_back(1);
149 //set10 164 //set10
150 temp10.push_back(1); 165 temp10.push_back(1);
151 temp10.push_back(0); 166 temp10.push_back(0);
152 temp10.push_back(1); 167 temp10.push_back(1);
153 //set11 168 //set11
154 temp11.push_back(1); 169 temp11.push_back(1);
155 temp11.push_back(1); 170 temp11.push_back(1);
156 temp11.push_back(1); 171 temp11.push_back(1);
157 //set12 172 //set12
158 temp12.push_back(1); 173 temp12.push_back(1);
159 temp12.push_back(0); 174 temp12.push_back(0);
160 temp12.push_back(1); 175 temp12.push_back(1);
161 //set13 176 //set13
162 temp13.push_back(1); 177 temp13.push_back(1);
163 temp13.push_back(0); 178 temp13.push_back(0);
164 temp13.push_back(1); 179 temp13.push_back(1);
165 //set14 180 //set14
166 temp14.push_back(1); 181 temp14.push_back(1);
167 temp14.push_back(1); 182 temp14.push_back(1);
168 temp14.push_back(0); 183 temp14.push_back(0);
169 //set15 184 //set15
170 temp15.push_back(0); 185 temp15.push_back(0);
171 temp15.push_back(1); 186 temp15.push_back(1);
172 temp15.push_back(1); 187 temp15.push_back(1);
173 //set16 188 //set16
174 temp16.push_back(1); 189 temp16.push_back(1);
175 temp16.push_back(0); 190 temp16.push_back(0);
176 temp16.push_back(1); 191 temp16.push_back(1);
177 std::cout << endl; 192 std::cout << endl;
178 //Filling up test DATASET 193 //Filling up test DATASET
179 DAT.push_back(temp); 194 DAT.push_back(temp);
180 DAT.push_back(temp2); 195 DAT.push_back(temp2);
181 DAT.push_back(temp3); 196 DAT.push_back(temp3);
182 DAT.push_back(temp4); 197 DAT.push_back(temp4);
183 DAT.push_back(temp5); 198 DAT.push_back(temp5);
184 DAT.push_back(temp6); 199 DAT.push_back(temp6);
185 DAT.push_back(temp7); 200 DAT.push_back(temp7);
186 DAT.push_back(temp8); 201 DAT.push_back(temp8);
187 DAT.push_back(temp9); 202 DAT.push_back(temp9);
188 DAT.push_back(temp10); 203 DAT.push_back(temp10);
189 DAT.push_back(temp11); 204 DAT.push_back(temp11);
190 DAT.push_back(temp12); 205 DAT.push_back(temp12);
191 DAT.push_back(temp13); 206 DAT.push_back(temp13);
192 DAT.push_back(temp14); 207 DAT.push_back(temp14);
193 DAT.push_back(temp15); 208 DAT.push_back(temp15);
194 DAT.push_back(temp16); 209 DAT.push_back(temp16);
195 std::cout << endl; 210 std::cout << endl;
196 size_t totvars = DAT[1].size(); 211 size_t totvars = DAT[1].size();
197 size_t tottuples = DAT.size(); 212 size_t tottuples = DAT.size();
198 */ 213 */
199 //Time before user input 214 //Time before user input
200 time_t now = time(0); 215 time_t now = time(0);
201 char* dt = ctime(&now); 216 char* dt = ctime(&now);
202 std::cout << "The local date and time is: " << dt << std::endl; 217 std::cout << "The local date and time is: " << dt << std::endl;
203 string FILENAME; 218 string FILENAME;
204 std::cout << "What is the name of your file?: "; 219 std::cout << "What is the name of your file?: ";
205 std::cin >> FILENAME; 220 std::cin >> FILENAME;
206 // Variable declarations 221 // Variable declarations
207 fstream file; 222 fstream file;
208 int COLS; 223 int COLS;
209 std::cout << "How many variables in your data file?: "; 224 std::cout << "How many variables in your data file?: ";
210 std::cin >> COLS; 225 std::cin >> COLS;
211 std::cout << endl; 226 std::cout << endl;
212 vector < vector <int> > DAT; // 2d array as a vector of vectors 227 vector < vector <int> > DAT; // 2d array as a vector of vectors
213 vector <int> rowVector(COLS); // vector to add into 'array' (represents a row) 228 vector <int> rowVector(COLS); // vector to add into 'array' (represents a row)
214 int row = 0; // Row counter 229 int row = 0; // Row counter
215 230
216 // Read file 231 // Read file
217 file.open(FILENAME.c_str(), ios::in); // Open file 232 file.open(FILENAME.c_str(), ios::in); // Open file
218 if (file.is_open()) { // If file has correctly opened... 233 if (file.is_open()) { // If file has correctly opened...
219 // Output debug message 234 // Output debug message
220 cout << "File correctly opened" << endl; 235 cout << "File correctly opened" << endl;
221 236
222 // Dynamically store data into array 237 // Dynamically store data into array
223 while (file.good()) { // ... and while there are no errors, 238 while (file.good()) { // ... and while there are no errors,
224 DAT.push_back(rowVector); // add a new row, 239 DAT.push_back(rowVector); // add a new row,
225 for (int col = 0; col<COLS; col++) { 240 for (int col = 0; col<COLS; col++) {
226 file >> DAT[row][col]; // fill the row with col elements 241 file >> DAT[row][col]; // fill the row with col elements
227 } 242 }
228 row++; // Keep track of actual row 243 row++; // Keep track of actual row
229 } 244 }
230 } 245 }
231 else cout << "Unable to open file" << endl; 246 else cout << "Unable to open file" << endl;
232 file.close(); 247 file.close();
248 DAT.erase(DAT.end());
233 249
234 size_t totvars = DAT[1].size(); //column number 250 size_t totvars = DAT[1].size(); //column number
235 size_t tottuples = DAT.size();//row number 251 size_t tottuples = DAT.size();//row number
236 252
237 253
238 254
239 255
240 //Ask User for the order they would like to test Order starts at 0 e.g. 0,1,2 is valid order for three variables 256 //Ask User for the order they would like to test Order starts at 0 e.g. 0,1,2 is valid order for three variables
241 int order_i; 257 int order_i;
242 vector <int> order; 258 vector <int> order;
243 std::cout << "Which is the order you would like to test?: "; 259 std::cout << "Which is the order you would like to test?: ";
244 for (int i = 0; i < totvars;++i) { 260 for (int i = 0; i < totvars;++i) {
245 std::cin >> order_i; 261 std::cin >> order_i;
246 order.push_back(order_i); 262 order.push_back(order_i);
247 } 263 }
248 std::cout << endl; 264 std::cout << endl;
249 //Ask user for max value and min value of each variable being inputed 265 //Ask user for max value and min value of each variable being inputed
250 vector <int> max_cat; 266 vector <int> max_cat;
251 vector <int> min_cat; 267 vector <int> min_cat;
252 int max_cat_input; 268 int max_cat_input;
253 int min_cat_input; 269 int min_cat_input;
254 270
255 std::cout << "What is the maximum categorical value of each variable?: "; 271 std::cout << "What is the maximum categorical value of each variable?: ";
256 for (int i = 0; i < totvars; ++i) { 272 for (int i = 0; i < totvars; ++i) {
257 std::cin >> max_cat_input; 273 std::cin >> max_cat_input;
258 max_cat.push_back(max_cat_input); 274 max_cat.push_back(max_cat_input);
259 } 275 }
260 for (size_t i = 0; i < max_cat.size();++i) { 276 for (size_t i = 0; i < max_cat.size();++i) {
261 std::cout << max_cat[i] << " "; 277 std::cout << max_cat[i] << " ";
262 } 278 }
263 std::cout << endl; 279 std::cout << endl;
264 std::cout << endl; 280 std::cout << endl;
265 std::cout << "What is the minimum categorical value of each variable?: "; 281 std::cout << "What is the minimum categorical value of each variable?: ";
266 for (int i = 0; i < totvars; ++i) { 282 for (int i = 0; i < totvars; ++i) {
267 std::cin >> min_cat_input; 283 std::cin >> min_cat_input;
268 min_cat.push_back(min_cat_input); 284 min_cat.push_back(min_cat_input);
269 } 285 }
270 for (size_t i = 0; i < min_cat.size();++i) { 286 for (size_t i = 0; i < min_cat.size();++i) {
271 std::cout << min_cat[i] << " "; 287 std::cout << min_cat[i] << " ";
272 } 288 }
273 std::cout << endl; 289 std::cout << endl;
274 //Set the maximum for the amount of parents for any given variable 290 //Set the maximum for the amount of parents for any given variable
275 unsigned int maxparents; 291 unsigned int maxparents;
276 std::cout << "What is the maximum amount of parents to be considered for any given variable?: "; 292 std::cout << "What is the maximum amount of parents to be considered for any given variable?: ";
277 std::cin >> maxparents; 293 std::cin >> maxparents;
278 std::cout << endl; 294 std::cout << endl;
279 295
280 double PERCENT; 296 double PERCENT;
281 std::cout << "What is the percentage (Please use the format of 98.99, not 0.9899)?: "; 297 std::cout << "What is the percentage (Please use the format of 98.99, not 0.9899)?: ";
282 std::cin >> PERCENT; 298 std::cin >> PERCENT;
283 std::cout << endl; 299 std::cout << endl;
284 300
301 string filepath;
302 std::cout << "What is the output file name?: ";
303 std::cin >> filepath;
304 std::cout << endl;
305
285 //Time that program begins 306 //Time that program begins
286 time_t start = time(0); 307 time_t start = time(0);
287 char* dt_start = ctime(&start); 308 char* dt_start = ctime(&start);
288 std::cout << "The local date and time is: " << dt_start << std::endl; 309 std::cout << "The local date and time is: " << dt_start << std::endl;
289 310
290 311
291 //Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category: 312 //Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category:
292 //categories in i 313 //categories in i
293 vector <int> catsi; 314 vector <int> catsi;
294 //total combinations of variables 315 //total combinations of variables
295 //int totcombos = 1; 316 //int totcombos = 1;
296 //how many catagori for every variable 317 //how many catagori for every variable
297 for ( int i = 0; i < totvars; ++i) { 318 for ( int i = 0; i < totvars; ++i) {
298 catsi.push_back((max_cat[i] - min_cat[i]) + 1); 319 catsi.push_back((max_cat[i] - min_cat[i]) + 1);
299 //totcombos = totcombos * ((max_cat[i] - min_cat[i]) + 1); 320 //totcombos = totcombos * ((max_cat[i] - min_cat[i]) + 1);
300 } 321 }
301 /*size_t mincatsi = 10000, mincatvar; 322 /*size_t mincatsi = 10000, mincatvar;
302 for (size_t i = 0; i < catsi.size(); ++i) { 323 for (size_t i = 0; i < catsi.size(); ++i) {
303 if (catsi[i] < mincatsi) { 324 if (catsi[i] < mincatsi) {
304 mincatsi = catsi[i]; 325 mincatsi = catsi[i];
305 mincatvar = i; 326 mincatvar = i;
306 } 327 }
307 }*/ 328 }*/
308 //print out catsi 329 //print out catsi
309 for (size_t i = 0; i < catsi.size();++i) { 330 for (size_t i = 0; i < catsi.size();++i) {
310 std::cout << catsi[i] << " "; 331 std::cout << catsi[i] << " ";
311 } 332 }
312 std::cout << endl; 333 std::cout << endl;
313 334
314 //Total Families Ui,alpha for a particular variable in the order 335 //Total Families Ui,alpha for a particular variable in the order
315 vector <unsigned long long> families;//is vector, first element is the number of parentset for the first variabel...... 336 vector <boost::multiprecision::cpp_int> families;
316 for (unsigned int i = 0; i < totvars; ++i) { 337 for (unsigned int i = 0; i < totvars; ++i) {
317 int numparents = i; 338 int numparents = i;
318 if (numparents == 0) { 339 if (numparents == 0) {
319 families.push_back(1); 340 families.push_back(1);
320 } 341 }
321 else { 342 else {
322 343
323 unsigned long long numfams = 0; 344 boost::multiprecision::cpp_int numfams = 0;
324 for (unsigned int j = 0; j <= i; ++j) { 345 for (unsigned int j = 0; j <= i; ++j) {
325 if (j <= maxparents) { 346 if (j <= maxparents) {
326 unsigned long long jFactorial = 1; 347 unsigned long long jFactorial = 1;
327 unsigned long long ijFactorial = 1; 348 unsigned long long ijFactorial = 1;
328 unsigned long long iFactorial = 1; 349 unsigned long long iFactorial = 1;
329 //Calculate j! 350 //Calculate j!
330 for (unsigned int g = 0; g <= j; ++g) { 351 for (unsigned int g = 0; g <= j; ++g) {
331 if (g != 0) { 352 if (g != 0) {
332 jFactorial *= g; 353 jFactorial *= g;
333 } 354 }
334 } 355 }
335 //Calculate i! 356 //Calculate i!
336 for (unsigned int g = 0; g <= i; ++g) { 357 for (unsigned int g = 0; g <= i; ++g) {
337 if (g != 0) { 358 if (g != 0) {
338 iFactorial *= g; 359 iFactorial *= g;
339 } 360 }
340 } 361 }
341 //Calculate (i-j)! 362 //Calculate (i-j)!
342 for (unsigned int g = 0; g <= (i - j); ++g) { 363 for (unsigned int g = 0; g <= (i - j); ++g) {
343 if (g != 0) { 364 if (g != 0) {
344 ijFactorial *= g; 365 ijFactorial *= g;
345 } 366 }
346 } 367 }
347 numfams += (iFactorial) / (jFactorial * ijFactorial); 368 numfams += BinomialCoefficient(i,j);
348 } 369 }
349 else { 370 else {
350 break; 371 break;
351 } 372 }
352 } 373 }
353 families.push_back(numfams); 374 families.push_back(numfams);
354 } 375 }
355 } 376 }
356 377
357 378
358 //How many parent combinations for each step? As well as there counts 379 //How many parent combinations for each step? As well as there counts
359 vector< vector <int> > ParentCombos; 380 vector< vector <int> > ParentCombos;
360 vector< vector <int> > fullNijkvector; 381 vector< vector <int> > fullNijkvector;
361 vector< vector <string> > indexofvar; // This is label or index for each variable. 382 vector< vector <string> > indexofvar; // This is label or index for each variable.
362 for (size_t i = 0; i < order.size(); ++i) { 383 for (size_t i = 0; i < order.size(); ++i) {
363 //i represents the order of the variable 384 //i represents the order of the variable
364 if (i == 0) { 385 if (i == 0) {
365 vector <int> tmp,Nijkovercombos1; 386 vector <int> tmp,Nijkovercombos1;
366 vector <string> tempstring; 387 vector <string> tempstring;
367 tempstring.push_back("[0]"); 388 // tempstring.push_back("[0]");
389 tempstring.push_back("[" + int_to_str(order[0])+ "]");
368 tmp.push_back(1); 390 tmp.push_back(1);
369 ParentCombos.push_back(tmp); 391 ParentCombos.push_back(tmp);
370 //counting the amount of times that a value of the first variable in the order occurs 392 //counting the amount of times that a value of the first variable in the order occurs
371 //this starts with the maximum value for that variable 393 //this starts with the maximum value for that variable
372 for (int hello = max_cat[order[0]];hello >= min_cat[order[0]];--hello) { 394 for (int hello = max_cat[order[0]];hello >= min_cat[order[0]];--hello) {
373 //hello cycles through the categories of the first variable in the order 395 //hello cycles through the categories of the first variable in the order
374 int Nijk1 = 0; 396 int Nijk1 = 0;
375 //green cycles through tuples 397 //green cycles through tuples
376 for (int green = 0; green < tottuples; ++green){ 398 for (int green = 0; green < tottuples; ++green){
377 if (DAT[green][order[0]] == hello) { 399 if (DAT[green][order[0]] == hello) {
378 Nijk1 += 1; 400 Nijk1 += 1;
379 } 401 }
380 } 402 }
381 Nijkovercombos1.push_back(Nijk1); 403 Nijkovercombos1.push_back(Nijk1);
382 } 404 }
383 fullNijkvector.push_back(Nijkovercombos1); 405 fullNijkvector.push_back(Nijkovercombos1);
384 indexofvar.push_back(tempstring); 406 indexofvar.push_back(tempstring);
385 407
386 } 408 }
387 else { 409 else {
388 vector <int> tmp,Nijkovercombos1; 410 vector <int> tmp,Nijkovercombos1;
389 vector <string> tempstring; 411 vector <string> tempstring;
390 string tempa = "[" + int_to_str(i) +"]"; 412 // string tempa = "[" + int_to_str(i) +"]";
413 string tempa = "[" + int_to_str(order[i]) +"]";
391 tempstring.push_back(tempa); 414 tempstring.push_back(tempa);
392 tmp.push_back(1); 415 tmp.push_back(1);
393 int numparnts = i; 416 int numparnts = i;
394 //counting the amount of times that a value of the last variable in the current order size occurs 417 //counting the amount of times that a value of the last variable in the current order size occurs
395 //this starts with the maximum value for that variable 418 //this starts with the maximum value for that variable
396 for (int hello = max_cat[order[numparnts]];hello >= min_cat[order[numparnts]];--hello) { 419 for (int hello = max_cat[order[numparnts]];hello >= min_cat[order[numparnts]];--hello) {
397 //hello cycles through the categories of the first variable in the order 420 //hello cycles through the categories of the first variable in the order
398 int Nijk1 = 0; 421 int Nijk1 = 0;
399 //green cycles through tuples 422 //green cycles through tuples
400 for (int green = 0; green < tottuples; ++green) { 423 for (int green = 0; green < tottuples; ++green) {
401 if (DAT[green][order[numparnts]] == hello) { 424 if (DAT[green][order[numparnts]] == hello) {
402 Nijk1 += 1; 425 Nijk1 += 1;
403 } 426 }
404 } 427 }
405 Nijkovercombos1.push_back(Nijk1); 428 Nijkovercombos1.push_back(Nijk1);
406 } 429 }
407 fullNijkvector.push_back(Nijkovercombos1); 430 fullNijkvector.push_back(Nijkovercombos1);
408 431
409 //j representing the number of parents 432 //j representing the number of parents
410 for (int it = 1; it <= numparnts; ++it) { 433 for (int it = 1; it <= numparnts; ++it) {
411 //(333)Creating a vector that uses the right combination 434 //(333)Creating a vector that uses the right combination
412 double Nloopy = 0, NcolFactorial = 1, iFactorial = 1, NiFactorial = 1; 435 boost::multiprecision::cpp_int Nloopy = 0;
436 //, NcolFactorial = 1, iFactorial = 1, NiFactorial = 1;
413 /*std::cout << "This is for " << numparnts << " choose " << it << endl; 437 /*std::cout << "This is for " << numparnts << " choose " << it << endl;
414 std::cout << "The iteration number is: " << it << endl;*/ 438 std::cout << "The iteration number is: " << it << endl;*/
415 //Accounting for the limit of parent quantity 439 //Accounting for the limit of parent quantity
416 if (it > maxparents) { 440 if (it > maxparents) {
417 break; 441 break;
418 } 442 }
419 else { 443 else {
420 vector <int> NewMat(numparnts, 0); 444 vector <int> NewMat(numparnts, 0);
421 for (int p = 0; p < it; ++p) { 445 for (int p = 0; p < it; ++p) {
422 NewMat[p] = 1; 446 NewMat[p] = 1;
423 } 447 }
424 for (int g = 2; g <= numparnts; ++g) { 448 /*for (int g = 2; g <= numparnts; ++g) {
425 NcolFactorial *= g; 449 NcolFactorial *= g;
426 } 450 }
427 for (int g = 2; g <= it; ++g) { 451 for (int g = 2; g <= it; ++g) {
428 iFactorial *= g; 452 iFactorial *= g;
429 } 453 }
430 for (int g = 2; g <= (numparnts - it); ++g) { 454 for (int g = 2; g <= (numparnts - it); ++g) {
431 NiFactorial *= g; 455 NiFactorial *= g;
432 } 456 }*/
433 //Nloopy represents the result of numparnts choose i e.g. numparnts choose 1 equals numparnts 457 //Nloopy represents the result of numparnts choose i e.g. numparnts choose 1 equals numparnts
434 Nloopy = NcolFactorial / (iFactorial * NiFactorial); 458 Nloopy = BinomialCoefficient(numparnts,it);
435 for (int iNloopy = 0; iNloopy < Nloopy; ++iNloopy) { 459 for (int iNloopy = 0; iNloopy < Nloopy; ++iNloopy) {
436 int combsparents = 1; 460 int combsparents = 1;
437 vector <int> parsetv; 461 vector <int> parsetv;
438 for (int par = 0; par < NewMat.size(); ++par){ 462 for (int par = 0; par < NewMat.size(); ++par){
439 if (NewMat[par] == 1){ 463 if (NewMat[par] == 1){
440 parsetv.push_back(par); 464 parsetv.push_back(par);
441 } 465 }
442 } 466 }
443 467
444 string tempstring2; 468 string tempstring2;
445 469
446 for (int par2 = 0; par2 < parsetv.size(); ++par2){ 470 for (int par2 = 0; par2 < parsetv.size(); ++par2){
447 if(par2+1==parsetv.size()){ 471 if(par2+1==parsetv.size()){
448 tempstring2 = tempstring2 + int_to_str(parsetv[par2]); 472 // tempstring2 = tempstring2 + int_to_str(parsetv[par2]);
473 tempstring2 = tempstring2 + int_to_str(order[parsetv[par2]]);
449 //tempstring2 = tempstring2 +","+"|" + int_to_str(i); 474 //tempstring2 = tempstring2 +","+"|" + int_to_str(i);
450 tempstring2 = "[" + int_to_str(i)+"|"+tempstring2 +"]"; 475 // tempstring2 = "[" + int_to_str(i)+"|"+tempstring2 +"]";
476 tempstring2 = "[" + int_to_str(order[i])+"|"+tempstring2 +"]";
451 } 477 }
452 else{ 478 else{
453 //tempstring2 = tempstring2 + int_to_str(parsetv[par2])+","; 479 //tempstring2 = tempstring2 + int_to_str(parsetv[par2])+",";
454 tempstring2 = tempstring2 + int_to_str(parsetv[par2])+":"; 480 // tempstring2 = tempstring2 + int_to_str(parsetv[par2])+":";
481 tempstring2 = tempstring2 + int_to_str(order[parsetv[par2]])+":";
455 } 482 }
456 483
457 484
458 } 485 }
459 tempstring.push_back(tempstring2); 486 tempstring.push_back(tempstring2);
460 487
461 /*std::cout << "This is iNloopy: " << (iNloopy + 1) << endl; 488 /*std::cout << "This is iNloopy: " << (iNloopy + 1) << endl;
462 std::cout << "Here comes the NewMat:" << endl;*/ 489 std::cout << "Here comes the NewMat:" << endl;*/
463 //(444)This sets up the process for changing 490 //(444)This sets up the process for changing
464 //PosOne tells me the position of the last one in the vector 491 //PosOne tells me the position of the last one in the vector
465 //We want to change when the position is the last position available in the vector 492 //We want to change when the position is the last position available in the vector
466 int SumOnes = 0, PosOne = 0, SumOnes2 = 0, PosOne2, NxtOne = 0, FrstOne = 0; 493 int SumOnes = 0, PosOne = 0, SumOnes2 = 0, PosOne2, NxtOne = 0, FrstOne = 0;
467 int SumOnes3 = 0, SumOnes4 = 0, SumY = 0; 494 int SumOnes3 = 0, SumOnes4 = 0, SumY = 0;
468 for (PosOne = (numparnts - 1); PosOne >= 0; --PosOne) { 495 for (PosOne = (numparnts - 1); PosOne >= 0; --PosOne) {
469 if (NewMat[PosOne] == 1) { 496 if (NewMat[PosOne] == 1) {
470 break; 497 break;
471 } 498 }
472 } 499 }
473 for (int y = (numparnts - 1); y >= (numparnts - it); --y) { 500 for (int y = (numparnts - 1); y >= (numparnts - it); --y) {
474 //SumOnes tells you the amount of ones in the last i columns 501 //SumOnes tells you the amount of ones in the last i columns
475 //These are the last columns being considered 502 //These are the last columns being considered
476 SumOnes += NewMat[y]; 503 SumOnes += NewMat[y];
477 } 504 }
478 for (PosOne2 = (numparnts - 1); PosOne2 >= 0; --PosOne2) { 505 for (PosOne2 = (numparnts - 1); PosOne2 >= 0; --PosOne2) {
479 //SumOnes2 tells you the amount of ones before you reach the next zero 506 //SumOnes2 tells you the amount of ones before you reach the next zero
480 //PosOne2 keeps track of the position of the coming zero 507 //PosOne2 keeps track of the position of the coming zero
481 SumOnes2 += NewMat[PosOne2]; 508 SumOnes2 += NewMat[PosOne2];
482 if ((SumOnes2 > 0) & (NewMat[PosOne2] == 0)) { 509 if ((SumOnes2 > 0) & (NewMat[PosOne2] == 0)) {
483 break; 510 break;
484 } 511 }
485 } 512 }
486 for (FrstOne = 0; FrstOne < numparnts; ++FrstOne) { 513 for (FrstOne = 0; FrstOne < numparnts; ++FrstOne) {
487 //FrstOne tells you the position of the first number 1 starting from the left hand side 514 //FrstOne tells you the position of the first number 1 starting from the left hand side
488 if (NewMat[FrstOne] == 1) { 515 if (NewMat[FrstOne] == 1) {
489 break; 516 break;
490 } 517 }
491 } 518 }
492 for (int x = (numparnts - 1); x >= (numparnts - it + 1); --x) { 519 for (int x = (numparnts - 1); x >= (numparnts - it + 1); --x) {
493 //SumOnes4 helps keep track of the sum of all ones located in the last i - 1 positions 520 //SumOnes4 helps keep track of the sum of all ones located in the last i - 1 positions
494 SumOnes4 += NewMat[x]; 521 SumOnes4 += NewMat[x];
495 } 522 }
496 //Prints out NewMat 523 //Prints out NewMat
497 /*for (int u = 0; u < numparnts; ++u) { 524 /*for (int u = 0; u < numparnts; ++u) {
498 std::cout << NewMat[u] << " "; 525 std::cout << NewMat[u] << " ";
499 } 526 }
500 std::cout << endl;*/ 527 std::cout << endl;*/
501 528
502 529
503 //Adding in the code that will allow counts parent combinations for this particular variable 530 //Adding in the code that will allow counts parent combinations for this particular variable
504 for (int q = 0; q < i; ++q) { 531 for (int q = 0; q < i; ++q) {
505 if (NewMat[q] == 1) { 532 if (NewMat[q] == 1) {
506 combsparents *= catsi[order[q]]; 533 combsparents *= catsi[order[q]];
507 } 534 }
508 } 535 }
509 tmp.push_back(combsparents);//made the whole parentset configure for a variable 536 tmp.push_back(combsparents);//made the whole parentset configure for a variable
510 /*std::cout << "This is combsparents: " << combsparents << endl; 537 /*std::cout << "This is combsparents: " << combsparents << endl;
511 std::cout << endl;*/ 538 std::cout << endl;*/
512 vector <int> hvect; 539 vector <int> hvect;
513 //hvect tells us which variables are being considered always the last variable is being considered 540 //hvect tells us which variables are being considered always the last variable is being considered
514 //e.g if ABC is our order and we are on i equals 1 then we are looking at relationships between A and B only 541 //e.g if ABC is our order and we are on i equals 1 then we are looking at relationships between A and B only
515 //continued: A is the only one that is either a parent or isn't a parent so hvect will be < 0 1 > 542 //continued: A is the only one that is either a parent or isn't a parent so hvect will be < 0 1 >
516 //for A C hvect will be < 0 2 > 543 //for A C hvect will be < 0 2 >
517 for (int h = 0; h < i; ++h) { 544 for (int h = 0; h < i; ++h) {
518 if (NewMat[h] == 1) { 545 if (NewMat[h] == 1) {
519 hvect.push_back(h); 546 hvect.push_back(h);
520 } 547 }
521 } 548 }
522 hvect.push_back(numparnts); 549 hvect.push_back(numparnts);
523 size_t shvect = hvect.size(); 550 size_t shvect = hvect.size();
524 //Prints out hvect 551 //Prints out hvect
525 /*for (int u = 0; u < shvect; ++u) { 552 /*for (int u = 0; u < shvect; ++u) {
526 std::cout << hvect[u] << " "; 553 std::cout << hvect[u] << " ";
527 }*/ 554 }*/
528 //std::cout << endl; 555 //std::cout << endl;
529 //Counting the amount of values in the data that have that particular parent combination 556 //Counting the amount of values in the data that have that particular parent combination
530 vector <int> Nijkovercombos; 557 vector <int> Nijkovercombos;
531 for (int last = min_cat[order[numparnts]]; last <= max_cat[order[numparnts]]; ++last) { 558 for (int last = min_cat[order[numparnts]]; last <= max_cat[order[numparnts]]; ++last) {
532 //(333)Creating a vector that uses the right combination 559 //(333)Creating a vector that uses the right combination
533 /*std::cout << "This is for " << i << " place in the order with value of variable equal to" << last << endl;*/ 560 /*std::cout << "This is for " << i << " place in the order with value of variable equal to" << last << endl;*/
534 vector <int> Test(shvect, last), maxtest; 561 vector <int> Test(shvect, last), maxtest;
535 for (int p = 0; p < (shvect-1); ++p) { 562 for (int p = 0; p < (shvect-1); ++p) {
536 Test[p] = max_cat[order[hvect[p]]]; 563 Test[p] = max_cat[order[hvect[p]]];
537 } 564 }
538 maxtest = Test; 565 maxtest = Test;
539 for (int i2Nloopy = 0; i2Nloopy < combsparents; ++i2Nloopy) { 566 for (int i2Nloopy = 0; i2Nloopy < combsparents; ++i2Nloopy) {
540 //std::cout << endl; 567 //std::cout << endl;
541 /*std::cout << endl; 568 /*std::cout << endl;
542 std::cout << "This is i2Nloopy: " << (i2Nloopy + 1) << endl; 569 std::cout << "This is i2Nloopy: " << (i2Nloopy + 1) << endl;
543 std::cout << "Here comes the Test:" << endl;*/ 570 std::cout << "Here comes the Test:" << endl;*/
544 //(444)This sets up the process for changing 571 //(444)This sets up the process for changing
545 //NMpos tells me the position of the last non minimum value in the vector 572 //NMpos tells me the position of the last non minimum value in the vector
546 //We want to change when the position is the last position available in the vector 573 //We want to change when the position is the last position available in the vector
547 int NMpos = 0, minpos = 0; 574 int NMpos = 0, minpos = 0;
548 for (NMpos = (shvect - 2); NMpos >= 0; --NMpos) { 575 for (NMpos = (shvect - 2); NMpos >= 0; --NMpos) {
549 if (Test[NMpos] != min_cat[order[hvect[NMpos]]]) { 576 if (Test[NMpos] != min_cat[order[hvect[NMpos]]]) {
550 break; 577 break;
551 } 578 }
552 } 579 }
553 for (minpos = (shvect - 2); minpos >= 0; --minpos) { 580 for (minpos = (shvect - 2); minpos >= 0; --minpos) {
554 //minpos tells you the position of the last minimum value 581 //minpos tells you the position of the last minimum value
555 if (Test[minpos] == min_cat[order[hvect[minpos]]]) { 582 if (Test[minpos] == min_cat[order[hvect[minpos]]]) {
556 break; 583 break;
557 } 584 }
558 } 585 }
559 //Prints out Test 586 //Prints out Test
560 /*for (int u = 0; u < shvect; ++u) { 587 /*for (int u = 0; u < shvect; ++u) {
561 std::cout << Test[u] << " "; 588 std::cout << Test[u] << " ";
562 } 589 }
563 std::cout << endl; 590 std::cout << endl;
564 std::cout << endl; 591 std::cout << endl;
565 std::cout << endl;*/ 592 std::cout << endl;*/
566 //Count how many occurrences of the value are present in the data 593 //Count how many occurrences of the value are present in the data
567 int Nijk = 0; 594 int Nijk = 0;
568 for (int num2size = 0; num2size < tottuples; ++num2size) { 595 for (int num2size = 0; num2size < tottuples; ++num2size) {
569 int countcorrect = 0; 596 int countcorrect = 0;
570 for (size_t g = 0; g < Test.size(); ++g) { 597 for (size_t g = 0; g < Test.size(); ++g) {
571 //num2size cycles through tuples 598 //num2size cycles through tuples
572 //order[hvect[g]] represents the variable in the order that we are considering as a parent 599 //order[hvect[g]] represents the variable in the order that we are considering as a parent
573 if (DAT[num2size][order[hvect[g]]] == Test[g]) { 600 if (DAT[num2size][order[hvect[g]]] == Test[g]) {
574 countcorrect += 1; 601 countcorrect += 1;
575 } 602 }
576 } 603 }
577 if (countcorrect == Test.size()) { 604 if (countcorrect == Test.size()) {
578 Nijk += 1; 605 Nijk += 1;
579 } 606 }
580 } 607 }
581 //Nijkovercombos displays data as follows 608 //Nijkovercombos displays data as follows
582 //it starts with the smallest value for the last variable in hvect 609 //it starts with the smallest value for the last variable in hvect
583 //and the largest values in the first n-1 variables in hvect 610 //and the largest values in the first n-1 variables in hvect
584 //max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0 611 //max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0
585 //count,count,count,count e.g. 13, 2, 2, 3, 4, 10 612 //count,count,count,count e.g. 13, 2, 2, 3, 4, 10
586 Nijkovercombos.push_back(Nijk); 613 Nijkovercombos.push_back(Nijk);
587 //(666)Now that the values have been calculated find out what the next combination of variables should be 614 //(666)Now that the values have been calculated find out what the next combination of variables should be
588 if ((NMpos == -1) & (minpos == (shvect - 2))) { 615 if ((NMpos == -1) & (minpos == (shvect - 2))) {
589 //break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000 616 //break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000
590 break; 617 break;
591 } 618 }
592 if (minpos < NMpos) { 619 if (minpos < NMpos) {
593 Test[NMpos] = Test[NMpos] - 1; 620 Test[NMpos] = Test[NMpos] - 1;
594 } 621 }
595 else if (NMpos < minpos) { 622 else if (NMpos < minpos) {
596 Test[NMpos] = Test[NMpos] - 1; 623 Test[NMpos] = Test[NMpos] - 1;
597 for (int filler = NMpos + 1; filler < (shvect - 1); ++filler) { 624 for (int filler = NMpos + 1; filler < (shvect - 1); ++filler) {
598 Test[filler] = maxtest[filler]; 625 Test[filler] = maxtest[filler];
599 } 626 }
600 } 627 }
601 } 628 }
602 } 629 }
603 fullNijkvector.push_back(Nijkovercombos); 630 fullNijkvector.push_back(Nijkovercombos);
604 //(666)Now that the unique values have been calculated find out what the next combination of variables should be 631 //(666)Now that the unique values have been calculated find out what the next combination of variables should be
605 if ((PosOne == (numparnts - 1)) & (SumOnes == it)) { 632 if ((PosOne == (numparnts - 1)) & (SumOnes == it)) {
606 break; 633 break;
607 } 634 }
608 else if ((PosOne == (numparnts - 1)) & (SumOnes != it)) { 635 else if ((PosOne == (numparnts - 1)) & (SumOnes != it)) {
609 for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) { 636 for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) {
610 //NxtOne tells you the position of the next closest number 1 that we would 637 //NxtOne tells you the position of the next closest number 1 that we would
611 //like to change the position of (we will call it the important number one) 638 //like to change the position of (we will call it the important number one)
612 //SumOnes3 helps keep track of the sum of all ones between now and the next important number one 639 //SumOnes3 helps keep track of the sum of all ones between now and the next important number one
613 SumOnes3 += NewMat[NxtOne]; 640 SumOnes3 += NewMat[NxtOne];
614 if (SumOnes3 == (SumOnes2 + 1)) { 641 if (SumOnes3 == (SumOnes2 + 1)) {
615 break; 642 break;
616 } 643 }
617 } 644 }
618 if (SumOnes4 == (it - 1)) { 645 if (SumOnes4 == (it - 1)) {
619 //If all except one of the 1's are found in the last it - 1 columns 646 //If all except one of the 1's are found in the last it - 1 columns
620 for (int x = 0; x < numparnts; ++x) { 647 for (int x = 0; x < numparnts; ++x) {
621 if (((x <= (NxtOne + SumOnes3)) & (x > NxtOne)) | (x == (FrstOne + 1))) { 648 if (((x <= (NxtOne + SumOnes3)) & (x > NxtOne)) | (x == (FrstOne + 1))) {
622 //If 649 //If
623 NewMat[x] = 1; 650 NewMat[x] = 1;
624 } 651 }
625 else { 652 else {
626 NewMat[x] = 0; 653 NewMat[x] = 0;
627 } 654 }
628 } 655 }
629 } 656 }
630 else { 657 else {
631 for (int x = 0; x < numparnts; ++x) { 658 for (int x = 0; x < numparnts; ++x) {
632 if (((x <= (NxtOne + SumOnes3)) & (x > NxtOne)) | (x == FrstOne)) { 659 if (((x <= (NxtOne + SumOnes3)) & (x > NxtOne)) | (x == FrstOne)) {
633 //If the position is that of the first 1 or it falls between the changed number one and the total 660 //If the position is that of the first 1 or it falls between the changed number one and the total
634 //amount of ones that are on that side of the zero 10111 661 //amount of ones that are on that side of the zero 10111
635 NewMat[x] = 1; 662 NewMat[x] = 1;
636 } 663 }
637 else if ((x != FrstOne) & (x != NxtOne) & (NewMat[x] == 1) & (x < PosOne2)) { 664 else if ((x != FrstOne) & (x != NxtOne) & (NewMat[x] == 1) & (x < PosOne2)) {
638 //If it is not the position of the first 1 and it is not the position of the 1 whose position we are interested in changing 665 //If it is not the position of the first 1 and it is not the position of the 1 whose position we are interested in changing
639 //and the previous value at this position was 1 and the postion is below the value of the first zero spotted from the right 666 //and the previous value at this position was 1 and the postion is below the value of the first zero spotted from the right
640 NewMat[x] = 1; 667 NewMat[x] = 1;
641 } 668 }
642 else { 669 else {
643 NewMat[x] = 0; 670 NewMat[x] = 0;
644 } 671 }
645 } 672 }
646 } 673 }
647 } 674 }
648 else if ((PosOne != (numparnts - 1)) & (SumOnes != it)) { 675 else if ((PosOne != (numparnts - 1)) & (SumOnes != it)) {
649 for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) { 676 for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) {
650 //NxtOne tells you the position of the next closest number 1 that we would 677 //NxtOne tells you the position of the next closest number 1 that we would
651 //like to change the position of (we will call it the important number one) 678 //like to change the position of (we will call it the important number one)
652 //SumOnes3 helps keep track of the sum of all ones between now and the next important number one 679 //SumOnes3 helps keep track of the sum of all ones between now and the next important number one
653 SumOnes3 += NewMat[NxtOne]; 680 SumOnes3 += NewMat[NxtOne];
654 if (SumOnes3 == 1) { 681 if (SumOnes3 == 1) {
655 break; 682 break;
656 } 683 }
657 } 684 }
658 if (it != 1) { 685 if (it != 1) {
659 for (int x = 0; x < numparnts; ++x) { 686 for (int x = 0; x < numparnts; ++x) {
660 if (x == (NxtOne + 1)) { 687 if (x == (NxtOne + 1)) {
661 NewMat[x] = 1; 688 NewMat[x] = 1;
662 } 689 }
663 else if (x == NxtOne) { 690 else if (x == NxtOne) {
664 NewMat[x] = 0; 691 NewMat[x] = 0;
665 } 692 }
666 else if ((NewMat[x] == 1) & (x != NxtOne)) { 693 else if ((NewMat[x] == 1) & (x != NxtOne)) {
667 NewMat[x] = 1; 694 NewMat[x] = 1;
668 } 695 }
669 else { 696 else {
670 NewMat[x] = 0; 697 NewMat[x] = 0;
671 } 698 }
672 } 699 }
673 } 700 }
674 else { 701 else {
675 for (int x = 0; x < numparnts; ++x) { 702 for (int x = 0; x < numparnts; ++x) {
676 if ((x == (NxtOne + 1))) { 703 if ((x == (NxtOne + 1))) {
677 NewMat[x] = 1; 704 NewMat[x] = 1;
678 } 705 }
679 else { 706 else {
680 NewMat[x] = 0; 707 NewMat[x] = 0;
681 } 708 }
682 } 709 }
683 } 710 }
684 } 711 }
685 } 712 }
686 713
687 } 714 }
688 } 715 }
689 ParentCombos.push_back(tmp); 716 ParentCombos.push_back(tmp);
690 indexofvar.push_back(tempstring); 717 indexofvar.push_back(tempstring);
691 } 718 }
692 719
693 } 720 }
694 //std::cout << endl; 721 //std::cout << endl;
695 std::cout << endl; 722 std::cout << endl;
696 //printing out the ParentCombos matrix just created above 723 //printing out the ParentCombos matrix just created above
697 /*for (size_t i = 0; i < ParentCombos.size(); ++i) { 724 /*for (size_t i = 0; i < ParentCombos.size(); ++i) {
698 for (size_t j = 0; j < ParentCombos[i].size(); ++j) { 725 for (size_t j = 0; j < ParentCombos[i].size(); ++j) {
699 std::cout << ParentCombos[i][j] << " "; 726 std::cout << ParentCombos[i][j] << " ";
700 } 727 }
701 std::cout << endl; 728 std::cout << endl;
702 }*/ 729 }*/
703 std::cout << endl; 730 std::cout << endl;
704 //printing out the fullNijkvector matrix just created above 731 //printing out the fullNijkvector matrix just created above
705 /*for (size_t i = 0; i < fullNijkvector.size(); ++i) { 732 /*for (size_t i = 0; i < fullNijkvector.size(); ++i) {
706 for (size_t j = 0; j < fullNijkvector[i].size(); ++j) { 733 for (size_t j = 0; j < fullNijkvector[i].size(); ++j) {
707 std::cout << fullNijkvector[i][j] << " "; 734 std::cout << fullNijkvector[i][j] << " ";
708 } 735 }
709 std::cout << endl; 736 std::cout << endl;
710 } 737 }
711 std::cout << endl;*/ 738 std::cout << endl;*/
712 //Print out Data 739 //Print out Data
713 /*for (int i = 0; i < DAT.size(); ++i) { 740 /*for (int i = 0; i < DAT.size(); ++i) {
714 for (int j = 0; j < DAT[i].size(); ++j) { 741 for (int j = 0; j < DAT[i].size(); ++j) {
715 std::cout << DAT[i][j] << " "; 742 std::cout << DAT[i][j] << " ";
716 } 743 }
717 std::cout << endl; 744 std::cout << endl;
718 }*/ 745 }*/
719 //Print out the families size 746 //Print out the families size
720 /*for (size_t i = 0; i < families.size(); ++i) { 747 /*for (size_t i = 0; i < families.size(); ++i) {
721 std::cout << families[i] << " "; 748 std::cout << families[i] << " ";
722 }*/ 749 }*/
723 750
724 //Obtaining the actual score from this information 751 //Obtaining the actual score from this information
725 //varinorder cycles through families (the amount of parent families that should be considered for the variable with a particular order starting 752 //varinorder cycles through families (the amount of parent families that should be considered for the variable with a particular order starting
726 //the first variable in the order) 753 //the first variable in the order)
727 //keeping track of the position within the fullNijkvector associated with the varinorder and the qi_Uialpha 754 //keeping track of the position within the fullNijkvector associated with the varinorder and the qi_Uialpha
728 int posinfull = 0; 755 int posinfull = 0;
729 //finlogscore is the final score in natural log format 756 //finlogscore is the final score in natural log format
730 long double finlogscore = 0.0; 757 long double finlogscore = 0.0;
731 vector< vector <double> > vecvarparset; 758 vector< vector <double> > vecvarparset;
732 for (size_t varinorder = 0; varinorder < families.size(); ++varinorder) { 759 for (size_t varinorder = 0; varinorder < families.size(); ++varinorder) {
733 //sumovUialpha is the the sum over all parent sets for a particular variable 760 //sumovUialpha is the the sum over all parent sets for a particular variable
734 long double sumovUialpha = 0.0; 761 long double sumovUialpha = 0.0;
735 //vector of all values of seclastgamma 762 //vector of all values of seclastgamma
736 vector <double> vec2ndlastgamma; 763 vector <double> vec2ndlastgamma;
737 long double maxseclastgamma; 764 long double maxseclastgamma;
738 //Uialpha cycles through all the parent sets for a particular family 765 //Uialpha cycles through all the parent sets for a particular family
739 for (int Uialpha = 0; Uialpha < families[varinorder]; ++Uialpha) { 766 for (int Uialpha = 0; Uialpha < families[varinorder]; ++Uialpha) {
740 // nijkprime represents the value of 1/(ri * qi) 767 // nijkprime represents the value of 1/(ri * qi)
741 long double nijkprime, nijprime; 768 long double nijkprime, nijprime;
742 long double rij = catsi[order[varinorder]], PCs = ParentCombos[varinorder][Uialpha]; 769 long double rij = catsi[order[varinorder]], PCs = ParentCombos[varinorder][Uialpha];
743 770
744 nijprime = 1.0 / (PCs); 771 nijprime = 1.0 / (PCs);
745 nijkprime = 1.0 / (rij * PCs); 772 nijkprime = 1.0 / (rij * PCs);
746 //seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic 773 //seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic
747 long double seclastgamma = 0.0; 774 long double seclastgamma = 0.0;
748 //qi_Uialpha cycles through the combinations for the parents in a set 775 //qi_Uialpha cycles through the combinations for the parents in a set
749 for (int qi_Uialpha = 0; qi_Uialpha < ParentCombos[varinorder][Uialpha];++qi_Uialpha) { 776 for (int qi_Uialpha = 0; qi_Uialpha < ParentCombos[varinorder][Uialpha];++qi_Uialpha) {
750 long double lastgamma = 0.0; 777 long double lastgamma = 0.0;
751 long double nij = 0.0; 778 long double nij = 0.0;
752 //countijk cycles through the categories of the variable with a particular order 779 //countijk cycles through the categories of the variable with a particular order
753 //catsi is in the order that data is input and so one must use the order[varinorder] to first obtain the variable that we are referring to 780 //catsi is in the order that data is input and so one must use the order[varinorder] to first obtain the variable that we are referring to
754 //and then find the categories for it 781 //and then find the categories for it
755 for (int countijk = 0; countijk < catsi[order[varinorder]]; ++countijk) { 782 for (int countijk = 0; countijk < catsi[order[varinorder]]; ++countijk) {
756 long double topy; 783 long double topy;
757 //rightcol lets you find the right column/position of the value that you need for a particular category within the 784 //rightcol lets you find the right column/position of the value that you need for a particular category within the
758 int rightcol = qi_Uialpha + (countijk * ParentCombos[varinorder][Uialpha]); 785 int rightcol = qi_Uialpha + (countijk * ParentCombos[varinorder][Uialpha]);
759 nij += fullNijkvector[posinfull][rightcol]; 786 nij += fullNijkvector[posinfull][rightcol];
760 topy = (nijkprime + fullNijkvector[posinfull][rightcol]); 787 topy = (nijkprime + fullNijkvector[posinfull][rightcol]);
761 788
762 //Using boost lgamma function for the product over categories and parent combinations 789 //Using boost lgamma function for the product over categories and parent combinations
763 lastgamma += boost::math::lgamma(topy) - boost::math::lgamma(nijkprime); 790 lastgamma += boost::math::lgamma(topy) - boost::math::lgamma(nijkprime);
764 791
765 } 792 }
766 long double boty = nij + nijprime; 793 long double boty = nij + nijprime;
767 seclastgamma += lastgamma + boost::math::lgamma(nijprime) - boost::math::lgamma(boty); 794 seclastgamma += lastgamma + boost::math::lgamma(nijprime) - boost::math::lgamma(boty);
768 795
769 } 796 }
770 vec2ndlastgamma.push_back(seclastgamma); 797 vec2ndlastgamma.push_back(seclastgamma);
771 798
772 799
773 800
774 //Calculate sumovUialpha based on the logsumexp concept 801 //Calculate sumovUialpha based on the logsumexp concept
775 if (Uialpha + 1 == families[varinorder]) { 802 if (Uialpha + 1 == families[varinorder]) {
776 803
777 for (size_t que = 0; que < vec2ndlastgamma.size(); ++que) { 804 for (size_t que = 0; que < vec2ndlastgamma.size(); ++que) {
778 //change the value of maxseclastgamma if new value is larger than the previous value 805 //change the value of maxseclastgamma if new value is larger than the previous value
779 if (que == 0) { 806 if (que == 0) {
780 maxseclastgamma = vec2ndlastgamma[0]; 807 maxseclastgamma = vec2ndlastgamma[0];
781 } 808 }
782 else { 809 else {
783 if (maxseclastgamma < vec2ndlastgamma[que]) { 810 if (maxseclastgamma < vec2ndlastgamma[que]) {
784 maxseclastgamma = vec2ndlastgamma[que]; 811 maxseclastgamma = vec2ndlastgamma[que];
785 } 812 }
786 } 813 }
787 } 814 }
788 for (size_t what = 0; what < vec2ndlastgamma.size(); ++what) { 815 for (size_t what = 0; what < vec2ndlastgamma.size(); ++what) {
789 sumovUialpha += exp(vec2ndlastgamma[what] - maxseclastgamma); 816 sumovUialpha += exp(vec2ndlastgamma[what] - maxseclastgamma);
790 } 817 }
791 //add info on parent set scores for each variable to this vector of vectors 818 //add info on parent set scores for each variable to this vector of vectors
792 vecvarparset.push_back(vec2ndlastgamma); 819 vecvarparset.push_back(vec2ndlastgamma);
793 } 820 }
794 821
795 /*std::cout << endl; 822 /*std::cout << endl;
796 std::cout << seclastgamma; 823 std::cout << seclastgamma;
797 std::cout << endl;*/ 824 std::cout << endl;*/
798 posinfull += 1; 825 posinfull += 1;
799 //std::cout << posinfull << endl; 826 //std::cout << posinfull << endl;
800 } 827 }
801 finlogscore += log(sumovUialpha) + maxseclastgamma; 828 finlogscore += log(sumovUialpha) + maxseclastgamma;
802 } 829 }
803 830
804 vector < vector<string> > parSet; 831 vector < vector<string> > parSet;
805 vector< map <double, string, greater <double> > > parSetScoreSorted; 832 vector< map <double, string, greater <double> > > parSetScoreSorted;
806 vector< map <double, string, greater <double> > > strucScore; 833 vector< map <double, string, greater <double> > > strucScore;
807 834
808 835
809 //Below is another way to match the index or label sets with the scores sets, and store the (score, label) into a map vector. And for each vector element the map is a sorted map. 836 //Below is another way to match the index or label sets with the scores sets, and store the (score, label) into a map vector. And for each vector element the map is a sorted map.
810 837
811 for (unsigned i = 0; i < indexofvar.size(); ++i){ 838 for (unsigned i = 0; i < indexofvar.size(); ++i){
812 839
813 map <double, string, greater <double> > tempMap; 840 map <double, string, greater <double> > tempMap;
814 for (unsigned j=0; j< indexofvar[i].size(); ++j){ 841 for (unsigned j=0; j< indexofvar[i].size(); ++j){
815 tempMap.insert(make_pair(vecvarparset[i][j], indexofvar[i][j])); 842 tempMap.insert(make_pair(vecvarparset[i][j], indexofvar[i][j]));
816 } 843 }
817 parSetScoreSorted.push_back(tempMap); 844 parSetScoreSorted.push_back(tempMap);
818 845
819 } 846 }
820 847
821 848
822 pair <double, string> bestStrScore; 849 pair <double, string> bestStrScore;
823 double bestScore = 0; 850 double bestScore = 0;
824 string bestLable; 851 string bestLable;
825 852
826 853
827 for (unsigned i = 0; i < parSetScoreSorted.size(); ++i){ 854 for (unsigned i = 0; i < parSetScoreSorted.size(); ++i){
828 map <double, string> :: iterator itr; 855 map <double, string> :: iterator itr;
829 itr = parSetScoreSorted[i].begin(); 856 itr = parSetScoreSorted[i].begin();
830 bestScore = bestScore + (itr->first); 857 bestScore = bestScore + (itr->first);
831 bestLable = bestLable +(itr->second); 858 bestLable = bestLable +(itr->second);
832 } 859 }
833 860
834 bestStrScore = make_pair(bestScore, bestLable);//This is the best score. 861 bestStrScore = make_pair(bestScore, bestLable);//This is the best score.
835 862
836 vector < pair <double, string> > sortedStru;//This store all the structures in the percentage. 863 vector < pair <double, string> > sortedStru;//This store all the structures in the percentage.
837 vector < vector < pair <double, string > > > deltaC; 864 vector < vector < pair <double, string > > > deltaC;
838 865
839 for (unsigned l = 1; l< parSetScoreSorted.size(); ++l){ 866 for (unsigned l = 1; l< parSetScoreSorted.size(); ++l){
840 map <double, string> :: iterator itr0, itr1; 867 map <double, string> :: iterator itr0, itr1;
841 vector < pair <double, string > > tempDelta; 868 vector < pair <double, string > > tempDelta;
842 double tempDeltaS; 869 double tempDeltaS;
843 string tempDeltaL; 870 string tempDeltaL;
844 itr0 = parSetScoreSorted[l].begin(); 871 itr0 = parSetScoreSorted[l].begin();
845 itr1 = parSetScoreSorted[l].begin(); 872 itr1 = parSetScoreSorted[l].begin();
846 double tem1=exp(itr1->first), tem2 = exp(itr1->first); 873 double tem1=exp(itr1->first), tem2 = exp(itr1->first);
847 for (unsigned m = 1; m< parSetScoreSorted[l].size(); ++m){ 874 for (unsigned m = 1; m< parSetScoreSorted[l].size(); ++m){
848 tem1 = tem2; 875 tem1 = tem2;
849 itr1 = ++itr1; 876 itr1 = ++itr1;
850 tem2 = tem1 + exp(itr1->first); 877 tem2 = tem1 + exp(itr1->first);
851 double tem = (tem1/tem2)*100; 878 double tem = (tem1/tem2)*100;
852 if(tem <= PERCENT){ 879 if(tem <= PERCENT){
853 tempDeltaS = (itr1->first)-(itr0->first); 880 tempDeltaS = (itr1->first)-(itr0->first);
854 tempDeltaL = itr1->second; 881 tempDeltaL = itr1->second;
855 tempDelta.push_back(make_pair(tempDeltaS, tempDeltaL)); 882 tempDelta.push_back(make_pair(tempDeltaS, tempDeltaL));
856 } 883 }
857 } 884 }
858 885
859 deltaC.push_back( tempDelta); 886 deltaC.push_back( tempDelta);
860 } 887 }
861 888
862 889
863 for(unsigned i=0; i< deltaC.size(); ++i){ 890 for(unsigned i=0; i< deltaC.size(); ++i){
864 for (unsigned j=0; j< deltaC[i].size(); ++j) 891 for (unsigned j=0; j< deltaC[i].size(); ++j)
865 { 892 {
866 double score = bestStrScore.first + deltaC[i][j].first; 893 double score = bestStrScore.first + deltaC[i][j].first;
867 string lab = bestStrScore.second; 894 string lab = bestStrScore.second;
868 findAndReplaceAll(lab, parSetScoreSorted[i+1].begin()->second, deltaC[i][j].second); 895 findAndReplaceAll(lab, parSetScoreSorted[i+1].begin()->second, deltaC[i][j].second);
869 sortedStru.push_back(make_pair(score, lab)); 896 sortedStru.push_back(make_pair(score, lab));
870 } 897 }
871 } 898 }
872 899
873 sort(sortedStru.begin(),sortedStru.end(),compareDe); 900 sort(sortedStru.begin(),sortedStru.end(),compareDe);
874 901
902 vector <double> strP;
903 double s = exp(bestScore);
875 904
876 905 for(unsigned i=0; i< sortedStru.size(); ++i ){
877 /* for (unsigned l = 1; l< parSetScoreSorted.size(); ++l){ //l means each variable 906 double s0 = s;
878 for (unsigned m = contl; m < parSetScoreSorted[l].size(); ++m){ //m means the number of parents sets 907 s = s + exp(sortedStru[i].first);
879 contl =m; 908 double te = s0/s*100;
880 vector < pair <string, double > > delta; // store the different of the highest score and the second highest score 909 strP.push_back(te);
881 for (unsigned i = l; i < parSetScoreSorted.size(); ++i){
882 map <double, string> :: iterator itr0, itr1;
883 double tempDeltaS;
884 string tempDeltaL;
885 itr0 = parSetScoreSorted[i].begin();
886 itr1 = itr0;
887 for (int kr = 0; kr < contl; ++kr){
888 itr1 = itr0++;
889 }
890 tempDeltaS = (itr1->first)-(itr0->first);
891 tempDeltaL = int_to_str(i);
892 delta.push_back(make_pair(tempDeltaL, tempDeltaS));
893 }
894
895 sort(delta.begin(),delta.end(),compareI);
896
897 double const stopLimit = PERCENT;
898 double tempbeforeS = bestStrScore.first, temafterS;
899
900 for (unsigned i = 0; i < delta.size(); ++i){
901 pair <double, string > temppair;
902 double tempLime;
903 int ind = str_to_int(delta[i].first);
904 string s = parSetScoreSorted[ind].begin()->second;
905 string sRep = (++parSetScoreSorted[ind].begin())->second;
906 temppair = bestStrScore;
907 findAndReplaceAll(temppair.second, s, sRep);
908 temppair.first = temppair.first - delta[i].second;
909
910 temafterS = logAB(tempbeforeS, temppair.first);
911
912 tempLime = persentageXofY(temafterS, tempbeforeS) ;
913
914 tempbeforeS = temafterS;
915
916 if(tempLime > stopLimit){
917 goto finish;
918 }
919
920 sortedStru.push_back(temppair);
921
922 // std::cout << ind << " "<< s << " "<< sRep << " "<< temppair.first << " " << temppair.second <<endl;
923 // std::cout << temppair.first << " " << temppair.second <<endl;
924 // std::cout << std::endl;
925
926 }
927
928 }
929 910
930 } 911 }
931 912
932 finish:*/ 913 /* double fi = s/exp(finlogscore)*100;
933 914 strP.push_back(fi);*/
934 915
935 std::cout << std::endl; 916 std::cout << std::endl;
936 std::cout << "Total Score: "<< boost::lexical_cast<string>(finlogscore) << std::endl; 917 std::cout << "Total Score: "<< boost::lexical_cast<string>(finlogscore) << std::endl;
937 std::cout << std::endl; 918 std::cout << std::endl;
938 919
939 920
921 vector < pair <double, string> > bestorders;
922 double struPer = 0;
923
940 cout << "Best several structures are:" << endl; 924 cout << "Best several structures are:" << endl;
941 cout << bestScore << " : " << bestLable << endl; 925 cout << bestScore << " : " << bestLable << endl;
942 926
927 bestorders.push_back(make_pair(bestScore, bestLable));
928 struPer = struPer + exp(bestScore);
929
943 for (unsigned i=0; i< sortedStru.size(); ++i ) 930 for (unsigned i=0; i< sortedStru.size(); ++i )
944 { 931 {
945 cout << sortedStru[i].first << " : " << sortedStru[i].second << endl; 932 cout << sortedStru[i].first << " : " << sortedStru[i].second << endl;
933 bestorders.push_back(make_pair(sortedStru[i].first, sortedStru[i].second));
934 struPer = struPer + exp(sortedStru[i].first);
935 if(strP[i+1] > PERCENT){
936 break;
937 }
938 }
939
940 /* cout << "Best several structures are:" << endl;
941 cout << bestScore << " : " << bestLable << " " << endl;
942
943 for (unsigned i=0; i< sortedStru.size(); ++i )
944 {
945 cout << sortedStru[i].first << " : " << sortedStru[i].second << " "<< endl;
946 } 946 }
947 947
948 std::cout << std::endl; 948 std::cout << std::endl;
949 std::cout << std::endl; 949 std::cout << std::endl;*/
950 950
951 for(unsigned i=0; i< parSetScoreSorted.size(); ++i){ 951 /* for(unsigned i=0; i< parSetScoreSorted.size(); ++i){
952 map <double, string> :: iterator itr; 952 map <double, string> :: iterator itr;
953 std::cout << "For variable "<< i <<":"<< endl; 953 std::cout << "For variable "<< order[i] <<":"<< endl;
954 for (itr = parSetScoreSorted[i].begin(); itr != parSetScoreSorted[i].end(); ++itr) 954 for (itr = parSetScoreSorted[i].begin(); itr != parSetScoreSorted[i].end(); ++itr)
955 { 955 {
956 cout << itr->second << " : " << itr->first << "; "; 956 cout << itr->second << " : " << itr->first << "; ";
957 } 957 }
958 cout << endl; 958 cout << endl;
959 cout << endl; 959 cout << endl;
960 960
961 }*/
962
963
964 string orderstring;
965 for (unsigned i=0; i< order.size(); ++i){
966 orderstring.append(int_to_str(order[i])+" ");
961 } 967 }
962 968
969 ofstream myfile;
970 myfile.open( filepath.c_str());
963 971
964 /* ofstream myfile;
965 myfile.open("/home/zgong001/Documents/SprinklerDataset/bestStructure.txt");
966 if (myfile.is_open()) 972 if (myfile.is_open())
967 { 973 {