LCOV - code coverage report
Current view: top level - src - pdt_plan.cxx (source / functions) Hit Total Coverage
Test: code_coverage_filter.info Lines: 526 1159 45.4 %
Date: 2024-03-28 16:04:17 Functions: 40 96 41.7 %
Branches: 565 1912 29.6 %

           Branch data     Line data    Source code
       1                 :            : /** 
       2                 :            :  * \file pdt_plan.cxx \brief PDT Plan Implementation
       3                 :            :  *
       4                 :            :  * This source files implements the different methods for the pdt_plan class
       5                 :            :  *
       6                 :            :  *  @author: Abdul-Amir Yassine
       7                 :            :  *  @date: March 6th, 2018
       8                 :            :  */
       9                 :            : 
      10                 :            : #include "pdt_plan.h"
      11                 :            : #include "read_mesh.h"
      12                 :            : #include "src_placer.h"
      13                 :            : #include "fixed_sources_convex_program.h"
      14                 :            : #include "run_fullmonte.h"
      15                 :            : #include "rt_feedback.h"
      16                 :            : #include "sa_engine.h"
      17                 :            : #include "file_parser.h"
      18                 :            : #include "light_source.h"
      19                 :            : 
      20                 :            : //! Default Constructor
      21         [ #  # ]:          0 : pdt_plan::pdt_plan()
      22                 :            : {
      23                 :          0 :     _placer = nullptr;
      24                 :          0 :     _data = nullptr;
      25                 :          0 :     _my_lp = nullptr;
      26                 :            : 
      27                 :          0 :     _rt_feedback_system = nullptr;
      28                 :            :     
      29         [ #  # ]:          0 :     _wavelength = "";
      30         [ #  # ]:          0 :     _source_type = "";
      31         [ #  # ]:          0 :     _placement_type = "";
      32                 :          0 :     _tumor_weight = 0; 
      33                 :            : 
      34                 :          0 :     _source_tailored = false;
      35                 :            : 
      36                 :          0 :     _num_packets = 0;
      37                 :          0 :     _total_energy = 0;
      38                 :            : 
      39                 :          0 :     _tumor_region_number.clear();
      40                 :          0 :     _mesh_size = 0;
      41                 :          0 :     _num_sources = 0;
      42                 :            : 
      43                 :          0 :     _tumor_volume = 0;
      44                 :          0 :     _total_mesh_volume = 0;
      45                 :          0 :     _max_region_index = 0;
      46                 :            : 
      47                 :          0 :     _ignore_unnecessary_tetra = false;
      48                 :          0 :     _actual_mesh_size = 0; 
      49                 :            : 
      50         [ #  # ]:          0 :     _params_file = "";
      51                 :          0 :     _read_data_from_vtk = false;
      52                 :            : 
      53                 :          0 :     _guardband_thresholds_called = false;
      54                 :          0 :     _is_p_max_set = false;
      55                 :          0 :     _is_unnecessary_removed = false;
      56                 :          0 :     _is_tetra_data_called = false; 
      57                 :          0 :     _is_lp_solved = false;
      58                 :            : 
      59                 :          0 :     _num_non_zero_sources = 0;
      60                 :          0 :     _total_power = 0;
      61                 :            : 
      62                 :          0 :     _vary_tumor_weight = false;
      63                 :            : 
      64                 :          0 :     _beta_cvar = 0.0;
      65                 :          0 :     _op_cvar_approach = false;
      66                 :          0 :     _op_variance_constraints = false;
      67                 :            : 
      68                 :          0 :     _fullmonte_engine = nullptr;
      69                 :            : 
      70         [ #  # ]:          0 :         _sa_engine_progress_placement_file = "";
      71                 :          0 :     _sa_engine_beta = 0.02; // Default value 
      72                 :          0 :         _sa_engine_lambda = 0.3; // Default value
      73                 :          0 :     _sa_engine_moves_per_temp = 0.0;
      74                 :          0 :         _run_sa = false;
      75                 :          0 :     _sa_engine_early_terminate = false;
      76                 :            : 
      77         [ #  # ]:          0 :     _final_distribution_output_file = "";
      78                 :            : 
      79                 :          0 :     _run_rt_feedback = false; 
      80                 :            : 
      81         [ #  # ]:          0 :     _final_distribution_output_file = "";
      82                 :            : 
      83                 :          0 :     _run_rt_feedback = false; 
      84                 :            : 
      85                 :          0 :         _target_tumor_v100 = 0.0;
      86                 :          0 :         _upper_tumor_v100_bound = 0.0;
      87                 :          0 :         _lower_tumor_v100_bound = 0.0; 
      88                 :            : 
      89                 :          0 :         _max_num_lp_iterations = std::numeric_limits<int>::max(); 
      90                 :            : 
      91                 :          0 :         _lp_constraints_defined = false;
      92                 :          0 :         _lp_objective_defined = false; 
      93                 :            :     
      94                 :          0 :     _rt_feedback_rand_seed = 1;
      95                 :          0 :     _rt_feedback_detector_radius = 0.5f; 
      96                 :          0 :     _rt_feedback_detector_na = 0.5f; 
      97                 :          0 :     _rt_feedback_lut_size = 50;
      98         [ #  # ]:          0 :     _rt_feedback_lut_file_name = "";
      99                 :            : 
     100                 :          0 :     _power_variance_optimization = false;
     101                 :          0 :     _sigma_power = 0.0;
     102                 :          0 : }
     103                 :            : 
     104                 :            : //!
     105                 :            : //! Regular Constructor
     106                 :            : //! @param params_file: [in] the path to the parameters file
     107                 :            : //!
     108         [ +  - ]:         72 : pdt_plan::pdt_plan(string params_file)
     109                 :            : {
     110                 :         36 :     _sa_engine_beta = 0.02; // default value in case it is not supplied
     111                 :         36 :         _sa_engine_lambda = 0.3; // default value in case it is not supplied
     112                 :         36 :     _sa_engine_moves_per_temp = 0.0;
     113                 :         36 :         _run_sa = false;
     114                 :            : 
     115                 :         36 :     _run_rt_feedback = false; 
     116                 :            : 
     117                 :         36 :     _tumor_volume = 0;
     118                 :         36 :     _total_mesh_volume = 0;
     119                 :         36 :     _max_region_index = 0;
     120                 :            : 
     121 [ +  - ][ +  - ]:         36 :     _placer = new src_placer();
     122                 :         36 :     _my_lp = nullptr;
     123                 :            : 
     124                 :         36 :     _rt_feedback_system = nullptr;
     125                 :            : 
     126                 :         36 :     _source_tailored = false;
     127                 :         36 :     _sa_engine_early_terminate = false;
     128                 :            : 
     129                 :            :     // Reading the parameters file
     130         [ +  - ]:         36 :     _params_file = params_file;
     131                 :            : 
     132         [ +  - ]:         72 :     string file_type = _params_file.substr(_params_file.size() - 3);
     133 [ +  - ][ +  - ]:         36 :     if (file_type == "xml") {
     134                 :            :         // read xml file
     135 [ +  - ][ +  - ]:         36 :         _fparser = new Parser (_params_file);
                 [ +  - ]
     136                 :            :     } else {
     137         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmUnsupported file type: %s. Exiting...\033[0m\n", 31, file_type);
     138                 :          0 :         exit(-1);
     139                 :            :     }
     140                 :            : 
     141         [ +  - ]:         36 :     read_params_file();
     142                 :            : 
     143                 :            :     // in read_params_file
     144         [ +  - ]:         36 :     _placer->set_num_packets(_num_packets);
     145         [ +  - ]:         36 :     _placer->set_total_energy(_total_energy);
     146                 :            : 
     147         [ +  - ]:         36 :     _fparser->get_injection_point_options(_constrain_injection_point/*, _placement_mode*/);
     148         [ +  - ]:         36 :     _num_injection_points = _fparser->get_num_injection_points();
     149                 :            : 
     150                 :            :     // Append project directory path to files specified in input parameters for testing
     151         [ +  - ]:         36 :     if (_file_path_override) {
     152         [ +  - ]:         36 :         override_file_paths();
     153                 :            :     }
     154                 :            : 
     155                 :            :     // Read the mesh 
     156         [ +  - ]:         36 :     read_mesh(); 
     157                 :            :     
     158                 :            :     // Read tissue properties to set _d_max and _d_min vectors
     159         [ +  - ]:         36 :     read_tissue_properties();
     160                 :            : 
     161                 :            :     // Read the fluence matrix from FullMonte output vtk
     162 [ +  + ][ +  - ]:         36 :     if (_running_tests && !_run_sa) {
     163         [ +  - ]:          5 :                 read_fullmonte_output_matrix(); 
     164                 :            :         }
     165                 :            : 
     166 [ +  - ][ +  - ]:         36 :     _fullmonte_engine = new WrapperFullMonteSW(_mesh_file, _fparser, _read_data_from_vtk, _use_cuda); 
                 [ +  - ]
     167                 :            : 
     168                 :            :     // ignoring unnecessary tetra
     169                 :         36 :     _ignore_unnecessary_tetra = true;
     170                 :         36 :     _actual_mesh_size = _mesh_size;
     171         [ +  - ]:         36 :     _tetra_to_consider.resize(_mesh_size);
     172         [ +  - ]:         36 :     fill(_tetra_to_consider.begin(), _tetra_to_consider.end(), true);
     173                 :            : 
     174                 :         36 :     _guardband_thresholds_called = false;
     175                 :         36 :     _is_p_max_set = false;
     176                 :         36 :     _is_unnecessary_removed = false;
     177                 :         36 :     _is_tetra_data_called = false; 
     178                 :         36 :     _is_lp_solved = false;
     179                 :            : 
     180                 :         36 :     _num_non_zero_sources = 0;
     181                 :         36 :     _total_power = 0;
     182                 :            : 
     183                 :         36 :         _lower_tumor_v100_bound = _target_tumor_v100 - 1e-1; // Default value
     184                 :         36 :         _upper_tumor_v100_bound = _target_tumor_v100 + 5e-2; // Default value
     185                 :            : 
     186                 :         36 :         _max_num_lp_iterations = std::numeric_limits<int>::max(); 
     187                 :         36 :         _num_lp_iterations = 0; 
     188                 :            :     
     189                 :         36 :     _beta_cvar = 0.0;
     190                 :         36 :     _op_cvar_approach = false;
     191                 :            : 
     192                 :         36 :     _op_variance_constraints = false;
     193                 :            :         
     194                 :         36 :         _lp_constraints_defined = false;
     195                 :         36 :         _lp_objective_defined = false; 
     196                 :            :     
     197                 :         36 :     _power_variance_optimization = false;
     198                 :         36 :     _sigma_power = 0.0;
     199                 :         36 : }
     200                 :            : 
     201                 :            : //! Destructor
     202                 :         62 : pdt_plan::~pdt_plan()
     203                 :            : {
     204         [ +  - ]:         31 :     if (_placer) 
     205         [ +  - ]:         31 :         delete _placer;
     206         [ +  - ]:         31 :     if (_data)
     207         [ +  - ]:         31 :         delete _data;
     208         [ +  + ]:         31 :     if (_my_lp)
     209         [ +  - ]:          5 :         delete _my_lp;
     210         [ +  - ]:         31 :     if (_fullmonte_engine)
     211         [ +  - ]:         31 :         delete _fullmonte_engine;
     212         [ -  + ]:         31 :     if (_rt_feedback_system)
     213         [ #  # ]:          0 :         delete _rt_feedback_system; 
     214         [ +  - ]:         31 :     if (_fparser)
     215         [ +  - ]:         31 :         delete _fparser;
     216                 :            : 
     217         [ -  + ]:         31 :     for (unsigned int i = 0; i < _placers.size(); i++)
     218         [ #  # ]:          0 :         if (_placers[i])
     219         [ #  # ]:          0 :             delete _placers[i];
     220                 :         31 :     _placers.clear();
     221                 :         31 : }
     222                 :            : 
     223                 :            : // General methods
     224                 :            : //! This method reads the mesh and stores it inside the _data pointer
     225                 :         58 : void pdt_plan::read_mesh()
     226                 :            : {
     227         [ +  - ]:         58 :     fprintf(stderr, "\n\033[1;%dmReading Mesh %s.\033[0m\n", 35, _mesh_file.c_str());
     228                 :            : 
     229 [ +  - ][ +  - ]:         58 :     _data = new mesh(_mesh_file, _fparser, _read_data_from_vtk);
                 [ +  - ]
     230                 :            : 
     231         [ +  - ]:         58 :     _fparser->get_tumor_regions(_tumor_region_number);
     232                 :            : 
     233 [ +  - ][ +  - ]:         58 :     cout << "Number of Different Materials is: " << _data->get_number_of_materials() << endl;
         [ +  - ][ +  - ]
     234                 :            : 
     235 [ +  - ][ +  - ]:         58 :     cout << "Number of points is: " << _data->get_number_of_points() << endl;
         [ +  - ][ +  - ]
     236                 :            : 
     237 [ +  - ][ +  - ]:         58 :     cout << "Number of Tetrahedra is: " << _data->get_mesh_size() << endl;
         [ +  - ][ +  - ]
     238                 :            : 
     239                 :            :     // Finding max and min coords
     240                 :         58 :     float max_x = _NEG_INFINITY_;
     241                 :         58 :     float max_y = _NEG_INFINITY_;
     242                 :         58 :     float max_z = _NEG_INFINITY_;
     243                 :         58 :     float min_x = _INFINITY_;
     244                 :         58 :     float min_y = _INFINITY_;
     245                 :         58 :     float min_z = _INFINITY_;
     246 [ +  - ][ +  + ]:    4073224 :     for (unsigned int i = 0; i < _data->get_list_of_points().size(); i++)
     247                 :            :     {
     248 [ +  - ][ +  - ]:    4073166 :         min_x = min(min_x, _data->get_list_of_points()[i]->get_xcoord());
     249 [ +  - ][ +  - ]:    4073166 :         min_y = min(min_y, _data->get_list_of_points()[i]->get_ycoord());
     250 [ +  - ][ +  - ]:    4073166 :         min_z = min(min_z, _data->get_list_of_points()[i]->get_zcoord());
     251 [ +  - ][ +  - ]:    4073166 :         max_x = max(max_x, _data->get_list_of_points()[i]->get_xcoord());
     252 [ +  - ][ +  - ]:    4073166 :         max_y = max(max_y, _data->get_list_of_points()[i]->get_ycoord());
     253 [ +  - ][ +  - ]:    4073166 :         max_z = max(max_z, _data->get_list_of_points()[i]->get_zcoord());
     254                 :            :     }
     255 [ +  - ][ +  - ]:         58 :     cout << "max_x: " << max_x << " min_x: " << min_x << endl;
         [ +  - ][ +  - ]
                 [ +  - ]
     256 [ +  - ][ +  - ]:         58 :     cout << "max_y: " << max_y << " min_y: " << min_y << endl;
         [ +  - ][ +  - ]
                 [ +  - ]
     257 [ +  - ][ +  - ]:         58 :     cout << "max_z: " << max_z << " min_z: " << min_z << endl;
         [ +  - ][ +  - ]
                 [ +  - ]
     258                 :            : 
     259                 :            : 
     260                 :            :     /**************************************************************************************/
     261                 :            :     
     262         [ +  - ]:         58 :     _mesh_size = _data->get_mesh_size();
     263                 :            :   
     264                 :            :     //cout << "Tumor Region index is: " << _tumor_region_number << endl;
     265                 :         58 : }
     266                 :            : 
     267                 :            : //! This method reads the tissue properties file to set d_max and d_min vectors
     268                 :         36 : void pdt_plan::read_tissue_properties()
     269                 :            : {
     270                 :         36 :     fprintf(stderr, "\n\033[1;%dmReading Tissue Properties.\033[0m\n", 33);
     271                 :            : 
     272                 :         36 :     _fparser->get_tissue_properties(_d_min_tissues, _d_max_tissues, _guardband_tissues, _tissue_names_ids);
     273                 :            : 
     274                 :         36 :     double scale_factor_packets = _num_packets/_total_energy;
     275         [ +  + ]:        216 :     for (unsigned int i = 0; i < _data->get_number_of_materials(); i++)
     276                 :            :     {
     277         [ +  + ]:        180 :         if (_d_max_tissues[i] != _INFINITY_)
     278                 :        144 :             _d_max_tissues[i] *= scale_factor_packets;
     279                 :        180 :         _d_min_tissues[i] *= scale_factor_packets;
     280                 :        180 :         cout << i << " dmax " << _d_max_tissues[i] << " dmin " << _d_min_tissues[i] << " guardband " << _guardband_tissues[i] << endl;
     281                 :            :     } 
     282                 :            : 
     283                 :         36 : }
     284                 :            : 
     285                 :            : //! This method reads the Fluence matrix G from FullMonte output VTK file
     286                 :          5 : void pdt_plan::read_fullmonte_output_matrix()
     287                 :            : {
     288 [ +  - ][ +  - ]:         10 :     profiler fullmonte_read("Reading Fluence Matrix");
     289                 :            : 
     290         [ +  - ]:         10 :     stringstream tcl_path;
     291         [ +  - ]:          5 :     tcl_path << _fullmonte_output_file; 
     292                 :            : 
     293         [ +  - ]:         10 :     string vtk_file_path = tcl_path.str();
     294                 :            : 
     295                 :          5 :     Util::BEGIN_TIME = chrono::steady_clock::now();
     296 [ +  - ][ +  - ]:          5 :     _placer->read_vtk_file(vtk_file_path, _mesh_size);
     297                 :          5 :     Util::END_TIME = chrono::steady_clock::now();
     298                 :            : 
     299                 :            :     double vtk_read_time = (double)chrono::duration_cast<chrono::seconds>
     300                 :          5 :                     (Util::END_TIME - Util::BEGIN_TIME).count(); 
     301                 :            : 
     302 [ +  - ][ +  - ]:          5 :     cout << "Time to read VTK file is: " << vtk_read_time << " secs" << endl;
         [ +  - ][ +  - ]
     303                 :            :     
     304         [ +  - ]:          5 :     _num_sources = (unsigned int)_placer->get_source_element_fluence_matrix().size();
     305                 :            : 
     306 [ +  - ][ +  - ]:          5 :     cout << "Number of candidate sources is: " << _num_sources << endl;
                 [ +  - ]
     307                 :          5 : }
     308                 :            : 
     309                 :            : //! Helper method: generates the initial solution
     310                 :          0 : void pdt_plan::read_initial_placement()
     311                 :            : {
     312                 :            :     // read file
     313                 :          0 :     _fparser->get_sources(_initial_placement);
     314                 :            : 
     315         [ #  # ]:          0 :     for (unsigned i = 0; i < _placement_results.size(); i++) {
     316         [ #  # ]:          0 :         for (unsigned j = 0; j < _placement_results[i].size(); j++) {
     317 [ #  # ][ #  # ]:          0 :             if (_placement_results[i][j]) delete _placement_results[i][j];
     318                 :            :         }
     319                 :            :     }
     320                 :          0 :     _placement_results.clear();
     321                 :          0 :     _placement_results.push_back(_initial_placement);
     322                 :          0 : }
     323                 :            : 
     324                 :            : //! This method returns the thresholds of the tissues
     325                 :          0 : vector<double> pdt_plan::get_thresholds() {
     326         [ #  # ]:          0 :     vector<double> thresholds(_d_max_tissues.size());
     327                 :            : 
     328         [ #  # ]:          0 :     for (unsigned i = 0; i < thresholds.size(); i++) {
     329 [ #  # ][ #  # ]:          0 :         if (_tumor_region_number.find(i+1) == _tumor_region_number.end()) {
     330                 :          0 :             thresholds[i] = _d_max_tissues[i];
     331                 :            :         } else {
     332                 :          0 :             thresholds[i] = _d_min_tissues[i];
     333                 :            :         }
     334                 :            :     }
     335                 :            : 
     336                 :          0 :     return thresholds;
     337                 :            : }
     338                 :            : 
     339                 :            : //! This method returns the thresholds of the tissues multiplied by the specified guardband
     340                 :         10 : vector<double> pdt_plan::get_thresholds_with_guardband()
     341                 :            : {
     342         [ +  - ]:         10 :     vector<double> thresholds(_d_max_tissues.size());
     343                 :            : 
     344 [ +  - ][ +  - ]:         10 :     _slopes.resize(_data->get_number_of_materials()); // this is used in weighting the tetras 
     345                 :            :                                 // slopes are inversely proportional to the thresholds
     346                 :            : 
     347                 :         10 :     double first_threshold = -1; // used in setting the slopes
     348                 :            : 
     349         [ +  + ]:         60 :     for (unsigned i = 0; i < thresholds.size(); i++) {
     350 [ +  - ][ +  + ]:         50 :         if (_tumor_region_number.find(i+1) == _tumor_region_number.end()) {
     351                 :         40 :             thresholds[i] = _guardband_tissues[i]*_d_max_tissues[i];
     352                 :            :         } else {
     353                 :         10 :             thresholds[i] = _guardband_tissues[i]*_d_min_tissues[i];
     354                 :            :         }
     355         [ +  + ]:         50 :         if (i == 0) {
     356                 :         10 :             first_threshold = thresholds[i];
     357                 :         10 :             _slopes[i] = 1.0;
     358                 :            :         } else {
     359                 :         40 :             _slopes[i] = (first_threshold)/thresholds[i];
     360                 :            :         }
     361                 :         50 :         _slopes[i] *= _num_packets/_total_energy;
     362                 :            :     }
     363                 :            : 
     364                 :         10 :     _guardband_thresholds_called = true;
     365                 :            :     
     366                 :         10 :     return thresholds;
     367                 :            : }
     368                 :            : 
     369                 :            : //! This method removes the unnecessary tetra that do not affect the optimizations based on the thresholds given
     370                 :         16 : void pdt_plan::remove_unnecessary_tetra(const vector<double> & guardbanded_thresholds)
     371                 :            : {
     372                 :         16 :         _d_max.clear(); 
     373                 :         16 :         _d_min.clear(); 
     374                 :         16 :         _tetra_volumes.clear(); 
     375                 :         16 :         _tetra_weights.clear(); 
     376                 :         16 :         _tet_material_ids.clear(); 
     377                 :         16 :         _tumor_indices.clear(); 
     378                 :         16 :         _tumor_indices.resize(0); 
     379                 :            : 
     380         [ +  - ]:         16 :     if (_ignore_unnecessary_tetra)
     381                 :         16 :         _placer->get_ignored_tetras(_tumor_region_number, _data, guardbanded_thresholds, _actual_mesh_size,
     382                 :         16 :              _tetra_to_consider);
     383                 :            : 
     384                 :         16 :     cout << "Number of Tetra to consider: " << _actual_mesh_size << endl;
     385                 :            : 
     386                 :         16 :     _d_max.resize(_actual_mesh_size);
     387                 :         16 :     _d_min.resize(_actual_mesh_size);
     388                 :         16 :     _tetra_volumes.resize(_actual_mesh_size);
     389                 :         16 :     _tetra_weights.resize(_actual_mesh_size);
     390                 :         16 :     _tet_material_ids.resize(_actual_mesh_size);
     391                 :            : 
     392                 :         16 :     _is_unnecessary_removed = true;
     393                 :            : 
     394         [ -  + ]:         16 :     if (_op_variance_constraints)
     395                 :          0 :         compute_tetra_dose_variance();
     396                 :         16 : }
     397                 :            : 
     398                 :            : //!
     399                 :            : //! This method computes the data of all tetra. It updates the following vectors: _original_tet_material_ids, 
     400                 :            : //! _original_tetra_volumes, _original_tumor_indices, _d_max_orig, _d_min_orig
     401                 :            : //!
     402                 :         10 : void pdt_plan::compute_original_tetra_data(const vector<double> & thresholds)
     403                 :            : {
     404                 :            : 
     405                 :            :     tetrahedron *t;
     406                 :         10 :     int tumor_count = 0; 
     407                 :         10 :     _total_mesh_volume = 0.0; _tumor_volume = 0.0;
     408                 :         20 :     unordered_map<unsigned int, unsigned int> mat_ids; // stores the material ids
     409                 :            : 
     410                 :            :     // Initialize the size of all important vectors 
     411                 :         10 :     _original_tetra_volumes.clear();
     412         [ +  - ]:         10 :     _original_tetra_volumes.resize(_mesh_size);
     413                 :            : 
     414                 :         10 :     _original_tet_material_ids.clear();
     415         [ +  - ]:         10 :     _original_tet_material_ids.resize(_mesh_size);
     416                 :            :     
     417                 :         10 :     _d_max_orig.clear(); _d_min_orig.clear(); 
     418         [ +  - ]:         10 :     _d_max_orig.resize(_mesh_size);
     419         [ +  - ]:         10 :     _d_min_orig.resize(_mesh_size);
     420                 :            : 
     421                 :            :     // Initialize coordinates to compute center of tumor
     422                 :         10 :     float max_x = _NEG_INFINITY_, max_y = _NEG_INFINITY_, max_z = _NEG_INFINITY_, 
     423                 :         10 :             min_x = _INFINITY_, min_y = _INFINITY_, min_z = _INFINITY_;
     424                 :            : 
     425         [ +  + ]:    4233770 :     for (unsigned int i = 0; i < _mesh_size; i++)
     426                 :            :     {
     427         [ +  - ]:    4233760 :         t = _data->get_list_of_tetrahedra()[i];
     428                 :            : 
     429         [ +  - ]:    4233760 :         _original_tetra_volumes[i] = _data->compute_tet_volume(t);
     430                 :    4233760 :         _total_mesh_volume += _original_tetra_volumes[i]; 
     431                 :            : 
     432         [ +  - ]:    4233760 :         _original_tet_material_ids[i] = t->get_material_id();
     433                 :            : 
     434 [ +  - ][ +  - ]:    4233760 :         if (_material_volumes.find(t->get_material_id()) == _material_volumes.end())
                 [ +  + ]
     435 [ +  - ][ +  - ]:         60 :             _material_volumes[t->get_material_id()] = _original_tetra_volumes[i];
     436                 :            :         else
     437 [ +  - ][ +  - ]:    4233700 :             _material_volumes[t->get_material_id()] += _original_tetra_volumes[i];
     438                 :            :         
     439 [ +  - ][ +  - ]:    4233760 :         if (mat_ids.find(t->get_material_id()) == mat_ids.end())
                 [ +  + ]
     440 [ +  - ][ +  - ]:         60 :             mat_ids[t->get_material_id()] = 1;
     441                 :            :         else
     442 [ +  - ][ +  - ]:    4233700 :             mat_ids[t->get_material_id()] += 1;
     443                 :            : 
     444 [ +  - ][ +  - ]:    4233760 :         if (_tumor_region_number.find((unsigned) t->get_material_id()) != _tumor_region_number.end()) 
                 [ +  + ]
     445                 :            :         {
     446                 :      58890 :             tumor_count++;
     447 [ +  - ][ +  - ]:      58890 :             _original_tumor_indices.push_back(t->get_id());
     448                 :            : 
     449 [ +  - ][ +  - ]:      58890 :             max_x = max(max_x, _data->get_list_of_points()[t->get_p1_id()-1]->get_xcoord());
                 [ +  - ]
     450 [ +  - ][ +  - ]:      58890 :             max_x = max(max_x, _data->get_list_of_points()[t->get_p2_id()-1]->get_xcoord());
                 [ +  - ]
     451 [ +  - ][ +  - ]:      58890 :             max_x = max(max_x, _data->get_list_of_points()[t->get_p3_id()-1]->get_xcoord());
                 [ +  - ]
     452 [ +  - ][ +  - ]:      58890 :             max_x = max(max_x, _data->get_list_of_points()[t->get_p4_id()-1]->get_xcoord());
                 [ +  - ]
     453 [ +  - ][ +  - ]:      58890 :             min_x = min(min_x, _data->get_list_of_points()[t->get_p1_id()-1]->get_xcoord());
                 [ +  - ]
     454 [ +  - ][ +  - ]:      58890 :             min_x = min(min_x, _data->get_list_of_points()[t->get_p2_id()-1]->get_xcoord());
                 [ +  - ]
     455 [ +  - ][ +  - ]:      58890 :             min_x = min(min_x, _data->get_list_of_points()[t->get_p3_id()-1]->get_xcoord());
                 [ +  - ]
     456 [ +  - ][ +  - ]:      58890 :             min_x = min(min_x, _data->get_list_of_points()[t->get_p4_id()-1]->get_xcoord());
                 [ +  - ]
     457                 :            : 
     458 [ +  - ][ +  - ]:      58890 :             max_y = max(max_y, _data->get_list_of_points()[t->get_p1_id()-1]->get_ycoord());
                 [ +  - ]
     459 [ +  - ][ +  - ]:      58890 :             max_y = max(max_y, _data->get_list_of_points()[t->get_p2_id()-1]->get_ycoord());
                 [ +  - ]
     460 [ +  - ][ +  - ]:      58890 :             max_y = max(max_y, _data->get_list_of_points()[t->get_p3_id()-1]->get_ycoord());
                 [ +  - ]
     461 [ +  - ][ +  - ]:      58890 :             max_y = max(max_y, _data->get_list_of_points()[t->get_p4_id()-1]->get_ycoord());
                 [ +  - ]
     462 [ +  - ][ +  - ]:      58890 :             min_y = min(min_y, _data->get_list_of_points()[t->get_p1_id()-1]->get_ycoord());
                 [ +  - ]
     463 [ +  - ][ +  - ]:      58890 :             min_y = min(min_y, _data->get_list_of_points()[t->get_p2_id()-1]->get_ycoord());
                 [ +  - ]
     464 [ +  - ][ +  - ]:      58890 :             min_y = min(min_y, _data->get_list_of_points()[t->get_p3_id()-1]->get_ycoord());
                 [ +  - ]
     465 [ +  - ][ +  - ]:      58890 :             min_y = min(min_y, _data->get_list_of_points()[t->get_p4_id()-1]->get_ycoord());
                 [ +  - ]
     466                 :            : 
     467 [ +  - ][ +  - ]:      58890 :             max_z = max(max_z, _data->get_list_of_points()[t->get_p1_id()-1]->get_zcoord());
                 [ +  - ]
     468 [ +  - ][ +  - ]:      58890 :             max_z = max(max_z, _data->get_list_of_points()[t->get_p2_id()-1]->get_zcoord());
                 [ +  - ]
     469 [ +  - ][ +  - ]:      58890 :             max_z = max(max_z, _data->get_list_of_points()[t->get_p3_id()-1]->get_zcoord());
                 [ +  - ]
     470 [ +  - ][ +  - ]:      58890 :             max_z = max(max_z, _data->get_list_of_points()[t->get_p4_id()-1]->get_zcoord());
                 [ +  - ]
     471 [ +  - ][ +  - ]:      58890 :             min_z = min(min_z, _data->get_list_of_points()[t->get_p1_id()-1]->get_zcoord());
                 [ +  - ]
     472 [ +  - ][ +  - ]:      58890 :             min_z = min(min_z, _data->get_list_of_points()[t->get_p2_id()-1]->get_zcoord());
                 [ +  - ]
     473 [ +  - ][ +  - ]:      58890 :             min_z = min(min_z, _data->get_list_of_points()[t->get_p3_id()-1]->get_zcoord());
                 [ +  - ]
     474 [ +  - ][ +  - ]:      58890 :             min_z = min(min_z, _data->get_list_of_points()[t->get_p4_id()-1]->get_zcoord());
                 [ +  - ]
     475                 :            :         }
     476                 :            : 
     477 [ +  - ][ +  + ]:    4233760 :         if (t->get_material_id() == 0) // only one element has material id of 0
     478                 :            :         {
     479                 :         10 :             _d_max_orig[i] = 0; _d_min_orig[i] = 0; 
     480                 :            :         }
     481                 :            :         else
     482                 :            :         {
     483 [ +  - ][ +  - ]:    4233750 :                 if (_tumor_region_number.find((unsigned) t->get_material_id()) == _tumor_region_number.end()) 
                 [ +  + ]
     484                 :            :             {
     485         [ +  - ]:    4174860 :                 _d_max_orig[i] = thresholds[t->get_material_id() - 1];
     486                 :    4174860 :                 _d_min_orig[i] = 0;
     487                 :            : 
     488                 :            :             }
     489                 :            :             else
     490                 :            :             {
     491                 :      58890 :                 _tumor_volume += _original_tetra_volumes[i];
     492                 :      58890 :                 _d_max_orig[i] = _INFINITY_; 
     493         [ +  - ]:      58890 :                 _d_min_orig[i] = thresholds[t->get_material_id() - 1];
     494                 :            :             }
     495                 :            :         }
     496                 :            :     }
     497                 :            : 
     498                 :            :     // Compute center of tumor
     499         [ +  - ]:         10 :     _tumor_center = SP_Point((max_x + min_x)/2, (max_y + min_y)/2, (max_z + min_z)/2);
     500                 :            : 
     501 [ +  - ][ +  - ]:         10 :     cout << endl << "Tumor size = " << tumor_count << " tetrahedra, with volume = " 
         [ +  - ][ +  - ]
     502 [ +  - ][ +  - ]:         10 :                             << _tumor_volume << " unit_length^3" << endl;
                 [ +  - ]
     503 [ +  - ][ +  - ]:         10 :     cout << "NON-Tumor size = " << _mesh_size - tumor_count << " tetrahedra, with volume = " 
                 [ +  - ]
     504 [ +  - ][ +  - ]:         10 :                         << _total_mesh_volume - _tumor_volume << " unit_length^3" << endl;
                 [ +  - ]
     505 [ +  - ][ +  - ]:         10 :     cout << "Tumor center = " << _tumor_center << endl << endl;
         [ +  - ][ +  - ]
     506                 :            : 
     507 [ +  - ][ +  + ]:         60 :     for (unsigned int current_idx = 1; current_idx <= _data->get_number_of_materials(); current_idx++)
     508 [ +  - ][ +  - ]:         50 :         cout << "Number of elements in Region " << current_idx << " is " << mat_ids[current_idx]
         [ +  - ][ +  - ]
                 [ +  - ]
     509 [ +  - ][ +  - ]:        100 :             << "     with volume " << _material_volumes[current_idx]/1000 << " cm^3" << endl; 
         [ +  - ][ +  - ]
                 [ +  - ]
     510         [ +  - ]:         10 :     cout << endl;
     511                 :         10 : }
     512                 :            : 
     513                 :            : //!
     514                 :            : //! This method computes the data of all tetra. It updates the following vectors: _tetra_weights, _tetra_volumes, 
     515                 :            : //! _tet_material_ids, _d_max, _d_min
     516                 :            : //!
     517                 :         16 : void pdt_plan::compute_tetra_data(const vector<double> & thresholds)
     518                 :            : {
     519                 :            : 
     520         [ -  + ]:         16 :     if (!_is_unnecessary_removed)
     521                 :            :     {
     522                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::compute_tetra_data(): remove_unnecessary_tetra() was not called. Exiting!\033[0m\n", 31);
     523                 :          0 :         exit(-1);
     524                 :            : 
     525                 :            :     }
     526                 :            : 
     527                 :            :     tetrahedron *t;
     528                 :         16 :     int tumor_count = 0; 
     529                 :         16 :     double total_volume = 0.0, tumor_volume = 0.0;
     530                 :            : 
     531                 :         16 :     unsigned int index_to_consider = 0;
     532         [ +  + ]:    6774032 :     for (unsigned int i = 0; i < _mesh_size; i++)
     533                 :            :     {
     534                 :    6774016 :         t = _data->get_list_of_tetrahedra()[i];
     535                 :            : 
     536         [ +  + ]:    6774016 :         if (!_tetra_to_consider[i]) // ignored tetra
     537                 :    4408627 :             continue;
     538                 :            : 
     539                 :    2365389 :         _tetra_volumes[index_to_consider] = _data->compute_tet_volume(t);
     540                 :    2365389 :         total_volume += _tetra_volumes[index_to_consider]; 
     541                 :            :         
     542 [ +  - ][ +  - ]:    2365389 :         if (_tumor_region_number.find((unsigned) t->get_material_id()) != _tumor_region_number.end()) 
                 [ +  + ]
     543                 :            :         {
     544                 :            : 
     545         [ +  - ]:      94224 :             _tumor_indices.push_back(index_to_consider+1);
     546                 :            : 
     547                 :      94224 :             tumor_count++; 
     548                 :            : 
     549                 :            :             // computing volume of the tumor
     550                 :      94224 :             tumor_volume += _tetra_volumes[index_to_consider];
     551                 :            :         }
     552                 :            : 
     553         [ -  + ]:    2365389 :         if (t->get_material_id() == 0) // only one element has material id of 0
     554                 :            :         {
     555                 :          0 :             _d_max[index_to_consider] = 0; _d_min[index_to_consider] = 0; _tetra_weights[index_to_consider] = 0;
     556                 :            :         }
     557                 :            :         else
     558                 :            :         {
     559 [ +  - ][ +  - ]:    2365389 :                         if (_tumor_region_number.find((unsigned) t->get_material_id()) == _tumor_region_number.end()) 
                 [ +  + ]
     560                 :            :             {
     561                 :    2271165 :                 _d_max[index_to_consider] = thresholds[t->get_material_id() - 1];
     562                 :    2271165 :                 _d_min[index_to_consider] = 0;
     563                 :            :             }
     564                 :            :             else
     565                 :            :             {
     566                 :      94224 :                 _d_max[index_to_consider] = _INFINITY_; 
     567                 :      94224 :                 _d_min[index_to_consider] = thresholds[t->get_material_id() - 1];
     568                 :            :             }
     569                 :    2365389 :             _tetra_weights[index_to_consider] = _slopes[t->get_material_id() - 1]*_tetra_volumes[index_to_consider];
     570                 :            : 
     571 [ +  - ][ +  - ]:    2365389 :                         if (_tumor_region_number.find((unsigned) t->get_material_id()) != _tumor_region_number.end()) 
                 [ +  + ]
     572                 :      94224 :                 _tetra_weights[index_to_consider] *= _tumor_weight;
     573                 :            :         }
     574                 :            : 
     575                 :    2365389 :         _tet_material_ids[index_to_consider] = t->get_material_id();
     576                 :            : 
     577                 :    2365389 :         index_to_consider++;
     578                 :            :     }
     579                 :            : 
     580                 :         16 :     cout << endl << "After removing unnecessary tetra: " << endl;
     581                 :         16 :     cout << "   Tumor size = " << tumor_count << " tetrahedra, with volume = " 
     582                 :         16 :                             << tumor_volume << " unit_length^3" << endl;
     583                 :         16 :     cout << "   NON-Tumor size = " << _actual_mesh_size - tumor_count << " tetrahedra, with volume = " 
     584                 :         16 :                         << total_volume - tumor_volume << " unit_length^3" << endl << endl;
     585                 :            : 
     586                 :         16 :     _is_tetra_data_called = true;
     587                 :         16 : }
     588                 :            : 
     589                 :            : //! This method initializes the LP with its required parameters. The instance initialized assumes the default number of fixed sources.
     590                 :         16 : void pdt_plan::initialize_lp()
     591                 :            : {
     592                 :            :     // delete previous instance of LP
     593         [ +  + ]:         16 :     if (_my_lp != nullptr)
     594         [ +  - ]:         10 :         delete _my_lp;
     595                 :            : 
     596                 :            :     // checking if _p_max is set
     597         [ -  + ]:         16 :     if (_p_max.size() != _num_sources)
     598                 :            :     {
     599                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::initialize_lp(): _p_max vector size does not match number of sources. Exiting!\033[0m\n", 31);
     600                 :          0 :         exit(-1);
     601                 :            :     }
     602                 :            : 
     603 [ +  - ][ +  - ]:         16 :     if (!(_guardband_thresholds_called && _is_unnecessary_removed && _is_tetra_data_called))
                 [ -  + ]
     604                 :            :     {
     605                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::initialize_lp(): Dependent functions were not called. Exiting!\033[0m\n", 31);
     606                 :          0 :         exit(-1);
     607                 :            :     }
     608                 :            :     
     609                 :            :     // Creating a new instance of the fixed source optimization lp
     610                 :            :     // TODO: placer->get_num_tumor_elements() returns a 0 but the fixed_sources_cp program does not need it 
     611                 :         16 :     _my_lp = new fixed_sources_cp(_actual_mesh_size, _num_sources, _placer->get_num_tumor_elements(),
     612                 :            :                     _d_max, _d_min, 
     613                 :         16 :                     _p_max, _placer->get_source_element_fluence_matrix(),
     614                 :            :                     _tetra_weights,
     615         [ +  - ]:         16 :                     _tetra_to_consider);
     616                 :            : 
     617 [ +  - ][ -  + ]:         16 :     if (_op_cvar_approach || _op_variance_constraints)
     618                 :          0 :         _my_lp->initialize_cvar_parameters(_beta_cvar, _tet_material_ids); //_original_tet_material_ids);
     619                 :            : 
     620                 :            :         // Defining constraints and objective
     621                 :         16 :         _lp_constraints_defined = _my_lp->define_mosek_constraints();
     622         [ +  - ]:         16 :         if (_lp_constraints_defined)
     623                 :         16 :                 _lp_objective_defined = _my_lp->define_mosek_objective(); 
     624                 :         16 : }
     625                 :            : 
     626                 :            : 
     627                 :            : //! This method initializes the LP with its required parameters. The instance initialized assumes the given number of sources.
     628                 :          0 : void pdt_plan::initialize_lp(unsigned int num_sources, double p_max_value)
     629                 :            : {
     630                 :            :     // delete previous instance of LP
     631         [ #  # ]:          0 :     if (_my_lp != nullptr)
     632         [ #  # ]:          0 :         delete _my_lp;
     633                 :            : 
     634                 :            :     // checking if the specified number of sources is greater than the maximum number of sources found in FullMonte output
     635         [ #  # ]:          0 :     if (num_sources > _num_sources)
     636                 :            :     {
     637                 :            :         fprintf(stderr, "\n\033[1;%dmpdt_plan::initialize_lp(num_sources, p_max_value): number of sources specified "
     638                 :          0 :                             "is greater than the maximum number of sources from the light simulater output. Exiting!\033[0m\n", 31);
     639                 :          0 :         exit(-1);
     640                 :            :     }
     641                 :            : 
     642                 :            :     // checking if _p_max is set
     643         [ #  # ]:          0 :     if (_p_max.size() != num_sources)
     644                 :            :     {
     645                 :          0 :         set_p_max_vector(p_max_value); // TODO this won't work
     646                 :            :     }
     647                 :            : 
     648 [ #  # ][ #  # ]:          0 :     if (!(_guardband_thresholds_called && _is_unnecessary_removed && _is_tetra_data_called))
                 [ #  # ]
     649                 :            :     {
     650                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::initialize_lp(num_sources): Dependent functions were not called. Exiting!\033[0m\n", 31);
     651                 :          0 :         exit(-1);
     652                 :            :     }
     653                 :            :     
     654                 :            :     // Creating a new instance of the fixed source optimization lp
     655                 :          0 :     _my_lp = new fixed_sources_cp(_actual_mesh_size, num_sources, _placer->get_num_tumor_elements(),
     656                 :            :                     _d_max, _d_min, 
     657                 :          0 :                     _p_max, _placer->get_source_element_fluence_matrix(),
     658                 :            :                     _tetra_weights,
     659         [ #  # ]:          0 :                     _tetra_to_consider);
     660                 :            : 
     661 [ #  # ][ #  # ]:          0 :     if (_op_cvar_approach || _op_variance_constraints)
     662                 :          0 :         _my_lp->initialize_cvar_parameters(_beta_cvar, _tet_material_ids);//_original_tet_material_ids);
     663                 :            : 
     664                 :            :         // Defining constraints and objective
     665                 :          0 :         _lp_constraints_defined = _my_lp->define_mosek_constraints();
     666         [ #  # ]:          0 :         if (_lp_constraints_defined)
     667                 :          0 :                 _lp_objective_defined = _my_lp->define_mosek_objective(); 
     668                 :          0 : }
     669                 :            : 
     670                 :            : //!
     671                 :            : //! This method initializes the LP with its required parameters. The instance initialized assumes the given number of sources
     672                 :            : //! and the given fluence matrix. This method is used in the virtual source reduction algorithm 
     673                 :            : //!
     674                 :          0 : void pdt_plan::initialize_lp(unsigned int num_sources, double p_max_value, const vector<vector<float>> & fluence_matrix)
     675                 :            : {
     676                 :            :     // delete previous instance of LP
     677         [ #  # ]:          0 :     if (_my_lp != nullptr)
     678         [ #  # ]:          0 :         delete _my_lp;
     679                 :            : 
     680                 :            :     // checking if _p_max is set
     681         [ #  # ]:          0 :     if (_p_max.size() != num_sources)
     682                 :            :     {
     683                 :          0 :                 _p_max.clear(); 
     684                 :          0 :                 _p_max.resize(num_sources);
     685                 :            :     }
     686                 :          0 :         fill(_p_max.begin(), _p_max.end(), p_max_value); 
     687                 :            : 
     688 [ #  # ][ #  # ]:          0 :     if (!(_guardband_thresholds_called && _is_unnecessary_removed && _is_tetra_data_called))
                 [ #  # ]
     689                 :            :     {
     690                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::initialize_lp(num_sources): Dependent functions were not called. Exiting!\033[0m\n", 31);
     691                 :          0 :         exit(-1);
     692                 :            :     }
     693                 :            :     
     694                 :            :     // Creating a new instance of the fixed source optimization lp
     695                 :          0 :     _my_lp = new fixed_sources_cp(_actual_mesh_size, num_sources, _placer->get_num_tumor_elements(),
     696                 :            :                     _d_max, _d_min, 
     697                 :            :                     _p_max, fluence_matrix,
     698                 :            :                     _tetra_weights,
     699         [ #  # ]:          0 :                     _tetra_to_consider);
     700                 :            : 
     701 [ #  # ][ #  # ]:          0 :     if (_op_cvar_approach || _op_variance_constraints)
     702                 :          0 :         _my_lp->initialize_cvar_parameters(_beta_cvar, _tet_material_ids);
     703                 :            : 
     704                 :            :         // Defining constraints and objective
     705                 :          0 :         _lp_constraints_defined = _my_lp->define_mosek_constraints();
     706         [ #  # ]:          0 :         if (_lp_constraints_defined)
     707                 :          0 :                 _lp_objective_defined = _my_lp->define_mosek_objective(); 
     708                 :          0 : }
     709                 :            : 
     710                 :            : //! This method solves the LP by calling the suitable MOSEK function. source_tailored is true if we are optimizing over line sources with tailored profiles.
     711                 :         16 : void pdt_plan::solve_lp(const vector<double> & guardbanded_thresholds)
     712                 :            : {
     713 [ +  - ][ +  - ]:         32 :     profiler opt_time("Optimization time"); 
     714                 :            : 
     715         [ +  - ]:         16 :     fprintf(stderr, "\n\033[1;%dmSolving the linear program...\033[0m\n", 35);
     716                 :            :     
     717                 :         16 :     double optimizer_time = 0;
     718                 :         16 :     double max_tumor_weight = -1; // initialize to -1 (indicates INF)
     719                 :         16 :     double min_tumor_weight = 0;  // initialize to 0
     720                 :            :     
     721         [ +  + ]:         16 :     if (_vary_tumor_weight) // keep varying the weight of the tumor to get a v100 of 98%
     722                 :            :     {
     723                 :          2 :         double prev_weight = _tumor_weight; // initial previous weight
     724                 :          2 :                 double prev_tumor_v100 = 0.0; // used to terminate solve if tumor v100 dearch when it should have increased (mosek error) 
     725                 :          2 :         double alpha = 1.0;
     726                 :          2 :                 _num_lp_iterations = 0; 
     727                 :            : 
     728                 :          2 :         bool first_trial = true;
     729                 :            :             
     730                 :          2 :                 bool lp_constraints_updated = true;
     731                 :          3 :         while (true) // keep looping if it didn't converge
     732                 :            :         { 
     733                 :          5 :                         _num_lp_iterations++;
     734 [ +  - ][ +  - ]:          5 :             cout << endl << "Current Tumor weight is: " << _tumor_weight << endl;
         [ +  - ][ +  - ]
     735                 :            : 
     736                 :            :             // change the tetra weights accordingly 
     737                 :            : 
     738         [ +  + ]:      29450 :             for (unsigned int i = 0; i < _tumor_indices.size(); i++)
     739                 :            :             {
     740                 :      29445 :                 _tetra_weights[_tumor_indices[i]-1] *= (_tumor_weight/prev_weight);
     741                 :            :             }
     742                 :            : 
     743                 :            :             // create a new lp object
     744         [ -  + ]:          5 :                         if (_my_lp == nullptr) {
     745         [ #  # ]:          0 :                 initialize_lp();  
     746         [ +  + ]:          5 :                         } else if (!first_trial) { 
     747         [ +  - ]:          3 :                                 _my_lp->set_tetra_weights(_tetra_weights); 
     748         [ +  - ]:          3 :                                 lp_constraints_updated = _my_lp->update_tumor_mosek_constraints(_tumor_weight);
     749                 :            :                         }
     750                 :            : 
     751 [ +  - ][ +  - ]:          5 :                         if (!_lp_constraints_defined || !_lp_objective_defined || !lp_constraints_updated)
                 [ -  + ]
     752                 :            :                         {
     753                 :            :                                 fprintf(stderr, "\n\033[1;%dmpdt_plan::solve_lp::constraints and/or objectives are not defined. Exiting..."
     754         [ #  # ]:          0 :                                                 "\033[0m\n", 31); 
     755                 :          0 :                                 _num_lp_iterations = _max_num_lp_iterations + 1; 
     756                 :          2 :                                 return;
     757                 :            :                                 //std::exit(-1);
     758                 :            :                         }
     759                 :            : 
     760                 :            :             // solve the lp
     761         [ -  + ]:          5 :             if (_source_tailored) {
     762         [ #  # ]:          0 :                 _final_source_powers = _my_lp->run_mosek_optimization_with_tailored_profiles(_placer->
     763         [ #  # ]:          0 :                                                 get_fiber_elements(), optimizer_time, _num_non_zero_sources, _total_power);
     764                 :            :                         }
     765         [ -  + ]:          5 :             else if (_op_cvar_approach) {
     766         [ #  # ]:          0 :                 _final_source_powers = _my_lp->run_mosek_optimization_with_cvar_constraints(optimizer_time, _num_non_zero_sources, 
     767                 :          0 :                                     _total_power, _tumor_region_number, (unsigned int)_d_max_tissues.size(), 
     768                 :          0 :                                                                         _cvar_expected_dose_matrix);
     769                 :            :                         }
     770         [ -  + ]:          5 :             else if (_op_variance_constraints) {
     771         [ #  # ]:          0 :                 _final_source_powers = _my_lp->run_mosek_optimization_with_variance_constraints(optimizer_time, _num_non_zero_sources, 
     772                 :          0 :                                     _total_power, _tumor_region_number, (unsigned int)_d_max_tissues.size(), 
     773                 :          0 :                                                                         _source_element_dose_variance);
     774                 :            :                         }
     775         [ -  + ]:          5 :             else if (_power_variance_optimization) {
     776         [ #  # ]:          0 :                 _final_source_powers = _my_lp->run_mosek_optimization_with_power_variation(optimizer_time, _num_non_zero_sources, _total_power, _sigma_power, _power_variance_weight);
     777                 :            :             }
     778                 :            :             else {
     779         [ +  - ]:          5 :                                 _final_source_powers = _my_lp->run_mosek_optimization(optimizer_time, _num_non_zero_sources, _total_power);
     780                 :            :                         }
     781         [ -  + ]:          5 :                         if (_final_source_powers.empty()) {
     782                 :          0 :                                 _num_lp_iterations = _max_num_lp_iterations + 1; 
     783                 :          0 :                                 return;
     784                 :            :                         }
     785                 :            : 
     786                 :          5 :             _is_lp_solved = true; 
     787                 :            : 
     788                 :            :                         // clear old fluence values
     789                 :          5 :                         _final_fluence.clear();
     790                 :            : 
     791         [ +  - ]:          8 :             vector<double> v100_tissues = compute_v_alpha_tissues(guardbanded_thresholds, alpha);
     792                 :            :         
     793         [ +  - ]:          5 :             print_v100s(v100_tissues);
     794                 :            : 
     795         [ +  + ]:          5 :             double temp_prev = first_trial ? 1.0 : prev_weight; // store previous weight
     796                 :          5 :             prev_weight = _tumor_weight;
     797                 :            : 
     798                 :          5 :             double lower_v100_bound = _lower_tumor_v100_bound;
     799                 :          5 :                         double upper_v100_bound = _upper_tumor_v100_bound;
     800                 :            : 
     801                 :            :                         // Compute v100 of tumor using all tumor regions 
     802                 :          5 :                         double tumor_v100 = 0.0; 
     803                 :          5 :                         double total_tumor_volume = 0.0; 
     804         [ +  + ]:         10 :                         for (auto tumor_region : _tumor_region_number) {
     805         [ +  - ]:          5 :                                 total_tumor_volume += _material_volumes[tumor_region];
     806                 :            : 
     807         [ +  - ]:          5 :                                 tumor_v100 += ((v100_tissues[tumor_region - 1]*_material_volumes[tumor_region])/100);
     808                 :            :                         }
     809                 :          5 :                         tumor_v100 = (tumor_v100/total_tumor_volume)*100;
     810                 :            :            
     811 [ +  + ][ +  + ]:          5 :             if ((tumor_v100/*v100_tissues[*(_tumor_region_number.begin()) - 1]*/ > lower_v100_bound && 
     812         [ -  + ]:          3 :                                 tumor_v100/*v100_tissues[*(_tumor_region_number.begin()) - 1]*/ < upper_v100_bound) || 
     813                 :          3 :                                 _num_lp_iterations > _max_num_lp_iterations/* || _power_variance_optimization*/) // within target
     814                 :            :             {
     815                 :            :                 break;
     816                 :            :             }
     817 [ +  + ][ -  + ]:          3 :                         else if (tumor_v100/*v100_tissues[*(_tumor_region_number.begin()) - 1]*/ < prev_tumor_v100 && 
     818                 :          1 :                                         temp_prev < _tumor_weight) 
     819                 :            :                         {
     820 [ #  # ][ #  # ]:          0 :                                 std::cerr << "Error during optimization" << std::endl;
     821                 :          0 :                                 _num_lp_iterations = _max_num_lp_iterations + 1; 
     822                 :            :                                 //return; 
     823                 :          0 :                 break;
     824                 :            :                         }
     825         [ +  + ]:          3 :             else if (first_trial)
     826                 :            :             {
     827                 :          1 :                 first_trial = false;
     828                 :            :                 
     829         [ -  + ]:          1 :                 if (tumor_v100/*v100_tissues[*(_tumor_region_number.begin()) - 1]*/ > upper_v100_bound)
     830                 :            :                 {
     831                 :          0 :                     max_tumor_weight = _tumor_weight;
     832                 :          0 :                     _tumor_weight /= 2;
     833                 :            :                 }
     834                 :            :                 else 
     835                 :            :                 {
     836                 :          1 :                     min_tumor_weight = _tumor_weight;
     837                 :          1 :                     _tumor_weight *= 2;
     838                 :            :                 }
     839                 :            : 
     840                 :            :             }
     841                 :            :             else
     842                 :            :             {
     843         [ +  - ]:          2 :                 if (tumor_v100/*v100_tissues[*(_tumor_region_number.begin()) - 1]*/ > upper_v100_bound)
     844                 :            :                 {
     845                 :          2 :                     max_tumor_weight = _tumor_weight;
     846         [ -  + ]:          2 :                     if (min_tumor_weight == 0) {
     847                 :          0 :                         _tumor_weight /= 10; // set to 10 to quickly restore tumor weight
     848                 :            :                     } else {
     849                 :          2 :                         _tumor_weight = (min_tumor_weight + max_tumor_weight) / 2;
     850                 :            :                     }
     851                 :            :                     
     852         [ -  + ]:          2 :                     if (_tumor_weight >= prev_weight) {
     853 [ #  # ][ #  # ]:          0 :                         std::cerr << "Tumor coverage cannot be lowered further." << std::endl;
     854                 :          0 :                         break;
     855                 :            :                     } 
     856                 :            :                 }
     857                 :            :                 else 
     858                 :            :                 {
     859                 :          0 :                     min_tumor_weight = _tumor_weight;
     860         [ #  # ]:          0 :                     if (max_tumor_weight == -1) {
     861                 :          0 :                         _tumor_weight *= 10; // set to 10 to achieve larger tumor weight faster
     862         [ #  # ]:          0 :                         if (_tumor_weight >= 1e11) {
     863 [ #  # ][ #  # ]:          0 :                             std::cerr << "Not able to achieve tumor coverage." << std::endl;
     864                 :          0 :                             _num_lp_iterations = _max_num_lp_iterations + 1; 
     865                 :            :                             //return; 
     866                 :          0 :                             break;
     867                 :            :                         }
     868                 :            :                     } else {
     869                 :          0 :                         _tumor_weight = (min_tumor_weight + max_tumor_weight) / 2;
     870                 :            :                     }
     871                 :            :                 }
     872                 :            :             }
     873                 :            : 
     874         [ +  + ]:          5 :                         prev_tumor_v100 = tumor_v100;//v100_tissues[*(_tumor_region_number.begin()) - 1]; 
     875                 :            : 
     876                 :            :                         //delete _my_lp;
     877                 :            :                         //_my_lp = nullptr;
     878                 :            :         }
     879                 :            :     }   
     880                 :            :     else
     881                 :            :     {
     882                 :            : 
     883         [ -  + ]:         14 :         if (_my_lp == nullptr) // LP is not initialized
     884                 :            :         {
     885         [ #  # ]:          0 :             fprintf(stderr, "\n\033[1;%dmpdt_plan::solve_lp(): Linear program is not initialized. Exiting!\033[0m\n", 31);
     886                 :          0 :             exit(-1);
     887                 :            :         }
     888 [ +  - ][ -  + ]:         14 :                 if (!_lp_constraints_defined || !_lp_objective_defined)
     889                 :            :                 {
     890                 :            :                         fprintf(stderr, "\n\033[1;%dmpdt_plan::solve_lp::constraints and or objectives are not defined. Exiting..."
     891         [ #  # ]:          0 :                                         "\033[0m\n", 31); 
     892                 :          0 :                         std::exit(-1);
     893                 :            :                 }
     894                 :            :    
     895                 :            :         // running the optimization
     896         [ +  + ]:         14 :         if (_source_tailored)
     897 [ +  - ][ +  - ]:          2 :             _final_source_powers = _my_lp->run_mosek_optimization_with_tailored_profiles(_placer->get_fiber_elements(), optimizer_time,
     898                 :          1 :                                     _num_non_zero_sources, _total_power);
     899         [ -  + ]:         13 :         else if (_op_cvar_approach)
     900         [ #  # ]:          0 :             _final_source_powers = _my_lp->run_mosek_optimization_with_cvar_constraints(optimizer_time, _num_non_zero_sources, 
     901                 :          0 :                                 _total_power, _tumor_region_number, (unsigned int)_d_max_tissues.size(), _cvar_expected_dose_matrix);
     902         [ -  + ]:         13 :         else if (_op_variance_constraints)
     903         [ #  # ]:          0 :             _final_source_powers = _my_lp->run_mosek_optimization_with_variance_constraints(optimizer_time, _num_non_zero_sources, 
     904                 :          0 :                                 _total_power, _tumor_region_number, (unsigned int)_d_max_tissues.size(), _source_element_dose_variance);
     905         [ -  + ]:         13 :         else if (_power_variance_optimization) {
     906         [ #  # ]:          0 :             _final_source_powers = _my_lp->run_mosek_optimization_with_power_variation(optimizer_time, _num_non_zero_sources, _total_power, _sigma_power, _power_variance_weight);
     907                 :            :         }
     908                 :            :         else
     909         [ +  - ]:         13 :             _final_source_powers = _my_lp->run_mosek_optimization(optimizer_time, _num_non_zero_sources, _total_power);
     910                 :            : 
     911                 :            :         // setiing source powers to original source powers from population average of OP
     912                 :            :         {
     913                 :            : //             // tumor 4
     914                 :            : //             _final_source_powers.resize(3);
     915                 :            : //             _final_source_powers[0] = 0.00215877;//1.94776e-5;
     916                 :            : //             _final_source_powers[1] = 0.00983784;//8.78593e-5;
     917                 :            : //             _final_source_powers[2] = 0.00506754;//9.70509e-5;
     918                 :            :         }
     919                 :            : 
     920                 :            :         {
     921                 :            :             // tumor 1
     922                 :            : //             _final_source_powers.resize(19);
     923                 :            : //             _final_source_powers[0] = 7.67176e-10;
     924                 :            : //             _final_source_powers[1] = 8.45444e-08;
     925                 :            : //             _final_source_powers[2] = 0;
     926                 :            : //             _final_source_powers[3] = 4.94303e-06;
     927                 :            : //             _final_source_powers[4] = 5.95508e-07;
     928                 :            : //             _final_source_powers[5] = 1.27716e-05;
     929                 :            : //             _final_source_powers[6] = 2.52222e-07;
     930                 :            : //             _final_source_powers[7] = 3.76555e-08;
     931                 :            : //             _final_source_powers[8] = 1.26749e-06;
     932                 :            : //             _final_source_powers[9] = 3.65716e-07;
     933                 :            : //             _final_source_powers[10] = 1.3793e-06;
     934                 :            : //             _final_source_powers[11] = 5.57241e-07;
     935                 :            : //             _final_source_powers[12] = 1.57735e-06;
     936                 :            : //             _final_source_powers[13] = 3.22146e-08;
     937                 :            : //             _final_source_powers[14] = 6.72433e-07;
     938                 :            : //             _final_source_powers[15] = 1.3181e-06;
     939                 :            : //             _final_source_powers[16] = 1.61248e-06;
     940                 :            : //             _final_source_powers[17] = 5.04861e-08;
     941                 :            : //             _final_source_powers[18] = 1.30405e-09;
     942                 :            :         }
     943                 :            :         {
     944                 :            :             // tumor 6
     945                 :            : //             _final_source_powers.resize(14);
     946                 :            : //             _final_source_powers[0] = 8.32786e-08;
     947                 :            : //             _final_source_powers[1] = 3.95238e-05;
     948                 :            : //             _final_source_powers[2] = 3.3434e-05;
     949                 :            : //             _final_source_powers[3] = 8.56056e-07;
     950                 :            : //             _final_source_powers[4] = 1.42737e-06;
     951                 :            : //             _final_source_powers[5] = 0.000127729;
     952                 :            : //             _final_source_powers[6] = 9.52769e-07;
     953                 :            : //             _final_source_powers[7] = 5.8649e-06;
     954                 :            : //             _final_source_powers[8] = 5.52204e-07;
     955                 :            : //             _final_source_powers[9] = 2.03246e-07;
     956                 :            : //             _final_source_powers[10] = 0.000304194;
     957                 :            : //             _final_source_powers[11] = 8.242e-05;
     958                 :            : //             _final_source_powers[12] = 1.07207e-05;
     959                 :            : //             _final_source_powers[13] = 4.97289e-07;
     960                 :            :         }
     961                 :            :     
     962                 :            :     }
     963                 :            :     
     964                 :            : //             _final_source_powers.resize(3);
     965                 :            : //             _final_source_powers[0] = 3.15699e-5;//0.00215877;//1.94776e-5;
     966                 :            : //             _final_source_powers[1] = 0.000410393;//0.00983784;//8.78593e-5;
     967                 :            : //             _final_source_powers[2] = 0.00010072;//0.00506754;//9.70509e-5;
     968                 :            :     // setting the optimizer time vector
     969         [ +  - ]:         16 :     _lp_runtime.push_back(optimizer_time);
     970                 :            : 
     971                 :            :     // setting the cost vector
     972 [ +  - ][ +  - ]:         16 :     _lp_cost.push_back(_my_lp->get_cost());
     973                 :            : 
     974         [ +  - ]:         16 :     _is_lp_solved = true; 
     975                 :            : }
     976                 :            : 
     977                 :            : //! This method computes the cost of the last optimization
     978                 :         22 : double pdt_plan::get_last_cost()
     979                 :            : {
     980                 :         22 :     return _my_lp->get_cost();
     981                 :            : }
     982                 :            : 
     983                 :            : //! This method computes the fluence distribution after optimization and stores the result in _final_fluence
     984                 :         21 : void pdt_plan::compute_post_opt_fluence()
     985                 :            : {
     986         [ -  + ]:         21 :     if (!_is_lp_solved)
     987                 :            :     {
     988                 :          0 :         fprintf(stderr, "\n\033[1;%dmpdt_plan::compute_post_opt_fluence(): LP is not solved yet. Exiting!\033[0m\n", 31);
     989                 :          0 :         exit(-1);
     990                 :            :     }
     991                 :            : 
     992                 :         21 :     _placer->compute_current_fluence(_mesh_size, _final_source_powers, _final_fluence);
     993                 :         21 : }
     994                 :            : 
     995                 :            : //! This method computes the v_alpha of all tissues after optimization at the given thresholds. alpha = 1 ==> v_100, alpha = 0.9 ==>v90, ...etc.
     996                 :         21 : vector<double> pdt_plan::compute_v_alpha_tissues(const vector<double> & guardbanded_thresholds, double alpha)
     997                 :            : {
     998         [ +  + ]:         21 :     if (_final_fluence.size() == 0) // did not compute the final fluence distribution
     999                 :          5 :         compute_post_opt_fluence();
    1000                 :            : 
    1001                 :            : 
    1002                 :         21 :     vector<double> v_alpha = _placer->compute_v_alpha_tissues(_final_fluence, _original_tet_material_ids,
    1003                 :         21 :                     _original_tetra_volumes, guardbanded_thresholds, _tumor_region_number, alpha);
    1004                 :            : 
    1005                 :         21 :     return v_alpha;
    1006                 :            : }
    1007                 :            : 
    1008                 :            : //!
    1009                 :            : //! This method evaluates the placer with a different number of packets
    1010                 :            : //! @param vtk_matrix_file [in]: this specifies the file and name of the vtk file from FullMonte that has the fluence distribution of the
    1011                 :            : //!                       sources according to the new number of packets
    1012                 :            : //! @param guardband [in]: guardband on the healthy tissue thresholds. Range: (0,1]
    1013                 :            : //!
    1014                 :            : // void pdt_plan::evaluate_placer(double new_num_packets, string vtk_matrix_file, double guardband)
    1015                 :            : // {
    1016                 :            : 
    1017                 :            : //     // TODO Depricate this 
    1018                 :            : 
    1019                 :            : //     if (!_is_lp_solved)
    1020                 :            : //     {
    1021                 :            : //         fprintf(stderr, "\n\033[1;%dmpdt_plan::evaluate_placer(): LP is not solved yet. Exiting!\033[0m\n", 31);
    1022                 :            : //         exit(-1);
    1023                 :            : //     }
    1024                 :            : 
    1025                 :            : //     // Creating a new src_placer object
    1026                 :            : //     src_placer * evaluation_placer = new src_placer();
    1027                 :            : //     evaluation_placer->set_num_packets(new_num_packets);
    1028                 :            : //     evaluation_placer->set_total_energy(_total_energy);
    1029                 :            : 
    1030                 :            : //     // Reading VTK file
    1031                 :            : //     evaluation_placer->read_vtk_file(vtk_matrix_file, _mesh_size);
    1032                 :            :     
    1033                 :            : //     // computing the fluence using the placer function since it takes into account the ignored tetra
    1034                 :            : //     vector<double> final_fluence;
    1035                 :            : //     evaluation_placer->compute_current_fluence(_mesh_size, _final_source_powers, final_fluence);
    1036                 :            : 
    1037                 :            : //     double original_scale_factor_packets = _num_packets/_total_energy;
    1038                 :            : //     double new_scale_factor_packets = new_num_packets/_total_energy;
    1039                 :            : //     for (unsigned int i = 0; i < _d_max_tissues.size(); i++)
    1040                 :            : //     {
    1041                 :            : //         if (_d_max_tissues[i] != _INFINITY_)
    1042                 :            : //             _d_max_tissues[i] *= (new_scale_factor_packets/original_scale_factor_packets);
    1043                 :            : //         _d_min_tissues[i] *= (new_scale_factor_packets/original_scale_factor_packets);
    1044                 :            : //     }
    1045                 :            :    
    1046                 :            : //     vector<double> new_thresholds(5);
    1047                 :            : //     new_thresholds[0] = _d_min_tissues[4];
    1048                 :            : //     new_thresholds[4] = _d_min_tissues[4];
    1049                 :            : //     new_thresholds[2] = guardband*(_d_max_tissues[2]/0.1);
    1050                 :            : //     new_thresholds[3] = guardband*(_d_max_tissues[3]/0.1); // dividing by 0.1 since given data is already guardbanded by 10%
    1051                 :            : //     new_thresholds[1] = 10*max(_d_min_tissues[4], max(new_thresholds[2], new_thresholds[3]));
    1052                 :            : 
    1053                 :            : 
    1054                 :            : //     // computing v90 of the different tissues
    1055                 :            : //     double alpha = 0.9; // TODO get as input
    1056                 :            : //     vector<double> v_alpha_tissues_evaluation = evaluation_placer->compute_v_alpha_tissues(final_fluence, _tet_material_ids,
    1057                 :            : //                                                                 _tetra_volumes, new_thresholds,
    1058                 :            : //                                                                 _tumor_region_number, alpha);
    1059                 :            : 
    1060                 :            : //     cout << endl << "Evaluation Grey Matter V90 is:  " << v_alpha_tissues_evaluation[2] << " cm^3" << endl;
    1061                 :            : //     cout <<         "Evaluation White Matter V90 is: " << v_alpha_tissues_evaluation[3] << " cm^3" << endl;
    1062                 :            : //     cout <<         "Evaluation Tumor V90 is:        " << v_alpha_tissues_evaluation[4] << " %" << endl << endl;
    1063                 :            : 
    1064                 :            : //     //*
    1065                 :            : //     bpt::ptime local_time = Util::get_local_time();
    1066                 :            : //     ofstream results_out;
    1067                 :            : //     stringstream results_out_ss;
    1068                 :            : //     results_out_ss << "Results/line_sources/fixed_length/photons10e5/post_opt/" << _wavelength << ".txt";  // TODO: specify as input
    1069                 :            : //     results_out.open(results_out_ss.str().c_str(), ofstream::app);
    1070                 :            : //     results_out << local_time << "\t" << _data_name << "\t" << _source_type << "\t" << _placement_type << "\t" << 
    1071                 :            : //             "interior point\t" << _num_non_zero_sources << "\t" << _tumor_weight << "\t" << _my_lp->get_cost() << 
    1072                 :            : //             "\t" << v_alpha_tissues_evaluation[4] << "\t" << v_alpha_tissues_evaluation[2] << "\t" << v_alpha_tissues_evaluation[3] << "\t" << _lp_runtime[0] 
    1073                 :            : //             << "\t" << _total_power << "\t" << _total_energy << endl;
    1074                 :            : //     results_out.close();
    1075                 :            : //     // */
    1076                 :            : 
    1077                 :            : //     delete evaluation_placer;
    1078                 :            : // }
    1079                 :            : 
    1080                 :            : //!
    1081                 :            : //! This method runs the conservative approach of the OP variation by reading all num_op_cases fullmonte output files, and 
    1082                 :            : //! builds the matrix by assigning each healthy tissue element its maximum fluence, and each tumor element its minimum element
    1083                 :            : //!
    1084                 :          0 : void pdt_plan::conservative_approach_op_variation(unsigned int num_op_cases)
    1085                 :            : {
    1086 [ #  # ][ #  # ]:          0 :     profiler opt_time("Conservative OP Variation"); 
    1087                 :            : 
    1088         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmGenerating dose matrix for conservative OP variation approach...\033[0m\n", 35);
    1089                 :            : 
    1090                 :            :     // get the file path and partial name of the fullmonte output files (convention is <file_name>_<file_number>.vtk where 0 < file_number < num_op_cases
    1091         [ #  # ]:          0 :     string vtk_file_path = _fullmonte_output_file.substr(0, _fullmonte_output_file.size()-5);
    1092                 :            : 
    1093                 :            : //     // create a vector of placers of size num_op_cases
    1094                 :            : //     _placers.resize(num_op_cases);
    1095                 :            : // 
    1096                 :            : //     // initialize the placers 
    1097                 :            : //     for (unsigned int i = 0; i < num_op_cases; i++)
    1098                 :            : //         _placers[i] = new src_placer(); 
    1099                 :            : // 
    1100                 :            : //     // read the num_op_cases fullmonte outputs
    1101                 :            : //      Util::BEGIN_TIME = chrono::steady_clock::now();
    1102                 :            : //     for (unsigned int i = 0; i < num_op_cases; i++)
    1103                 :            : //     {
    1104                 :            : //         stringstream vtk_file_ss; 
    1105                 :            : //         vtk_file_ss << vtk_file_path << i << ".vtk";
    1106                 :            : // 
    1107                 :            : //         _placers[i]->read_vtk_file(vtk_file_ss.str(), _mesh_size);
    1108                 :            : //     }
    1109                 :            : //     Util::END_TIME = chrono::steady_clock::now();
    1110                 :            : 
    1111                 :            :     double vtk_read_time = (double)chrono::duration_cast<chrono::seconds>
    1112                 :          0 :                     (Util::END_TIME - Util::BEGIN_TIME).count(); 
    1113                 :            : 
    1114 [ #  # ][ #  # ]:          0 :     cout << "Time to read all VTK files is: " << vtk_read_time << " secs" << endl;
         [ #  # ][ #  # ]
    1115                 :            :     
    1116 [ #  # ][ #  # ]:          0 :     read_fullmonte_output_diff_op(vtk_file_path, num_op_cases);
    1117                 :            : 
    1118                 :            :     // safety check
    1119 [ #  # ][ #  # ]:          0 :     if (_data->get_list_of_tetrahedra().size() != _mesh_size)
    1120                 :            :     {
    1121         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmconservative_approach_op_variation::Size doesn't match. Exiting...\033[0m\n", 31);
    1122                 :          0 :         exit(-1);
    1123                 :            :     }
    1124                 :            :         
    1125                 :            :     // checking over all matrices to find max fluence for each healthy element and min fluence for each tumor element
    1126         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1127                 :            :     {
    1128         [ #  # ]:          0 :         for (unsigned int j = 0; j < _num_sources; j++)
    1129                 :            :         {
    1130         [ #  # ]:          0 :             for (unsigned int k = 0; k < _mesh_size; k++)
    1131                 :            :             {
    1132 [ #  # ][ #  # ]:          0 :                 if (_tumor_region_number.find((unsigned)_data->get_list_of_tetrahedra()[k]->get_material_id()) 
         [ #  # ][ #  # ]
    1133                 :          0 :                                                 != _tumor_region_number.end())
    1134                 :            :                 {
    1135 [ #  # ][ #  # ]:          0 :                     _placer->get_source_element_fluence_matrix()[j][k] = min(_placer->get_source_element_fluence_matrix()[j][k],
    1136         [ #  # ]:          0 :                                                                         _placers[i]->get_source_element_fluence_matrix()[j][k]);
    1137                 :            :                 }
    1138                 :            :                 else
    1139                 :            :                 {
    1140 [ #  # ][ #  # ]:          0 :                     _placer->get_source_element_fluence_matrix()[j][k] = max(_placer->get_source_element_fluence_matrix()[j][k],
    1141         [ #  # ]:          0 :                                                                         _placers[i]->get_source_element_fluence_matrix()[j][k]);
    1142                 :            :                 }
    1143                 :            :             }
    1144                 :            :         }
    1145                 :            :     }
    1146         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmDone Generating dose matrix for conservative OP variation approach...\033[0m\n", 35);
    1147                 :            : 
    1148                 :          0 : }
    1149                 :            : 
    1150                 :            : //!
    1151                 :            : //! method that reads all fullmonte output cases for tumor op variation, 
    1152                 :            : //! assigns each healthy tissue element fluence to the <percentile> of that element 
    1153                 :            : //! (value that is largest than <percentile> of the elements),
    1154                 :            : //! and each tumor element fluence to the 1-percentile of that element    
    1155                 :            : //! 0 < percentile < 1
    1156                 :            : //!
    1157                 :          0 : void pdt_plan::percentile_approach_op_variation(string vtk_files_path, unsigned int num_op_cases, double percentile)
    1158                 :            : {
    1159 [ #  # ][ #  # ]:          0 :     profiler opt_time("Percentile OP Variation"); 
    1160                 :            : 
    1161         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmGenerating dose matrix for Percentile OP variation approach...\033[0m\n", 35);
    1162                 :            : 
    1163                 :            :     
    1164 [ #  # ][ #  # ]:          0 :     if (percentile < 0 || percentile > 1)
    1165                 :            :     {
    1166         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmpercentile_approach_op_variation::percentile should be between 0 and 1. Exiting...\033[0m\n", 31);
    1167                 :          0 :         exit(-1);
    1168                 :            :     }
    1169                 :            :     
    1170                 :            :     // get the file path and partial name of the fullmonte output files (convention is <file_name>_<file_number>.vtk where 0 < file_number < num_op_cases
    1171                 :            : //     string vtk_file_path = _fullmonte_output_file.substr(0, _fullmonte_output_file.size()-5);
    1172         [ #  # ]:          0 :     string vtk_file_path = vtk_files_path.substr(0, vtk_files_path.size()-5);
    1173                 :            : 
    1174                 :            : //     // create a vector of placers of size num_op_cases
    1175                 :            : //     _placers.resize(num_op_cases);
    1176                 :            : // 
    1177                 :            : //     // initialize the placers 
    1178                 :            : //     for (unsigned int i = 0; i < num_op_cases; i++)
    1179                 :            : //         _placers[i] = new src_placer(); 
    1180                 :            : // 
    1181                 :            : //     // read the num_op_cases fullmonte outputs
    1182                 :            : //      Util::BEGIN_TIME = chrono::steady_clock::now();
    1183                 :            : //     for (unsigned int i = 0; i < num_op_cases; i++)
    1184                 :            : //     {
    1185                 :            : //         stringstream vtk_file_ss; 
    1186                 :            : //         vtk_file_ss << vtk_file_path << i << ".vtk";
    1187                 :            : // 
    1188                 :            : //         _placers[i]->read_vtk_file(vtk_file_ss.str(), _mesh_size);
    1189                 :            : //     }
    1190                 :            : //     Util::END_TIME = chrono::steady_clock::now();
    1191                 :            : // 
    1192                 :            : //     double vtk_read_time = (double)chrono::duration_cast<chrono::seconds>
    1193                 :            : //                     (Util::END_TIME - Util::BEGIN_TIME).count(); 
    1194                 :            : // 
    1195                 :            : //     cout << "Time to read all VTK files is: " << vtk_read_time << " secs" << endl;
    1196                 :            : 
    1197 [ #  # ][ #  # ]:          0 :     read_fullmonte_output_diff_op(vtk_file_path, num_op_cases);
    1198                 :            :     
    1199                 :            :     // safety check
    1200 [ #  # ][ #  # ]:          0 :     if (_data->get_list_of_tetrahedra().size() != _mesh_size)
    1201                 :            :     {
    1202         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmpercentile_approach_op_variation::Size doesn't match. Exiting...\033[0m\n", 31);
    1203                 :          0 :         exit(-1);
    1204                 :            :     }
    1205                 :            :         
    1206                 :            :     // checking over all matrices to find <percentile> fluence for each healthy element and <1-percentile> fluence for each tumor element
    1207                 :            :     
    1208                 :            :     // creating an array for each element of size ceil(percentile*mesh_size) to store the max values for healthy tissue elements (and taking the minimum) 
    1209                 :            :     // or the minimum values for tumor elements (and taking the maximum)
    1210                 :            :     // The temporary arrays will be filled in ascending order for OAR and descending order for tumor elements
    1211                 :            :     
    1212                 :          0 :     unsigned int perc_array_size = 1 + (unsigned int)((1-percentile)*num_op_cases);
    1213         [ #  # ]:          0 :     vector<vector<vector<float>>> percentile_array(_num_sources);
    1214         [ #  # ]:          0 :     vector<vector<vector<unsigned int>>> percentile_placer_index(_num_sources); // stores the index of the case where the percentile was taken from
    1215                 :            : 
    1216         [ #  # ]:          0 :     for (unsigned int i = 0; i < _num_sources; i++)
    1217                 :            :     {
    1218         [ #  # ]:          0 :         percentile_array[i].resize(_mesh_size);
    1219         [ #  # ]:          0 :         percentile_placer_index[i].resize(_mesh_size);
    1220         [ #  # ]:          0 :         for (unsigned int j = 0; j < _mesh_size; j++)
    1221                 :            :         {
    1222         [ #  # ]:          0 :             percentile_array[i][j].resize(perc_array_size);
    1223         [ #  # ]:          0 :             percentile_placer_index[i][j].resize(perc_array_size);
    1224                 :            : 
    1225 [ #  # ][ #  # ]:          0 :             if (_tumor_region_number.find((unsigned)_data->get_list_of_tetrahedra()[j]->get_material_id())
         [ #  # ][ #  # ]
    1226                 :          0 :                                         != _tumor_region_number.end())
    1227         [ #  # ]:          0 :                 fill(percentile_array[i][j].begin(), percentile_array[i][j].end(), _INFINITY_);
    1228                 :            :             else
    1229         [ #  # ]:          0 :                 fill(percentile_array[i][j].begin(), percentile_array[i][j].end(), 0.0);
    1230                 :            : 
    1231         [ #  # ]:          0 :             fill(percentile_placer_index[i][j].begin(), percentile_placer_index[i][j].end(), 0);
    1232                 :            :         }
    1233                 :            :     }
    1234                 :            : 
    1235                 :            :     // looping over all num_op_cases
    1236         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1237                 :            :     {
    1238         [ #  # ]:          0 :         for (unsigned int j = 0; j < _num_sources; j++)
    1239                 :            :         {
    1240         [ #  # ]:          0 :             for (unsigned int k = 0; k < _mesh_size; k++)
    1241                 :            :             {
    1242         [ #  # ]:          0 :                 float fluence_value = _placers[i]->get_source_element_fluence_matrix()[j][k];
    1243 [ #  # ][ #  # ]:          0 :                 if (_tumor_region_number.find((unsigned)_data->get_list_of_tetrahedra()[k]->get_material_id())
         [ #  # ][ #  # ]
    1244                 :          0 :                                                         != _tumor_region_number.end())
    1245                 :            :                 {
    1246         [ #  # ]:          0 :                     for (unsigned int p = 0; p < perc_array_size; p++)
    1247                 :            :                     {
    1248         [ #  # ]:          0 :                         if (fluence_value < percentile_array[j][k][p]) // if the new value is less than one of the elements then shift up one and replace
    1249                 :            :                         {
    1250         [ #  # ]:          0 :                             for (unsigned int z = perc_array_size - 1; z > p; z--)
    1251                 :            :                             {
    1252                 :          0 :                                 percentile_array[j][k][z] = percentile_array[j][k][z-1];
    1253                 :          0 :                                 percentile_placer_index[j][k][z] = percentile_placer_index[j][k][z-1];
    1254                 :            :                             }
    1255                 :          0 :                             percentile_array[j][k][p] = fluence_value;
    1256                 :          0 :                             percentile_placer_index[j][k][p] = i;
    1257                 :          0 :                             break;
    1258                 :            :                         }
    1259                 :            :                     }
    1260                 :            :                 }
    1261                 :            :                 else
    1262                 :            :                 {
    1263         [ #  # ]:          0 :                     for (int p = perc_array_size - 1; p >= 0; p--)
    1264                 :            :                     {
    1265         [ #  # ]:          0 :                         if (fluence_value > percentile_array[j][k][p]) // if the new value is greater than one of the elements then shift down one and replace
    1266                 :            :                         {
    1267         [ #  # ]:          0 :                             for (int z = 0; z < p; z++)
    1268                 :            :                             {
    1269                 :          0 :                                 percentile_array[j][k][z] = percentile_array[j][k][z+1];
    1270                 :          0 :                                 percentile_placer_index[j][k][z] = percentile_placer_index[j][k][z+1];
    1271                 :            :                             }
    1272                 :          0 :                             percentile_array[j][k][p] = fluence_value;
    1273                 :          0 :                             percentile_placer_index[j][k][p] = i;
    1274                 :          0 :                             break;
    1275                 :            :                         }
    1276                 :            :                     }
    1277                 :            :                 }
    1278                 :            :             }
    1279                 :            :         }
    1280                 :            :     }
    1281                 :            :     
    1282                 :            :     // getting the smallest value in the percentile array for the OAR elements and the largest for the tumor elements
    1283         [ #  # ]:          0 :     for (unsigned int j = 0; j < _num_sources; j++)
    1284                 :            :     {
    1285                 :          0 :         float max_oar_percentile = 0.0;
    1286                 :            :         //int index_max = -1;
    1287         [ #  # ]:          0 :         for (unsigned int k = 0; k < _mesh_size; k++)
    1288                 :            :         {
    1289                 :            :             // for each source check the highest percentile value for each OAR tetra and get the fluence matrix of 
    1290                 :            :             // the placer that gives the highest percentile among all tetra (can make it for a specific tissue type)
    1291 [ #  # ][ #  # ]:          0 :             if (_tumor_region_number.find((unsigned)_data->get_list_of_tetrahedra()[k]->get_material_id())
         [ #  # ][ #  # ]
    1292                 :          0 :                                         == _tumor_region_number.end())
    1293                 :            :             {
    1294         [ #  # ]:          0 :                 if (percentile_array[j][k][perc_array_size - 1] > max_oar_percentile)
    1295                 :            :                 {
    1296                 :            :                     //index_max = percentile_placer_index[j][k][perc_array_size - 1];
    1297                 :          0 :                     max_oar_percentile = percentile_array[j][k][perc_array_size - 1];
    1298                 :            :                 }
    1299                 :            :             }
    1300                 :            :         }
    1301                 :            : 
    1302         [ #  # ]:          0 :         for (unsigned int k = 0; k < _mesh_size; k++)
    1303                 :            :         {
    1304         [ #  # ]:          0 :             _placer->get_source_element_fluence_matrix()[j][k] =  percentile_array[j][k][perc_array_size - 1];
    1305                 :            : //              _placer->get_source_element_fluence_matrix()[j][k] =  _placers[index_max]->get_source_element_fluence_matrix()[j][k];
    1306                 :            :         }
    1307                 :            :     }
    1308                 :            : 
    1309                 :            : 
    1310         [ #  # ]:          0 :     for (unsigned int i = 0; i < percentile_array.size(); i++)
    1311                 :            :     {
    1312         [ #  # ]:          0 :         for (unsigned int j = 0; j < percentile_array[i].size(); j++)
    1313                 :            :         {
    1314                 :          0 :             percentile_array[i][j].clear();
    1315                 :          0 :             percentile_placer_index[i][j].clear();
    1316                 :            :         }
    1317                 :          0 :         percentile_array[i].clear();
    1318                 :          0 :         percentile_placer_index[i].clear();
    1319                 :            :     }
    1320                 :          0 :     percentile_array.clear();
    1321                 :          0 :     percentile_placer_index.clear();
    1322                 :            : 
    1323         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmDone Generating dose matrix for conservative OP variation approach...\033[0m\n", 35);
    1324                 :          0 : }
    1325                 :            : 
    1326                 :            : //! This method returns a set of v100 for different tissues with different sets of OP
    1327                 :          0 : vector<vector<double>> pdt_plan::compute_v_alpha_tissues_percentile_approach(const vector<double> & guardbanded_thresholds, double alpha, unsigned int num_op_cases)
    1328                 :            : {
    1329         [ #  # ]:          0 :     if (_placers.size() == 0)
    1330                 :            :     {
    1331                 :          0 :         fprintf(stderr, "\033[1;%dmcompute_v_alpha_percentile_approach::FullMonte output for different cases was not read before. Exiting...\033[0m\n", 31);
    1332                 :          0 :         exit(-1);
    1333                 :            :     }
    1334                 :            : 
    1335         [ #  # ]:          0 :     vector<vector<double>> v_alphas(num_op_cases);
    1336                 :            : 
    1337         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1338                 :            :     {
    1339                 :          0 :         vector<double> post_opt_fluence;
    1340         [ #  # ]:          0 :         _placers[i]->compute_current_fluence(_mesh_size, _final_source_powers, post_opt_fluence);
    1341                 :            : 
    1342         [ #  # ]:          0 :         v_alphas[i] = _placers[i]->compute_v_alpha_tissues(post_opt_fluence, _original_tet_material_ids,
    1343                 :          0 :                     _original_tetra_volumes, guardbanded_thresholds, _tumor_region_number, alpha);
    1344                 :            :     }
    1345                 :            :     
    1346                 :          0 :     return v_alphas;
    1347                 :            : }
    1348                 :            : 
    1349                 :            : //!
    1350                 :            : //! This method creates a feasibility problem that guarantees that no tetra exceeds its threshold dose 
    1351                 :            : //! with a probability of 1-delta when the optical properties are varying
    1352                 :            : //!
    1353                 :          0 : void pdt_plan::op_variation_feasibility_problem(const vector<double> & /*guardbanded_thresholds*/, 
    1354                 :            :                                                                                 unsigned int num_op_cases, double delta)
    1355                 :            : {
    1356 [ #  # ][ #  # ]:          0 :     profiler feasibility("OP Variation Feasibility Study"); 
    1357                 :            :     /*
    1358                 :            :      *  The method follows the following steps:
    1359                 :            :      *  1- Set the mean value for each tetrahedron to the dose it gets from the average optical properties
    1360                 :            :      *  2- Read all Fullmonte outputs for different sets of optical properties.
    1361                 :            :      *  3- Fit a normal distribution for the fluence at each tetrahedron due to each source by assigning a certain standard deviation
    1362                 :            :      *  4- Compute the covariance matrix for each tetrahedron between the different sources 
    1363                 :            :      *  5- build the feasibility problem and solve it
    1364                 :            :      *  6- If a solution is found compute the probablity of tissues v100s and output the dose-expected volume histogram
    1365                 :            :      */
    1366                 :            : 
    1367         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmStarting the feasibility study of the OP variation...\033[0m\n", 35);
    1368                 :            : 
    1369 [ #  # ][ #  # ]:          0 :     vector<vector<float>> mean_tetra_fluence = _placer->get_source_element_fluence_matrix();
    1370                 :            : 
    1371                 :            :     
    1372 [ #  # ][ #  # ]:          0 :     if (delta < 0 || delta > 1)
    1373                 :            :     {
    1374         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmop_variation_feasibility_problem::delta should be between 0 and 1. Exiting...\033[0m\n", 31);
    1375                 :          0 :         exit(-1);
    1376                 :            :     }
    1377                 :            :     
    1378                 :            :     // get the file path and partial name of the fullmonte output files (convention is <file_name>_<file_number>.vtk where 0 < file_number < num_op_cases
    1379                 :            :     string vtk_file_path = "/media/yassineab/DATAPART1/PhD_Research_temp_data/variation_optical_props/line_sources/flat_profiles/"
    1380         [ #  # ]:          0 :                             "normal_distribution/vary_mu_a_s/Colin27_tagged_tumor_4_line_sources_tet_fluence_v1_635nm_10p_volume_tumor_";
    1381                 :            : 
    1382 [ #  # ][ #  # ]:          0 :     read_fullmonte_output_diff_op(vtk_file_path, num_op_cases);
    1383                 :            : 
    1384                 :            : 
    1385                 :            :     // TODO: Continue this approach if time permits
    1386                 :          0 : }
    1387                 :            : 
    1388                 :            : //!
    1389                 :            : //! This method minimizes the expected cost while constraining the expected dose
    1390                 :            : //! received by the highest 100*(1-beta)% of OAR tissues and lowest 100*(1-beta)% of the tumore 
    1391                 :            : //!
    1392                 :          0 : void pdt_plan::cvar_optimization_op_problem(unsigned int num_op_cases, double beta)
    1393                 :            : {
    1394 [ #  # ][ #  # ]:          0 :     profiler cvar_opt("CVAR OPT problem");
    1395                 :            :     
    1396         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmStarting the CVAR optimization study of the OP variation...\033[0m\n", 35);
    1397                 :            : 
    1398 [ #  # ][ #  # ]:          0 :     if (beta < 0 || beta > 1)
    1399                 :            :     {
    1400         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmcvar_optimization_op_problem::beta should be between 0 and 1. Exiting...\033[0m\n", 31);
    1401                 :          0 :         exit(-1);
    1402                 :            :     }
    1403                 :            : 
    1404                 :            :     // get the file path and partial name of the fullmonte output files (convention is <file_name>_<file_number>.vtk where 0 < file_number < num_op_cases
    1405                 :            :     string vtk_file_path = "/media/yassineab/DATAPART1/PhD_Research_temp_data/variation_optical_props/line_sources/flat_profiles/"
    1406         [ #  # ]:          0 :                             "normal_distribution/vary_mu_a_s/Colin27_tagged_tumor_6_line_sources_tet_fluence_v1_635nm_10p_volume_tumor_1.vtk";
    1407                 :            : //     string vtk_file_path = "/media/yassineab/DATAPART1/PhD_Research_temp_data/variation_optical_props/line_sources/flat_profiles/"
    1408                 :            : //                             "normal_distribution/vary_mu_a_s/with_sampling_probabilities/Colin27_tagged_tumor_4_line_sources_tet_fluence_v1_635nm_10p_volume_tumor_1.vtk";
    1409                 :            : 
    1410                 :            : //     read_fullmonte_output_diff_op(vtk_file_path, num_op_cases);
    1411                 :            : 
    1412                 :            :     // TODO: build the dose response function from the average fluence (assign probabilities to bins)
    1413                 :            : //     percentile_approach_op_variation(vtk_file_path, num_op_cases, 0.5);
    1414                 :            : 
    1415                 :            : 
    1416                 :            :     // reading discretized probabilities of the rand OPs
    1417                 :            :     string probs_file = "/home/yassineab/PhD_Research/src_placement/TCL_scripts/fixed_sources_vtk/line_sources/fixed_length/variation_optical_props/635nm/"
    1418         [ #  # ]:          0 :                             "normal_distribution/vary_mu_a_s/with_sampling_probabilities/mu_eff_probs.txt";
    1419                 :            : 
    1420         [ #  # ]:          0 :     ifstream probs_in(probs_file);
    1421 [ #  # ][ #  # ]:          0 :     if (!probs_in)
    1422                 :            :     {
    1423                 :          0 :         fprintf(stderr, "\033[1;%dmpdt_plan::cvar_optimization_op_problem: Could not read "
    1424         [ #  # ]:          0 :                         "Probabilities file %s\033[0m\n", 31, probs_file.c_str());
    1425                 :          0 :         exit(-1);
    1426                 :            :     }
    1427                 :            : 
    1428         [ #  # ]:          0 :     vector<double> mu_eff_probs(num_op_cases, 0);
    1429         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1430                 :            :     {
    1431                 :          0 :         string line;
    1432         [ #  # ]:          0 :         getline(probs_in, line);
    1433         [ #  # ]:          0 :         stringstream str(line);
    1434         [ #  # ]:          0 :         str >> mu_eff_probs[i];
    1435                 :            :     }
    1436                 :            : 
    1437                 :          0 :     double sum_probs = 0; 
    1438         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1439                 :          0 :         sum_probs += mu_eff_probs[i];
    1440                 :            : 
    1441                 :            :     // normalizing the probabilities just in case we used a different number of cases
    1442         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1443                 :            :     {
    1444                 :          0 :         mu_eff_probs[i] = mu_eff_probs[i]/sum_probs;
    1445                 :            :     }
    1446                 :            : 
    1447                 :            :     // Computing the expected dose matrix from the probabilities and the _placers vector
    1448         [ #  # ]:          0 :     _cvar_expected_dose_matrix.resize(_num_sources);
    1449         [ #  # ]:          0 :     for (unsigned int i = 0; i < _num_sources; i++)
    1450                 :            :     {
    1451         [ #  # ]:          0 :         _cvar_expected_dose_matrix[i].resize(_mesh_size);
    1452         [ #  # ]:          0 :         fill(_cvar_expected_dose_matrix[i].begin(), _cvar_expected_dose_matrix[i].end(), 0);
    1453                 :            :         
    1454         [ #  # ]:          0 :         for (unsigned int j = 0; j < _mesh_size; j++)
    1455                 :            :         {
    1456         [ #  # ]:          0 :             for (unsigned int k = 0; k < num_op_cases; k++)
    1457                 :            :             {
    1458                 :            :        //         _cvar_expected_dose_matrix[i][j] += (float)mu_eff_probs[k]*_placers[k]->get_source_element_fluence_matrix()[i][j];
    1459                 :            :             }
    1460                 :            :         }
    1461                 :            :     }
    1462                 :            :      
    1463                 :            : 
    1464                 :          0 :     _op_cvar_approach = true;
    1465                 :          0 :     _beta_cvar = beta; 
    1466                 :            : 
    1467                 :            : //     _op_variance_constraints = true;
    1468                 :            :     //compute_tetra_dose_variance(); // TODO: doesn't work if we are calling ignore_unnecessary_tetra
    1469                 :            : 
    1470                 :            :     //_ignore_unnecessary_tetra = false; // TODO: doesn't work if we are using cvar_expected_dose_matrix
    1471                 :          0 : }
    1472                 :            : 
    1473                 :            : //! This method reads different outputs of FullMonte corresponding to different sets of optical properties, and places the results in _placers vector
    1474                 :          0 : void pdt_plan::read_fullmonte_output_diff_op(string vtk_files_path, unsigned int num_op_cases)
    1475                 :            : {
    1476 [ #  # ][ #  # ]:          0 :     profiler op_cases("Reading Fullmonte output diff OP");
    1477                 :            : 
    1478                 :            :     // create a vector of placers of size num_op_cases
    1479         [ #  # ]:          0 :     for (unsigned int i = 0; i < _placers.size(); i++)
    1480         [ #  # ]:          0 :         if (_placers[i] != nullptr)
    1481         [ #  # ]:          0 :             delete _placers[i];
    1482                 :            : 
    1483                 :          0 :     _placers.clear();
    1484         [ #  # ]:          0 :     _placers.resize(num_op_cases);
    1485                 :            : 
    1486                 :            :     // initialize the placers 
    1487         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1488 [ #  # ][ #  # ]:          0 :         _placers[i] = new src_placer(); 
    1489                 :            : 
    1490                 :            :     // read the num_op_cases fullmonte outputs
    1491                 :          0 :         Util::BEGIN_TIME = chrono::steady_clock::now();
    1492         [ #  # ]:          0 :     for (unsigned int i = 0; i < num_op_cases; i++)
    1493                 :            :     {
    1494         [ #  # ]:          0 :         stringstream vtk_file_ss; 
    1495 [ #  # ][ #  # ]:          0 :         vtk_file_ss << vtk_files_path << i << ".vtk";
                 [ #  # ]
    1496                 :            : 
    1497 [ #  # ][ #  # ]:          0 :         _placers[i]->read_vtk_file(vtk_file_ss.str(), _mesh_size);
    1498                 :            :     }
    1499                 :          0 :     Util::END_TIME = chrono::steady_clock::now();
    1500                 :            : 
    1501                 :            :     double vtk_read_time = (double)chrono::duration_cast<chrono::seconds>
    1502                 :          0 :                     (Util::END_TIME - Util::BEGIN_TIME).count(); 
    1503                 :            : 
    1504 [ #  # ][ #  # ]:          0 :     cout << "Time to read all VTK files is: " << vtk_read_time << " secs" << endl;
         [ #  # ][ #  # ]
    1505                 :          0 : }
    1506                 :            : 
    1507                 :            : //! This method computes the variance of the dose at each tetrahedron from each source with unit power due to variation in OP
    1508                 :          0 : void pdt_plan::compute_tetra_dose_variance()
    1509                 :            : {
    1510 [ #  # ][ #  # ]:          0 :         profiler variance_compu("Dose variance computation");
    1511                 :            : 
    1512         [ #  # ]:          0 :     fprintf(stderr, "\n\033[1;%dmComputing the variance of the dose at each tetrahedron...\033[0m\n", 35);
    1513                 :            : 
    1514         [ #  # ]:          0 :         if (_placers.size() == 0)
    1515                 :            :         {
    1516         [ #  # ]:          0 :                 fprintf(stderr, "\n\033[1;%dmpdt_plan::compute_tetra_dose_variance()::No variation cases have been read! Exiting...\033[0m\n", 31);
    1517                 :          0 :                 exit(-1);       
    1518                 :            :         }
    1519                 :            : 
    1520                 :          0 :     unsigned int num_op_cases = (unsigned int)_placers.size();
    1521                 :            : 
    1522                 :            :     // Clearing previously created variance matrix
    1523         [ #  # ]:          0 :     for (unsigned int i = 0; i < _source_element_dose_variance.size(); i++)
    1524                 :            :     {
    1525         [ #  # ]:          0 :         if (_source_element_dose_variance[i].size() != 0)
    1526                 :            :         {
    1527                 :          0 :             _source_element_dose_variance[i].clear();
    1528         [ #  # ]:          0 :             _source_element_dose_variance[i].resize(0);
    1529                 :            :         }
    1530                 :            :     }
    1531                 :          0 :     _source_element_dose_variance.clear();
    1532                 :            :     
    1533                 :            :     // Creating the variance matrix
    1534         [ #  # ]:          0 :     _source_element_dose_variance.resize(_num_sources);
    1535         [ #  # ]:          0 :     for (unsigned int i = 0; i < _num_sources; i++)
    1536                 :            :     {
    1537         [ #  # ]:          0 :         _source_element_dose_variance[i].resize(_actual_mesh_size);
    1538                 :            :         
    1539                 :          0 :         int actual_tet_idx = -1;
    1540         [ #  # ]:          0 :         for (unsigned int j = 0; j < _mesh_size; j++)
    1541                 :            :         {
    1542         [ #  # ]:          0 :             if (!_tetra_to_consider[j]) // ignored tetra
    1543                 :          0 :                 continue; 
    1544                 :            :             
    1545                 :          0 :             actual_tet_idx++; 
    1546                 :            : 
    1547                 :            :             // computing the expected value of the dose first
    1548                 :          0 :             double tetra_expected_value = 0; 
    1549         [ #  # ]:          0 :             for (unsigned int k = 0; k < num_op_cases; k++)
    1550                 :            :             {
    1551         [ #  # ]:          0 :                 tetra_expected_value += _placers[k]->get_source_element_fluence_matrix()[i][j];
    1552                 :            :             }
    1553                 :          0 :             tetra_expected_value = tetra_expected_value/num_op_cases; 
    1554                 :            : 
    1555                 :            :             // Computing the variance
    1556                 :          0 :             double tetra_variance = 0.0;
    1557         [ #  # ]:          0 :             for (unsigned int k = 0; k < num_op_cases; k++)
    1558                 :            :             {
    1559         [ #  # ]:          0 :                 tetra_variance += pow(_placers[k]->get_source_element_fluence_matrix()[i][j] - tetra_expected_value, 2);
    1560                 :            :             }
    1561                 :          0 :             tetra_variance = tetra_variance/(num_op_cases - 1); // N-1 degrees of freedom 
    1562                 :            : 
    1563                 :          0 :             _source_element_dose_variance[i][actual_tet_idx] = tetra_variance;
    1564                 :            :         }
    1565         [ #  # ]:          0 :         if ((int)_actual_mesh_size != actual_tet_idx+1)
    1566                 :            :         {
    1567 [ #  # ][ #  # ]:          0 :             cerr << endl << "ERROR in computing dose variance matrix. Exiting...." << endl;
                 [ #  # ]
    1568                 :          0 :             exit (-1);
    1569                 :            :         }
    1570                 :            :     }
    1571                 :          0 : }
    1572                 :            : 
    1573                 :          0 : vector<vector<Src_Abstract*> > & pdt_plan::get_placement_results() {
    1574                 :          0 :     return _placement_results;
    1575                 :            : }
    1576                 :            : 
    1577                 :            : //! This method creates and runs an SA engine for source placement
    1578                 :          1 : void pdt_plan::run_sa_engine()
    1579                 :            : {
    1580 [ +  - ][ +  - ]:          2 :     profiler SA("SA optimization"); 
    1581         [ -  + ]:          1 :     for (unsigned i = 0; i < _placement_results.size(); i++) {
    1582         [ #  # ]:          0 :         for (unsigned j = 0; j < _placement_results[i].size(); j++) {
    1583 [ #  # ][ #  # ]:          0 :             if (_placement_results[i][j]) delete _placement_results[i][j];
    1584                 :            :         }
    1585                 :            :     }
    1586                 :          1 :     _placement_results.clear();
    1587                 :            :     
    1588                 :            :         sa_engine* annealer = new sa_engine(_fparser, _fullmonte_engine, _data, 
    1589                 :            :                         _d_min_orig, _d_max_orig, _original_tumor_indices, _tumor_volume, 
    1590 [ +  - ][ +  - ]:          1 :                         _file_path_override, _sa_engine_early_terminate, this);
    1591         [ -  + ]:          1 :     if (!_constrain_injection_point) {
    1592         [ #  # ]:          0 :         annealer->run_sa_engine(-1);
    1593                 :            : 
    1594                 :            :         // TODO: update to save up to 5 best placements
    1595 [ #  # ][ #  # ]:          0 :         cout << "Finished running SA." << endl;
    1596                 :            : 
    1597         [ #  # ]:          0 :         _initial_placement = annealer->get_placement_solution(); 
    1598         [ #  # ]:          0 :         run_fullmonte_multiple_sources(_initial_placement); 
    1599                 :            : 
    1600         [ #  # ]:          0 :         _placement_results.push_back(_initial_placement);
    1601                 :            :     } else {
    1602         [ +  + ]:          2 :         for (unsigned j = 0; j < _num_injection_points; j++) {
    1603         [ +  - ]:          1 :             annealer->run_sa_engine(j);
    1604                 :            : 
    1605                 :            :             // TODO: update this for multiple injection points
    1606 [ +  - ][ +  - ]:          1 :             cout << "Finished running SA." << endl;
    1607                 :            : 
    1608         [ +  - ]:          1 :             _initial_placement = annealer->get_placement_solution(); 
    1609         [ +  - ]:          1 :             run_fullmonte_multiple_sources(_initial_placement); 
    1610                 :            : 
    1611         [ +  - ]:          1 :             _placement_results.push_back(_initial_placement);
    1612                 :            :         }
    1613                 :            :     } 
    1614                 :            : 
    1615         [ +  - ]:          1 :     delete annealer;
    1616                 :          1 : }
    1617                 :            : 
    1618                 :            : //! This method returns an annealer for unit testing purposes
    1619                 :          3 : sa_engine* pdt_plan::get_annealer_for_testing() {
    1620                 :            :     sa_engine* annealer = new sa_engine(_fparser, _fullmonte_engine, _data, 
    1621                 :            :                         _d_min_orig, _d_max_orig, _original_tumor_indices, _tumor_volume, 
    1622         [ +  - ]:          3 :                         _file_path_override, _sa_engine_early_terminate, this);
    1623                 :            :     
    1624                 :          3 :     return annealer;
    1625                 :            : }
    1626                 :            : 
    1627                 :            : 
    1628                 :          0 : bool pdt_plan::get_constrain_injection_point(){
    1629                 :          0 :     return _constrain_injection_point;
    1630                 :            : }
    1631                 :          0 : unsigned pdt_plan::get_num_injection_points(){
    1632                 :          0 :     return _num_injection_points;
    1633                 :            : }
    1634                 :            : 
    1635                 :            : //! This method either runs SA or reads the initial placement and then runs FM to generate dose matrix
    1636                 :          0 : void pdt_plan::get_diffuser_placement_and_dose_matrix() 
    1637                 :            : {
    1638         [ #  # ]:          0 :     if (_running_tests) {
    1639                 :          0 :         return; 
    1640                 :            :     }
    1641                 :            : 
    1642         [ #  # ]:          0 :         if (_run_sa) {
    1643                 :          0 :                 run_sa_engine();
    1644                 :            :         }
    1645                 :            :         else {
    1646                 :          0 :                 read_initial_placement(); 
    1647                 :            :         
    1648                 :            : //         run_rt_feedback(); 
    1649                 :            : 
    1650         [ #  # ]:          0 :         if (_fullmonte_output_file != "") {
    1651                 :          0 :             read_fullmonte_output_matrix(); 
    1652                 :            :         }
    1653                 :            :         else {
    1654                 :          0 :                     run_fullmonte_multiple_sources(get_initial_placement()); 
    1655                 :            :             }
    1656                 :            :     }
    1657                 :            : }
    1658                 :            : 
    1659                 :            : //!
    1660                 :            : //! This method runs a large number of virtual sources, allocates the power across them
    1661                 :            : //!     then it starts eliminating least power sources and re-optimizes.
    1662                 :            : //! The algorithm generates multiple plans the trades-off number of sources to plan quality
    1663                 :            : //! @param num_plans [in]: number of desired plans
    1664                 :            : //! @param thresholds [in]: vector of dose thresholds
    1665                 :            : //! @return Returns a vector of vectors of size equal to num_plans. Each vector is the placement of sources of the corresponding plan.
    1666                 :            : //!
    1667                 :          0 : std::vector<std::vector<Src_Abstract*>> pdt_plan::run_virtual_source_reduction(unsigned num_plans, 
    1668                 :            :                                         const std::vector<double> & thresholds)
    1669                 :            : {
    1670 [ #  # ][ #  # ]:          0 :         profiler virtual_profile("Running virtual source reduction"); 
    1671                 :            : 
    1672 [ #  # ][ #  # ]:          0 :         if (_placement_type != "virtual") {
    1673                 :            :         fprintf(stderr, "\033[1;%dmpdt_plan::run_virtual_source_reduction::Running light source reduction algorithm while "
    1674         [ #  # ]:          0 :                         "placement_type is not virtual. Exiting...\033[0m\n", 31);
    1675                 :          0 :         exit(-1);
    1676                 :            :         }
    1677                 :            : 
    1678                 :            :         /* The method assumes that you already run Fullmonte with the virtual sources */
    1679                 :            : 
    1680                 :          0 :         unsigned num_sources_chosen = 51; // We first choose the highest 51 power sources and then eliminate one at a time. 
    1681                 :            : 
    1682 [ #  # ][ #  # ]:          0 :         if (num_sources_chosen > _num_sources || num_plans > num_sources_chosen)
    1683                 :            :         {
    1684                 :            :                 fprintf(stderr, "\033[1;%dmpdt_plan::run_virtual_source_reduction::Number of chosen sources is greater than "
    1685         [ #  # ]:          0 :                                         "the number of virtual sources. Exiting.....\033[0m\n", 31);
    1686                 :          0 :                 exit(-1);
    1687                 :            :         }
    1688                 :            :         
    1689                 :            :         // We don't need to reallocate
    1690                 :          0 :         bool old_vary_tumor_weight = _vary_tumor_weight;
    1691                 :          0 :         _vary_tumor_weight = false;
    1692                 :            : 
    1693                 :            :         // 1- Run Power allocation on virtual sources
    1694         [ #  # ]:          0 :         initialize_lp();
    1695         [ #  # ]:          0 :         solve_lp(thresholds); 
    1696                 :            : 
    1697                 :            :         
    1698                 :            :         // 2- Copy the original fluence matrix and sources
    1699         [ #  # ]:          0 :         vector<vector<float>> old_source_fluence_matrix(_num_sources); 
    1700         [ #  # ]:          0 :         for (unsigned i = 0; i < _num_sources; i++) {
    1701 [ #  # ][ #  # ]:          0 :                 old_source_fluence_matrix[i].resize(_placer->get_source_element_fluence_matrix()[i].size());
    1702 [ #  # ][ #  # ]:          0 :                 for (unsigned j = 0; j < _placer->get_source_element_fluence_matrix()[i].size(); j++) {
    1703         [ #  # ]:          0 :                         old_source_fluence_matrix[i][j] = _placer->get_source_element_fluence_matrix()[i][j]; 
    1704                 :            :                 }
    1705                 :            :         }
    1706         [ #  # ]:          0 :         vector<Src_Abstract*> old_source_placement = _initial_placement;
    1707                 :            : 
    1708                 :            :         // 3- Start elimination loop
    1709         [ #  # ]:          0 :         vector<vector<Src_Abstract*>> new_plans(num_plans+1);
    1710         [ #  # ]:          0 :     new_plans[0] = old_source_placement;
    1711                 :          0 :         unsigned curr_plan_idx = 1;
    1712         [ #  # ]:          0 :         for (unsigned new_num_sources = num_sources_chosen; new_num_sources >= 1; new_num_sources--)
    1713                 :            :         {
    1714 [ #  # ][ #  # ]:          0 :                 vector<double> curr_source_powers = get_final_source_powers(); 
    1715                 :            : 
    1716                 :            :                 // sort the source powers 
    1717         [ #  # ]:          0 :                 std::vector<std::pair<unsigned, double>> power_idx_pairs(curr_source_powers.size()); 
    1718         [ #  # ]:          0 :                 for (unsigned i = 0; i < curr_source_powers.size(); i++) {
    1719                 :          0 :                         power_idx_pairs[i].first = i;
    1720                 :          0 :                         power_idx_pairs[i].second = curr_source_powers[i];
    1721                 :            :                 }
    1722                 :          0 :                 std::sort(power_idx_pairs.begin(), power_idx_pairs.end(), [] (const std::pair<unsigned, double> & p1, 
    1723         [ #  # ]:          0 :                                         const std::pair<unsigned, double> & p2) -> bool { return p1.second > p2.second; });
    1724                 :            : 
    1725                 :            :                 
    1726                 :            :                 // building the new matrix based on the new sources chosen
    1727                 :          0 :                 vector<vector<float>> new_source_fluence_matrix; 
    1728                 :          0 :                 vector<Src_Abstract*> new_source_placement;
    1729         [ #  # ]:          0 :                 for (unsigned i = 0; i < new_num_sources; i++) {
    1730         [ #  # ]:          0 :                         new_source_fluence_matrix.push_back(old_source_fluence_matrix[power_idx_pairs[i].first]);
    1731                 :            : 
    1732         [ #  # ]:          0 :                         new_source_placement.push_back(old_source_placement[power_idx_pairs[i].first]);
    1733                 :            :                 }
    1734                 :            :                         
    1735         [ #  # ]:          0 :                 if (new_num_sources <= num_plans)
    1736                 :            :                 {
    1737         [ #  # ]:          0 :                         new_plans[curr_plan_idx] = new_source_placement;
    1738                 :          0 :                         curr_plan_idx++; 
    1739                 :            :                 }
    1740                 :            : 
    1741 [ #  # ][ #  # ]:          0 :                 cout << "Number of sources for re-optimization is: " << new_num_sources << endl;
                 [ #  # ]
    1742                 :            : 
    1743                 :            :                 // initialize new lp
    1744         [ #  # ]:          0 :                 initialize_lp(new_num_sources, _INFINITY_, new_source_fluence_matrix); 
    1745                 :            : 
    1746                 :            :                 // Solve lp 
    1747         [ #  # ]:          0 :                 solve_lp(thresholds); 
    1748                 :            : 
    1749                 :            :                 // Re-initialize the source fluence matrix and sources vector
    1750                 :          0 :                 old_source_placement.clear(); 
    1751         [ #  # ]:          0 :                 old_source_placement = new_source_placement;
    1752                 :            : 
    1753         [ #  # ]:          0 :                 for (unsigned i = 0; i < old_source_fluence_matrix.size(); i++) {
    1754                 :          0 :                         old_source_fluence_matrix[i].clear();
    1755                 :            :                 }
    1756                 :          0 :                 old_source_fluence_matrix.clear(); 
    1757                 :            : 
    1758         [ #  # ]:          0 :                 old_source_fluence_matrix.resize(new_num_sources); 
    1759         [ #  # ]:          0 :                 for (unsigned i = 0; i < new_num_sources; i++) 
    1760                 :            :                 {
    1761         [ #  # ]:          0 :                         old_source_fluence_matrix[i] = new_source_fluence_matrix[i]; 
    1762                 :            :                 }
    1763                 :            :         } // end elimination loop
    1764                 :            : 
    1765                 :          0 :         _vary_tumor_weight = old_vary_tumor_weight;
    1766                 :            : 
    1767                 :          0 :         return new_plans;
    1768                 :            : }
    1769                 :            : 
    1770                 :            : //! This method runs the problem with power variation
    1771                 :          0 : void pdt_plan::power_variation_problem(double weight, double uncertainty, vector<double>& thresholds)
    1772                 :            : {
    1773                 :          0 :     printf("\nRunning power variation...\n");
    1774                 :          0 :     _num_packets *= 2;
    1775                 :          0 :     _power_variance_optimization = true;
    1776                 :          0 :     _power_variance_weight = weight;
    1777                 :          0 :     _sigma_power = uncertainty/2.0;
    1778                 :            : 
    1779                 :            :     // approach 2: run lp
    1780                 :            :     //initialize_lp();
    1781                 :          0 :     solve_lp(thresholds);
    1782                 :            : 
    1783                 :          0 :     _power_variance_optimization = false;
    1784                 :          0 :     _num_packets /= 2;
    1785                 :          0 : }
    1786                 :            : 
    1787                 :            : // Getters
    1788                 :            : 
    1789                 :            : //! This method returns the final source powers after optimizations
    1790                 :          4 : vector<double> & pdt_plan::get_final_source_powers()
    1791                 :            : {
    1792                 :          4 :     return _final_source_powers;
    1793                 :            : }
    1794                 :            : 
    1795                 :            : //! This method returns the final fluence distribution of the mesh
    1796                 :          8 : vector<double> & pdt_plan::get_final_fluence()
    1797                 :            : {
    1798                 :          8 :     return _final_fluence;
    1799                 :            : }
    1800                 :            : 
    1801                 :            : //! This method returns the vector that stores the lp costs after optimization (usually it is of size 1 unless we are running the LSR algorithm)
    1802                 :          0 : vector<double> & pdt_plan::get_lp_cost()
    1803                 :            : {
    1804                 :          0 :     return _lp_cost;
    1805                 :            : }
    1806                 :            : 
    1807                 :            : //! This method returns the vector that stores the lp runtimes after optimization (usually it is of size 1 unless we are running the LSR algorithm)
    1808                 :          0 : vector<double> & pdt_plan::get_lp_runtime()
    1809                 :            : {
    1810                 :          0 :     return _lp_runtime;
    1811                 :            : }
    1812                 :            : 
    1813                 :            : //! This method returns the vector of max thresholds for all elements
    1814                 :          0 : vector<double> & pdt_plan::get_d_max_vec()
    1815                 :            : {
    1816                 :          0 :     return _d_max;
    1817                 :            : }
    1818                 :            : 
    1819                 :            : //! This method returns the vector of min thresholds for all elements
    1820                 :          0 : vector<double> & pdt_plan::get_d_min_vec()
    1821                 :            : {
    1822                 :          0 :     return _d_min;
    1823                 :            : }
    1824                 :            : 
    1825                 :            : //! This method returns the slopes vector
    1826                 :          0 : vector<double> & pdt_plan::get_slopes()
    1827                 :            : {
    1828                 :          0 :     return _slopes;
    1829                 :            : }
    1830                 :            : 
    1831                 :            : //! This method returns a vector of max thresholds per tissue type (size is the same as number of materials)
    1832                 :          0 : vector<double> & pdt_plan::get_d_max_tissues()
    1833                 :            : {
    1834                 :          0 :     return _d_max_tissues;
    1835                 :            : }
    1836                 :            : 
    1837                 :            : //! This method returns a vector of min thresholds per tissue type (size is the same as number of materials)
    1838                 :          0 : vector<double> & pdt_plan::get_d_min_tissues()
    1839                 :            : {
    1840                 :          0 :     return _d_min_tissues;
    1841                 :            : }
    1842                 :            : 
    1843                 :            : //! This method returns a vector of safety multipliers per tissue type (size is the same as number of materials)
    1844                 :          0 : vector<double> & pdt_plan::get_guardband_tissues()
    1845                 :            : {
    1846                 :          0 :     return _guardband_tissues;
    1847                 :            : }
    1848                 :            : 
    1849                 :            : //! This method returns the volume of each tetra before ignoring the unnecessary tetra
    1850                 :          0 : vector<double> & pdt_plan::get_original_tetra_volumes()
    1851                 :            : {
    1852                 :          0 :     return _original_tetra_volumes;
    1853                 :            : }
    1854                 :            : 
    1855                 :            : //! This method returns the volume of each tetra after ignoring the unnecessary tetra
    1856                 :          0 : vector<double> & pdt_plan::get_tetra_volumes()
    1857                 :            : {
    1858                 :          0 :     return _tetra_volumes;
    1859                 :            : }
    1860                 :            : 
    1861                 :            : //! This method returns the weight of each tetra considered in the optimization
    1862                 :          0 : vector<double> & pdt_plan::get_tetra_weights()
    1863                 :            : {
    1864                 :          0 :     return _tetra_weights;
    1865                 :            : }
    1866                 :            : 
    1867                 :            : //! This method returns a vector of the material ids of all tetras in the mesh
    1868                 :          0 : vector<unsigned int> & pdt_plan::get_original_tetra_material_ids()
    1869                 :            : {
    1870                 :          0 :     return _original_tet_material_ids;
    1871                 :            : }    
    1872                 :            : 
    1873                 :            : //! This method returns a vector of the material ids of the tetras considered in the optimization
    1874                 :          0 : vector<unsigned int> & pdt_plan::get_tetra_material_ids()
    1875                 :            : {
    1876                 :          0 :     return _tet_material_ids;
    1877                 :            : }
    1878                 :            : 
    1879                 :            : //! This method returns a vector of size _mesh_size that indicates if a tetra is considered in the optimization or not
    1880                 :          0 : vector<bool> & pdt_plan::get_tetra_to_consider()
    1881                 :            : {
    1882                 :          0 :     return _tetra_to_consider;
    1883                 :            : }
    1884                 :            : 
    1885                 :            : //! This method returns a vector of the initial placement read from file 
    1886                 :          0 : vector<Src_Abstract*> & pdt_plan::get_initial_placement()
    1887                 :            : {
    1888                 :          0 :         return _initial_placement;
    1889                 :            : }
    1890                 :            : 
    1891                 :            : // //! This method returns a vector of the cut-end fibers initial directions read from a file
    1892                 :            : // vector<SP_Point> & pdt_plan::get_cut_end_dirs()
    1893                 :            : // {
    1894                 :            : //      return _cut_end_dirs;
    1895                 :            : // }
    1896                 :            : 
    1897                 :            : // //! This method returns a vector of the cut-end fibers radii read from file 
    1898                 :            : // vector<double> & pdt_plan::get_cut_end_radii()
    1899                 :            : // {
    1900                 :            : //      return _cut_end_radii;
    1901                 :            : // }
    1902                 :            : 
    1903                 :            : // //! This method returns a vector of the cut-end fibers NA read from file 
    1904                 :            : // vector<double> & pdt_plan::get_cut_end_na()
    1905                 :            : // {
    1906                 :            : //      return _cut_end_na;
    1907                 :            : // }
    1908                 :            : 
    1909                 :            : //! This method returns a vector of indices of tumor elements in the mesh after ignoring unnecessary tetra
    1910                 :          0 : vector<unsigned int> & pdt_plan::get_tumor_indices()
    1911                 :            : {
    1912                 :          0 :     return _tumor_indices;
    1913                 :            : }
    1914                 :            : 
    1915                 :            : //! This method returns the number of photon packets used in FullMonte
    1916                 :         10 : double pdt_plan::get_num_packets()
    1917                 :            : {
    1918                 :         10 :     return _num_packets;
    1919                 :            : }
    1920                 :            : 
    1921                 :            : //! This method returns the total energy used which is read from the params file
    1922                 :         22 : double pdt_plan::get_total_energy()
    1923                 :            : {
    1924                 :         22 :     return _total_energy;
    1925                 :            : }
    1926                 :            : 
    1927                 :            : //! This method returns the index of the tumor region (usually it is the highest index of the material indices in the mesh)
    1928                 :          4 : unordered_set<unsigned> pdt_plan::get_tumor_region_number()
    1929                 :            : {
    1930                 :          4 :     return _tumor_region_number;
    1931                 :            : }
    1932                 :            : 
    1933                 :            : //! This method returns the mesh size in terms of number of tetrahedra 
    1934                 :          0 : unsigned int pdt_plan::get_mesh_size()
    1935                 :            : {
    1936                 :          0 :     return _mesh_size;
    1937                 :            : }
    1938                 :            : 
    1939                 :            : //! This method returns the mesh size after removing the unnecessary tetra
    1940                 :          0 : unsigned int pdt_plan::get_actual_mesh_size()
    1941                 :            : {
    1942                 :          0 :     return _actual_mesh_size;
    1943                 :            : }
    1944                 :            : 
    1945                 :            : //! This method returns the number of sources to place
    1946                 :          0 : unsigned int pdt_plan::get_num_sources()
    1947                 :            : {
    1948                 :          0 :     return _num_sources;
    1949                 :            : }
    1950                 :            : 
    1951                 :            : //! This method returns a pointer to the placer object
    1952                 :          0 : src_placer* pdt_plan::get_placer()
    1953                 :            : {
    1954                 :          0 :     return _placer;
    1955                 :            : }
    1956                 :            : 
    1957                 :            : //! This method returns a pointer to the mesh 
    1958                 :         23 : mesh* pdt_plan::get_mesh()
    1959                 :            : {
    1960                 :         23 :     return _data;
    1961                 :            : }
    1962                 :            : 
    1963                 :            : //! This method returns the tumor volume. It should be called after the compute_original_tetra_data function
    1964                 :          0 : double pdt_plan::get_tumor_volume()
    1965                 :            : {
    1966                 :          0 :     return _tumor_volume;
    1967                 :            : }
    1968                 :            : 
    1969                 :            : //! This method returns the tumor weight. 
    1970                 :         15 : double pdt_plan::get_tumor_weight()
    1971                 :            : {
    1972                 :         15 :         return _tumor_weight; 
    1973                 :            : }
    1974                 :            : 
    1975                 :            : //! This method returns the mesh volume. It should be called after the compute_original_tetra_data function
    1976                 :          0 : double pdt_plan::get_total_mesh_volume()
    1977                 :            : {
    1978                 :          0 :     return _total_mesh_volume;
    1979                 :            : }
    1980                 :            : 
    1981                 :            : //! This method returns the center of the rectangular volume containing the tumor volume (NOT the centroid)
    1982                 :          0 : SP_Point pdt_plan::get_tumor_center() 
    1983                 :            : {
    1984                 :          0 :     return _tumor_center;
    1985                 :            : }
    1986                 :            : 
    1987                 :            : //! This method returns the total power from the optimizer
    1988                 :          0 : double pdt_plan::get_total_power()
    1989                 :            : {
    1990                 :          0 :     return _total_power;
    1991                 :            : }
    1992                 :            : 
    1993                 :            : //! This method returns the a map from tissue names to their IDs (1-based index)
    1994                 :          0 : unordered_map<unsigned, string>& pdt_plan::get_tissue_names() 
    1995                 :            : {
    1996                 :          0 :     return _tissue_names_ids;
    1997                 :            : }
    1998                 :            : 
    1999                 :            : //! This method returns the number of sources with non-zero power after optimization
    2000                 :          0 : unsigned int pdt_plan::get_num_non_zero_sources()
    2001                 :            : {
    2002                 :          0 :     return _num_non_zero_sources;
    2003                 :            : }
    2004                 :            : 
    2005                 :            : //! This method returns the source type inputted by the user
    2006                 :          0 : string pdt_plan::get_source_type()
    2007                 :            : {
    2008                 :          0 :     return _source_type;
    2009                 :            : }
    2010                 :            : 
    2011                 :            : //! This method returns the placement type inputted by the user
    2012                 :          5 : string pdt_plan::get_placement_type()
    2013                 :            : {
    2014                 :          5 :     return _placement_type;
    2015                 :            : }
    2016                 :            : 
    2017                 :            : //! This method returns the file path and name of the fullmonte vtk output that is read from the input parameters to pdt-space
    2018                 :          0 : string pdt_plan::get_fullmonte_output_file()
    2019                 :            : {
    2020                 :          0 :     return _fullmonte_output_file;
    2021                 :            : }
    2022                 :            : 
    2023                 :            : //! This method returns the wavelength inputted by the user 
    2024                 :          0 : string pdt_plan::get_wavelength()
    2025                 :            : {
    2026                 :          0 :     return _wavelength;
    2027                 :            : }
    2028                 :            : 
    2029                 :            : //! This method returns a boolean variable indicating whether the sources emission profiles are tailored or not depending on the user input
    2030                 :          0 : bool pdt_plan::is_source_tailored()
    2031                 :            : {
    2032                 :          0 :     return _source_tailored;
    2033                 :            : }
    2034                 :            : 
    2035                 :            : //! This method returns the number of iterations needed to converge to target tumor v100
    2036                 :         11 : int pdt_plan::get_num_iterations()
    2037                 :            : {
    2038                 :         11 :         return _num_lp_iterations; 
    2039                 :            : }
    2040                 :            : 
    2041                 :            : //! This method indicates if we are running regression tests or not
    2042                 :          0 : bool pdt_plan::is_running_tests() 
    2043                 :            : {
    2044                 :          0 :         return _running_tests;
    2045                 :            : }
    2046                 :            : 
    2047                 :            : //! This method returns the parameters file
    2048                 :          4 : string pdt_plan::get_params_file()
    2049                 :            : {
    2050                 :          4 :     return _params_file;
    2051                 :            : }
    2052                 :            : 
    2053                 :            : // Setters
    2054                 :            : //! This method resizes the _p_max vector to _num_sources and initializes it with p_value
    2055                 :         16 : void pdt_plan::set_p_max_vector(double p_value)
    2056                 :            : {
    2057                 :         16 :     _p_max.clear();
    2058                 :         16 :     _p_max.resize(_num_sources);
    2059                 :         16 :     fill(_p_max.begin(), _p_max.end(), p_value);
    2060                 :            : 
    2061                 :         16 :     _is_p_max_set = true;
    2062                 :         16 : }
    2063                 :            : 
    2064                 :            : //! This method sets the number of packets used to run FM
    2065                 :          4 : void pdt_plan::set_num_packets(double num_packets)
    2066                 :            : {
    2067                 :          4 :     _num_packets = num_packets;
    2068                 :          4 : }
    2069                 :            : 
    2070                 :            : //! This method sets the tumor weight to a new weight 
    2071                 :         11 : void pdt_plan::set_tumor_weight(double weight)
    2072                 :            : {
    2073                 :         11 :         _tumor_weight = weight; 
    2074                 :         11 : }
    2075                 :            : 
    2076                 :            : //! This method enables or disables varying the tumor weight to converge to 98%
    2077                 :          2 : void pdt_plan::vary_tumor_weight(bool vary)
    2078                 :            : {
    2079                 :          2 :         _vary_tumor_weight = vary;
    2080                 :          2 : }
    2081                 :            : 
    2082                 :            : //! This method sets the target tumor v100
    2083                 :          4 : void pdt_plan::set_target_tumor_v100(double target)
    2084                 :            : {
    2085                 :          4 :         _target_tumor_v100 = target;
    2086                 :          4 : }
    2087                 :            : 
    2088                 :            : //! This method sets the upper and lower bounds on the target tumor v100 desired 
    2089                 :          4 : void pdt_plan::set_tumor_v100_bounds(double lower, double upper)
    2090                 :            : {
    2091                 :          4 :         _lower_tumor_v100_bound = lower;
    2092                 :          4 :         _upper_tumor_v100_bound = upper;
    2093                 :          4 : }
    2094                 :            : 
    2095                 :            : //! This method sets the maximum number of lp iterations to terminate early
    2096                 :          4 : void pdt_plan::set_max_num_lp_iterations(int num)
    2097                 :            : {
    2098                 :          4 :         _max_num_lp_iterations = num; 
    2099                 :          4 : }
    2100                 :            : 
    2101                 :            : //! This method sets whether early termination of SA engine is required
    2102                 :          1 : void pdt_plan::set_sa_early_termination()
    2103                 :            : {
    2104                 :          1 :     _sa_engine_early_terminate = true;
    2105                 :          1 : }
    2106                 :            : 
    2107                 :            : 
    2108                 :            : // Helper Methods
    2109                 :            : 
    2110                 :            : //! This method reads the parameters file and initializes the corresponding parameters map
    2111                 :         36 : void pdt_plan::read_params_file()
    2112                 :            : {
    2113 [ +  - ][ +  - ]:         36 :     _params_to_values_map = _fparser->get_params();
    2114                 :            : 
    2115                 :            :     // Reading if data is read from the vtk file
    2116 [ +  - ][ +  - ]:         72 :     string read_vtk = _params_to_values_map["read_vtk"];
                 [ +  - ]
    2117                 :         36 :     transform(read_vtk.begin(), read_vtk.end(), read_vtk.begin(), ::tolower);
    2118 [ +  - ][ +  + ]:         36 :     if (read_vtk == "true") 
    2119                 :         13 :         _read_data_from_vtk = true;
    2120                 :            :     else 
    2121                 :         23 :         _read_data_from_vtk = false;
    2122                 :            : 
    2123                 :            :     // Setting the wavelength
    2124 [ +  - ][ +  - ]:         72 :     string wavelength = _params_to_values_map["wavelength"];
                 [ +  - ]
    2125                 :            : 
    2126         [ +  + ]:        144 :     for (char const & c : wavelength) {
    2127         [ -  + ]:        108 :         if (std::isdigit(c) == 0) {
    2128         [ #  # ]:          0 :             fprintf(stderr, "\033[1;%dmInvalid wavelength: %s. Exiting...\033[0m\n", 31, wavelength);
    2129                 :          0 :             exit(-1);
    2130                 :            :         }
    2131                 :            :     }
    2132         [ +  - ]:         72 :     stringstream wavelength_ss;
    2133 [ +  - ][ +  - ]:         36 :     wavelength_ss << wavelength << "nm";
    2134         [ +  - ]:         36 :     _wavelength = wavelength_ss.str();
    2135                 :            : 
    2136                 :            :     // Reading the source and placement types
    2137         [ +  - ]:         36 :     _source_type = _fparser->get_source_type();
    2138                 :         36 :     transform(_source_type.begin(), _source_type.end(), _source_type.begin(), ::tolower);
    2139 [ +  - ][ -  + ]:         36 :     if (_source_type == "mixed") {
    2140         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmWarning: Mixed source types.\033[0m\n", 33);
    2141 [ +  - ][ +  + ]:         36 :     } else if (_source_type != "point" && _source_type != "line" && _source_type != "tailored" && _source_type != "cut-end") {
         [ +  - ][ +  + ]
         [ +  - ][ -  + ]
         [ #  # ][ #  # ]
                 [ -  + ]
    2142         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmUnsupported source type: %s. Exiting...\033[0m\n", 31, _source_type);
    2143                 :          0 :         exit(-1);
    2144                 :            :     }
    2145                 :            : 
    2146 [ +  - ][ +  - ]:         36 :     _placement_type = _params_to_values_map["placement_type"]; // tells whether we are using fixed or virtual sources
                 [ +  - ]
    2147                 :         36 :     transform(_placement_type.begin(), _placement_type.end(), _placement_type.begin(), ::tolower);
    2148 [ +  - ][ -  + ]:         36 :     if (_placement_type != "fixed" && _placement_type != "virtual" && _placement_type != "sa") {
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ -  + ]
    2149                 :            :     
    2150         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmUnsupported placement type: %s. Exiting...\033[0m\n", 31, _placement_type);
    2151                 :          0 :         exit(-1);
    2152                 :            :     }
    2153                 :            : 
    2154                 :            :     // set if tailored
    2155                 :         36 :     _source_tailored = _fparser->get_source_tailored();
    2156                 :            : 
    2157                 :            :     // TODO: check if empty 
    2158                 :            :     // _data_name = _params_to_values_map["DATA_NAME"];
    2159                 :            :     // _data_dir_path = _params_to_values_map["DATA_DIR"];
    2160                 :            :         // _optical_prop_file = _params_to_values_map["OPTICAL_FILE"]; 
    2161 [ +  - ][ +  - ]:         36 :     _mesh_file = _params_to_values_map["mesh_file"];
                 [ +  - ]
    2162                 :            : 
    2163                 :            :     // TODO: check if empty 
    2164 [ +  - ][ +  - ]:        113 :     _fullmonte_output_file = _params_to_values_map.find("fm_output_file") != _params_to_values_map.end() ?
         [ +  + ][ +  - ]
           [ +  -  +  + ]
                 [ #  # ]
    2165 [ +  - ][ +  - ]:         77 :                         _params_to_values_map["fm_output_file"] : "";
         [ +  + ][ +  + ]
           [ #  #  #  # ]
    2166                 :            :     
    2167                 :            :     // TODO: check if empty 
    2168 [ +  - ][ +  - ]:        108 :     _final_distribution_output_file = _params_to_values_map.find("final_dist_file") != _params_to_values_map.end() ?
         [ -  + ][ #  # ]
           [ +  -  +  - ]
                 [ #  # ]
    2169 [ #  # ][ #  # ]:         72 :                         _params_to_values_map["final_dist_file"] : "";
         [ -  + ][ -  + ]
           [ #  #  #  # ]
    2170                 :            :    
    2171                 :            :     // TODO: check if these are not numbers 
    2172 [ +  - ][ +  - ]:         36 :     _num_packets = atof(_params_to_values_map["num_packets"].c_str());
    2173 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("pnf") != _params_to_values_map.end()) {
                 [ +  - ]
    2174 [ +  - ][ +  - ]:         36 :         _total_energy = atof(_params_to_values_map["pnf"].c_str());
    2175                 :            :     } else {
    2176                 :          0 :         _total_energy = 4e11;
    2177                 :            :     }
    2178 [ +  - ][ +  - ]:         36 :     _tumor_weight = atof(_params_to_values_map["tumor_weight"].c_str());
    2179                 :            : 
    2180         [ +  - ]:         36 :     _placer->set_num_packets(_num_packets);
    2181         [ +  - ]:         36 :     _placer->set_total_energy(_total_energy);
    2182                 :            : 
    2183 [ +  - ][ +  - ]:         36 :         if (_params_to_values_map.find("run_tests") != _params_to_values_map.end()) {
                 [ +  - ]
    2184 [ +  - ][ +  - ]:         72 :                 string run_tests = _params_to_values_map["run_tests"];
                 [ +  - ]
    2185                 :         36 :         transform(run_tests.begin(), run_tests.end(), run_tests.begin(), ::tolower);
    2186                 :            : 
    2187         [ +  - ]:         36 :                 _running_tests = run_tests == "true" ? true : false;
    2188                 :            :         }
    2189                 :            :         else {
    2190                 :          0 :                 _running_tests = false;
    2191                 :            :         }
    2192                 :            : 
    2193 [ +  - ][ +  - ]:         36 :         if (_params_to_values_map.find("USE_CUDA") != _params_to_values_map.end()) {
                 [ -  + ]
    2194 [ #  # ][ #  # ]:          0 :                 string use_cuda = _params_to_values_map["USE_CUDA"];
                 [ #  # ]
    2195                 :          0 :         transform(use_cuda.begin(), use_cuda.end(), use_cuda.begin(), ::tolower);
    2196                 :            : 
    2197         [ #  # ]:          0 :                 _use_cuda = use_cuda == "true" ? true : false;
    2198                 :            : 
    2199         [ #  # ]:          0 :         if (!USE_CUDA && _use_cuda) {
    2200         [ #  # ]:          0 :             fprintf(stderr, "\033[1;%dmCuda not supported, using CPU instead...\033[0m\n", 31);
    2201                 :          0 :             _use_cuda = false;
    2202                 :            :         }
    2203                 :            :         }
    2204                 :            :         else {
    2205                 :         36 :                 _use_cuda = false;
    2206                 :            :         }
    2207                 :            : 
    2208         [ +  + ]:         36 :     if (_running_tests) _file_path_override = true;
    2209                 :         31 :     else _file_path_override = false;
    2210                 :            : 
    2211 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("override_file_path") != _params_to_values_map.end()) {
                 [ +  - ]
    2212 [ +  - ][ +  - ]:         72 :                 string file_path_override = _params_to_values_map["override_file_path"];
                 [ +  - ]
    2213                 :         36 :         transform(file_path_override.begin(), file_path_override.end(), file_path_override.begin(), ::tolower);
    2214                 :            : 
    2215         [ +  - ]:         36 :                 _file_path_override = file_path_override == "true" ? true : false;
    2216                 :            :         }
    2217                 :            : 
    2218                 :         36 :     _run_sa = _fparser->get_run_sa();
    2219                 :            : 
    2220 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("run_rt_feedback") != _params_to_values_map.end()) {
                 [ +  - ]
    2221 [ +  - ][ +  - ]:         72 :                 string run_rt_feedback_s = _params_to_values_map["run_rt_feedback"];
                 [ +  - ]
    2222                 :         36 :         transform(run_rt_feedback_s.begin(), run_rt_feedback_s.end(), run_rt_feedback_s.begin(), ::tolower);
    2223                 :            : 
    2224         [ +  - ]:         36 :                 _run_rt_feedback = run_rt_feedback_s == "true" ? true : false;
    2225                 :            :         }
    2226                 :            :         else {
    2227                 :          0 :                 _run_rt_feedback = false;
    2228                 :            :         }
    2229                 :            : 
    2230 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("vary_tumor_weight") != _params_to_values_map.end()) {
                 [ +  - ]
    2231 [ +  - ][ +  - ]:         72 :         string vary = _params_to_values_map["vary_tumor_weight"];
                 [ +  - ]
    2232                 :         36 :         transform(vary.begin(), vary.end(), vary.begin(), ::tolower);
    2233                 :            : 
    2234         [ +  - ]:         36 :         _vary_tumor_weight = (vary == "true") ? true : false;
    2235                 :            :     } else { // default
    2236                 :          0 :         _vary_tumor_weight = true;
    2237                 :            :     }
    2238                 :            : 
    2239                 :            :     // TODO: check if these are numbers
    2240 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("target_tumor_cov") != _params_to_values_map.end()) {
                 [ +  - ]
    2241 [ +  - ][ +  - ]:         36 :         _target_tumor_v100 = std::abs(atof(_params_to_values_map["target_tumor_cov"].c_str()));
    2242                 :            :     } else { // default 
    2243                 :          0 :         _target_tumor_v100 = 98.0;
    2244                 :            :     }
    2245                 :            : 
    2246 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("rt_feedback_seed") != _params_to_values_map.end()) {
                 [ -  + ]
    2247 [ #  # ][ #  # ]:          0 :         _rt_feedback_rand_seed = atoi(_params_to_values_map["rt_feedback_seed"].c_str());
    2248                 :            :     } else {
    2249                 :         36 :         _rt_feedback_rand_seed = 5; // Default value
    2250                 :            :     }
    2251                 :            : 
    2252 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("rt_feedback_det_rad") != _params_to_values_map.end()) {
                 [ -  + ]
    2253 [ #  # ][ #  # ]:          0 :         _rt_feedback_detector_radius = static_cast<float>(atof(_params_to_values_map["rt_feedback_det_rad"].c_str()));
    2254                 :            :     } else {
    2255                 :         36 :         _rt_feedback_detector_radius = 0.5f; // Default value
    2256                 :            :     }
    2257                 :            :     
    2258 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("rt_feedback_det_na") != _params_to_values_map.end()) {
                 [ -  + ]
    2259 [ #  # ][ #  # ]:          0 :         _rt_feedback_detector_na = static_cast<float>(atof(_params_to_values_map["rt_feedback_det_na"].c_str()));
    2260                 :            :     } else {
    2261                 :         36 :         _rt_feedback_detector_na = 0.5f; // Default value
    2262                 :            :     }
    2263                 :            : 
    2264 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("rt_feedback_lut_file") != _params_to_values_map.end()) {
                 [ -  + ]
    2265 [ #  # ][ #  # ]:          0 :         _rt_feedback_lut_file_name = _params_to_values_map["rt_feedback_lut_file"]; 
                 [ #  # ]
    2266                 :            :     } else {
    2267         [ +  - ]:         36 :         _rt_feedback_lut_file_name = "";
    2268                 :            :     }
    2269                 :            : 
    2270                 :            :     // check if a number
    2271 [ +  - ][ +  - ]:         36 :     if (_params_to_values_map.find("rt_feedback_lut_size") != _params_to_values_map.end()) {
                 [ -  + ]
    2272 [ #  # ][ #  # ]:          0 :         _rt_feedback_lut_size = atoi(_params_to_values_map["rt_feedback_lut_size"].c_str()); 
    2273                 :            :     } else {
    2274                 :         36 :         _rt_feedback_lut_size = 50; 
    2275                 :            :     }
    2276                 :         36 : }
    2277                 :            : 
    2278                 :            : /**
    2279                 :            :  * This method runs FullMonteSW with one source at a time with unit power to generate the dose matrix
    2280                 :            :  * @param source [in]: The sources to run FullMonte with 
    2281                 :            :  */
    2282                 :         12 : void pdt_plan::run_fullmonte_multiple_sources(const vector<Src_Abstract*> & sources)
    2283                 :            : {
    2284 [ +  - ][ +  - ]:         24 :     profiler fm_engine("Running FM one source at a time"); 
    2285                 :            :     // create a vector of sources 
    2286         [ +  - ]:         24 :     std::vector<Source::Abstract*> srcs(sources.size()); 
    2287         [ +  + ]:         48 :     for (unsigned i = 0; i < sources.size(); i++) {
    2288         [ +  - ]:         36 :         srcs[i] = sources[i]->get_fm_source();
    2289                 :            :     }
    2290                 :            : 
    2291 [ +  - ][ +  + ]:         42 :     for (unsigned i = 0; i < _placer->get_source_element_fluence_matrix().size(); i++)
    2292         [ +  - ]:         30 :         _placer->get_source_element_fluence_matrix()[i].clear(); 
    2293         [ +  - ]:         12 :     _placer->get_source_element_fluence_matrix().clear(); 
    2294                 :            : 
    2295                 :            :     // default 10^6 packets
    2296 [ +  - ][ +  - ]:         24 :     _fullmonte_engine->run_multiple_fullmonte(srcs, _placer->get_source_element_fluence_matrix(), false, "", 
    2297         [ +  - ]:         12 :                         static_cast<unsigned>(_num_packets)); 
    2298         [ +  - ]:         12 :         _num_sources = static_cast<unsigned>(_placer->get_source_element_fluence_matrix().size()); 
    2299                 :            : 
    2300                 :            :     // memory cleanup
    2301         [ +  + ]:         48 :     for (unsigned i = 0; i < srcs.size(); i++) {
    2302 [ +  - ][ +  - ]:         36 :         if (srcs[i]) delete srcs[i];
    2303                 :            :     }
    2304                 :         12 : }
    2305                 :            : 
    2306                 :            : //! This method runs a composite source of all sources in _initial_placement vector
    2307                 :          0 : void pdt_plan::run_fullmonte_composite_source(const vector<float> & src_powers)
    2308                 :            : {
    2309 [ #  # ][ #  # ]:          0 :     profiler fm_engine("Running FM composite source"); 
    2310                 :            : 
    2311 [ #  # ][ #  # ]:          0 :     if (_final_distribution_output_file == "") {
    2312                 :          0 :         return;
    2313                 :            :     }
    2314                 :            : 
    2315                 :            :     // create a vector of sources 
    2316 [ #  # ][ #  # ]:          0 :     std::vector<Source::Abstract*> srcs(_initial_placement.size()); 
    2317         [ #  # ]:          0 :     for (unsigned i = 0; i < _initial_placement.size(); i++) {
    2318         [ #  # ]:          0 :         srcs[i] = _initial_placement[i]->get_fm_source();
    2319                 :          0 :         srcs[i]->power(src_powers[i]);
    2320                 :            :     }
    2321                 :            : 
    2322                 :            :     // default 10^6 packets
    2323                 :          0 :         vector<float> fluence; 
    2324         [ #  # ]:          0 :     _fullmonte_engine->run_composite_fullmonte(srcs, fluence, true, _final_distribution_output_file, 
    2325         [ #  # ]:          0 :                         static_cast<unsigned>(_num_packets*static_cast<double>(_initial_placement.size()))); 
    2326                 :            : 
    2327                 :            :     // memory cleanup
    2328         [ #  # ]:          0 :     for (unsigned i = 0; i < srcs.size(); i++) {
    2329 [ #  # ][ #  # ]:          0 :         if (srcs[i]) delete srcs[i];
    2330                 :            :     }
    2331                 :            : }
    2332                 :            : 
    2333                 :            : //!
    2334                 :            : //! This method generates a vector of vectors of tetra IDs that enclose the points through which each source passes through,
    2335                 :            : //! Then it tags the tetra of each source with a new material and generates a new vtk mesh
    2336                 :            : //!
    2337                 :          0 : void pdt_plan::mesh_sources() 
    2338                 :            : {
    2339                 :            :  // TODO This is not used anymore
    2340                 :            :     /* 
    2341                 :            :      * Step 1: Discretize source
    2342                 :            :      *         Each source is discretized into 0.2mm elements, and the enclosing tetra is found for each point
    2343                 :            :      *
    2344                 :            :      */
    2345                 :          0 :     std::vector<std::vector<unsigned>> enclosing_tetra_ids;
    2346         [ #  # ]:          0 :     for (auto src : _initial_placement) 
    2347                 :            :     {
    2348                 :          0 :         std::vector<unsigned> src_enclosing_tetra;
    2349                 :            : 
    2350         [ #  # ]:          0 :         if (src->get_type() == Src_Abstract::Point) // point or cut-end source
    2351                 :            :         {
    2352                 :          0 :             Src_Point* src_p = (Src_Point*)src;
    2353 [ #  # ][ #  # ]:          0 :             src_enclosing_tetra.push_back(_fullmonte_engine->get_tetra_enclosing_point(src_p->get_xcoord(), 
         [ #  # ][ #  # ]
    2354         [ #  # ]:          0 :                     src_p->get_ycoord(), src_p->get_zcoord()));
    2355                 :            :         }
    2356         [ #  # ]:          0 :         else if (src->get_type() == Src_Abstract::CutEnd) // point or cut-end source
    2357                 :            :         {
    2358                 :          0 :             SP_Point endp = ((Src_CutEnd*)src)->get_endpoint();
    2359 [ #  # ][ #  # ]:          0 :             src_enclosing_tetra.push_back(_fullmonte_engine->get_tetra_enclosing_point(endp.get_xcoord(), 
         [ #  # ][ #  # ]
    2360         [ #  # ]:          0 :                     endp.get_ycoord(), endp.get_zcoord()));
    2361                 :            :         }
    2362                 :            :         else // line source
    2363                 :            :         {
    2364                 :          0 :             Src_Line* src_l = (Src_Line*)src;
    2365         [ #  # ]:          0 :             double src_length = src_l->get_src_length();
    2366                 :            : 
    2367                 :            :             // Get direction vector
    2368 [ #  # ][ #  # ]:          0 :             SP_Point src_dir = (src_l->end1 - src_l->end0)/src_length;
    2369                 :            : 
    2370                 :          0 :             double elem_length = 0.2; // mm
    2371                 :            :             
    2372                 :          0 :             unsigned num_pts = static_cast<unsigned>(floor(src_length/elem_length) + 1);
    2373                 :            :             
    2374                 :            :             
    2375                 :            :             // Get num_rand_perp random perpendicular vectors and store them
    2376                 :          0 :             unsigned num_rand_perp = 10; 
    2377         [ #  # ]:          0 :             std::vector<SP_Point> rand_perp(num_rand_perp); 
    2378         [ #  # ]:          0 :             for (unsigned i = 0; i < num_rand_perp; i++) 
    2379                 :            :             {
    2380         [ #  # ]:          0 :                 rand_perp[i] = src_l->get_rand_perpendicular_vector(); 
    2381                 :            :             }
    2382                 :            : 
    2383                 :            :             // Get number of points in perpendicular direction
    2384                 :          0 :             double cylinder_radius = 2.5; // assuming 5mm cylindrical diffuser diameter
    2385                 :          0 :             unsigned num_pts_diameter = static_cast<unsigned>(floor(cylinder_radius/elem_length) + 1);
    2386                 :            : 
    2387                 :            :             // find enclosing tetra
    2388         [ #  # ]:          0 :             for (unsigned i = 0; i < num_pts; i++) 
    2389                 :            :             {
    2390 [ #  # ][ #  # ]:          0 :                 SP_Point curr_pt = src_l->end0 + (elem_length*i*src_dir); 
    2391                 :            :                
    2392 [ #  # ][ #  # ]:          0 :                 src_enclosing_tetra.push_back(_fullmonte_engine->get_tetra_enclosing_point(curr_pt.get_xcoord(), 
         [ #  # ][ #  # ]
    2393         [ #  # ]:          0 :                        curr_pt.get_ycoord(), curr_pt.get_zcoord()));
    2394                 :            : 
    2395                 :            :                 // Get num_rand_perp random perpendicular vectors and find enclosing tetra
    2396         [ #  # ]:          0 :                 for (unsigned j = 0; j < num_rand_perp; j++) 
    2397                 :            :                 {
    2398         [ #  # ]:          0 :                     for (unsigned k = 0; k < num_pts_diameter; k++)
    2399                 :            :                     {
    2400 [ #  # ][ #  # ]:          0 :                         SP_Point curr_pt_polar = curr_pt + (elem_length*k*rand_perp[j]); 
    2401                 :            :                        
    2402 [ #  # ][ #  # ]:          0 :                         src_enclosing_tetra.push_back(_fullmonte_engine->get_tetra_enclosing_point(curr_pt_polar.get_xcoord(), 
         [ #  # ][ #  # ]
    2403         [ #  # ]:          0 :                                curr_pt_polar.get_ycoord(), curr_pt_polar.get_zcoord()));
    2404                 :            :                     }
    2405                 :            :                 }
    2406                 :            :             }
    2407                 :            :         }
    2408         [ #  # ]:          0 :         enclosing_tetra_ids.push_back(src_enclosing_tetra); 
    2409                 :            :     }
    2410                 :            : 
    2411                 :            : 
    2412                 :            :     /*
    2413                 :            :      * 
    2414                 :            :      * Step 2: Tag each tetra in enclosing_tetra_ids with new region IDs (only for FullMonte, as it does not affect PDT-SPACE) 
    2415                 :            :      */
    2416                 :            :     // Get FullMonte Regions 
    2417         [ #  # ]:          0 :     WrapperFullMonteSW* fullmonte_engine_clone = _fullmonte_engine->clone();
    2418         [ #  # ]:          0 :     Partition* fm_regions = fullmonte_engine_clone->get_tetra_mesh()->regions(); 
    2419                 :            : 
    2420 [ #  # ][ #  # ]:          0 :     cout << "Number of regions is: " << fm_regions->count() << endl;
         [ #  # ][ #  # ]
    2421                 :            : 
    2422 [ #  # ][ #  # ]:          0 :     cout << "Number of sources: " << _initial_placement.size() << endl;
                 [ #  # ]
    2423         [ #  # ]:          0 :     for (unsigned i = 0; i < enclosing_tetra_ids.size(); i++)
    2424                 :            :     {
    2425         [ #  # ]:          0 :         for (unsigned j = 0; j < enclosing_tetra_ids[i].size(); j++) 
    2426                 :            :         {
    2427         [ #  # ]:          0 :             fm_regions->assign(enclosing_tetra_ids[i][j], _max_region_index + 1);//+ i);
    2428                 :            :         }
    2429         [ #  # ]:          0 :         fm_regions->recountRegions();   
    2430                 :            :     }
    2431 [ #  # ][ #  # ]:          0 :     cout << "Number of regions after: " << fm_regions->count() << endl;
         [ #  # ][ #  # ]
    2432 [ #  # ][ #  # ]:          0 :     cout << "Number of materials: " << fullmonte_engine_clone->get_material_set()->size() << endl;
                 [ #  # ]
    2433                 :            : 
    2434 [ #  # ][ #  # ]:          0 :     fullmonte_engine_clone->write_mesh("logs/rt_feedback_results/test_src_mesh.vtk");
    2435                 :          0 :     std::exit(0);
    2436                 :            : }
    2437                 :            : 
    2438                 :            : //! This method prints the tissue names with their v100
    2439                 :          5 : void pdt_plan::print_v100s(const std::vector<double>& v100s)
    2440                 :            : {
    2441                 :          5 :     std::cout << std::endl;
    2442                 :          5 :     double tumor_volume = 0;
    2443                 :          5 :     double coverage_volume = 0;
    2444         [ +  + ]:         30 :     for (unsigned i = 1; i <= _d_max_tissues.size(); i++) 
    2445                 :            :     {
    2446 [ +  - ][ +  + ]:         25 :         if (_tumor_region_number.find(i) != _tumor_region_number.end()) {
    2447         [ +  - ]:         10 :             stringstream tumor;
    2448 [ +  - ][ +  - ]:          5 :             tumor << _tissue_names_ids[i];
    2449 [ +  - ][ +  - ]:          5 :             std::cout << std::left << std::setw(20) << std::setfill(' ') << tumor.str();
         [ +  - ][ +  - ]
                 [ +  - ]
    2450         [ +  - ]:          5 :             tumor_volume += _material_volumes[i];
    2451         [ +  - ]:          5 :             coverage_volume += v100s[i - 1] * _material_volumes[i];
    2452 [ +  - ][ +  - ]:          5 :             std::cout << std::fixed << std::setprecision(5) << " V100 is: " << v100s[i - 1] << " ";
         [ +  - ][ +  - ]
                 [ +  - ]
    2453         [ +  - ]:          5 :             std::cout << "%";
    2454                 :            :         } else {
    2455 [ +  - ][ +  - ]:         20 :             std::cout << std::left << std::setw(20) << std::setfill(' ') << _tissue_names_ids[i];
         [ +  - ][ +  - ]
                 [ +  - ]
    2456 [ +  - ][ +  - ]:         20 :             std::cout << std::fixed << std::setprecision(5) << " V100 is: " << v100s[i - 1] << " ";
         [ +  - ][ +  - ]
                 [ +  - ]
    2457         [ +  - ]:         20 :             std::cout << "unit_length^3";
    2458                 :            :         }
    2459                 :            : 
    2460         [ +  - ]:         25 :         std::cout << std::endl;
    2461                 :            :     }
    2462                 :            : 
    2463                 :          5 :     std::cout << std::left << std::setw(20) << std::setfill(' ') << "Total Tumor" 
    2464                 :         10 :                 << std::fixed << std::setprecision(5) << " V100 is: " << coverage_volume/tumor_volume << " %\n";
    2465                 :            : 
    2466                 :          5 :     std::cout << std::endl;
    2467                 :          5 : }
    2468                 :            : 
    2469                 :            : //! This method appends all directory paths with PDT_SPACE_SRC_DIR
    2470                 :         36 : void pdt_plan::override_file_paths () {
    2471         [ +  - ]:         36 :     if (_mesh_file != "")
    2472 [ +  - ][ +  - ]:         36 :         _mesh_file = string(PDT_SPACE_SRC_DIR) + _mesh_file;
    2473         [ +  + ]:         36 :     if (_fullmonte_output_file != "")
    2474 [ +  - ][ +  - ]:          5 :         _fullmonte_output_file = string(PDT_SPACE_SRC_DIR) + _fullmonte_output_file;
    2475         [ -  + ]:         36 :     if (_final_distribution_output_file != "")
    2476 [ #  # ][ #  # ]:          0 :         _final_distribution_output_file = string(PDT_SPACE_SRC_DIR) + _final_distribution_output_file; 
    2477         [ -  + ]:         36 :     if (_sa_engine_progress_placement_file != "")
    2478 [ #  # ][ #  # ]:          0 :         _sa_engine_progress_placement_file = string(PDT_SPACE_SRC_DIR) + _sa_engine_progress_placement_file;
    2479                 :         36 : }
    2480                 :            : 
    2481                 :            : // //!
    2482                 :            : // //! This method sets the parameters used to read mesh files
    2483                 :            : // //! @param data_dir_path: [in] Value for _data_dir_path
    2484                 :            : // //! @param data_name: [in] Name of the mesh file with no extension
    2485                 :            : // //! @param optical_prop_file: [in] Value for _optical_prop_file
    2486                 :            : // //! @param read_data_from_vtk: [in] Value for _read_data_from_vtk
    2487                 :            : // //! @param wavelength: [in] Value for _wavelength
    2488                 :            : // //!
    2489                 :            : // void pdt_plan::set_parameters (string data_dir_path, string data_name, 
    2490                 :            : //         string optical_prop_file, bool read_data_from_vtk, string wavelength) {
    2491                 :            : //     _data_dir_path = data_dir_path;
    2492                 :            : //     _data_name = data_name;
    2493                 :            : //     _optical_prop_file = optical_prop_file;
    2494                 :            : //     _read_data_from_vtk = read_data_from_vtk;
    2495                 :            : //     _wavelength = wavelength;
    2496                 :            : // }
    2497                 :            : 
    2498                 :            : //! 
    2499                 :            : //! This method acts as a setter for the _params_file and reads it
    2500                 :            : //! @param params_file: [in] Value for _params_file
    2501                 :            : //!
    2502                 :          0 : void pdt_plan::set_params_file (string params_file) {
    2503                 :          0 :     _params_file = params_file;
    2504 [ #  # ][ #  # ]:          0 :     if (_placer) delete _placer;
    2505         [ #  # ]:          0 :     _placer = new src_placer();
    2506 [ #  # ][ #  # ]:          0 :     if (_fparser) delete _fparser;
    2507 [ #  # ][ #  # ]:          0 :     _fparser = new Parser(params_file);
    2508                 :          0 :     read_params_file();
    2509                 :          0 : }
    2510                 :            : 
    2511                 :            : //!
    2512                 :            : //! This method sets the initial placement vector. This setter is not managing any memory operations on sources
    2513                 :            : //! @param placement: [in] vector of Src_Abstract* of the source placement to be set to
    2514                 :            : //!
    2515                 :          0 : void pdt_plan::set_initial_placement(vector<Src_Abstract*> placement) {
    2516                 :          0 :     _initial_placement = placement;
    2517                 :          0 : }
    2518                 :            : 
    2519                 :            : //!
    2520                 :            : //! This method sets the final source power vector. Use with caution as it sets the result vector from conv opt
    2521                 :            : //! @param source_powers: [in] vector of double of the source powers to be set to
    2522                 :            : //!
    2523                 :          0 : void pdt_plan::set_final_source_powers(vector<double> source_powers) {
    2524                 :          0 :     _final_source_powers = source_powers;
    2525                 :          0 :     _is_lp_solved = true;
    2526                 :          0 : }
    2527                 :            : 
    2528                 :            : /**
    2529                 :            :  * This method initializes an RtFeedback object and runs a realtime feedback system to recover optical properties 
    2530                 :            :  */
    2531                 :          0 : void pdt_plan::run_rt_feedback()
    2532                 :            : {
    2533         [ #  # ]:          0 :     if (!_run_rt_feedback) {
    2534                 :          0 :         return;
    2535                 :            :     }
    2536                 :            : 
    2537         [ #  # ]:          0 :     if (_tumor_region_number.size() > 1) {
    2538         [ #  # ]:          0 :         fprintf(stderr, "\033[1;%dmpdt_plan::RtFeedback does not support more than one tumor region ID. Exiting...\033[0m\n", 31); 
    2539                 :          0 :         std::exit(-1); 
    2540                 :            :     }
    2541                 :            : 
    2542                 :            :     // 1- Initialize RtFeedback object
    2543                 :            :     {
    2544 [ #  # ][ #  # ]:          0 :         profiler init_rt("Initializing RtFeedback");
    2545                 :          0 :         _rt_feedback_system = new RtFeedback(_fullmonte_engine, _initial_placement, static_cast<unsigned long long>(_num_packets), 
    2546                 :          0 :                             *(_tumor_region_number.begin()), _tumor_volume, _rt_feedback_rand_seed, _rt_feedback_detector_radius, 
    2547 [ #  # ][ #  # ]:          0 :                                     _rt_feedback_detector_na, _rt_feedback_lut_size); 
    2548                 :            :     }
    2549                 :            :     
    2550         [ #  # ]:          0 :     std::vector<float> powers(_initial_placement.size(), 1.0f); 
    2551                 :            :     { 
    2552 [ #  # ][ #  # ]:          0 :         profiler LUT_rt("Building mu_eff lookup table"); 
    2553 [ #  # ][ #  # ]:          0 :         if (_rt_feedback_lut_file_name == "") {
    2554         [ #  # ]:          0 :             _rt_feedback_system->build_lut(_initial_placement, powers); // TODO: make it an option
    2555                 :            :         } else {
    2556 [ #  # ][ #  # ]:          0 :             _rt_feedback_system->build_lut_from_file(_rt_feedback_lut_file_name, _initial_placement.size()); 
    2557                 :            :         }
    2558                 :            :     }
    2559                 :            : 
    2560                 :            : //         _rt_feedback_system->prepare_training_data(_initial_placement); // TODO: make it an option 
    2561                 :            :     {
    2562 [ #  # ][ #  # ]:          0 :         profiler prof_run_rt("Runnning RT Feedback System"); 
    2563                 :          0 :         unsigned num_simulations = 25; // 25 recoveries per tumor 
    2564                 :            : 
    2565                 :            :         // 2- Define loop criterion 
    2566         [ #  # ]:          0 :         vector<double> thresholds = get_thresholds_with_guardband(); 
    2567         [ #  # ]:          0 :         for (unsigned i = 0; i < num_simulations; i++) {
    2568                 :            : 
    2569                 :            :             // 3- While criterion not met, call recover_tumor_optical_properties, and re-optimize powers
    2570                 :            : 
    2571                 :            :             // 3- While criterion not met, call recover_tumor_optical_properties, and re-optimize powers
    2572                 :            : 
    2573                 :            : 
    2574                 :            :         
    2575                 :            :             float new_mu_a, new_mu_s; 
    2576         [ #  # ]:          0 :             _rt_feedback_system->recover_tumor_optical_properties(_initial_placement, powers, new_mu_a, new_mu_s); 
    2577                 :            : 
    2578         [ #  # ]:          0 :             run_fullmonte_multiple_sources(_initial_placement);
    2579         [ #  # ]:          0 :             remove_unnecessary_tetra(thresholds); 
    2580         [ #  # ]:          0 :             compute_tetra_data(thresholds); 
    2581         [ #  # ]:          0 :             initialize_lp(); 
    2582         [ #  # ]:          0 :             solve_lp(thresholds); 
    2583         [ #  # ]:          0 :             compute_post_opt_fluence();
    2584         [ #  # ]:          0 :             vector<double> v_alpha = compute_v_alpha_tissues(thresholds, 1.0); 
    2585                 :            : 
    2586 [ #  # ][ #  # ]:          0 :             std::cout << "Simulation: " << i << " V100s post op recovery:" << std::endl;
         [ #  # ][ #  # ]
    2587         [ #  # ]:          0 :             print_v100s(v_alpha); 
    2588 [ #  # ][ #  # ]:          0 :             std::cout << "Simulation: "  << i << " Printing source powers after optimization post op recovery" << std::endl;
         [ #  # ][ #  # ]
    2589         [ #  # ]:          0 :             for (unsigned int j = 0; j < _final_source_powers.size(); j++) {
    2590 [ #  # ][ #  # ]:          0 :                 std::cout << "x[" << j << "] = " << std::fixed << std::setprecision(8) << _final_source_powers[j] << std::endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
    2591                 :            :             }
    2592                 :            :         }
    2593                 :            :     }
    2594 [ +  - ][ +  - ]:          4 : }

Generated by: LCOV version 1.12