69 | | |
70 | | int row_count = 0; |
71 | | for(set<vertex*>::iterator vert_ref = (++simplex.begin()); vert_ref!=simplex.end();vert_ref++) |
72 | | { |
73 | | vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); |
74 | | |
75 | | // cout << (*vert_ref)->get_coordinates() << endl; |
76 | | // cout << relative_coords << endl; |
77 | | |
78 | | jacobian.set_row(row_count,relative_coords); |
79 | | |
80 | | double new_a = relative_coords*cur_condition; |
81 | | |
82 | | pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); |
83 | | if(returned.second == false) |
84 | | { |
85 | | (*returned.first).second++; |
| 72 | vertex* base_vertex; |
| 73 | set<vertex*>::iterator base_ref = simplex.begin(); |
| 74 | bool order_correct; |
| 75 | |
| 76 | |
| 77 | do |
| 78 | { |
| 79 | order_correct = true; |
| 80 | as.clear(); |
| 81 | |
| 82 | base_vertex = (*base_ref); |
| 83 | |
| 84 | cout << "Base coords:" << base_vertex->get_coordinates() << endl; |
| 85 | |
| 86 | a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition[0]; |
| 87 | |
| 88 | int row_count = 0; |
| 89 | for(set<vertex*>::iterator vert_ref = simplex.begin(); vert_ref!=simplex.end();vert_ref++) |
| 90 | { |
| 91 | if(vert_ref != base_ref) |
| 92 | { |
| 93 | vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); |
| 94 | |
| 95 | |
| 96 | |
| 97 | jacobian.set_row(row_count,relative_coords); |
| 98 | |
| 99 | double new_a = relative_coords*cur_condition; |
| 100 | |
| 101 | if(new_a + a_0 == 0) |
| 102 | { |
| 103 | base_ref++; |
| 104 | |
| 105 | if(base_ref == simplex.end()) |
| 106 | { |
| 107 | throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!"); |
| 108 | } |
| 109 | |
| 110 | order_correct = false; |
| 111 | |
| 112 | break; |
| 113 | } |
| 114 | |
| 115 | cout << "Absolute coords:(V" << row_count << ")" << (*vert_ref)->get_coordinates() << endl; |
| 116 | //cout << "Relative coords:(V" << row_count << ")" << relative_coords << endl; |
| 117 | |
| 118 | pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); |
| 119 | if(returned.second == false) |
| 120 | { |
| 121 | (*returned.first).second++; |
| 122 | } |
| 123 | /* |
| 124 | else |
| 125 | { |
| 126 | cout << "a[" << row_count << "] = " << new_a << endl; |
| 127 | } |
| 128 | */ |
| 129 | |
| 130 | //as.ins(as.size(),new_a); |
| 131 | |
| 132 | row_count++; |
| 133 | } |
109 | | double multiplier = det_jacobian; |
110 | | |
111 | | if(a_0!=(*as_ref).first) |
112 | | { |
113 | | int as1_order = (*as_ref).second; |
| 147 | double multiplier = det_jacobian; |
| 148 | |
| 149 | int as1_order = (*as_ref).second; |
| 150 | |
| 151 | correction_term /= -pow((*as_ref).first,as1_order); |
| 152 | |
| 153 | current_emlig->set_correction_factors(as1_order); |
| 154 | |
| 155 | vector<double> factors; |
| 156 | int number_of_factors = 0; |
| 157 | for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) |
| 158 | { |
117 | | current_emlig->set_correction_factors(as1_order); |
118 | | |
119 | | vector<double> factors; |
120 | | int number_of_factors = 0; |
121 | | for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) |
122 | | { |
123 | | |
124 | | if(as2_ref!=as_ref) |
125 | | { |
126 | | for(int k = 0;k<(*as2_ref).second;k++) |
127 | | { |
128 | | factors.push_back((*as_ref).first-(*as2_ref).first); |
129 | | } |
130 | | |
131 | | multiplier /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second); |
132 | | number_of_factors += (*as2_ref).second; |
133 | | } |
134 | | else |
135 | | { |
136 | | factors.push_back((*as_ref).first); |
137 | | |
138 | | multiplier /= (*as_ref).first; |
139 | | number_of_factors += 1; |
140 | | } |
141 | | |
142 | | } |
143 | | |
144 | | double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow(a_0+(*as_ref).first,condition_order-number_of_factors+1); |
145 | | for(int k = 0;k < as1_order-1;k++) |
146 | | { |
147 | | double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow(a_0+(*as_ref).first,condition_order-number_of_factors-k); |
148 | | |
149 | | 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++) |
150 | | { |
151 | | double bracket_combination = 1; |
152 | | for(int j = 0;j<(*combi_ref).size();j++) |
153 | | { |
154 | | bracket_combination /= factors[(*combi_ref)[j]]; |
155 | | } |
156 | | |
157 | | bracket+=bracket_factor*bracket_combination; |
158 | | } |
159 | | } |
160 | | |
| 178 | } |
| 179 | |
| 180 | double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow(a_0+(*as_ref).first,condition_order-number_of_factors); |
| 181 | for(int k = 0;k < as1_order-1;k++) |
| 182 | { |
| 183 | double bracket_factor = pow((double)-1,k+1)*fact(condition_order-number_of_factors-k)/fact(as1_order-2-k)/pow(a_0+(*as_ref).first,condition_order-number_of_factors-k); |