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 : }
|