root/applications/robust/robustlib.cpp @ 1361

Revision 1357, 8.9 kB (checked in by sindj, 13 years ago)

Opravena chyba s uvolnovanim pameti pri ruseni ve for_merging poli a dalsi drobne nedostatky. JS

Line 
1#include "robustlib.h"
2
3
4double polyhedron::triangulate(bool should_integrate)
5        {
6                for(set<simplex*>::iterator t_ref = triangulation.begin();t_ref!=triangulation.end();t_ref++)
7                {
8                        delete (*t_ref);
9                }               
10                triangulation.clear();
11
12                if(should_integrate)
13                {
14                        ((toprow *)this)->probability = 0.0;
15                }               
16               
17                if(vertices.size()==1)
18                {
19                        simplex* vert_simplex = new simplex((vertex*)this);                     
20
21                        triangulation.insert(vert_simplex);
22                }
23
24                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
25                {                       
26                        for(set<simplex*>::iterator s_ref = (*child_ref)->triangulation.begin();s_ref!=(*child_ref)->triangulation.end();s_ref++)
27                        {
28                                simplex* new_simplex = new simplex((*s_ref)->vertices);                                                         
29
30                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex->vertices.insert(*vertices.begin());
31
32                                if(ret_val.second == true)
33                                {
34                                        double cur_prob = 0;
35
36                                        /*
37                                        if(should_integrate&&new_simplex->vertices.size()!=3)
38                                        {
39                                                cout << "Error: Wrong vertex count for integration!";
40                                        }
41                                        */
42
43                                        if(should_integrate)
44                                        {
45                                                cur_prob = ((toprow *)this)->integrate_simplex(new_simplex, 'S');
46
47                                                ((toprow *)this)->probability += cur_prob;
48                                        }
49
50                                        triangulation.insert(new_simplex);
51                                }
52                                else
53                                {
54                                        delete new_simplex;
55                                }
56                        }       
57                }       
58
59                if(should_integrate)
60                {
61                        return ((toprow *)this)->probability;
62                }
63                else
64                {
65                        return 0.0;
66                }
67               
68                /*
69                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
70                {                       
71                        for(list<set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++)
72                        {
73                                set<vertex*> new_simplex;
74                                new_simplex.insert((*t_ref).begin(),(*t_ref).end());                   
75
76                                for(set<vertex*>::iterator suitable_vert_ref = vertices.begin();*suitable_vert_ref<*(*t_ref).begin();suitable_vert_ref++)
77                                {
78                                        set<vertex*> suitable_simplex;
79                                        suitable_simplex.insert(new_simplex.begin(),new_simplex.end());
80                                       
81                                        suitable_simplex.insert(*suitable_vert_ref);
82
83                                        triangulation.push_back(suitable_simplex);
84
85                                        if(should_integrate)
86                                        {
87                                                ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(suitable_simplex, 'S');
88                                        }
89                                }                               
90                        }       
91                }*/             
92        }
93
94
95        double toprow::integrate_simplex(simplex* simplex, char c)
96        {
97                // cout << ((toprow*)this)->condition << endl;         
98
99                int condition_order = ((toprow*)this)->condition_order-2; // -2 by bylo, pokud chceme uniformni apriorno
100
101                int sigma_order = condition_order-my_emlig->number_of_parameters;
102
103                // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
104                // pause(0.1);
105
106                if(sigma_order >= 0)
107                {
108
109                        //cout << endl;
110                        //cout << ((toprow*)this)->condition << endl;
111                        //cout << "C:" << condition_order+2 << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters+2 << endl;
112                       
113
114                        emlig* current_emlig;
115                        simplex->clear_gammas();
116
117                        if(this->my_emlig!=NULL)
118                        {
119                                current_emlig = this->my_emlig;
120                        }
121                        else
122                        {
123                                throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!");
124                        }                                               
125
126                        toprow* as_toprow = (toprow*)this;
127
128                        vec cur_condition = as_toprow->condition_sum.get(1,as_toprow->condition_sum.size()-1);
129
130                        // cout << as_toprow->condition << endl;                       
131
132                        int dimension = simplex->vertices.size()-1;
133
134                        mat jacobian(dimension,dimension);                     
135
136                        double a_0;
137                        map<double,int> as;
138                        vertex* base_vertex;
139                        set<vertex*>::iterator base_ref = simplex->vertices.begin();
140                        bool order_correct;
141                                               
142
143                        do
144                        {
145                                order_correct = true;
146                                as.clear();
147
148                                base_vertex = (*base_ref);
149
150                       
151                                // cout << endl << "Base coords:" << base_vertex->get_coordinates() << endl;
152
153                                a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition_sum[0];
154                               
155
156                                int row_count = 0;
157                                for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++)
158                                {
159                                        if(vert_ref != base_ref)
160                                        {
161                                                vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();                                           
162
163                                                jacobian.set_row(row_count,relative_coords);
164
165                                                double new_a = -relative_coords*cur_condition;
166
167                                                if(new_a + a_0 == 0 || new_a == 0)
168                                                {
169                                                        base_ref++;
170
171                                                        if(base_ref == simplex->vertices.end()&&simplex->vertices.size()!=2)
172                                                        {
173                                                                throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!");
174                                                        }
175                                                        /*
176                                                        else if(implex->vertices.size()==2)
177                                                        {
178                                                                return relative_coords[0]*exp
179                                                        }
180                                                        */
181                                                       
182                                                        order_correct = false;                                         
183
184                                                        break;
185                                                }                                               
186                                               
187                                               
188
189                                                if(a_0<current_emlig->min_ll)
190                                                {
191                                                        current_emlig->minimal_vertex = base_vertex;
192                                                        current_emlig->min_ll = a_0;
193                                                }
194
195                                                double a_m = (*vert_ref)->get_coordinates()*cur_condition-as_toprow->condition_sum[0];
196                                                if(a_m<current_emlig->min_ll)
197                                                {
198                                                        current_emlig->minimal_vertex = (*vert_ref);
199                                                        current_emlig->min_ll = a_m;                                           
200                                                }
201
202                                                //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl;
203
204                                                // cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl;
205                                                //cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl;
206
207                                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
208                                                if(returned.second == false)
209                                                {
210                                                        (*returned.first).second++;
211                                                }
212                                                /*
213                                                else
214                                                {
215                                                        cout << "a[" << row_count << "] = " << new_a << endl;
216                                                }
217                                                */
218                                               
219                                                //as.ins(as.size(),new_a);                                                     
220                                               
221                                                row_count++;
222                                        }
223                                }
224                        }
225                        while(!order_correct);
226                       
227                        /*
228                        cout << "a_0: " << a_0 << "    ";
229                        int as_count = 1;
230                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
231                        {
232                                cout << "a_" << as_count << ": " << (*as_ref).first << "    ";
233                                as_count++;
234                        }*/
235
236                        double int_value = 0;                   
237
238                        // cout << jacobian << endl;
239
240                        double det_jacobian    = abs(det(jacobian));
241                        double correction_term = det_jacobian;                 
242                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
243                        {
244                                double multiplier = det_jacobian;                               
245                               
246                                int as1_order = (*as_ref).second;
247                               
248                                correction_term /= pow(-(*as_ref).first,as1_order);                                     
249                               
250                                current_emlig->set_correction_factors(as1_order);
251
252                                vector<double> factors;
253                                int number_of_factors = 0;                                                             
254                                for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
255                                {
256                                       
257                                        if(as2_ref!=as_ref)
258                                        {                                                                               
259                                                for(int k = 0;k<(*as2_ref).second;k++)
260                                                {
261                                                        factors.push_back((*as_ref).first-(*as2_ref).first);
262                                                }
263
264                                                multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second);
265                                                number_of_factors += (*as2_ref).second;
266                                        }
267                                        else
268                                        {
269                                                factors.push_back((*as_ref).first);
270
271                                                multiplier        /= (*as_ref).first;
272                                                number_of_factors += 1;
273                                        }
274                                       
275                                }
276                               
277                                                               
278
279                                double bracket = fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+1);
280                                                               
281                                simplex->insert_gamma(0,bracket*multiplier,a_0-(*as_ref).first);                                                               
282
283                                // bracket *= fact(sigma_order);
284
285                                for(int k = 0;k < as1_order-1;k++)
286                                {
287                                        double local_bracket = 0;
288                                        double bracket_factor = -pow((double)-1,k+1)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order-k);
289                                       
290                                        ivec control_vec = ivec();
291                                        control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1);                                 
292                                       
293                                        for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].upper_bound(my_ivec(control_vec));combi_ref++)
294                                        {
295                                                double bracket_combination = 1;
296                                                for(int j = 0;j<(*combi_ref).size();j++)
297                                                {
298                                                        //cout << "Factor vector:" << (*combi_ref) << endl;
299                                                       
300                                                        bracket_combination /= factors[(*combi_ref)[j]-1];
301                                                }                                               
302                                                                                               
303                                                local_bracket += bracket_factor*bracket_combination;                                                                   
304                                        }
305
306                                        simplex->insert_gamma(k+1,local_bracket*multiplier,a_0-(*as_ref).first);
307
308                                        int division_factor = 1;
309                                        for(int s = 0;s<k;s++)
310                                        {
311                                                division_factor *= k-s;
312                                        }
313                                       
314                                        bracket += local_bracket/division_factor; //local_bracket*fact(sigma_order-k);
315                                }                               
316
317                                int_value += multiplier*bracket;
318                                                                                                                                               
319                        }
320
321                        double correction_term_base = correction_term/pow(a_0,sigma_order+1);
322
323                        simplex->insert_gamma(0,correction_term_base,a_0);                     
324
325                        correction_term = correction_term_base;//fact(sigma_order)*correction_term_base;
326
327                        //cout << c << int_value << endl;
328
329                        int_value += correction_term;
330               
331
332                        //cout << "Probability:" << int_value << endl;
333                        //pause(0.100);
334
335                        simplex->probability = int_value;
336                       
337                        return int_value;       
338                       
339                }
340                else
341                {
342                        return 0.0;
343                }
344        }
Note: See TracBrowser for help on using the browser.