| 197 | } |
| 198 | |
| 199 | double bracket_combination = 0; |
| 200 | |
| 201 | for(int k = 0;k < a_order;k++) |
| 202 | { |
| 203 | if(k==gamma_facs.size()) |
| 204 | { |
| 205 | gamma_facs.push_back(0.0); |
| 206 | } |
| 207 | |
| 208 | double factor_multiplier = pow((double)-1,k)/pow((*a_ref).first,sigma_order-k);//pow((double)-1,k+1) |
| 209 | |
| 210 | if(k!=0) |
| 211 | { |
| 212 | ivec control_vec = ivec(); |
| 213 | control_vec.ins(0,my_emlig->number_of_parameters-a_order+1); |
| 214 | for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k-1].begin();combi_ref!=this->my_emlig->correction_factors[k-1].upper_bound(my_ivec(control_vec));combi_ref++) |
| 215 | { |
| 216 | double bracket_factor = 1/fact(a_order-1-k); |
| 217 | |
| 218 | for(int j = 0;j<(*combi_ref).size();j++) |
| 219 | { |
| 220 | //cout << "Factor vector:" << (*combi_ref) << endl; |
| 221 | |
| 222 | bracket_factor /= factors[(*combi_ref)[j]-1]; |
| 223 | } |
| 224 | |
| 225 | double value = bracket_factor*factor_multiplier*k_multiplier; |
| 226 | |
| 227 | simplex->insert_gamma(k,value,(*a_ref).first); |
| 228 | gamma_facs[k] += value; |
| 229 | } |
| 230 | } |
| 231 | else |
| 232 | { |
| 233 | double value = k_multiplier*factor_multiplier/fact(a_order-1); |
| 234 | |
| 235 | simplex->insert_gamma(0,value,(*a_ref).first); |
| 236 | gamma_facs[k] += value; |
| 237 | } |
| 238 | } |
| 239 | } |
| 240 | |
| 241 | double int_value = 0; |
| 242 | int gamma_size = (int)gamma_facs.size(); |
| 243 | if(sigma_order-gamma_size>=0) |
| 244 | { |
| 245 | for(int k = 0;k<gamma_size;k++) |
| 246 | { |
| 247 | int_value += det_jacobian*fact(sigma_order-gamma_size+k)/fact(jacobian.rows())*gamma_facs[k]/pow(2.0,sigma_order); |
199 | | |
200 | | double cur_as_ref = (*a_ref).first; |
201 | | |
202 | | double gamma_multiplier = -cur_as_ref-a_0; |
203 | | |
204 | | double bracket = fact(as1_order-1)/pow(gamma_multiplier,sigma_order); |
205 | | |
206 | | simplex->insert_gamma(0,bracket*multiplier,gamma_multiplier); |
207 | | |
208 | | // bracket *= fact(sigma_order); |
209 | | |
210 | | for(int k = 0;k < as1_order-1;k++) |
211 | | { |
212 | | double local_bracket = 0; |
213 | | double bracket_factor = 1/fact(as1_order-1-k)/pow(gamma_multiplier,sigma_order-k-1);//pow((double)-1,k+1) |
214 | | |
215 | | ivec control_vec = ivec(); |
216 | | control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1); |
217 | | |
218 | | 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++) |
219 | | { |
220 | | double bracket_combination = 1; |
221 | | for(int j = 0;j<(*combi_ref).size();j++) |
222 | | { |
223 | | //cout << "Factor vector:" << (*combi_ref) << endl; |
224 | | |
225 | | bracket_combination /= factors[(*combi_ref)[j]-1]; |
226 | | } |
227 | | |
228 | | local_bracket += bracket_factor*bracket_combination; |
229 | | } |
230 | | |
231 | | simplex->insert_gamma(k+1,local_bracket*multiplier,gamma_multiplier); |
232 | | |
233 | | int division_factor = 1; |
234 | | for(int s = 0;s<k;s++) |
235 | | { |
236 | | division_factor *= k-s; |
237 | | } |
238 | | |
239 | | bracket += local_bracket/division_factor; //local_bracket*fact(sigma_order-k); |
240 | | } |
241 | | |
242 | | int_value += multiplier*bracket; |
243 | | |
244 | | |
245 | | |
246 | | } |
247 | | |
248 | | double correction_term_base = correction_term/pow(-a_0,sigma_order); |
249 | | |
250 | | simplex->insert_gamma(0,correction_term_base,-a_0); |
251 | | |
252 | | correction_term = correction_term_base;//fact(sigma_order)*correction_term_base; |
253 | | |
254 | | //cout << c << int_value << endl; |
255 | | |
256 | | int_value += correction_term; |
257 | | |
| 249 | } |
| 250 | else |
| 251 | { |
| 252 | int_value = -1; |
| 253 | } |
| 254 | |
| 255 | simplex->probability = int_value; |