Changeset 1376 for applications/robust/robustlib.cpp
- Timestamp:
- 06/21/11 19:09:56 (13 years ago)
- Files:
-
- 1 modified
-
applications/robust/robustlib.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/robustlib.cpp
r1367 r1376 121 121 } 122 122 123 toprow* as_toprow = (toprow*)this; 124 125 vec cur_condition = as_toprow->condition_sum.get(1,as_toprow->condition_sum.size()-1); 126 127 // cout << as_toprow->condition << endl; 123 toprow* as_toprow = (toprow*)this; 128 124 129 125 int dimension = simplex->vertices.size()-1; … … 131 127 mat jacobian(dimension,dimension); 132 128 133 double a_0;134 129 map<double,int> as; 135 vertex* base_vertex; 136 set<vertex*>::iterator base_ref = simplex->vertices.begin(); 137 bool order_correct; 138 139 140 do 141 { 142 order_correct = true; 143 as.clear(); 144 145 base_vertex = (*base_ref); 146 147 /* 148 cout << "Base coords: " << base_vertex->get_coordinates() << endl; 149 cout << "Condition: " << as_toprow->condition_sum << endl; 150 pause(0.1); 151 */ 152 153 a_0 = -base_vertex->get_coordinates()*cur_condition+as_toprow->condition_sum[0]; 130 vertex* base_vertex = simplex->vertices.begin(); 131 154 132 155 int row_count = 0; 156 for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 157 { 158 if(vert_ref != base_ref) 159 { 160 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 161 162 jacobian.set_row(row_count,relative_coords); 163 164 double new_a = -relative_coords*cur_condition; 165 166 if(new_a + a_0 == 0 || new_a == 0) 167 { 168 base_ref++; 169 170 if(base_ref == simplex->vertices.end()&&simplex->vertices.size()!=2) 171 { 172 throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!"); 173 } 174 /* 175 else if(implex->vertices.size()==2) 176 { 177 return relative_coords[0]*exp 178 } 179 */ 180 181 order_correct = false; 182 183 break; 184 } 185 186 187 188 if(a_0<current_emlig->min_ll) 189 { 190 current_emlig->minimal_vertex = base_vertex; 191 current_emlig->min_ll = a_0; 192 } 193 194 double a_m = -(*vert_ref)->get_coordinates()*cur_condition+as_toprow->condition_sum[0]; 195 if(a_m<current_emlig->min_ll) 196 { 197 current_emlig->minimal_vertex = (*vert_ref); 198 current_emlig->min_ll = a_m; 199 } 200 201 //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl; 202 203 // cout << "Absolute coords:(V" << row_count << ")" << (*vert_ref)->get_coordinates() << endl; 204 //cout << "Relative coords:(V" << row_count << ")" << relative_coords << endl; 205 206 pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); 207 if(returned.second == false) 208 { 209 (*returned.first).second++; 210 } 211 /* 212 else 213 { 214 cout << "a[" << row_count << "] = " << new_a << endl; 215 } 216 */ 217 218 //as.ins(as.size(),new_a); 219 220 row_count++; 221 } 222 } 223 } 224 while(!order_correct); 133 for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 134 { 135 if(vert_ref!=simplex->vertices.begin()) 136 { 137 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 138 jacobian.set_row(row_count,relative_coords); 139 } 140 141 double a = -((*vert_ref)->get_coordinates().ins(0,-1)*cur_condition); 142 if(a<current_emlig->min_ll) 143 { 144 current_emlig->minimal_vertex = (*vert_ref); 145 current_emlig->min_ll = a; 146 } 147 148 //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl; 149 150 // cout << "Absolute coords:(V" << row_count << ")" << (*vert_ref)->get_coordinates() << endl; 151 //cout << "Relative coords:(V" << row_count << ")" << relative_coords << endl; 152 153 pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(a,1)); 154 if(returned.second == false) 155 { 156 (*returned.first).second++; 157 } 158 } 159 225 160 226 161 /* … … 239 174 240 175 double det_jacobian = abs(det(jacobian)); 241 double correction_term = det_jacobian;242 for(map<double,int>::iterator a s_ref = as.begin();as_ref!=as.end();as_ref++)243 { 244 double multiplier = det_jacobian ;176 double correction_term; 177 for(map<double,int>::iterator a_ref = as.begin();as_ref!=as.end();as_ref++) 178 { 179 double multiplier = det_jacobian/fact(jacobian.rows()); 245 180 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); 181 int a_order = (*a_ref).second; 182 current_emlig->set_correction_factors(a_order); 251 183 252 184 vector<double> factors; 253 185 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) 186 for(map<double,int>::iterator a2_ref = as.begin();a2_ref!=as.end();a2_ref++) 187 { 188 if(a2_ref!=a_ref) 258 189 { 259 for(int k = 0;k<(*a s2_ref).second;k++)190 for(int k = 0;k<(*a2_ref).second;k++) 260 191 { 261 factors.push_back((*a s_ref).first-(*as2_ref).first);192 factors.push_back((*a_ref).first-(*a2_ref).first); 262 193 } 263 194 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 } 195 multiplier /= pow((*a_ref).first-(*a2_ref).first,(*a2_ref).second); 196 number_of_factors += (*a2_ref).second; 197 } 274 198 } 275 199 276 double cur_as_ref = (*a s_ref).first;200 double cur_as_ref = (*a_ref).first; 277 201 278 202 double gamma_multiplier = -cur_as_ref-a_0;
