Show
Ignore:
Timestamp:
01/10/11 19:11:25 (13 years ago)
Author:
sindj
Message:

Juchuu, pocita to. Odstraneny nektere chyby. Zbyva odstranit chybu s negativnimi pravdepodobnostmi. JS

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • applications/robust/robustlib.cpp

    r1270 r1271  
    33void polyhedron::triangulate(bool should_integrate) 
    44        {                
     5                if(should_integrate) 
     6                { 
     7                        ((toprow *)this)->probability = 0.0; 
     8                }                
     9                 
    510                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++) 
    6                 { 
    7                         set<double> simplex_integrals; 
    8                          
     11                {                        
    912                        for(list<set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++) 
    1013                        { 
    … …  
    1821                                        triangulation.push_back(new_simplex); 
    1922 
    20                                         int condition_order = ((toprow*)this)->condition_order-1; 
     23                                        if(should_integrate) 
     24                                        { 
     25                                                ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(new_simplex, 'S'); 
     26                                        } 
     27                                }  
     28                        }        
     29                }                
     30        } 
    2131 
    22                                         // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl; 
    23                                         // pause(0.1); 
    24                                         if(should_integrate && condition_order-my_emlig->number_of_parameters >= 0) 
     32 
     33        double toprow::integrate_simplex(set<vertex*> simplex, char c) 
     34        { 
     35                int condition_order = ((toprow*)this)->condition_order-1; 
     36 
     37                // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl; 
     38                // pause(0.1); 
     39 
     40                if(condition_order-my_emlig->number_of_parameters >= 0) 
     41                { 
     42                        emlig* current_emlig; 
     43 
     44                        if(this->my_emlig!=NULL) 
     45                        { 
     46                                current_emlig = this->my_emlig; 
     47                        } 
     48                        else 
     49                        { 
     50                                throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!"); 
     51                        }                                                
     52 
     53                        toprow* as_toprow = (toprow*)this; 
     54 
     55                        vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1); 
     56 
     57                        vertex* base_vertex = (*simplex.begin()); 
     58 
     59                        int dimension = simplex.size()-1; 
     60 
     61                        mat jacobian(dimension,dimension); 
     62 
     63                        double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0]; 
     64                        map<double,int> as; 
     65 
     66                        int row_count = 0; 
     67                        for(set<vertex*>::iterator vert_ref = (++simplex.begin()); vert_ref!=simplex.end();vert_ref++) 
     68                        { 
     69                                vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 
     70 
     71                                // cout << relative_coords << endl; 
     72 
     73                                jacobian.set_row(row_count,relative_coords); 
     74 
     75                                double new_a = relative_coords*cur_condition;                                                    
     76                                 
     77                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); 
     78                                if(returned.second == false) 
     79                                { 
     80                                        (*returned.first).second++; 
     81                                } 
     82                                else 
     83                                { 
     84                                        cout << "a[" << row_count << "] = " << new_a << endl; 
     85                                } 
     86                                //as.ins(as.size(),new_a);                                                       
     87                                 
     88                                row_count++; 
     89                        } 
     90 
     91                        // cout << endl; 
     92 
     93                        double int_value = 0; 
     94 
     95                        // cout << jacobian << endl; 
     96 
     97                        double det_jacobian = det(jacobian); 
     98                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) 
     99                        { 
     100                                double multiplier = det_jacobian; 
     101                                if(a_0!=(*as_ref).first) 
     102                                {                                                                
     103                                        int as1_order = (*as_ref).second; 
     104 
     105                                        current_emlig->set_correction_factors(as1_order); 
     106 
     107                                        vector<double> factors; 
     108                                        int number_of_factors = 0;                                                               
     109                                        for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) 
    25110                                        { 
    26                                                 emlig* current_emlig; 
     111                                                 
     112                                                if(as2_ref!=as_ref) 
     113                                                {                                                                                
     114                                                        for(int k = 0;k<(*as2_ref).second;k++) 
     115                                                        { 
     116                                                                factors.push_back((*as_ref).first-(*as2_ref).first); 
     117                                                        } 
    27118 
    28                                                 if(this->my_emlig!=NULL) 
    29                                                 { 
    30                                                         current_emlig = this->my_emlig; 
     119                                                        multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second); 
     120                                                        number_of_factors += (*as2_ref).second; 
    31121                                                } 
    32122                                                else 
    33123                                                { 
    34                                                         throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!"); 
    35                                                 }                                                
     124                                                        factors.push_back((*as_ref).first); 
    36125 
    37                                                 toprow* as_toprow = (toprow*)this; 
     126                                                        multiplier        /= (*as_ref).first; 
     127                                                        number_of_factors += 1; 
     128                                                }                                                                        
     129                                        }                                        
    38130 
    39                                                 vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1); 
     131                                        double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors); 
     132                                        for(int k = 0;k < as1_order-1;k++) 
     133                                        { 
     134                                                double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k); 
     135                                                 
     136                                                for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++) 
     137                                                { 
     138                                                        double bracket_combination = 1; 
     139                                                        for(int j = 0;j<(*combi_ref).size();j++) 
     140                                                        { 
     141                                                                bracket_combination /= factors[(*combi_ref)[j]]; 
     142                                                        } 
    40143 
    41                                                 vertex* base_vertex = (*new_simplex.begin()); 
     144                                                        bracket+=bracket_factor*bracket_combination;                                                                     
     145                                                }                                                                        
     146                                        }                                                                
    42147 
    43                                                 int dimension = new_simplex.size()-1; 
     148                                        int_value += multiplier*bracket;         
    44149 
    45                                                 mat jacobian(dimension,dimension); 
     150                                } 
     151                                else 
     152                                { 
     153                                        throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)"); 
     154                                } 
    46155 
    47                                                 double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0]; 
    48                                                 map<double,int> as; 
     156                                // cout << c << int_value << endl; 
    49157 
    50                                                 int row_count = 0; 
    51  
    52                                                 for(set<vertex*>::iterator vert_ref = (++new_simplex.begin()); vert_ref!=new_simplex.end();vert_ref++) 
    53                                                 { 
    54                                                         vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 
    55                                                         jacobian.set_row(row_count,relative_coords); 
    56  
    57                                                         double new_a = relative_coords*cur_condition;                                                    
    58                                                          
    59                                                         pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); 
    60                                                         if(returned.second == false) 
    61                                                         { 
    62                                                                 (*returned.first).second++; 
    63                                                         } 
    64                                                         //as.ins(as.size(),new_a);                                                       
    65                                                          
    66                                                         row_count++; 
    67                                                 } 
    68  
    69                                                 double int_value = 0; 
    70                                                 double det_jacobian = det(jacobian); 
    71                                                 for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) 
    72                                                 { 
    73                                                         double multiplier = det_jacobian; 
    74                                                         if(a_0!=(*as_ref).first) 
    75                                                         {                                                                
    76                                                                 int as1_order = (*as_ref).second; 
    77  
    78                                                                 current_emlig->set_correction_factors(as1_order); 
    79  
    80                                                                 vector<double> factors; 
    81                                                                 int number_of_factors = 0;                                                               
    82                                                                 for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) 
    83                                                                 { 
    84                                                                          
    85                                                                         if(as2_ref!=as_ref) 
    86                                                                         {                                                                                
    87                                                                                 for(int k = 0;k<(*as2_ref).second;k++) 
    88                                                                                 { 
    89                                                                                         factors.push_back((*as_ref).first-(*as2_ref).first); 
    90                                                                                 } 
    91  
    92                                                                                 multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second); 
    93                                                                                 number_of_factors += (*as2_ref).second; 
    94                                                                         } 
    95                                                                         else 
    96                                                                         { 
    97                                                                                 factors.push_back((*as_ref).first); 
    98  
    99                                                                                 multiplier        /= (*as_ref).first; 
    100                                                                                 number_of_factors += 1; 
    101                                                                         }                                                                        
    102                                                                 } 
    103  
    104                                                                 double facti = fact(0); 
    105  
    106                                                                 double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors); 
    107                                                                 for(int k = 0;k < as1_order-1;k++) 
    108                                                                 { 
    109                                                                         double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k); 
    110                                                                          
    111                                                                         for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++) 
    112                                                                         { 
    113                                                                                 double bracket_combination = 1; 
    114                                                                                 for(int j = 0;j<(*combi_ref).size();j++) 
    115                                                                                 { 
    116                                                                                         bracket_combination /= factors[(*combi_ref)[j]]; 
    117                                                                                 } 
    118  
    119                                                                                 bracket+=bracket_factor*bracket_combination;                                                                     
    120                                                                         }                                                                        
    121                                                                 }                                                                
    122  
    123                                                                 int_value += multiplier*bracket;         
    124  
    125                                                         } 
    126                                                         else 
    127                                                         { 
    128                                                                 throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)"); 
    129                                                         }                                                
    130                                                                                                          
    131                                                 } 
    132  
    133                                                 cout << int_value << endl; 
    134                                                 pause(0.1); 
    135  
    136  
    137                                                 simplex_integrals.insert(int_value); 
    138                                         } 
    139                                 }  
     158                                return int_value;                                                                                
    140159                        } 
    141  
    142                         if(should_integrate) 
    143                         { 
    144                                 ((toprow*)this)->probability = 0.0; 
    145  
    146                                 for(set<double>::iterator integ_ref = simplex_integrals.begin();integ_ref!=simplex_integrals.end();integ_ref++) 
    147                                 { 
    148                                         ((toprow*)this)->probability += (*integ_ref); 
    149                                 } 
    150                         } 
    151                 }                
     160                         
     161                } 
     162                else 
     163                { 
     164                        return 0.0; 
     165                } 
    152166        } 
    153167