Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
metric.cpp
Go to the documentation of this file.
1 #include "metric.h"
2 
3 
4 /* metric implemenation --------------------------------------------------*/
5 
6 metric::metric(std::vector<points*> datasets1, std::vector<points*> datasets2,
7  int var_ds1, int var_ds2,
8  ofstream &outfile, string outprefix) : outfile(outfile)
9 {
10 
11  this ->restrictToBoxRange = false;
12 
13  this->outprefix = outprefix;
14 
15  this->datasets1 = datasets1;
16  this->datasets2 = datasets2;
17 
18  this->var_ds1 = var_ds1;
19  this->var_ds2 = var_ds2;
20 }
21 
22 void metric::set_datasets1(std::vector<points*> datasets1){
23  this->datasets1 = datasets1;
24 }
25 
26 void metric::set_datasets2(std::vector<points*> datasets2){
27  this->datasets2 = datasets2;
28 }
29 
30 string metric::get_res(){
31  return string("No Metric Performed");
32 }
33 
34 string metric::get_metric_name(){
35  return string("Metric");
36 }
37 
38 void metric::error_out(string err_ps){
39 
40  string out;
41  out = "ERROR | Metric \"" + this->get_metric_name() + "\": " + err_ps;
42 
43 
44  outfile << out << std::endl;
45  cout << out << std::endl;
46 }
47 
48 
49 
50 /* mean square error -----------------------------------------------------*/
51 /*
52 mse::mse(std::vector<points*> datasets1, std::vector<points*> datasets2,
53  int var_ds1, int var_ds2,
54  ofstream &outfile, string outprefix)
55  : metric(datasets1, datasets2, var_ds1, var_ds2, outfile, outprefix){
56 
57  return;
58 }
59 
60 string mse::get_res(){
61  string ret;
62  std::stringstream out;
63  out.precision(10);
64  out << compare();
65  return out.str();
66 }
67 
68 long double mse::compare(){
69 
70  if(!check_card_match())
71  return -std::numeric_limits<long double>::infinity();
72 
73  long double total_sq_err = 0;
74 
75  for(int i=0; i<dataset1->get_num_points(); i++)
76  {
77 
78  pnt p1 = dataset1->get_point(i);
79  pnt p2 = dataset2->get_point(i);
80 
81  if( this ->restrictToBoxRange &&
82  this ->range ->hasPoint( p1.vals[0], p1.vals[1], 0.0 ) &&
83  this ->range ->hasPoint( p2.vals[0], p2.vals[1], 0.0 ) )
84  {
85 
86  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) ){
87  total_sq_err += pow( p1.vals[var_ds1] - p2.vals[var_ds2] , 2 );
88  }
89 
90  }
91  else if( ! this ->restrictToBoxRange )
92  {
93 
94  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) ){
95  total_sq_err += pow( p1.vals[var_ds1] - p2.vals[var_ds2] , 2 );
96  }
97 
98  }
99 
100  } // End for all points
101 
102  long double res = total_sq_err / ((long double)dataset1->get_num_points());
103 
104  return res;
105 }
106 
107 bool mse::check_card_match(){
108 
109  if(dataset1->get_num_points() == dataset2->get_num_points())
110  return true;
111 
112  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
113  err_ps += itoa(dataset1->get_num_points());
114  err_ps += ", Dataset 2 = ";
115  err_ps += itoa(dataset2->get_num_points());
116 
117  error_out(err_ps);
118 
119  return false;
120 }
121 
122 string mse::get_metric_name(){
123  return string("Mean Square Error");
124 }
125 */
126 
127 
128 /* root mean square error ------------------------------------------------*/
129 /*
130 rmse::rmse(std::vector<points*> datasets1, points *dataset2, int var_ds1, int var_ds2,
131  ofstream &outfile, string outprefix)
132  : metric(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
133 
134  return;
135 }
136 
137 string rmse::get_res(){
138  string ret;
139  std::stringstream out;
140  out.precision(10);
141  out << compare();
142  return out.str();
143 }
144 
145 long double rmse::compare(){
146 
147  if(!check_card_match())
148  return -std::numeric_limits<long double>::infinity();
149 
150  long double total_sq_err = 0;
151 
152  for(int i=0; i<dataset1->get_num_points(); i++)
153  {
154 
155  pnt p1 = dataset1->get_point(i);
156  pnt p2 = dataset2->get_point(i);
157 
158  if( this ->restrictToBoxRange &&
159  this ->range ->hasPoint( p1.vals[0], p1.vals[1], 0.0 ) &&
160  this ->range ->hasPoint( p2.vals[0], p2.vals[1], 0.0 ) )
161  {
162 
163  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) )
164  total_sq_err += pow( p1.vals[var_ds1] - p2.vals[var_ds2], 2 );
165  }
166  else if( ! this ->restrictToBoxRange )
167  {
168  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) )
169  total_sq_err += pow( p1.vals[var_ds1] - p2.vals[var_ds2], 2 );
170  }
171 
172  } // End for all points
173 
174  long double res = sqrt(total_sq_err / ((long double)dataset1->get_num_points()));
175 
176  return res;
177 
178 }
179 
180 bool rmse::check_card_match(){
181 
182  if(dataset1->get_num_points() == dataset2->get_num_points())
183  return true;
184 
185  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
186  err_ps += itoa(dataset1->get_num_points());
187  err_ps += ", Dataset 2 = ";
188  err_ps += itoa(dataset2->get_num_points());
189 
190  error_out(err_ps);
191 
192  return false;
193 }
194 
195 string rmse::get_metric_name(){
196  return string("Root Mean Square Error");
197 }
198 */
199 
200 
201 /* simple cross correlation / template matching---------------------------*/
202 /*
203 scc::scc(std::vector<points*> datasets1, points *dataset2, int var_ds1, int var_ds2,
204  ofstream &outfile, string outprefix)
205  : metric(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
206 
207  return;
208 }
209 
210 string scc::get_res(){
211  string ret;
212  std::stringstream out;
213  out.precision(10);
214  out << compare();
215  return out.str();
216 }
217 
218 long double scc::compare(){
219 
220  if(!check_card_match())
221  return -std::numeric_limits<long double>::infinity();
222 
223  long double total = 0;
224 
225  for(int i=0; i<dataset1->get_num_points(); i++)
226  {
227 
228  pnt p1 = dataset1->get_point(i);
229  pnt p2 = dataset2->get_point(i);
230 
231  if( this ->restrictToBoxRange &&
232  this ->range ->hasPoint( p1.vals[0], p1.vals[1], 0.0 ) &&
233  this ->range ->hasPoint( p2.vals[0], p2.vals[1], 0.0 ) )
234  {
235  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) )
236  total += p1.vals[var_ds1] * p2.vals[var_ds2];
237  }
238  else if( ! this ->restrictToBoxRange )
239  {
240  if( !isnan(p1.vals[var_ds1]) && !isnan(p2.vals[var_ds2]) )
241  total += p1.vals[var_ds1] * p2.vals[var_ds2];
242  }
243 
244  } // End for all points
245 
246  return total;
247 
248 }
249 
250 bool scc::check_card_match(){
251 
252  if(dataset1->get_num_points() == dataset2->get_num_points())
253  return true;
254 
255  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
256  err_ps += itoa(dataset1->get_num_points());
257  err_ps += ", Dataset 2 = ";
258  err_ps += itoa(dataset2->get_num_points());
259 
260  error_out(err_ps);
261 
262  return false;
263 }
264 
265 string scc::get_metric_name(){
266  return string("Cross-Correlation / Template Matching");
267 }
268 */
269 
270 
271 /* map -------------------------------------------------------------------*/
272 /*
273 map::map(std::vector<points*> datasets1, points *dataset2, int var_ds1, int var_ds2,
274  ofstream &outfile, string outprefix)
275  : metric(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
276 
277  return;
278 }
279 
280 rgb_color map::get_color_single_sided(long double pt_val, long double range, long double min){
281  int val = (int) (((pt_val - min) / range) * ((long double)(255)));
282 
283  if(isnan(pt_val)){
284 
285  rgb_color nanret;
286  nanret.r = 0;
287  nanret.g = 255;
288  nanret.b = 0;
289 
290  return nanret;
291  }
292 
293  rgb_color ret;
294  ret.r = gHeatedObject[val][0];
295  ret.g = gHeatedObject[val][1];
296  ret.b = gHeatedObject[val][2];
297 
298  return ret;
299 
300 }
301 
302 rgb_color map::get_color_double_sided(long double pt_val, long double range, long double min,
303  long double mid_point){
304 
305  if( mid_point >= min && mid_point <= (min+range) ){
306 
307  long double top_range = (range+min) - mid_point;
308  long double bot_range = mid_point - min;
309 
310  long double scale_factor;
311 
312  if(top_range > bot_range){
313  scale_factor = top_range / 127;
314  }
315  else if(bot_range >= top_range){
316  scale_factor = bot_range / 127;
317  }
318 
319  int gc_val = (int)(((pt_val-mid_point) / scale_factor) + (long double)(127));
320 
321  rgb_color ret;
322  ret.r = gBlueToCyan[gc_val][0];
323  ret.g = gBlueToCyan[gc_val][1];
324  ret.b = gBlueToCyan[gc_val][2];
325 
326  return ret;
327 
328  }
329 
330 
331  return get_color_single_sided(pt_val, range, min);
332 
333 }
334 */
335 
336 
337 /* difference map --------------------------------------------------------*/
338 /*
339 difmap::difmap(std::vector<indexed_points*> datasets1, indexed_points *dataset2, int var_ds1, int var_ds2,
340  ofstream &outfile, string outprefix)
341  : map(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
342 
343  return;
344 }
345 
346 string difmap::get_res(){
347  return compare();
348 }
349 
350 string difmap::compare(){
351 
352  //Difference map requires indexed data
353  //Try to convert the datasets to indexed_points datasets
354  indexed_points *ds1_ip = dynamic_cast<indexed_points*>(dataset1);
355  indexed_points *ds2_ip = dynamic_cast<indexed_points*>(dataset2);
356 
357  if(ds1_ip == NULL || ds2_ip == NULL)
358  return string("Improper Datasets");
359 
360  string ret("");
361 
362  if(!check_card_match()){
363  ret = "No Output";
364  return ret;
365  }
366 
367  //Create a new dataset to store the results
368  indexed_points* write_points = ds2_ip->clone();
369 
370  for(int x=0; x<ds1_ip->get_dim()[0]; x++)
371  {
372  for(int y=0; y<ds1_ip->get_dim()[1]; y++)
373  {
374  for(int z=0; z<ds1_ip->get_dim()[2]; z++)
375  {
376 
377  //TODO: loop over every dep_var
378 
379  std::vector<int> loc;
380  loc.push_back( x );
381  loc.push_back( y );
382  loc.push_back( z );
383 
384  pnt p1 = ds1_ip ->get_indexed_point( loc );
385  pnt p2 = ds2_ip ->get_indexed_point( loc );
386 
387  if( this ->restrictToBoxRange &&
388  this ->range ->hasPoint( p1.vals[0],p1.vals[1],0.0 ) &&
389  this ->range ->hasPoint( p2.vals[0],p2.vals[1],0.0 ) )
390  {
391 
392 
393  //Get the difference at each point
394  long double pt_val = fabs( ds1_ip->get_indexed_point(loc).vals[var_ds1]
395  - ds2_ip->get_indexed_point(loc).vals[var_ds2] );
396 
397  int index = x * write_points->get_dim()[1]
398  + y * write_points->get_dim()[2]
399  + z;
400 
401  //Store difference to dataset
402  write_points->set_point_val(index, var_ds1, pt_val);
403  }
404  else if( ! this ->restrictToBoxRange )
405  {
406  std::vector<int> loc;
407  loc.push_back(x);
408  loc.push_back(y);
409  loc.push_back(z);
410 
411  //Get the difference at each point
412  long double pt_val = fabs( ds1_ip->get_indexed_point(loc).vals[var_ds1]
413  - ds2_ip->get_indexed_point(loc).vals[var_ds2] );
414 
415  int index = x * write_points->get_dim()[1]
416  + y * write_points->get_dim()[2]
417  + z;
418 
419  //Store difference to dataset
420  write_points->set_point_val(index, var_ds1, pt_val);
421  }
422 
423  }
424  }
425  }
426 
427  //Construct file name
428  string output_name(outprefix);
429  output_name += "_difmap_";
430  output_name += "var1_";
431  output_name += itoa(var_ds1);
432  output_name += "_var2_";
433  output_name += itoa(var_ds2);
434  output_name += ".dat";
435 
436  //Write to tecplot file
437  TecplotOrderedWriter t_writer;
438  t_writer.setOutputFilename(output_name);
439 
440  if(!t_writer.init())
441  return string("Output file could not be opened");
442 
443  //Create variable names
444  std::vector<string> variables;
445  variables.push_back("Dep1");
446  variables.push_back("Dep2");
447  variables.push_back("Comparison");
448 
449  t_writer.addOrderedPartition(string("Partition 1"), write_points, var_ds1);
450  t_writer.writeFile(string("Comparison"), variables);
451 
452  delete write_points;
453 
454  return output_name;
455 }
456 
457 bool difmap::check_card_match(){
458  if(dataset1->get_num_points() != dataset2->get_num_points()){
459 
460  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
461  err_ps += itoa(dataset1->get_num_points());
462  err_ps += ", Dataset 2 = ";
463  err_ps += itoa(dataset2->get_num_points());
464 
465  error_out(err_ps);
466 
467  return false;
468  }
469 
470  //TODO
471  //Check that dimensions match
472 
473  return true;
474 }
475 
476 string difmap::get_metric_name(){
477  return string("Difference Map");
478 }
479 
480 long double difmap::get_max_value(int ds){
481 
482  long double max = -std::numeric_limits<long double>::infinity();
483 
484  int size = dataset1->get_num_points();
485 
486  for(int i=0; i<size; i++){
487 
488  if(ds == 1 && dataset1->get_point(i).vals[var_ds1] > max){
489  max = dataset1->get_point(i).vals[var_ds1];
490  }
491  if(ds == 2 && dataset2->get_point(i).vals[var_ds2] > max){
492  max = dataset2->get_point(i).vals[var_ds2];
493  }
494  }
495 
496  return max;
497 }
498 
499 long double difmap::get_min_value(int ds){
500 
501  long double min = std::numeric_limits<long double>::infinity();
502 
503  int size = dataset1->get_num_points();
504 
505  for(int i=0; i<size; i++){
506 
507  if(ds == 1 && dataset1->get_point(i).vals[var_ds1] < min){
508  min = dataset1->get_point(i).vals[var_ds1];
509  }
510  if(ds == 2 && dataset2->get_point(i).vals[var_ds2] < min){
511  min = dataset2->get_point(i).vals[var_ds2];
512  }
513 
514  }
515 
516  return min;
517 }
518 
519 
520 long double difmap::get_max_dif(){
521 
522  long double maxdif = -std::numeric_limits<long double>::infinity();
523 
524  int size = dataset1->get_num_points();
525 
526  for(int i=0; i<size; i++){
527 
528  long double dif = fabs( dataset1->get_point(i).vals[var_ds1]
529  - dataset2->get_point(i).vals[var_ds2] );
530 
531  if(dif > maxdif)
532  maxdif = dif;
533  }
534 
535  return maxdif;
536 }
537 
538 long double difmap::get_min_dif(){
539 
540  long double mindif = std::numeric_limits<long double>::infinity();
541 
542  int size = dataset1->get_num_points();
543 
544  for(int i=0; i<size; i++){
545 
546  long double dif = fabs( dataset1->get_point(i).vals[var_ds1]
547  - dataset2->get_point(i).vals[var_ds2] );
548  if(dif < mindif)
549  mindif = dif;
550  }
551 
552  return mindif;
553 }
554 
555 void difmap::output_info(string filename, long double max, long double min, int large_ds){
556  ofstream infowrite;
557  infowrite.open(filename.c_str(), ofstream::out);
558 
559  infowrite << "FILE 1 MAX: " << get_max_value(1) << endl;
560  infowrite << "FILE 1 MIN: " << get_min_value(1) << endl;
561 
562  infowrite << "FILE 2 MAX: " << get_max_value(2) << endl;
563  infowrite << "FILE 2 MIN: " << get_min_value(2) << endl;
564 
565  infowrite << "Max Difference: " << max << endl;
566  infowrite << "Min Difference: " << min << endl;
567  infowrite << "Range: " << (max-min) << endl;
568 
569  if(large_ds != -1)
570  infowrite << "Dataset with largest values: " << large_ds << endl;
571 
572  infowrite << "\n" << "Light green points: Values that could not be interpolated" << endl;
573  infowrite << "Dark green points: Divider" << endl;
574 }
575 */
576 
577 
578 /* difmap with key -------------------------------------------------------*/
579 /*
580 difmap_wkey::difmap_wkey(std::vector<indexed_points*> datasets1, indexed_points *dataset2,
581  int var_ds1, int var_ds2,
582  ofstream &outfile, string outprefix)
583  : difmap(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
584 
585  border_width = 4;
586  key_width = 20;
587 
588  border_color.r = 0;
589  border_color.g = 128;
590  border_color.b = 0;
591 
592  return;
593 }
594 
595 string difmap_wkey::get_metric_name(){
596  return string("Difference Map with Key");
597 }
598 
599 string difmap_wkey::compare(){
600 
601  indexed_points *ds1_ip = dynamic_cast<indexed_points*>(dataset1);
602  indexed_points *ds2_ip = dynamic_cast<indexed_points*>(dataset2);
603 
604  if(ds1_ip == NULL || ds2_ip == NULL)
605  return string("Improper Datasets");
606 
607  string ret("");
608 
609  if(!check_card_match()){
610  ret = "No Output";
611  return ret;
612  }
613 
614  long double max = get_max_dif();
615  long double min = get_min_dif();
616 
617  long double range = fabs( max - min );
618 
619  CImg<unsigned char> imgOut( ds1_ip->get_dim()[0] + border_width + key_width,
620  ds1_ip->get_dim()[1],
621  ds1_ip->get_dim()[2],
622  3, 0 );
623 
624  //Write difference map to image
625  for(int x=0; x<ds1_ip->get_dim()[0]; x++)
626  {
627  for(int y=0; y<ds1_ip->get_dim()[1]; y++)
628  {
629  for(int z=0; z<ds1_ip->get_dim()[2]; z++)
630  {
631 
632  //TODO: loop over every dep_var
633  std::vector<int> loc;
634  loc.push_back( x );
635  loc.push_back( y );
636  loc.push_back( z );
637 
638  pnt p1 = ds1_ip ->get_indexed_point( loc );
639  pnt p2 = ds2_ip ->get_indexed_point( loc );
640 
641  if( this ->restrictToBoxRange &&
642  this ->range ->hasPoint( p1.vals[0],p1.vals[1],0.0 ) &&
643  this ->range ->hasPoint( p2.vals[0],p2.vals[1],0.0 ) )
644  {
645  std::vector<int> loc;
646  loc.push_back(x);
647  loc.push_back(y);
648  loc.push_back(z);
649 
650  long double pt_val = fabs( ds1_ip->get_indexed_point(loc).vals[var_ds1]
651  - ds2_ip->get_indexed_point(loc).vals[var_ds2] );
652 
653  rgb_color col = get_color_double_sided(pt_val, range, min, 0);
654 
655  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 0) = col.r;
656  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 1) = col.g;
657  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 2) = col.b;
658  }
659  else if( ! this ->restrictToBoxRange )
660  {
661  std::vector<int> loc;
662  loc.push_back(x);
663  loc.push_back(y);
664  loc.push_back(z);
665 
666  long double pt_val = fabs( ds1_ip->get_indexed_point(loc).vals[var_ds1]
667  - ds2_ip->get_indexed_point(loc).vals[var_ds2] );
668 
669  rgb_color col = get_color_double_sided(pt_val, range, min, 0);
670 
671  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 0) = col.r;
672  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 1) = col.g;
673  imgOut(x, ds1_ip->get_dim()[1] - y - 1, z, 2) = col.b;
674  }
675 
676  }
677  }
678  }
679 
680 
681  //Determine if one dataset is always larger
682  int large_ds;
683  if( 0 >= min && 0 <= max )
684  large_ds = -1;
685  else{
686 
687  std::vector<int> loc;
688  loc.push_back(0);
689  loc.push_back(0);
690  loc.push_back(0);
691 
692  if(ds1_ip->get_indexed_point(loc).vals[var_ds1]
693  > ds2_ip->get_indexed_point(loc).vals[var_ds2]){
694 
695  large_ds = 1;
696  }
697  else
698  large_ds = 2;
699  }
700 
701  //Write border
702  int border_start = ds1_ip->get_dim()[0];
703 
704  for(int y=0; y < ds1_ip->get_dim()[1]; y++){
705  for(int x=border_start; x < border_start + border_width; x++){
706  for(int z=0; z<ds1_ip->get_dim()[2]; z++){
707 
708  imgOut( x, y, z, 0 ) = border_color.r;
709  imgOut( x, y, z, 1 ) = border_color.g;
710  imgOut( x, y, z, 2 ) = border_color.b;
711  }
712  }
713  }
714 
715  //Write map key
716  int x_start = ds1_ip->get_dim()[0] + border_width;
717 
718  long double k_height = ds1_ip->get_dim()[1];
719  long double unit = range / (k_height-1.0);
720 
721  for(int y=0; y<ds1_ip->get_dim()[1]; y++){
722  for(int x=x_start; x < x_start + key_width; x++){
723  for(int z=0; z<ds1_ip->get_dim()[2]; z++){
724  rgb_color col = get_color_double_sided(min+unit*y, range, min, 0);
725 
726  imgOut( x, ds1_ip->get_dim()[1] - y - 1, z, 0 ) = col.r;
727  imgOut( x, ds1_ip->get_dim()[1] - y - 1, z, 1 ) = col.g;
728  imgOut( x, ds1_ip->get_dim()[1] - y - 1, z, 2 ) = col.b;
729  }
730  }
731  }
732 
733  //Construct file name
734  string output_name(outprefix);
735  output_name += "_difmap_wkey_";
736  output_name += "var1_";
737  output_name += itoa(var_ds1);
738  output_name += "_var2_";
739  output_name += itoa(var_ds2);
740 
741  //Write info to txt file
742  output_info(output_name + ".txt", max, min, large_ds);
743 
744  //Write image
745  imgOut.save_bmp((output_name + ".bmp").c_str());
746 
747  ret = output_name;
748  return ret;
749 }
750 */
751 
752 
753 /* color map -------------------------------------------------------------*/
754 /*
755 colmap::colmap(std::vector<indexed_points*> datasets1, indexed_points *dataset2,
756  int var_ds1, int var_ds2,
757  ofstream &outfile, string outprefix)
758  : map(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
759 
760  return;
761 }
762 
763 string colmap::get_res(){
764  return compare();
765 }
766 
767 string colmap::compare(){
768 
769  string ret("");
770 
771  indexed_points *ds1_ip = dynamic_cast<indexed_points*>(dataset1);
772  indexed_points *ds2_ip = dynamic_cast<indexed_points*>(dataset2);
773 
774  if(ds1_ip == NULL || ds2_ip == NULL)
775  return string("Improper Datasets");
776 
777  for(int j=1; j<=2; j++){
778 
779  int var_num = -1;
780  int indep_vars = 0;
781  if(j == 1){
782  var_num = var_ds1;
783  indep_vars = ds1_ip->get_num_indep_vars();
784  }
785  else if(j == 2){
786  var_num = var_ds2;
787  indep_vars = ds2_ip->get_num_indep_vars();
788  }
789 
790  //Construct file name
791  string output_name(outprefix);
792  output_name += "_colmap";
793 
794  if(j == 1){
795  output_name += "_FIRST_";
796  output_name += "var1_";
797  output_name += itoa(var_ds1);
798  }
799  else if(j == 2){
800  output_name += "_SECOND_";
801  output_name += "var2_";
802  output_name += itoa(var_ds2);
803  }
804 
805  output_name += ".dat";
806 
807  ret += output_name + " ";
808 
809  //Write to tecplot file
810  TecplotOrderedWriter t_writer;
811  t_writer.setOutputFilename(output_name);
812 
813  if(!t_writer.init())
814  return string("Output file could not be opened");
815 
816  //Create variable names
817  std::vector<string> variables;
818 
819  //Construct variable names
820  for(int indep_var_num=1; indep_var_num <= indep_vars; indep_var_num++){
821  variables.push_back(string(itoa(indep_var_num)));
822  }
823  variables.push_back("Vals");
824 
825  if(j == 1)
826  t_writer.addOrderedPartition(string("Partition 1"), ds1_ip, var_ds1);
827  else if(j == 2)
828  t_writer.addOrderedPartition(string("Partition 1"), ds2_ip, var_ds2);
829 
830  t_writer.writeFile(string("Comparison"), variables);
831 
832  }
833 
834  return ret;
835 
836 }
837 
838 string colmap::get_metric_name(){
839  return string("Color Map");
840 }
841 
842 long double colmap::get_max_value(int ds, int var){
843 
844  long double max = -std::numeric_limits<long double>::infinity();
845 
846  int size = dataset1->get_num_points();
847 
848  for(int i=0; i<size; i++){
849 
850  if(ds == 1){
851  if(dataset1->get_point(i).vals[var] > max){
852  max = dataset1->get_point(i).vals[var];
853  }
854  }
855 
856  if(ds == 2){
857  if(dataset2->get_point(i).vals[var] > max){
858  max = dataset2->get_point(i).vals[var];
859  }
860  }
861  }
862 
863  return max;
864 }
865 
866 long double colmap::get_min_value(int ds, int var){
867 
868  long double min = std::numeric_limits<long double>::infinity();
869 
870  int size = dataset1->get_num_points();
871 
872  for(int i=0; i<size; i++){
873 
874  if(ds == 1){
875  if(dataset1->get_point(i).vals[var] < min){
876  min = dataset1->get_point(i).vals[var];
877  }
878  }
879 
880  if(ds == 2){
881  if(dataset2->get_point(i).vals[var] < min){
882  min = dataset2->get_point(i).vals[var];
883  }
884  }
885 
886  }
887 
888  return min;
889 }
890 */
891 
892 
893 /* sample correlation coefficient ----------------------------------------*/
894 /*
895 scorco::scorco(std::vector<points*> datasets1, points *dataset2, int var_ds1, int var_ds2,
896  ofstream &outfile, string outprefix)
897  : metric(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
898 
899  return;
900 }
901 
902 string scorco::get_res(){
903  string ret;
904  std::stringstream out;
905  out.precision(10);
906  out << compare();
907  return out.str();
908 }
909 
910 long double scorco::compare(){
911 
912  //Check that number of points matches
913  if(!check_card_match())
914  return -std::numeric_limits<long double>::infinity();
915 
916  //Get average of each dataset
917  long double f1_ave = get_average_f1();
918  long double f2_ave = get_average_f2();
919 
920  long double sum_f1 = 0;
921  long double sum_sq_f1 = 0;
922 
923  long double sum_f2 = 0;
924  long double sum_sq_f2 = 0;
925 
926  long double sum_mult = 0;
927 
928  for(int i=0; i<dataset1->get_num_points(); i++)
929  {
930 
931  pnt p1 = dataset1 ->get_point( i );
932  pnt p2 = dataset2 ->get_point( i );
933 
934  if( this ->restrictToBoxRange &&
935  this ->range ->hasPoint( p1.vals[0],p1.vals[1],0.0 ) &&
936  this ->range ->hasPoint( p2.vals[0],p2.vals[1],0.0 ) )
937  {
938  long double f1 = (long double)dataset1->get_point(i).vals[var_ds1];
939  long double f2 = (long double)dataset2->get_point(i).vals[var_ds2];
940 
941  if( !isnan(f1) && !isnan(f2) ){
942 
943  //Calculate the sum of dataset1
944  sum_f1 += f1;
945  //Calculate the square sum of dataset1
946  sum_sq_f1 += pow(f1, 2);
947 
948  //Calculate the sum of dataset2
949  sum_f2 += f2;
950  //Calculate the square sume of dataset2
951  sum_sq_f2 += pow(f2, 2);
952 
953  sum_mult += f1*f2;
954  }
955  }
956  else if( ! this ->restrictToBoxRange )
957  {
958  long double f1 = (long double)dataset1->get_point(i).vals[var_ds1];
959  long double f2 = (long double)dataset2->get_point(i).vals[var_ds2];
960 
961  if( !isnan(f1) && !isnan(f2) ){
962 
963  //Calculate the sum of dataset1
964  sum_f1 += f1;
965  //Calculate the square sum of dataset1
966  sum_sq_f1 += pow(f1, 2);
967 
968  //Calculate the sum of dataset2
969  sum_f2 += f2;
970  //Calculate the square sume of dataset2
971  sum_sq_f2 += pow(f2, 2);
972 
973  sum_mult += f1*f2;
974  }
975 
976  }
977 
978  }
979 
980  long double n = dataset1->get_num_points();
981 
982  long double res = 0;
983 
984  //Calculate sample correlation coefficient with calculated values
985  res = (n*sum_mult - (sum_f1*sum_f2)) /
986  ( sqrt(n*sum_sq_f1 - pow(sum_f1,2)) * sqrt(n*sum_sq_f2 - pow(sum_f2,2)) );
987 
988  return res;
989 
990 }
991 
992 long double scorco::get_average_f1(){
993 
994  long double total = 0;
995 
996  int num_indep = dataset1->get_num_indep_vars();
997 
998  for(int i=0; i<dataset1->get_num_points(); i++){
999  total += (long double)dataset1->get_point(i).vals[num_indep];
1000  }
1001 
1002  return (total / ((long double)dataset1->get_num_points()));
1003 }
1004 
1005 long double scorco::get_average_f2(){
1006 
1007  long double total = 0;
1008 
1009  int num_indep = dataset2->get_num_indep_vars();
1010 
1011  for(int i=0; i<dataset2->get_num_points(); i++){
1012  total += (long double)dataset2->get_point(i).vals[num_indep];
1013  }
1014 
1015  return (total / ((long double)dataset2->get_num_points()));
1016 }
1017 
1018 
1019 bool scorco::check_card_match(){
1020 
1021  if(dataset1->get_num_points() == dataset2->get_num_points())
1022  return true;
1023 
1024  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
1025  err_ps += itoa(dataset1->get_num_points());
1026  err_ps += ", Dataset 2 = ";
1027  err_ps += itoa(dataset2->get_num_points());
1028 
1029  error_out(err_ps);
1030 
1031  return false;
1032 }
1033 
1034 string scorco::get_metric_name(){
1035  return string("Sample Correlation Coefficient");
1036 }
1037 */
1038 
1039 
1040 /* modelling efficiency --------------------------------------------------*/
1041 /*
1042 modef::modef(std::vector<points*> datasets1, points *dataset2, int var_ds1, int var_ds2,
1043  ofstream &outfile, string outprefix)
1044  : metric(datasets1, dataset2, var_ds1, var_ds2, outfile, outprefix){
1045 
1046  return;
1047 }
1048 
1049 string modef::get_res(){
1050  string ret;
1051  std::stringstream out;
1052  out.precision(10);
1053  out << compare();
1054  return out.str();
1055 }
1056 
1057 long double modef::compare(){
1058 
1059  if(!check_card_match())
1060  return -std::numeric_limits<long double>::infinity();
1061 
1062  long double f1_ave = get_average_f1();
1063 
1064  long double sum_sq_dif = 0;
1065  long double sum_sq_mean_dif = 0;
1066 
1067  for(int i=0; i<dataset1->get_num_points(); i++)
1068  {
1069 
1070  pnt p1 = dataset1 ->get_point( i );
1071  pnt p2 = dataset2 ->get_point( i );
1072 
1073  if( this ->restrictToBoxRange &&
1074  this ->range ->hasPoint( p1.vals[0],p1.vals[1],0.0 ) &&
1075  this ->range ->hasPoint( p2.vals[0],p2.vals[1],0.0 ) )
1076  {
1077 
1078  long double f1 = (long double)dataset1->get_point(i).vals[var_ds1];
1079  long double f2 = (long double)dataset2->get_point(i).vals[var_ds2];
1080 
1081  if( !isnan(f1) && !isnan(f2) ){
1082  sum_sq_dif += pow(f1 - f2, 2);
1083  sum_sq_mean_dif += pow(f1 - f1_ave, 2);
1084  }
1085  }
1086  else if( ! this ->restrictToBoxRange )
1087  {
1088 
1089  long double f1 = (long double)dataset1->get_point(i).vals[var_ds1];
1090  long double f2 = (long double)dataset2->get_point(i).vals[var_ds2];
1091 
1092  if( !isnan(f1) && !isnan(f2) ){
1093  sum_sq_dif += pow(f1 - f2, 2);
1094  sum_sq_mean_dif += pow(f1 - f1_ave, 2);
1095  }
1096 
1097  }
1098 
1099 
1100  }
1101 
1102  long double res = 0;
1103 
1104  res = 1 - ((sum_sq_dif) / (sum_sq_mean_dif));
1105 
1106  return res;
1107 
1108 }
1109 
1110 long double modef::get_average_f1(){
1111 
1112  long double total = 0;
1113 
1114  int num_indep = dataset1->get_num_indep_vars();
1115 
1116  for(int i=0; i<dataset1->get_num_points(); i++){
1117  total += (long double)dataset1->get_point(i).vals[num_indep];
1118  }
1119 
1120  return (total / ((long double)dataset1->get_num_points()));
1121 }
1122 
1123 long double modef::get_average_f2(){
1124 
1125  long double total = 0;
1126 
1127  int num_indep = dataset2->get_num_indep_vars();
1128 
1129  for(int i=0; i<dataset2->get_num_points(); i++){
1130  total += (long double)dataset2->get_point(i).vals[num_indep];
1131  }
1132 
1133  return (total / ((long double)dataset2->get_num_points()));
1134 }
1135 
1136 
1137 bool modef::check_card_match(){
1138 
1139  if(dataset1->get_num_points() == dataset2->get_num_points())
1140  return true;
1141 
1142  string err_ps("Mismatched number of coordinates: Dataset 1 = ");
1143  err_ps += itoa(dataset1->get_num_points());
1144  err_ps += ", Dataset 2 = ";
1145  err_ps += itoa(dataset2->get_num_points());
1146 
1147  error_out(err_ps);
1148 
1149  return false;
1150 }
1151 
1152 string modef::get_metric_name(){
1153  return string("Modelling Efficiency");
1154 }
1155 */
1156 
1157 
1158 
1159 /* metric selection ------------------------------------------------------*/
1160 
1161 void print_metrics(){
1162 
1163  int i=0;
1164 
1165  ofstream unset;
1166  string nullstr("");
1167 /*
1168  mse one(NULL, NULL, -1, -1, unset, nullstr);
1169  cout<<"\t"<<++i<<": "<<one.get_metric_name()<<"\n";
1170 
1171  rmse two(NULL, NULL, -1, -1, unset, nullstr);
1172  cout<<"\t"<<++i<<": "<<two.get_metric_name()<<"\n";
1173 
1174  scc three(NULL, NULL, -1, -1, unset, nullstr);
1175  cout<<"\t"<<++i<<": "<<three.get_metric_name()<<"\n";
1176 
1177  difmap four(NULL, NULL, -1, -1, unset, nullstr);
1178  cout<<"\t"<<++i<<": "<<four.get_metric_name()<<"\n";
1179 
1180  difmap_wkey five(NULL, NULL, -1, -1, unset, nullstr);
1181  cout<<"\t"<<++i<<": "<<five.get_metric_name()<<"\n";
1182 
1183  colmap six(NULL, NULL, -1, -1, unset, nullstr);
1184  cout<<"\t"<<++i<<": "<<six.get_metric_name()<<"\n";
1185 
1186  scorco seven(NULL, NULL, -1, -1, unset, nullstr);
1187  cout<<"\t"<<++i<<": "<<seven.get_metric_name()<<"\n";
1188 
1189  modef eight(NULL, NULL, -1, -1, unset, nullstr);
1190  cout<<"\t"<<++i<<": "<<eight.get_metric_name()<<"\n";
1191 */
1192 }
1193 
1194 metric *metric_select(int m, std::vector<points*> dss1, std::vector<points*> dss2,
1195  int var_ds1, int var_ds2,
1196  ofstream &outfile, string outprefix){
1197 
1198  metric *ret;
1199 
1200  //Indexed points conversion
1201  // - For metrics that require indexed_points
1202  //indexed_points *ds1_ip = dynamic_cast<indexed_points*>(ds1);
1203  //indexed_points *ds2_ip = dynamic_cast<indexed_points*>(ds2);
1204 
1205  switch (m){
1206 /*
1207  case 1:
1208  ret = new mse(ds1, ds2, var_ds1, var_ds2, outfile, outprefix);
1209  break;
1210 
1211  case 2:
1212  ret = new rmse(ds1, ds2, var_ds1, var_ds2, outfile, outprefix);
1213  break;
1214 
1215  case 3:
1216  ret = new scc(ds1, ds2, var_ds1, var_ds2, outfile, outprefix);
1217  break;
1218 
1219  case 4:
1220  ret = new difmap(ds1_ip, ds2_ip, var_ds1, var_ds2, outfile, outprefix);
1221  break;
1222 
1223  case 5:
1224  ret = new difmap_wkey(ds1_ip, ds2_ip, var_ds1, var_ds2, outfile, outprefix);
1225  break;
1226 
1227  case 6:
1228  ret = new colmap(ds1_ip, ds2_ip, var_ds1, var_ds2, outfile, outprefix);
1229  break;
1230 
1231  case 7:
1232  ret = new scorco(ds1, ds2, var_ds1, var_ds2, outfile, outprefix);
1233  break;
1234 
1235  case 8:
1236  ret = new modef(ds1, ds2, var_ds1, var_ds2, outfile, outprefix);
1237  break;
1238 */
1239  default:
1240  ret = NULL;
1241  break;
1242 
1243  }
1244 
1245  return ret;
1246 
1247 }
1248 
metric::set_datasets1
void set_datasets1(std::vector< points * > datasets1)
Definition: metric.cpp:22
metric::get_res
virtual string get_res()
Definition: metric.cpp:30
metric::set_datasets2
void set_datasets2(std::vector< points * > datasets2)
Definition: metric.cpp:26
metric::datasets1
std::vector< points * > datasets1
Definition: metric.h:68
metric::metric
metric(std::vector< points * > datasets1, std::vector< points * > datasets2, int var_ds1, int var_ds2, ofstream &outfile, string outprefix)
Definition: metric.cpp:6
print_metrics
void print_metrics()
Sample Correlation Coefficient.
Definition: metric.cpp:1161
metric::var_ds2
int var_ds2
Definition: metric.h:72
metric::get_metric_name
virtual string get_metric_name()
Definition: metric.cpp:34
i
blockLoc i
Definition: read.cpp:79
metric::var_ds1
int var_ds1
Definition: metric.h:71
metric_select
metric * metric_select(int m, std::vector< points * > dss1, std::vector< points * > dss2, int var_ds1, int var_ds2, ofstream &outfile, string outprefix)
Definition: metric.cpp:1194
metric::error_out
void error_out(string err_ps)
Definition: metric.cpp:38
metric::outprefix
string outprefix
Definition: metric.h:66
metric::outfile
ofstream & outfile
Definition: metric.h:74
metric::datasets2
std::vector< points * > datasets2
Definition: metric.h:69
metric
Definition: metric.h:35