root/applications/robust/robustlib.cpp @ 1338

Revision 1338, 8.6 kB (checked in by sindj, 13 years ago)

Zacatky experimentovani v robustlib. Zacleneni vypoctu polynomialni rovnice pro optimalni akci s vyuzitim funkce itpp::roots(). Odstranovani zbylych chyb. 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())
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                                                        order_correct = false;                                         
177
178                                                        break;
179                                                }                                               
180                                               
181                                               
182
183                                                if(a_0<current_emlig->min_ll)
184                                                {
185                                                        current_emlig->minimal_vertex = base_vertex;
186                                                        current_emlig->min_ll = a_0;
187                                                }
188
189                                                double a_m = (*vert_ref)->get_coordinates()*cur_condition-as_toprow->condition_sum[0];
190                                                if(a_m<current_emlig->min_ll)
191                                                {
192                                                        current_emlig->minimal_vertex = (*vert_ref);
193                                                        current_emlig->min_ll = a_m;                                           
194                                                }
195
196                                                //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl;
197
198                                                // cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl;
199                                                //cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl;
200
201                                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
202                                                if(returned.second == false)
203                                                {
204                                                        (*returned.first).second++;
205                                                }
206                                                /*
207                                                else
208                                                {
209                                                        cout << "a[" << row_count << "] = " << new_a << endl;
210                                                }
211                                                */
212                                               
213                                                //as.ins(as.size(),new_a);                                                     
214                                               
215                                                row_count++;
216                                        }
217                                }
218                        }
219                        while(!order_correct);
220                       
221                        /*
222                        cout << "a_0: " << a_0 << "    ";
223                        int as_count = 1;
224                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
225                        {
226                                cout << "a_" << as_count << ": " << (*as_ref).first << "    ";
227                                as_count++;
228                        }*/
229
230                        double int_value = 0;                   
231
232                        // cout << jacobian << endl;
233
234                        double det_jacobian    = abs(det(jacobian));
235                        double correction_term = det_jacobian;                 
236                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
237                        {
238                                double multiplier = det_jacobian;                               
239                               
240                                int as1_order = (*as_ref).second;
241                               
242                                correction_term /= pow(-(*as_ref).first,as1_order);                                     
243                               
244                                current_emlig->set_correction_factors(as1_order);
245
246                                vector<double> factors;
247                                int number_of_factors = 0;                                                             
248                                for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
249                                {
250                                       
251                                        if(as2_ref!=as_ref)
252                                        {                                                                               
253                                                for(int k = 0;k<(*as2_ref).second;k++)
254                                                {
255                                                        factors.push_back((*as_ref).first-(*as2_ref).first);
256                                                }
257
258                                                multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second);
259                                                number_of_factors += (*as2_ref).second;
260                                        }
261                                        else
262                                        {
263                                                factors.push_back((*as_ref).first);
264
265                                                multiplier        /= (*as_ref).first;
266                                                number_of_factors += 1;
267                                        }
268                                       
269                                }
270                               
271                                                               
272
273                                double bracket = fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+1);
274                                                               
275                                simplex->insert_gamma(0,bracket*multiplier,a_0-(*as_ref).first);                                                               
276
277                                bracket *= fact(sigma_order);
278
279                                for(int k = 0;k < as1_order-1;k++)
280                                {
281                                        double local_bracket = 0;
282                                        double bracket_factor = -pow((double)-1,k+1)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order-k);
283                                       
284                                        ivec control_vec = ivec();
285                                        control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1);                                 
286                                       
287                                        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++)
288                                        {
289                                                double bracket_combination = 1;
290                                                for(int j = 0;j<(*combi_ref).size();j++)
291                                                {
292                                                        //cout << "Factor vector:" << (*combi_ref) << endl;
293                                                       
294                                                        bracket_combination /= factors[(*combi_ref)[j]-1];
295                                                }                                               
296                                                                                               
297                                                local_bracket += bracket_factor*bracket_combination;                                                                   
298                                        }
299
300                                        simplex->insert_gamma(k+1,local_bracket*multiplier,a_0-(*as_ref).first);
301                                       
302                                        bracket += local_bracket*fact(sigma_order-k);
303                                }                               
304
305                                int_value += multiplier*bracket;
306                                                                                                                                               
307                        }
308
309                        double correction_term_base = correction_term/pow(a_0,sigma_order+1);
310
311                        simplex->insert_gamma(0,correction_term_base,a_0);                     
312
313                        correction_term = fact(sigma_order)*correction_term_base;
314
315                        //cout << c << int_value << endl;
316
317                        int_value += correction_term;
318               
319
320                        //cout << "Probability:" << int_value << endl;
321                        //pause(0.100);
322
323                        simplex->probability = int_value;
324                       
325                        return int_value;       
326                       
327                }
328                else
329                {
330                        return 0.0;
331                }
332        }
Note: See TracBrowser for help on using the browser.