1 | #include "madx.h" |
---|
2 | |
---|
3 | /*______________________________________________________________ |
---|
4 | matchptcknobs.c |
---|
5 | Piotr Skowronski (CERN) 2006 |
---|
6 | MAD-X module responsible for matching with help of PTC knobs (parameters) |
---|
7 | |
---|
8 | Feb.2007: added proper support for variable limits |
---|
9 | |
---|
10 | It implements MAD-X the command: |
---|
11 | match, useptcknobs=true; |
---|
12 | (...) |
---|
13 | endmatch; |
---|
14 | |
---|
15 | It is basically a macro that performs automatically a set of commands |
---|
16 | depicted above as (...) within a rather complicated loop (points 3-17 below). |
---|
17 | It enables fast matching of any order parameters with PTC. |
---|
18 | |
---|
19 | The algorithm: |
---|
20 | |
---|
21 | 1. Buffer the key commands (ptc_varyknob, constraint, ptc_setswitch, ptc_twiss or ptc_normal, etc) |
---|
22 | appearing between |
---|
23 | match, useptcknobs=true; |
---|
24 | and any of matching actions calls (migrad,lmdif,jacobian, etc) |
---|
25 | |
---|
26 | 2. When matching action appears, |
---|
27 | a) set "The Current Variables Values" (TCVV) to zero |
---|
28 | b) perform THE LOOP, i.e. points 3-17 |
---|
29 | |
---|
30 | 3. Prepare PTC environment (ptc_createuniverse, ptc_createlayout) |
---|
31 | |
---|
32 | 4. Set the user defined knobs (with ptc_knob). |
---|
33 | |
---|
34 | 5. Set TCVV using ptc_setfieldcomp command |
---|
35 | |
---|
36 | 6. Run a PTC command (twiss or normal) |
---|
37 | |
---|
38 | 7. Run a runtime created script that performs a standard matching |
---|
39 | All the user defined knobs are variables of this matching. |
---|
40 | |
---|
41 | 8. Evaluate constraints expressions to get the matching function vector (I) |
---|
42 | |
---|
43 | 9. Add the matched values to TCVV |
---|
44 | |
---|
45 | 10. End PTC session (run ptc_end) |
---|
46 | |
---|
47 | 11. If the matched values are are not close enough to zeroes then goto 3 |
---|
48 | |
---|
49 | 12. Prepare PTC environment (ptc_createuniverse, ptc_createlayout) |
---|
50 | |
---|
51 | 13. Set TCVV using ptc_setfieldcomp command |
---|
52 | ( --- please note that knobs are not set in this case -- ) |
---|
53 | 14. Run a PTC command (twiss or normal) |
---|
54 | |
---|
55 | 15. Evaluate constraints expressions to get the matching function vector (II) |
---|
56 | |
---|
57 | 16. Evaluate a penalty function that compares matching function vectors (I) and (II) |
---|
58 | See points 7 and 14 |
---|
59 | |
---|
60 | 17 If the matching function vectors are not similar to each other |
---|
61 | within requested precision then goto 3 |
---|
62 | |
---|
63 | 18. Print TCVV, which are the matched values. |
---|
64 | |
---|
65 | ______________________________________________________________*/ |
---|
66 | |
---|
67 | #define MAX_CONTRAINS 100 |
---|
68 | #define MAX_KNOBS 100 |
---|
69 | #define COMM_LENGTH 500 |
---|
70 | |
---|
71 | struct madx_mpk_variable |
---|
72 | { |
---|
73 | char* name; /*variable name, if initial parameter just copy with mpk_ prefix, if field comp then {elname}_{"kn"/"ks"}{number}*/ |
---|
74 | char* namecv; |
---|
75 | double upper; |
---|
76 | double lower; |
---|
77 | double trustrange; |
---|
78 | double step; |
---|
79 | int knobidx; /*index of the knob that makes paramter of this variable */ |
---|
80 | double currentvalue; |
---|
81 | double oldvalue; |
---|
82 | int kn; /*if a field component then different from -1*/ |
---|
83 | int ks; |
---|
84 | char IsIniCond; /**/ |
---|
85 | }; |
---|
86 | |
---|
87 | struct madx_mpk_knob |
---|
88 | { |
---|
89 | char* elname; /*if a field property element name */ |
---|
90 | char* initial; /*if initial condition, name of a variable */ |
---|
91 | int* KN; /*array with filed components*/ |
---|
92 | int nKN; |
---|
93 | int* KS; |
---|
94 | int nKS; |
---|
95 | int exactnamematch; |
---|
96 | }; |
---|
97 | |
---|
98 | // private variables |
---|
99 | |
---|
100 | // static int madx_mpk_debug; |
---|
101 | |
---|
102 | // constraints |
---|
103 | static int madx_mpk_Nconstraints; |
---|
104 | static char* madx_mpk_constraints[MAX_CONTRAINS]; |
---|
105 | |
---|
106 | // knobs |
---|
107 | static int madx_mpk_Nknobs; |
---|
108 | static struct madx_mpk_knob madx_mpk_knobs[MAX_KNOBS]; |
---|
109 | static char* madx_mpk_setknobs[MAX_KNOBS]; /*this is ptc_setknobs*/ |
---|
110 | // static struct in_cmd* madx_mpk_varyknob_cmds[MAX_KNOBS]; /*this is ptc_setknobs*/ |
---|
111 | |
---|
112 | // other |
---|
113 | static int madx_mpk_Nvariables; /*total number if matching variables >= number of knob commands (one knob command may define more then one knob/variable)*/ |
---|
114 | static struct madx_mpk_variable madx_mpk_variables[MAX_KNOBS]; |
---|
115 | |
---|
116 | static struct in_cmd* madx_mpk_comm_createuniverse; |
---|
117 | static struct in_cmd* madx_mpk_comm_createlayout; |
---|
118 | static struct in_cmd* madx_mpk_comm_setswitch; |
---|
119 | static struct in_cmd* madx_mpk_comm_calculate;/*ptc_twiss or ptc_normal*/ |
---|
120 | |
---|
121 | // static char twisscommand[]="ptc_twiss, table=ptc_twiss, icase=6, no=2, betx=10, alfx=.0, bety=10., alfy=0, betz=10., alfz=0;"; |
---|
122 | |
---|
123 | // private functions |
---|
124 | |
---|
125 | static int |
---|
126 | factorial(int v) |
---|
127 | { |
---|
128 | if (v <= 1 ) |
---|
129 | { |
---|
130 | return 1; |
---|
131 | } |
---|
132 | else |
---|
133 | { |
---|
134 | return v*factorial(v-1); |
---|
135 | } |
---|
136 | } |
---|
137 | |
---|
138 | static int |
---|
139 | findsetknob(char* ename, int exactnamematch, char* initialpar) |
---|
140 | { |
---|
141 | /* |
---|
142 | finds a knob that corresponds to this name, and detects eventual errors |
---|
143 | if not found returns fresh empty knob |
---|
144 | */ |
---|
145 | int i; |
---|
146 | int cmpres; |
---|
147 | int result = madx_mpk_Nknobs; /*a new empty one*/ |
---|
148 | char* bconta; |
---|
149 | char* acontb; |
---|
150 | |
---|
151 | if (ename) |
---|
152 | { |
---|
153 | for (i = 0; i < madx_mpk_Nknobs; i++) |
---|
154 | { |
---|
155 | if (madx_mpk_knobs[i].elname == 0x0) continue; |
---|
156 | cmpres = strcmp(ename,madx_mpk_knobs[i].elname); |
---|
157 | if (cmpres == 0) |
---|
158 | { |
---|
159 | if ( exactnamematch == madx_mpk_knobs[i].exactnamematch ) |
---|
160 | { |
---|
161 | result = i; /*the same element or family*/ |
---|
162 | continue; /*to find out if there is any setting problem*/ |
---|
163 | } |
---|
164 | else |
---|
165 | { |
---|
166 | error("findsetknob","A knob for such named element(s) found, but name matching flag does not agree."); |
---|
167 | return -i; |
---|
168 | } |
---|
169 | } |
---|
170 | |
---|
171 | /* |
---|
172 | if a not-exact and b exact, b can not contain a |
---|
173 | if a not-exact and b non-exact, b can not contain a and a can not contain b |
---|
174 | if a exact, and b not-exact, a can not contain b |
---|
175 | */ |
---|
176 | acontb = strstr(ename, madx_mpk_knobs[i].elname); |
---|
177 | bconta = strstr(madx_mpk_knobs[i].elname, ename); |
---|
178 | if ( (acontb && (exactnamematch == 0)) || (bconta && (madx_mpk_knobs[i].exactnamematch == 0)) ) |
---|
179 | { |
---|
180 | error("findsetknob", |
---|
181 | "This variable (name %s, exactmatch %d) can cause ambiguity with another already defined variable (name %s, exactmatch %d)", |
---|
182 | ename, exactnamematch, madx_mpk_knobs[i].elname, madx_mpk_knobs[i].exactnamematch); |
---|
183 | return -i; |
---|
184 | } |
---|
185 | } |
---|
186 | } |
---|
187 | |
---|
188 | if (initialpar) |
---|
189 | { |
---|
190 | for (i = 0; i < madx_mpk_Nknobs; i++) |
---|
191 | { |
---|
192 | if (madx_mpk_knobs[i].initial) |
---|
193 | { |
---|
194 | if (( strcmp(initialpar,madx_mpk_knobs[i].initial) == 0 )) |
---|
195 | { |
---|
196 | |
---|
197 | error("findsetknob","Such initial parameter is already defined"); |
---|
198 | return -i; |
---|
199 | } |
---|
200 | } |
---|
201 | } |
---|
202 | } |
---|
203 | |
---|
204 | if (result == madx_mpk_Nknobs) |
---|
205 | { |
---|
206 | madx_mpk_Nknobs++; /*we have not found a setknob to add this vary knob*/ |
---|
207 | } |
---|
208 | |
---|
209 | return result; |
---|
210 | } |
---|
211 | |
---|
212 | static int |
---|
213 | mapptctomad(char* ptcname, char* madxname) |
---|
214 | { |
---|
215 | |
---|
216 | if( strcmp(ptcname,"beta11") == 0 ) |
---|
217 | { |
---|
218 | strcpy(madxname,"betx"); |
---|
219 | return 0; |
---|
220 | } |
---|
221 | |
---|
222 | if( strcmp(ptcname,"beta22") == 0 ) |
---|
223 | { |
---|
224 | strcpy(madxname,"bety"); |
---|
225 | return 0; |
---|
226 | } |
---|
227 | |
---|
228 | if( strcmp(ptcname,"beta33") == 0 ) |
---|
229 | { |
---|
230 | strcpy(madxname,"betz"); |
---|
231 | return 0; |
---|
232 | } |
---|
233 | |
---|
234 | if( strcmp(ptcname,"alfa11") == 0 ) |
---|
235 | { |
---|
236 | strcpy(madxname,"alfx"); |
---|
237 | return 0; |
---|
238 | } |
---|
239 | |
---|
240 | if( strcmp(ptcname,"alfa22") == 0 ) |
---|
241 | { |
---|
242 | strcpy(madxname,"alfy"); |
---|
243 | return 0; |
---|
244 | } |
---|
245 | |
---|
246 | if( strcmp(ptcname,"alfa33") == 0 ) |
---|
247 | { |
---|
248 | strcpy(madxname,"alfz"); |
---|
249 | return 0; |
---|
250 | } |
---|
251 | |
---|
252 | |
---|
253 | if( strcmp(ptcname,"disp1") == 0 ) |
---|
254 | { |
---|
255 | strcpy(madxname,"dx"); |
---|
256 | return 0; |
---|
257 | } |
---|
258 | |
---|
259 | if( strcmp(ptcname,"disp2") == 0 ) |
---|
260 | { |
---|
261 | strcpy(madxname,"dpz"); |
---|
262 | return 0; |
---|
263 | } |
---|
264 | |
---|
265 | if( strcmp(ptcname,"disp3") == 0 ) |
---|
266 | { |
---|
267 | strcpy(madxname,"dy"); |
---|
268 | return 0; |
---|
269 | } |
---|
270 | |
---|
271 | if( strcmp(ptcname,"disp4") == 0 ) |
---|
272 | { |
---|
273 | strcpy(madxname,"dpy"); |
---|
274 | return 0; |
---|
275 | } |
---|
276 | |
---|
277 | |
---|
278 | return 1; |
---|
279 | } |
---|
280 | |
---|
281 | static int |
---|
282 | readstartvalues(void) |
---|
283 | { |
---|
284 | /* |
---|
285 | reads the starting values of variable to initializ currentvalue |
---|
286 | */ |
---|
287 | |
---|
288 | int i; |
---|
289 | int n; |
---|
290 | int ncomp; |
---|
291 | double nvalue; |
---|
292 | struct node* node = 0x0; |
---|
293 | struct madx_mpk_variable* v = 0x0; |
---|
294 | struct madx_mpk_knob* kn = 0x0; |
---|
295 | char buff[50]; |
---|
296 | char* p; |
---|
297 | |
---|
298 | if (debuglevel) printf("\n\n\n READING INITIAL VALUES \n\n\n"); |
---|
299 | |
---|
300 | for (i = 0; i < madx_mpk_Nvariables; i++) |
---|
301 | { |
---|
302 | v = &(madx_mpk_variables[i]); |
---|
303 | kn = &(madx_mpk_knobs[v->knobidx]); |
---|
304 | |
---|
305 | if (v->IsIniCond) |
---|
306 | { |
---|
307 | mapptctomad(kn->initial,buff); |
---|
308 | v->currentvalue = command_par_value(buff, madx_mpk_comm_calculate->clone); |
---|
309 | if (debuglevel) printf("Initialized current value for %s to %f\n", kn->initial,v->currentvalue); |
---|
310 | } |
---|
311 | else if(kn->exactnamematch == 0) |
---|
312 | { |
---|
313 | if (debuglevel) printf("Family here\n"); |
---|
314 | n = 0; |
---|
315 | node = current_sequ->range_start; |
---|
316 | |
---|
317 | while (node != 0x0) |
---|
318 | { |
---|
319 | strcpy(buff,node->name); |
---|
320 | |
---|
321 | p = strstr(buff,kn->elname); |
---|
322 | if ( p == buff ) |
---|
323 | { |
---|
324 | break; |
---|
325 | } |
---|
326 | |
---|
327 | n++; |
---|
328 | node = node->next; |
---|
329 | |
---|
330 | if ( node == current_sequ->range_start ) |
---|
331 | { |
---|
332 | fatal_error("readstartvalues: Can not find element: ",kn->elname); |
---|
333 | return 1; |
---|
334 | } |
---|
335 | } |
---|
336 | |
---|
337 | if (v->kn >=0) |
---|
338 | { |
---|
339 | ncomp = v->kn; |
---|
340 | w_ptc_getnfieldcomp_(&n,&ncomp, &nvalue); |
---|
341 | v->currentvalue = nvalue; |
---|
342 | } |
---|
343 | else |
---|
344 | { |
---|
345 | ncomp = v->ks; |
---|
346 | w_ptc_getsfieldcomp_(&n,&ncomp, &nvalue); |
---|
347 | v->currentvalue = nvalue; |
---|
348 | } |
---|
349 | if (debuglevel) printf("Got first element %s of family %s, field is %f\n",kn->elname,buff,v->currentvalue); |
---|
350 | |
---|
351 | n++; |
---|
352 | node = node->next; |
---|
353 | |
---|
354 | while (node != 0x0) |
---|
355 | { |
---|
356 | strcpy(buff,node->name); |
---|
357 | |
---|
358 | p = strstr(buff,kn->elname); |
---|
359 | if ( p == buff ) |
---|
360 | { |
---|
361 | if (debuglevel) printf("Got another element %s of the family %s\n",node->name,kn->elname); |
---|
362 | |
---|
363 | if (v->kn >=0) |
---|
364 | { |
---|
365 | ncomp = v->kn; |
---|
366 | w_ptc_getnfieldcomp_(&n,&ncomp, &nvalue); |
---|
367 | } |
---|
368 | else |
---|
369 | { |
---|
370 | ncomp = v->ks; |
---|
371 | w_ptc_getsfieldcomp_(&n,&ncomp, &nvalue); |
---|
372 | } |
---|
373 | |
---|
374 | if (v->currentvalue != nvalue) |
---|
375 | { |
---|
376 | warningnew("matchptcknobs", |
---|
377 | "Element %s has incoherent field %f strngth with its family %f.\n", |
---|
378 | node->name, nvalue,v->currentvalue); |
---|
379 | } |
---|
380 | } |
---|
381 | |
---|
382 | n++; |
---|
383 | node = node->next; |
---|
384 | |
---|
385 | if ( node == current_sequ->range_start ) |
---|
386 | { |
---|
387 | break; |
---|
388 | } |
---|
389 | } |
---|
390 | } |
---|
391 | else |
---|
392 | { |
---|
393 | n = 0; |
---|
394 | node = current_sequ->range_start; |
---|
395 | |
---|
396 | while (node != 0x0) |
---|
397 | { |
---|
398 | strcpy(buff,node->name); |
---|
399 | p = strstr(buff,":"); |
---|
400 | if (p) *p = 0; |
---|
401 | |
---|
402 | if ( strcmp(buff,kn->elname) == 0 ) |
---|
403 | { |
---|
404 | break; |
---|
405 | } |
---|
406 | |
---|
407 | n++; |
---|
408 | |
---|
409 | node = node->next; |
---|
410 | |
---|
411 | if ( node == current_sequ->range_start ) |
---|
412 | { |
---|
413 | fatal_error("readstartvalues: Can not find element: ",kn->elname); |
---|
414 | return 1; |
---|
415 | } |
---|
416 | } |
---|
417 | |
---|
418 | if (v->kn >=0) |
---|
419 | { |
---|
420 | ncomp = v->kn; |
---|
421 | w_ptc_getnfieldcomp_(&n,&ncomp, &nvalue); |
---|
422 | v->currentvalue = nvalue; |
---|
423 | } |
---|
424 | else |
---|
425 | { |
---|
426 | ncomp = v->ks; |
---|
427 | w_ptc_getsfieldcomp_(&n,&ncomp, &nvalue); |
---|
428 | v->currentvalue = nvalue; |
---|
429 | } |
---|
430 | |
---|
431 | if (debuglevel) printf("Got %f for %s\n",v->currentvalue,kn->elname); |
---|
432 | |
---|
433 | } |
---|
434 | } |
---|
435 | |
---|
436 | return 0; |
---|
437 | } |
---|
438 | |
---|
439 | static int |
---|
440 | madx_mpk_scalelimits(int nv) |
---|
441 | { |
---|
442 | /*User specifies uppper and lower limits of the matched k-values |
---|
443 | here they are converted to the PTC units */ |
---|
444 | |
---|
445 | /* |
---|
446 | double l; |
---|
447 | struct element* el = 0x0; |
---|
448 | char* p; |
---|
449 | struct node* node = 0x0; |
---|
450 | */ |
---|
451 | int fieldorder;/*PTC nomenclature, 1 dipol, 2 quad ...*/ |
---|
452 | float fact = 0.0; |
---|
453 | |
---|
454 | struct madx_mpk_variable* v = 0x0; |
---|
455 | // struct madx_mpk_knob* kn = 0x0; // not used |
---|
456 | |
---|
457 | if ( (nv < 0) || (nv >= MAX_KNOBS) ) |
---|
458 | { |
---|
459 | error("madx_mpk_scalelimits","Passed variable out of range"); |
---|
460 | return 1; |
---|
461 | } |
---|
462 | |
---|
463 | |
---|
464 | v = &madx_mpk_variables[nv]; |
---|
465 | // kn = &(madx_mpk_knobs[v->knobidx]); // not used |
---|
466 | |
---|
467 | fieldorder = 1 + (v->kn >= 0)?v->kn:v->ks;/*PTC nomenclature, 1 dipol, 2 quad ...*/ |
---|
468 | |
---|
469 | fact = factorial(fieldorder); |
---|
470 | |
---|
471 | if ( ( v->kn == 0) || (v->ks == 0) ) |
---|
472 | { |
---|
473 | printf("madx_mpk_scalelimits: Dipol limits don't need to be scaled\n"); |
---|
474 | return 1; |
---|
475 | } |
---|
476 | |
---|
477 | /* If lenght of the magnet will be needed, thin elements ??? |
---|
478 | * if(kn->exactnamematch) |
---|
479 | * { |
---|
480 | * el = find_element(kn->elname, element_list); |
---|
481 | * } |
---|
482 | * else |
---|
483 | * { |
---|
484 | * node = current_sequ->range_start; |
---|
485 | * while (node != 0x0) |
---|
486 | * { |
---|
487 | * p = strstr(node->name,kn->elname); |
---|
488 | * if ( p == node->name ) |
---|
489 | * { |
---|
490 | * el = node->p_elem; |
---|
491 | * break; |
---|
492 | * } |
---|
493 | * |
---|
494 | * node = node->next; |
---|
495 | * if ( node == current_sequ->range_start ) |
---|
496 | * { |
---|
497 | * fatal_error("madx_mpk_scalelimits: Can not find element starting with: ",kn->elname); |
---|
498 | * return 1; |
---|
499 | * } |
---|
500 | * } |
---|
501 | * |
---|
502 | * |
---|
503 | * } |
---|
504 | * if (el == 0x0) |
---|
505 | * { |
---|
506 | * error("madx_mpk_scalelimits","Can not find element named %s in the current sequence",madx_mpk_knobs[v->knobidx].elname); |
---|
507 | * return 1; |
---|
508 | * } |
---|
509 | * |
---|
510 | * if (debuglevel) printf("Element %s\n",el->name); |
---|
511 | * |
---|
512 | * l = el_par_value("l",el); |
---|
513 | * if (l <= 0.0) |
---|
514 | * { |
---|
515 | * printf("Element %s has zero lenght\n",el->name); |
---|
516 | * l = 1.0; zero lentgh elements has field defined such, compatible with madx_ptc_module.f90:SUMM_MULTIPOLES_AND_ERRORS |
---|
517 | * |
---|
518 | */ |
---|
519 | |
---|
520 | v->upper = v->upper/fact; |
---|
521 | v->lower = v->lower/fact; |
---|
522 | |
---|
523 | if (debuglevel) printf("Set limits to field order (PTC) %d, fact=%f: %f %f\n", |
---|
524 | fieldorder, fact, v->lower, v->upper ); |
---|
525 | |
---|
526 | return 0; |
---|
527 | |
---|
528 | } |
---|
529 | |
---|
530 | static void |
---|
531 | madx_mpk_addfieldcomp(struct madx_mpk_knob* knob, int kn, int ks) |
---|
532 | { |
---|
533 | /* |
---|
534 | adds a new field component to a knob |
---|
535 | */ |
---|
536 | void* ptr = 0x0; |
---|
537 | |
---|
538 | if (knob == 0x0) |
---|
539 | { |
---|
540 | warning("madx_mpk_addfieldcomp","knob parameter is null"); |
---|
541 | return; |
---|
542 | } |
---|
543 | |
---|
544 | if (kn >= 0) |
---|
545 | { |
---|
546 | knob->nKN++; |
---|
547 | ptr = (void*)knob->KN; |
---|
548 | knob->KN = (int*)realloc(knob->KN, knob->nKN * sizeof(int)); |
---|
549 | knob->KN[knob->nKN - 1] = kn; |
---|
550 | |
---|
551 | if (ptr != knob->KN) |
---|
552 | { |
---|
553 | free(ptr); |
---|
554 | } |
---|
555 | } |
---|
556 | |
---|
557 | if (ks >= 0) |
---|
558 | { |
---|
559 | knob->nKS++; |
---|
560 | ptr = (void*)knob->KS; |
---|
561 | knob->KS = (int*)realloc(knob->KS, knob->nKS * sizeof(int)); |
---|
562 | knob->KS[knob->nKS - 1] = ks; |
---|
563 | |
---|
564 | if (ptr != knob->KS) |
---|
565 | { |
---|
566 | free(ptr); |
---|
567 | } |
---|
568 | } |
---|
569 | |
---|
570 | } |
---|
571 | |
---|
572 | static void |
---|
573 | makestdmatchfile(char* fname, char* matchactioncommand) |
---|
574 | { |
---|
575 | FILE* f = 0x0; |
---|
576 | struct madx_mpk_variable* v = 0x0; |
---|
577 | int i; |
---|
578 | unsigned int anumber = abs(time(0)*rand()); |
---|
579 | double lower, upper; |
---|
580 | |
---|
581 | if (debuglevel) printf("I am in makestdmatchfile, file name is >%s<\n", fname); |
---|
582 | |
---|
583 | while(f == 0x0) |
---|
584 | {/*repeat until generate unique file name*/ |
---|
585 | if (fname[0] == 0) |
---|
586 | { /*if string is empty*/ |
---|
587 | sprintf(fname,"/tmp/match_ptcknobs_%u.madx",anumber); |
---|
588 | } |
---|
589 | f = fopen(fname,"w"); |
---|
590 | if (f == 0x0) |
---|
591 | { |
---|
592 | warningnew("makestdmatchfile","Could not open file %s",fname); |
---|
593 | } |
---|
594 | } |
---|
595 | |
---|
596 | /*if (debuglevel < 2) fprintf(f,"assign, echo=/tmp/mpk_stdmatch.out;\n");*/ |
---|
597 | |
---|
598 | fprintf(f,"match, use_macro;\n"); |
---|
599 | |
---|
600 | for (i = 0; i<madx_mpk_Nvariables; i++) |
---|
601 | { |
---|
602 | |
---|
603 | v = &(madx_mpk_variables[i]); |
---|
604 | |
---|
605 | /*create vary command as string*/ |
---|
606 | if (debuglevel) printf("\ncurrentvalue=%f trustrange=%f lower=%f upper=%f\n", |
---|
607 | v->currentvalue, v->trustrange , v->lower ,v->upper); |
---|
608 | |
---|
609 | if ( (v->currentvalue - v->trustrange) < v->lower ) |
---|
610 | { |
---|
611 | lower = v->currentvalue - v->lower; |
---|
612 | } |
---|
613 | else |
---|
614 | { |
---|
615 | lower = - v->trustrange; |
---|
616 | } |
---|
617 | |
---|
618 | if ( (v->currentvalue + v->trustrange) > v->upper ) |
---|
619 | { |
---|
620 | upper = v->upper - v->currentvalue; |
---|
621 | } |
---|
622 | else |
---|
623 | { |
---|
624 | upper = v->trustrange; |
---|
625 | } |
---|
626 | |
---|
627 | fprintf(f," vary, name=%s, step= %g, lower= %g, upper = %g;\n", |
---|
628 | v->name, v->step, lower, upper); |
---|
629 | |
---|
630 | } |
---|
631 | fprintf(f," \n"); |
---|
632 | |
---|
633 | fprintf(f," m1: macro = \n"); |
---|
634 | fprintf(f," {\n"); |
---|
635 | |
---|
636 | for (i = 0; i<madx_mpk_Nvariables; i++) |
---|
637 | { |
---|
638 | if (madx_mpk_variables[i].IsIniCond ) |
---|
639 | { |
---|
640 | fprintf(f," ptc_setknobvalue ,element=%s, value=%s, refreshtables=false;\n", |
---|
641 | madx_mpk_knobs[ madx_mpk_variables[i].knobidx ].initial, |
---|
642 | madx_mpk_variables[i].name); |
---|
643 | } |
---|
644 | else |
---|
645 | { |
---|
646 | fprintf(f," ptc_setknobvalue ,element=%s, kn=%d ,ks=%d, value=%s, refreshtables=false;\n", |
---|
647 | madx_mpk_knobs[ madx_mpk_variables[i].knobidx ].elname, |
---|
648 | madx_mpk_variables[i].kn, |
---|
649 | madx_mpk_variables[i].ks, |
---|
650 | madx_mpk_variables[i].name); |
---|
651 | } |
---|
652 | if (debuglevel) fprintf(f," value , %s, %s;\n", madx_mpk_variables[i].name, madx_mpk_variables[i].namecv ); |
---|
653 | } |
---|
654 | |
---|
655 | fprintf(f," ptc_refreshpartables;\n"); |
---|
656 | |
---|
657 | fprintf(f," };\n"); |
---|
658 | |
---|
659 | |
---|
660 | for (i = 0; i<madx_mpk_Nconstraints; i++) |
---|
661 | { |
---|
662 | fprintf(f," %s;\n",madx_mpk_constraints[i]); |
---|
663 | } |
---|
664 | |
---|
665 | |
---|
666 | fprintf(f," %s\n",matchactioncommand); |
---|
667 | |
---|
668 | fprintf(f,"endmatch;\n"); |
---|
669 | |
---|
670 | /*if (debuglevel < 2) fprintf(f,"assign, echo=terminal;\n");*/ |
---|
671 | |
---|
672 | fclose(f); |
---|
673 | |
---|
674 | } |
---|
675 | |
---|
676 | static int |
---|
677 | run_ptccalculation(int setknobs, char* readstartval) |
---|
678 | { |
---|
679 | int i,n; |
---|
680 | char buff[500]; |
---|
681 | char comd[500]; |
---|
682 | char* iniparname; |
---|
683 | |
---|
684 | char **toks=madx_mpk_comm_calculate->tok_list->p; |
---|
685 | int ntoks = madx_mpk_comm_calculate->tok_list->curr; |
---|
686 | int start; |
---|
687 | |
---|
688 | struct madx_mpk_variable* v = 0x0; |
---|
689 | struct madx_mpk_knob* kn = 0x0; |
---|
690 | struct node* node = 0x0; |
---|
691 | char* p; |
---|
692 | |
---|
693 | |
---|
694 | this_cmd = madx_mpk_comm_createuniverse; |
---|
695 | current_command = madx_mpk_comm_createuniverse->clone; |
---|
696 | process(); |
---|
697 | |
---|
698 | this_cmd = madx_mpk_comm_createlayout; |
---|
699 | current_command = madx_mpk_comm_createlayout->clone; |
---|
700 | process(); |
---|
701 | |
---|
702 | if (madx_mpk_comm_setswitch) |
---|
703 | { |
---|
704 | this_cmd = madx_mpk_comm_createlayout; |
---|
705 | current_command = madx_mpk_comm_createlayout->clone; |
---|
706 | process(); |
---|
707 | } |
---|
708 | |
---|
709 | if (*readstartval == 0) |
---|
710 | { |
---|
711 | for(i=0;i<madx_mpk_Nvariables;i++) |
---|
712 | { |
---|
713 | |
---|
714 | v = &(madx_mpk_variables[i]); |
---|
715 | kn = &(madx_mpk_knobs[v->knobidx]); |
---|
716 | set_variable_(v->namecv, &(v->currentvalue)); |
---|
717 | |
---|
718 | if (v->IsIniCond) |
---|
719 | { /*Set the initial parameter in ptc_twiss*/ |
---|
720 | |
---|
721 | iniparname = kn->initial; |
---|
722 | mapptctomad(iniparname,buff); |
---|
723 | |
---|
724 | for(start=0; start<ntoks; start++) |
---|
725 | { |
---|
726 | if (strcmp(toks[start],buff)==0) |
---|
727 | { |
---|
728 | break; |
---|
729 | } |
---|
730 | } |
---|
731 | |
---|
732 | set_command_par_value( buff, madx_mpk_comm_calculate->clone ,v->currentvalue); |
---|
733 | |
---|
734 | if (debuglevel) printf("Setting Initial %s to CV %f, now it is %f\n", |
---|
735 | buff,v->currentvalue, |
---|
736 | command_par_value(buff, madx_mpk_comm_calculate->clone)); |
---|
737 | |
---|
738 | } |
---|
739 | else |
---|
740 | { |
---|
741 | |
---|
742 | if(kn->exactnamematch == 0) |
---|
743 | { |
---|
744 | n = 0; |
---|
745 | node = current_sequ->range_start; |
---|
746 | |
---|
747 | while (node != 0x0) |
---|
748 | { |
---|
749 | strcpy(buff,node->name); |
---|
750 | |
---|
751 | p = strstr(buff,kn->elname); |
---|
752 | if ( p == buff ) |
---|
753 | { |
---|
754 | p = strstr(buff,":"); |
---|
755 | if (p) *p = 0; |
---|
756 | |
---|
757 | |
---|
758 | sprintf(comd,"ptc_setfieldcomp, element=%s, kn=%d, ks=%d, value=%s;", |
---|
759 | buff, v->kn,v->ks,v->namecv); |
---|
760 | if (debuglevel) printf("%s\n",comd); |
---|
761 | pro_input_(comd); |
---|
762 | |
---|
763 | } |
---|
764 | |
---|
765 | n++; |
---|
766 | node = node->next; |
---|
767 | |
---|
768 | if ( node == current_sequ->range_start ) |
---|
769 | { |
---|
770 | break; |
---|
771 | } |
---|
772 | } |
---|
773 | |
---|
774 | } |
---|
775 | else |
---|
776 | { |
---|
777 | sprintf(comd,"ptc_setfieldcomp, element=%s, kn=%d, ks=%d, value=%s;", |
---|
778 | kn->elname, v->kn,v->ks,v->namecv); |
---|
779 | if (debuglevel) printf("%s\n",comd); |
---|
780 | pro_input_(comd); |
---|
781 | } |
---|
782 | } |
---|
783 | } |
---|
784 | } |
---|
785 | |
---|
786 | if(setknobs) |
---|
787 | { |
---|
788 | |
---|
789 | for(i=0;i<madx_mpk_Nknobs;i++) |
---|
790 | { |
---|
791 | /* printf("Setting knob %d: \n%s\n",i,madx_mpk_setknobs[i]);*/ |
---|
792 | pro_input_(madx_mpk_setknobs[i]); |
---|
793 | } |
---|
794 | } |
---|
795 | else |
---|
796 | { |
---|
797 | if (debuglevel) printf("Knob Setting Is not requested this time.\n"); |
---|
798 | } |
---|
799 | |
---|
800 | /* pro_input_(twisscommand);*/ |
---|
801 | |
---|
802 | if (debuglevel) printf("Running ptc_twiss or ptc_normal.\n"); |
---|
803 | |
---|
804 | this_cmd = madx_mpk_comm_calculate; |
---|
805 | current_command = madx_mpk_comm_calculate->clone; |
---|
806 | current_twiss = current_command; |
---|
807 | |
---|
808 | pro_ptc_twiss(); |
---|
809 | /* process();*/ |
---|
810 | |
---|
811 | if ( *readstartval ) |
---|
812 | { |
---|
813 | readstartvalues(); |
---|
814 | *readstartval = 0; |
---|
815 | } |
---|
816 | |
---|
817 | return 0; |
---|
818 | } |
---|
819 | |
---|
820 | static double |
---|
821 | match2_summary(void) |
---|
822 | { |
---|
823 | int i,j; |
---|
824 | |
---|
825 | double penalty; |
---|
826 | |
---|
827 | printf("\n"); |
---|
828 | printf("MATCH SUMMARY\n\n"); |
---|
829 | /* fprintf(prt_file, "Macro Constraint Value Penalty\n");*/ |
---|
830 | printf("--------------------------------------------------------------------\n"); |
---|
831 | penalty=0; |
---|
832 | for(i=0;i<MAX_MATCH_MACRO;i++) |
---|
833 | { |
---|
834 | if(match2_macro_name[i]==NULL) break; |
---|
835 | printf("macro: %-20s\n",match2_macro_name[i]); |
---|
836 | for(j=0;j<MAX_MATCH_CONS;j++) |
---|
837 | { |
---|
838 | if (match2_cons_name[i][j]==NULL) break; |
---|
839 | printf(" constraint: %-40s\n",match2_cons_name[i][j]); |
---|
840 | printf(" values: %+12.5e%c%+12.5e\n", |
---|
841 | match2_cons_value_lhs[i][j], |
---|
842 | match2_cons_sign[i][j], |
---|
843 | match2_cons_value_rhs[i][j]); |
---|
844 | printf(" weight: %+12.5e\n", match2_cons_weight[i][j]); |
---|
845 | printf(" penalty: %+12.5e\n\n",match2_cons_value[i][j]); |
---|
846 | penalty+=pow(match2_cons_value[i][j],2); |
---|
847 | } |
---|
848 | } |
---|
849 | |
---|
850 | return penalty; |
---|
851 | } |
---|
852 | |
---|
853 | static int |
---|
854 | madx_mpk_init(void) |
---|
855 | { |
---|
856 | /*performs initialization, |
---|
857 | prepares commands for exacution out of gethered user information*/ |
---|
858 | |
---|
859 | /*Create ptc_command: just copy paste the valid parameters - no error checking*/ |
---|
860 | |
---|
861 | char vary[COMM_LENGTH]; |
---|
862 | int i,k; |
---|
863 | |
---|
864 | for(i=0; i< madx_mpk_Nknobs; i++) |
---|
865 | { |
---|
866 | sprintf(vary,"ptc_knob, "); |
---|
867 | |
---|
868 | if (madx_mpk_knobs[i].elname) |
---|
869 | { |
---|
870 | sprintf(&(vary[strlen(vary)]),"element=%s, ",madx_mpk_knobs[i].elname); |
---|
871 | |
---|
872 | if (madx_mpk_knobs[i].nKN > 0) |
---|
873 | { |
---|
874 | |
---|
875 | sprintf(&(vary[strlen(vary)]),"kn="); |
---|
876 | |
---|
877 | for (k=0; k < madx_mpk_knobs[i].nKN; k++) |
---|
878 | { |
---|
879 | sprintf(&(vary[strlen(vary)]),"%d,",madx_mpk_knobs[i].KN[k]); |
---|
880 | } |
---|
881 | } |
---|
882 | |
---|
883 | if (madx_mpk_knobs[i].nKS > 0) |
---|
884 | { |
---|
885 | |
---|
886 | sprintf(&(vary[strlen(vary)]),"ks="); |
---|
887 | |
---|
888 | for (k=0; k < madx_mpk_knobs[i].nKS; k++) |
---|
889 | { |
---|
890 | sprintf(&(vary[strlen(vary)]),"%d,",madx_mpk_knobs[i].KS[k]); |
---|
891 | } |
---|
892 | } |
---|
893 | |
---|
894 | /*this is the last one anyway if ename is specified ; at the end*/ |
---|
895 | if (madx_mpk_knobs[i].exactnamematch) |
---|
896 | { |
---|
897 | sprintf(&(vary[strlen(vary)]),"exactmatch=%s; ","true"); |
---|
898 | } |
---|
899 | else |
---|
900 | { |
---|
901 | sprintf(&(vary[strlen(vary)]),"exactmatch=%s; ","false"); |
---|
902 | } |
---|
903 | |
---|
904 | } |
---|
905 | |
---|
906 | if (madx_mpk_knobs[i].initial) |
---|
907 | { |
---|
908 | sprintf(&(vary[strlen(vary)]),"initial=%s; ",madx_mpk_knobs[i].initial); |
---|
909 | } |
---|
910 | |
---|
911 | madx_mpk_setknobs[i] = (char*)mymalloc("madx_mpk_addvariable",1+strlen(vary)*sizeof(char)); |
---|
912 | strcpy(madx_mpk_setknobs[i],vary); |
---|
913 | |
---|
914 | if (debuglevel) printf("madx_mpk_setknobs[%d]= %s\n",i,madx_mpk_setknobs[i]); |
---|
915 | } |
---|
916 | |
---|
917 | return 0; |
---|
918 | } |
---|
919 | |
---|
920 | // public interface |
---|
921 | |
---|
922 | void |
---|
923 | madx_mpk_run(struct in_cmd* cmd) |
---|
924 | { |
---|
925 | /*the main routine of the module, called after at matching action (migrad, lmdif. etc...) */ |
---|
926 | char rout_name[] = "madx_mpk_run"; |
---|
927 | int i; |
---|
928 | double tolerance; |
---|
929 | int calls = 0, maxcalls; |
---|
930 | double ktar, penalty, penalty2, oldtar = 0.0, tar; |
---|
931 | double var; |
---|
932 | /* double* matchedvalues;*/ |
---|
933 | double* function_vector1 = 0x0; |
---|
934 | double* function_vector2 = 0x0; |
---|
935 | char ptcend[20]; |
---|
936 | char callmatchfile[200]; |
---|
937 | char matchfilename[300]; |
---|
938 | char matchactioncommand[400]; |
---|
939 | char maxNCallsExceeded = 0; |
---|
940 | char knobsatmatchedpoint = 0; /*flag indicating that matched values of knobs are with precision close to zero (expansion if valid only close to )*/ |
---|
941 | char matched2 = 0; |
---|
942 | char matched = 0; |
---|
943 | char readstartval = 1; |
---|
944 | // int retval; // not used |
---|
945 | int fieldorder; |
---|
946 | FILE* fdbg = fopen("currentvalues.txt","a+"); |
---|
947 | |
---|
948 | matchfilename[0] = 0; /*sign for makestdmatchfile to generate a new file name*/ |
---|
949 | |
---|
950 | tolerance = command_par_value("tolerance",cmd->clone); |
---|
951 | if (tolerance < 0) |
---|
952 | { |
---|
953 | warningnew("matchknobs.c: madx_mpk_run","Tolerance is less then 0. Command ignored."); |
---|
954 | return; |
---|
955 | } |
---|
956 | maxcalls = command_par_value("calls",cmd->clone); |
---|
957 | |
---|
958 | if(madx_mpk_comm_createuniverse == 0x0) |
---|
959 | { |
---|
960 | warningnew("matchknobs.c: madx_mpk_run","ptc_createuniverse command is missing."); |
---|
961 | return; |
---|
962 | } |
---|
963 | |
---|
964 | if(madx_mpk_comm_createlayout == 0x0) |
---|
965 | { |
---|
966 | warningnew("matchknobs.c: madx_mpk_run","ptc_createlayout is missing."); |
---|
967 | return; |
---|
968 | } |
---|
969 | |
---|
970 | if(madx_mpk_comm_calculate == 0x0) |
---|
971 | { |
---|
972 | warningnew("matchknobs.c: madx_mpk_run","Neither ptc_twiss nor ptc_normal seen since \"match, use_ptcknob;\""); |
---|
973 | return; |
---|
974 | } |
---|
975 | |
---|
976 | if (madx_mpk_Nknobs == 0) |
---|
977 | { |
---|
978 | warningnew("matchknobs.c: madx_mpk_run","no knobs seen yet."); |
---|
979 | return; |
---|
980 | } |
---|
981 | |
---|
982 | if (madx_mpk_Nvariables == 0) |
---|
983 | { |
---|
984 | warningnew("matchknobs.c: madx_mpk_run","no variables seen yet."); |
---|
985 | return; |
---|
986 | } |
---|
987 | |
---|
988 | |
---|
989 | sprintf(matchactioncommand,"%s,",cmd->tok_list->p[0]); /*there is no comma after the command name*/ |
---|
990 | i = 1; |
---|
991 | while(cmd->tok_list->p[i] != 0x0) |
---|
992 | { |
---|
993 | sprintf(&(matchactioncommand[strlen(matchactioncommand)])," %s",cmd->tok_list->p[i]); |
---|
994 | if (debuglevel) printf("%s",cmd->tok_list->p[i]); |
---|
995 | i++; |
---|
996 | } |
---|
997 | |
---|
998 | sprintf(&(matchactioncommand[strlen(matchactioncommand)]),";"); |
---|
999 | |
---|
1000 | madx_mpk_init(); // retval = not used |
---|
1001 | |
---|
1002 | |
---|
1003 | |
---|
1004 | strcpy(ptcend,"ptc_end;"); |
---|
1005 | /*check if ptc starting commands are fine */ |
---|
1006 | match_is_on = kMatch_NoMatch; |
---|
1007 | |
---|
1008 | match2_keepexpressions = 1; |
---|
1009 | |
---|
1010 | fprintf(fdbg,"\n###############################################\n"); |
---|
1011 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1012 | { |
---|
1013 | fprintf(fdbg,"%16s ",madx_mpk_variables[i].name); |
---|
1014 | } |
---|
1015 | fprintf(fdbg,"\n"); fflush(0); |
---|
1016 | |
---|
1017 | do |
---|
1018 | { |
---|
1019 | calls++; |
---|
1020 | |
---|
1021 | maxNCallsExceeded = ( (maxcalls > 0)&&(calls > maxcalls) ); |
---|
1022 | |
---|
1023 | printf("###########################################################################################\n"); |
---|
1024 | printf("\n\n\n Call %d\n",calls); |
---|
1025 | printf("###########################################################################################\n"); |
---|
1026 | |
---|
1027 | printf("\n CURRENT VALUES \n"); |
---|
1028 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1029 | { |
---|
1030 | printf("%16s ",madx_mpk_variables[i].name); |
---|
1031 | } |
---|
1032 | printf("\n"); |
---|
1033 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1034 | { |
---|
1035 | printf("%16f ",madx_mpk_variables[i].currentvalue); |
---|
1036 | } |
---|
1037 | printf("\n"); |
---|
1038 | |
---|
1039 | fflush(0); |
---|
1040 | |
---|
1041 | if (debuglevel) printf("\n\nCalling TWISS with KNOBS\n"); |
---|
1042 | run_ptccalculation(1,&readstartval); |
---|
1043 | |
---|
1044 | if (geterrorflag()) |
---|
1045 | { |
---|
1046 | error("Matching With Knobs","PTC calculation ended with an error. Check your setting and matching limits."); |
---|
1047 | pro_input_(ptcend); |
---|
1048 | goto cleaning; |
---|
1049 | } |
---|
1050 | |
---|
1051 | |
---|
1052 | /*set all variables to 0*/ |
---|
1053 | var = 0.0; |
---|
1054 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1055 | { |
---|
1056 | set_variable_(madx_mpk_variables[i].name, &var); |
---|
1057 | } |
---|
1058 | |
---|
1059 | if (debuglevel) printf("\n\nDoing Match on KNOBS\n"); |
---|
1060 | makestdmatchfile(matchfilename,matchactioncommand); |
---|
1061 | if (debuglevel) printf("matchfilename is %s\n",matchfilename); |
---|
1062 | |
---|
1063 | sprintf(callmatchfile,"call, file=\"%s\" ;",matchfilename); |
---|
1064 | pro_input_(callmatchfile); |
---|
1065 | |
---|
1066 | if (function_vector1 == 0x0) |
---|
1067 | { |
---|
1068 | function_vector1 = (double*)mymalloc("madx_mpk_run",(total_const+1)*sizeof(double)); |
---|
1069 | function_vector2 = (double*)mymalloc("madx_mpk_run",(total_const+1)*sizeof(double)); |
---|
1070 | } |
---|
1071 | |
---|
1072 | i = match2_evaluate_exressions(0,0,function_vector1); |
---|
1073 | if (debuglevel) printf("match2_evaluate_exressions returned fun_vector 1 of length %d\n",i); |
---|
1074 | |
---|
1075 | |
---|
1076 | pro_input_(ptcend); |
---|
1077 | |
---|
1078 | /* END OF STD MATCHING*/ |
---|
1079 | |
---|
1080 | /* CALCULATE KTAR*/ |
---|
1081 | |
---|
1082 | ktar = 0.0; |
---|
1083 | fprintf(fdbg,"\n%d V ",calls); |
---|
1084 | printf("MPK STEP %d",calls); |
---|
1085 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1086 | { |
---|
1087 | var = get_variable_(madx_mpk_variables[i].name); |
---|
1088 | fprintf(fdbg," %20.16f ", var); |
---|
1089 | ktar += var*var; |
---|
1090 | } |
---|
1091 | fprintf(fdbg,"\n"); |
---|
1092 | |
---|
1093 | |
---|
1094 | /* BIG KTAR*/ |
---|
1095 | /*BIG TAR*/ |
---|
1096 | /*SCALE AND BRING BACK OLD VALUES*/ |
---|
1097 | |
---|
1098 | tar = get_variable_("tar"); |
---|
1099 | |
---|
1100 | if ( (oldtar <= 0) || (tar < 10.0*oldtar)) /*correct within order of magnitude - we can not truest the oldtar */ |
---|
1101 | { /*because it also bears error of approximation*/ |
---|
1102 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1103 | { |
---|
1104 | madx_mpk_variables[i].oldvalue = madx_mpk_variables[i].currentvalue; |
---|
1105 | |
---|
1106 | var = get_variable_(madx_mpk_variables[i].name); |
---|
1107 | madx_mpk_variables[i].currentvalue += var; |
---|
1108 | } |
---|
1109 | } |
---|
1110 | else |
---|
1111 | { |
---|
1112 | fprintf(fdbg,"Narrowing trust ranges Bringing back the previous step values\n"); |
---|
1113 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1114 | { |
---|
1115 | madx_mpk_variables[i].trustrange /= 2.0; |
---|
1116 | madx_mpk_variables[i].step /= 2.0; |
---|
1117 | |
---|
1118 | madx_mpk_variables[i].currentvalue = madx_mpk_variables[i].oldvalue; |
---|
1119 | |
---|
1120 | } |
---|
1121 | |
---|
1122 | } |
---|
1123 | |
---|
1124 | /*SMALL TAR*/ |
---|
1125 | /*ADD DELTAS*/ |
---|
1126 | |
---|
1127 | /*CONTINUE*/ |
---|
1128 | |
---|
1129 | /* SMALL KTAR*/ |
---|
1130 | /*VERIFY*/ |
---|
1131 | |
---|
1132 | |
---|
1133 | fprintf(fdbg,"%d CV",calls); |
---|
1134 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1135 | { |
---|
1136 | fprintf(fdbg," %20.16f ", madx_mpk_variables[i].currentvalue); |
---|
1137 | } |
---|
1138 | fprintf(fdbg,"\n"); |
---|
1139 | |
---|
1140 | fprintf(fdbg," KTAR=%20.16E TAR=%20.16E OLDTAR=%20.16E \n",ktar, tar, oldtar); |
---|
1141 | printf(" KTAR=%20.16E TAR=%20.16E OLDTAR=%20.16E \n",ktar, tar, oldtar); |
---|
1142 | |
---|
1143 | if (debuglevel) printf("KTAR=%E \n", ktar ); |
---|
1144 | |
---|
1145 | fflush(0x0); |
---|
1146 | |
---|
1147 | |
---|
1148 | oldtar = tar; |
---|
1149 | |
---|
1150 | if (ktar < tolerance) |
---|
1151 | { |
---|
1152 | knobsatmatchedpoint = 1; |
---|
1153 | if (debuglevel) printf("Matched values of knobs are close to zeroes within the tolerance\n"); |
---|
1154 | } |
---|
1155 | else |
---|
1156 | { |
---|
1157 | |
---|
1158 | knobsatmatchedpoint = 0; |
---|
1159 | if (debuglevel) printf("Matched values of knobs are NOT close to zeroes within the tolerance\n"); |
---|
1160 | if (!maxNCallsExceeded) continue; |
---|
1161 | } |
---|
1162 | |
---|
1163 | |
---|
1164 | |
---|
1165 | if (debuglevel) printf("\n\nCalling TWISS with NO knobs\n"); |
---|
1166 | run_ptccalculation(0,&readstartval); |
---|
1167 | |
---|
1168 | i = match2_evaluate_exressions(0,0,function_vector2); |
---|
1169 | if (debuglevel) printf("match2_evaluate_exressions returned fun_vector 2 of length %d\n",i); |
---|
1170 | |
---|
1171 | pro_input_(ptcend); |
---|
1172 | |
---|
1173 | penalty = 0.0; |
---|
1174 | penalty2 = 0.0; |
---|
1175 | |
---|
1176 | for (i=0; i< total_const; i++) |
---|
1177 | { |
---|
1178 | penalty += function_vector2[i]*function_vector2[i]; |
---|
1179 | var = function_vector2[i] - function_vector1[i]; |
---|
1180 | /*if (debuglevel>1) */ |
---|
1181 | if (debuglevel) printf("Constraint%d: v1=%f v2=%f diff=%f \n",i,function_vector1[i],function_vector2[i],var); |
---|
1182 | penalty2 += var*var; |
---|
1183 | } |
---|
1184 | |
---|
1185 | if (debuglevel > 1) match2_summary(); |
---|
1186 | |
---|
1187 | fprintf(fdbg,"Verification: penalty=%E penalty2=%E \n",penalty, penalty2); |
---|
1188 | |
---|
1189 | if (debuglevel) printf("Closenes of the function vectors with and without knobs %e\n",penalty2); |
---|
1190 | |
---|
1191 | if (penalty2 < tolerance) |
---|
1192 | { |
---|
1193 | matched2 = 1; |
---|
1194 | if (debuglevel) printf("Constraints values with and without knobs within the tolerance\n"); |
---|
1195 | } |
---|
1196 | else |
---|
1197 | { |
---|
1198 | matched2 = 0; |
---|
1199 | if (debuglevel) printf("Constraints values with and without knobs NOT within the tolerance\n"); |
---|
1200 | } |
---|
1201 | |
---|
1202 | matched = matched2 && knobsatmatchedpoint; |
---|
1203 | |
---|
1204 | |
---|
1205 | if (debuglevel) printf("matched=%d, maxNCallsExceeded=%d\n",matched,maxNCallsExceeded); |
---|
1206 | |
---|
1207 | } |
---|
1208 | while( ( !(matched || maxNCallsExceeded )) && (ktar != 0.0) ); |
---|
1209 | /*if ktar is 0.0 we can not make better, because next iteration will be the same*/ |
---|
1210 | match2_keepexpressions = 0; |
---|
1211 | match_is_on = kMatch_PTCknobs; |
---|
1212 | |
---|
1213 | fflush(0); |
---|
1214 | tar = match2_summary(); |
---|
1215 | printf("\n\n"); |
---|
1216 | printf("================================================================================== \n"); |
---|
1217 | printf("== Matching finished successfully after %3.3d calls == \n",calls); |
---|
1218 | printf("== Final penalty function TAR = %E == \n", tar); |
---|
1219 | |
---|
1220 | printf("== Precision of Taylor expansion at the result point = %E ==\n",ktar); |
---|
1221 | printf("==------------------------------------------------------------------------------== \n"); |
---|
1222 | printf("== Variable Final Value == \n"); |
---|
1223 | printf("\n"); |
---|
1224 | |
---|
1225 | |
---|
1226 | for (i=0;i<madx_mpk_Nvariables;i++) |
---|
1227 | { |
---|
1228 | set_variable_(madx_mpk_variables[i].name,&(madx_mpk_variables[i].currentvalue)); |
---|
1229 | printf("== %2d %16s %+12.8E", |
---|
1230 | i,madx_mpk_variables[i].name, madx_mpk_variables[i].currentvalue); |
---|
1231 | |
---|
1232 | if (madx_mpk_variables[i].IsIniCond) |
---|
1233 | { |
---|
1234 | printf(" == \n"); |
---|
1235 | } |
---|
1236 | else |
---|
1237 | { |
---|
1238 | fieldorder = 0; |
---|
1239 | if (madx_mpk_variables[i].kn >= 0) |
---|
1240 | { |
---|
1241 | fieldorder = factorial(madx_mpk_variables[i].kn); |
---|
1242 | } |
---|
1243 | if (madx_mpk_variables[i].ks >= 0) |
---|
1244 | { |
---|
1245 | fieldorder = factorial(madx_mpk_variables[i].ks); |
---|
1246 | } |
---|
1247 | |
---|
1248 | if (fieldorder >= 0) |
---|
1249 | { |
---|
1250 | |
---|
1251 | printf(" PTC units -> MADX: %+12.8E == \n", |
---|
1252 | fieldorder*madx_mpk_variables[i].currentvalue); |
---|
1253 | |
---|
1254 | } |
---|
1255 | else |
---|
1256 | { |
---|
1257 | printf(" == \n"); |
---|
1258 | } |
---|
1259 | } |
---|
1260 | printf("\n"); |
---|
1261 | } |
---|
1262 | |
---|
1263 | printf("==========================================================================\n"); |
---|
1264 | |
---|
1265 | cleaning: |
---|
1266 | if (function_vector1) myfree(rout_name,function_vector1); |
---|
1267 | if (function_vector2) myfree(rout_name,function_vector2); |
---|
1268 | |
---|
1269 | function_vector1 = 0x0; |
---|
1270 | function_vector2 = 0x0; |
---|
1271 | |
---|
1272 | fclose(fdbg); |
---|
1273 | /* remove(matchfilename); */ |
---|
1274 | } |
---|
1275 | |
---|
1276 | void |
---|
1277 | madx_mpk_prepare(void) |
---|
1278 | { |
---|
1279 | /* prepares the internal variables*/ |
---|
1280 | int i; /**/ |
---|
1281 | |
---|
1282 | madx_mpk_Nconstraints = 0; |
---|
1283 | madx_mpk_Nknobs = 0; |
---|
1284 | madx_mpk_Nvariables = 0; |
---|
1285 | |
---|
1286 | madx_mpk_comm_createuniverse = 0x0; |
---|
1287 | madx_mpk_comm_createlayout = 0x0; |
---|
1288 | madx_mpk_comm_setswitch = 0x0; |
---|
1289 | madx_mpk_comm_calculate = 0x0; |
---|
1290 | |
---|
1291 | for (i = 0; i< MAX_KNOBS; i++) |
---|
1292 | { |
---|
1293 | |
---|
1294 | madx_mpk_knobs[i].elname = 0x0; |
---|
1295 | madx_mpk_knobs[i].initial = 0x0; |
---|
1296 | madx_mpk_knobs[i].KN = 0x0; |
---|
1297 | madx_mpk_knobs[i].nKN = 0; |
---|
1298 | madx_mpk_knobs[i].KS = 0x0; |
---|
1299 | madx_mpk_knobs[i].nKS = 0; |
---|
1300 | madx_mpk_knobs[i].exactnamematch = 1; |
---|
1301 | |
---|
1302 | madx_mpk_variables[i].name = 0x0; |
---|
1303 | madx_mpk_variables[i].namecv = 0x0; |
---|
1304 | madx_mpk_variables[i].upper = 0.0; |
---|
1305 | madx_mpk_variables[i].lower = 0.0; |
---|
1306 | madx_mpk_variables[i].trustrange = 0.0; |
---|
1307 | madx_mpk_variables[i].step = 0.0; |
---|
1308 | madx_mpk_variables[i].knobidx = -1; |
---|
1309 | madx_mpk_variables[i].currentvalue = 0.0; |
---|
1310 | |
---|
1311 | madx_mpk_variables[i].kn = -1; |
---|
1312 | madx_mpk_variables[i].ks = -1; |
---|
1313 | madx_mpk_variables[i].IsIniCond = -1; |
---|
1314 | } |
---|
1315 | } |
---|
1316 | |
---|
1317 | void |
---|
1318 | madx_mpk_addconstraint(const char* constr) |
---|
1319 | { |
---|
1320 | char* buff; |
---|
1321 | int l; |
---|
1322 | if (constr == 0x0) return; |
---|
1323 | l = strlen(constr); |
---|
1324 | if (l<1) return; |
---|
1325 | l++; |
---|
1326 | buff = (char*)mymalloc("madx_mpk_addconstraint",l*sizeof(char)); |
---|
1327 | strcpy(buff,constr); |
---|
1328 | madx_mpk_constraints[madx_mpk_Nconstraints++] = buff; |
---|
1329 | |
---|
1330 | } |
---|
1331 | |
---|
1332 | void |
---|
1333 | madx_mpk_addvariable(struct in_cmd* cmd) |
---|
1334 | { |
---|
1335 | char vary[COMM_LENGTH]; |
---|
1336 | int slen; |
---|
1337 | char* string; |
---|
1338 | char* ename; |
---|
1339 | char* initialpar; |
---|
1340 | int exactnamematch;/*0 for families with name starting with ename, 1 for for element having name ename*/ |
---|
1341 | int kn,ks; |
---|
1342 | /* int int_arr[100];*/ |
---|
1343 | /* int i;*/ |
---|
1344 | int knobidx; /*index of setknob command for this varyknob*/ |
---|
1345 | |
---|
1346 | struct madx_mpk_variable* v = 0x0; |
---|
1347 | |
---|
1348 | if (cmd == 0x0) return; |
---|
1349 | |
---|
1350 | /*get a name for the variable*/ |
---|
1351 | ename = command_par_string("element",cmd->clone); |
---|
1352 | initialpar = command_par_string("initial",cmd->clone); |
---|
1353 | if (( ename == 0x0 ) && ( initialpar == 0x0 )) |
---|
1354 | { |
---|
1355 | error("matchknobs.c: madx_mpk_addvariable", |
---|
1356 | "Neither element nor initial parameter specified. Command ignored!"); |
---|
1357 | return; |
---|
1358 | } |
---|
1359 | |
---|
1360 | if ( ename && initialpar) |
---|
1361 | { |
---|
1362 | error("matchknobs.c: madx_mpk_addvariable", |
---|
1363 | "Single command may define only one of two, field component or initial parameter. Command ignored!"); |
---|
1364 | return; |
---|
1365 | } |
---|
1366 | |
---|
1367 | kn = command_par_value("kn",cmd->clone); |
---|
1368 | ks = command_par_value("ks",cmd->clone); |
---|
1369 | |
---|
1370 | if ( ename && (kn >= 0) && (ks >=0) ) |
---|
1371 | { |
---|
1372 | error("matchknobs.c: madx_mpk_addvariable", |
---|
1373 | "Single command may define only one field component, not ks and kn together. Command ignored."); |
---|
1374 | return; |
---|
1375 | } |
---|
1376 | |
---|
1377 | exactnamematch = command_par_value("exactmatch",cmd->clone); |
---|
1378 | |
---|
1379 | knobidx = findsetknob(ename,exactnamematch,initialpar); |
---|
1380 | |
---|
1381 | if (knobidx < 0) |
---|
1382 | { |
---|
1383 | error("madx_mpk_addvariable","Error occured while adding this variable."); |
---|
1384 | return; |
---|
1385 | } |
---|
1386 | |
---|
1387 | /*so we make a new variable...*/ |
---|
1388 | |
---|
1389 | v = &(madx_mpk_variables[madx_mpk_Nvariables]); |
---|
1390 | |
---|
1391 | v->upper = command_par_value("upper",cmd->clone); |
---|
1392 | v->lower = command_par_value("lower",cmd->clone); |
---|
1393 | v->trustrange = command_par_value("trustrange",cmd->clone); |
---|
1394 | v->step = command_par_value("step",cmd->clone); |
---|
1395 | v->knobidx = knobidx; |
---|
1396 | |
---|
1397 | if (initialpar) |
---|
1398 | { |
---|
1399 | madx_mpk_knobs[knobidx].initial = (char*)mycalloc("madx_mpk_addvariable",1+strlen(initialpar),sizeof(char)); |
---|
1400 | strcpy(madx_mpk_knobs[knobidx].initial, initialpar); |
---|
1401 | |
---|
1402 | sprintf(vary,"mpk_%s",initialpar); |
---|
1403 | v->name = (char*)mycalloc("madx_mpk_addvariable",1+strlen(vary),sizeof(char)); |
---|
1404 | strcpy(v->name,vary); |
---|
1405 | |
---|
1406 | sprintf(vary,"mpk_%s_0",initialpar); |
---|
1407 | v->namecv = (char*)mycalloc("madx_mpk_addvariable",1+strlen(vary),sizeof(char)); |
---|
1408 | strcpy(v->namecv,vary); |
---|
1409 | v->IsIniCond = 1; |
---|
1410 | v->kn = -1; |
---|
1411 | v->ks = -1; |
---|
1412 | } |
---|
1413 | |
---|
1414 | if (ename) |
---|
1415 | { |
---|
1416 | madx_mpk_knobs[knobidx].elname = (char*)mymalloc("madx_mpk_addvariable",1+strlen(ename)*sizeof(char)); |
---|
1417 | strcpy(madx_mpk_knobs[knobidx].elname, ename); |
---|
1418 | |
---|
1419 | |
---|
1420 | if (exactnamematch != 0) |
---|
1421 | { |
---|
1422 | madx_mpk_knobs[knobidx].exactnamematch = 1; |
---|
1423 | } |
---|
1424 | else |
---|
1425 | { |
---|
1426 | madx_mpk_knobs[knobidx].exactnamematch = 0; |
---|
1427 | } |
---|
1428 | |
---|
1429 | if (kn >=0) |
---|
1430 | { |
---|
1431 | /*create variable name, use vary as a temp buffer */ |
---|
1432 | sprintf(vary,"%s_kn%d",ename,kn); |
---|
1433 | madx_mpk_addfieldcomp(&(madx_mpk_knobs[knobidx]), kn, -1); |
---|
1434 | } |
---|
1435 | |
---|
1436 | if (ks >=0) |
---|
1437 | { |
---|
1438 | /*create variable name, use vary as a temp buffer */ |
---|
1439 | sprintf(vary,"%s_ks%d",ename,kn); |
---|
1440 | madx_mpk_addfieldcomp(&(madx_mpk_knobs[knobidx]), -1, ks); |
---|
1441 | } |
---|
1442 | |
---|
1443 | |
---|
1444 | slen = strlen(vary); |
---|
1445 | string = (char*)mymalloc("madx_mpk_addvariable",(slen+1)*sizeof(char)); |
---|
1446 | strcpy(string,vary); |
---|
1447 | v->name = string; |
---|
1448 | |
---|
1449 | vary[slen ] = '_'; |
---|
1450 | vary[slen+1] = '0'; |
---|
1451 | vary[slen+2] = 0;/*end of string*/ |
---|
1452 | slen = slen+2; |
---|
1453 | |
---|
1454 | string = (char*)mymalloc("madx_mpk_addvariable",(slen+1)*sizeof(char)); |
---|
1455 | strcpy(string,vary); |
---|
1456 | v->namecv = string; |
---|
1457 | |
---|
1458 | v->kn = kn; |
---|
1459 | v->ks = ks; |
---|
1460 | v->IsIniCond = 0; |
---|
1461 | |
---|
1462 | madx_mpk_scalelimits(madx_mpk_Nvariables); |
---|
1463 | |
---|
1464 | } |
---|
1465 | |
---|
1466 | madx_mpk_Nvariables++; |
---|
1467 | |
---|
1468 | if (debuglevel) printf("Added new knobs: now there is\n knobs: %d\n variables: %d\n", |
---|
1469 | madx_mpk_Nknobs,madx_mpk_Nvariables); |
---|
1470 | } |
---|
1471 | |
---|
1472 | void |
---|
1473 | madx_mpk_end(void) |
---|
1474 | { |
---|
1475 | /* |
---|
1476 | remove_from_command_list(this_cmd->clone->name,stored_commands); |
---|
1477 | */ |
---|
1478 | match_is_on = kMatch_NoMatch; |
---|
1479 | } |
---|
1480 | |
---|
1481 | void |
---|
1482 | madx_mpk_setcreateuniverse(struct in_cmd* cmd) |
---|
1483 | { |
---|
1484 | cmd->clone_flag = 1; /* do not delete for the moment*/ |
---|
1485 | cmd->label = (char*)mymalloc("madx_mpk_setcreateuniverse",24*sizeof(char)); |
---|
1486 | strcpy(cmd->label,"matchptcknob_ptc_CU"); |
---|
1487 | |
---|
1488 | madx_mpk_comm_createuniverse = cmd; |
---|
1489 | } |
---|
1490 | |
---|
1491 | void |
---|
1492 | madx_mpk_setcreatelayout(struct in_cmd* cmd) |
---|
1493 | { |
---|
1494 | cmd->clone_flag = 1; /* do not delete for the moment*/ |
---|
1495 | cmd->label = (char*)mymalloc("madx_mpk_setcreatelayout",24*sizeof(char)); |
---|
1496 | strcpy(cmd->label,"matchptcknob_ptc_CL"); |
---|
1497 | madx_mpk_comm_createlayout = cmd; |
---|
1498 | } |
---|
1499 | |
---|
1500 | void |
---|
1501 | madx_mpk_setsetswitch(struct in_cmd* cmd) |
---|
1502 | { |
---|
1503 | cmd->clone_flag = 1; /* do not delete for the moment*/ |
---|
1504 | cmd->label = (char*)mymalloc("madx_mpk_setsetswitch",24*sizeof(char)); |
---|
1505 | strcpy(cmd->label,"matchptcknob_ptc_SSW"); |
---|
1506 | madx_mpk_comm_setswitch = cmd; |
---|
1507 | } |
---|
1508 | |
---|
1509 | void |
---|
1510 | madx_mpk_setcalc(struct in_cmd* cmd) |
---|
1511 | { |
---|
1512 | cmd->clone_flag = 1; /* do not delete for the moment*/ |
---|
1513 | cmd->label = (char*)mymalloc("madx_mpk_setcalc",24*sizeof(char)); |
---|
1514 | strcpy(cmd->label,"matchptcknob_ptc_CMD"); |
---|
1515 | madx_mpk_comm_calculate = cmd; |
---|
1516 | } |
---|
1517 | |
---|