Changeset 1271 for applications/robust

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

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

Location:
applications/robust
Files:
3 modified

Legend:

Unmodified
Added
Removed

  • applications/robust/main.cpp

    r1269 r1271  
    1313using namespace bdm; 
    1414 
    15 const int emlig_size = 2; 
     15const int emlig_size = 4; 
    1616 
    1717int main ( int argc, char* argv[] ) { 
    … …  
    4343        for(int i = 0;i<500;i++) 
    4444        { 
    45                 cout << "Step:" << i << endl; 
     45                cout << endl << "Step:" << i << endl; 
     46                 
     47                 
    4648 
    4749                double condition[emlig_size+1];          
    … …  
    5557                vec* condition_vec = new vec(condition,emlig_size+1); 
    5658                emlig1->add_condition(*condition_vec); 
     59 
     60                cout << emlig1->statistic_rowsize(emlig_size) << endl << endl; 
    5761         
     62                if(i-emlig1->number_of_parameters >= 0) 
     63                { 
     64                        pause(30); 
     65                } 
    5866 
    59                 emlig1->step_me(i); 
     67                // emlig1->step_me(i); 
    6068                 
    6169                vector<int> sizevector; 

  • 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 

  • applications/robust/robustlib.h

    r1270 r1271  
    250250                this->condition_order = condition_order; 
    251251        } 
     252 
     253        double integrate_simplex(set<vertex*> simplex, char c); 
    252254 
    253255}; 
    … …  
    726728                                                        else 
    727729                                                        { 
    728                                                                 ((toprow*)current_parent)->condition_order++; 
     730                                                                 
    729731 
    730732                                                                if(current_parent->negativechildren.size()>0) 
    … …  
    733735 
    734736                                                                        ((toprow*)current_parent)->condition-=toadd; 
     737 
     738                                                                         
    735739                                                                } 
    736740                                                                else if(current_parent->positivechildren.size()>0) 
    … …  
    743747                                                                { 
    744748                                                                        current_parent->raise_multiplicity();                                                            
     749                                                                } 
     750 
     751                                                                ((toprow*)current_parent)->condition_order++; 
     752 
     753                                                                if(level == number_of_parameters - 1) 
     754                                                                { 
     755                                                                        toprow* cur_par_toprow = ((toprow*)current_parent); 
     756                                                                        cur_par_toprow->probability = 0.0; 
     757                                                                         
     758                                                                        for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++) 
     759                                                                        { 
     760                                                                                cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C'); 
     761                                                                        }                                                                        
    745762                                                                } 
    746763 
    … …  
    955972                                        { 
    956973                                                vec coordinates1 = ((vertex*)(*(current_polyhedron->children.begin())))->get_coordinates();                                              
    957                                                 vec coordinates2 = ((vertex*)(*(current_polyhedron->children.begin()++)))->get_coordinates(); 
    958                                                 coordinates2.ins(0,-1.0); 
     974                                                vec coordinates2 = ((vertex*)(*(++current_polyhedron->children.begin())))->get_coordinates(); 
    959975                                                 
    960                                                 double t = (-toadd*coordinates2)/(toadd(1,toadd.size()-1)*coordinates1)+1; 
    961  
    962                                                 vec new_coordinates = coordinates1*t+(coordinates2(1,coordinates2.size()-1)-coordinates1);                                       
     976                                                vec extended_coord2 = coordinates2; 
     977                                                extended_coord2.ins(0,-1.0); 
     978 
     979                                                double t = (-toadd*extended_coord2)/((toadd(1,toadd.size()-1)*(coordinates1-coordinates2))); 
     980 
     981                                                vec new_coordinates = coordinates2+t*(coordinates1-coordinates2);                                        
     982 
     983                                                // cout << "c1:" << coordinates1 << endl << "c2:" << coordinates2 << endl << "nc:" << new_coordinates << endl; 
    963984 
    964985                                                vertex* neutral_vertex = new vertex(new_coordinates);