473,322 Members | 1,620 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,322 software developers and data experts.

a for loop that quits prematurely

I have a program which uses a for loop to iteratively fit a non-linear model to some data. The problem is, the loop always finishes after just one iteration (even though I might specify 1000) before the mimisation criteria have been met. The for loop is included below;
--------------------------------------------------------------------------------------------------------------------
Expand|Select|Wrap|Line Numbers
  1.  double a, b, c, d; 
  2. int n_max, n; //n_max = number of iterations & n = running index
  3.  
  4. cout << endl << "enter guess for a ... " << endl;
  5. cin >> a; 
  6.  
  7. cout << endl << "enter guess for b ... " << endl;
  8. cin >> b; 
  9.  
  10. cout << endl << "enter guess for c ... " << endl;
  11. cin >> c; 
  12.  
  13. cout << endl << "enter maximum number of iterations ... " << endl;
  14. cin >> n_max;
  15.  
  16. long double K_old[3][1], K_new[3][1], K_min[3][1];
  17. K_old [1][1] = a;
  18. K_old [2][1] = b;
  19. K_old [3][1] = c;
  20. double chi_min;
  21. chi_min = 2000;
  22. long double chi[n_max];
  23. for(n = 0; n <= n_max; n++)
  24. {
  25. long double Chi, Chi_prime[1][3], Chi_prime_prime[3][3][n_max];
  26. double X, i[9096];
  27. int x;
  28. Chi = 0;
  29. Chi_prime_prime[1][1][n] = 0; 
  30. Chi_prime_prime[1][2][n] = 0;
  31. Chi_prime_prime[1][3][n] = 0;
  32. Chi_prime_prime[2][1][n] = 0;
  33. Chi_prime_prime[2][2][n] = 0;
  34. Chi_prime_prime[2][3][n] = 0;
  35. Chi_prime_prime[3][1][n] = 0;
  36. Chi_prime_prime[3][2][n] = 0;
  37. Chi_prime_prime[3][3][n] = 0;
  38. X = 2 * S * S * -1; 
  39.  
  40. for (j = 0; j < m; j++) //calculating grad's and concavity
  41. {
  42. i[j] = Imax * exp((lambda[j] - average)/X) - (a * a * lambda[j] + b * lambda[j] + c); 
  43. Chi = Chi + pow(((I[j] - i[j])/S), 2);
  44. Chi_prime[1][1] = Chi_prime[1][1] + 2*((I[j] - i[j])/S)*((2 * a * lambda[j])/S); //derivative w.r.t. a
  45. Chi_prime[2][1] = Chi_prime[2][1] + 2*((I[j] - i[j])/S)*(lambda[j]/S); //derivative b
  46. Chi_prime[3][1] = Chi_prime[3][1] + 2*((I[j] - i[j])/S)*(1/S); //derivative c
  47. Chi_prime_prime[1][1][n] = Chi_prime_prime[1][1][n] + ((I[j] - i[j])/S)*((4 * lambda[j])/S) + ((4 * a * lambda[j])/S)*((2 * a * lambda[j])/S);
  48. Chi_prime_prime[1][2][n] = Chi_prime_prime[1][2][n] + (((4 * a * lambda[j])/S)*(lambda[j]/S));
  49. Chi_prime_prime[1][3][n] = Chi_prime_prime[1][3][n] + ((4 * a * lambda[j])/S)*(1/S);
  50. Chi_prime_prime[2][1][n] = Chi_prime_prime[2][1][n] + ((2 * lambda[j])/S)*((2 * a * lambda[j])/S);
  51. Chi_prime_prime[2][2][n] = Chi_prime_prime[2][2][n] + (((2 * lambda[j])/S)*(lambda[j]/S));
  52. Chi_prime_prime[2][3][n] = Chi_prime_prime[2][3][n] + (((2 * lambda[j])/S)*(1/S));
  53. Chi_prime_prime[3][1][n] = Chi_prime_prime[3][1][n] + (2/S)*((2 * a * lambda[j])/S);
  54. Chi_prime_prime[3][2][n] = Chi_prime_prime[3][2][n] + ((2/S)*(lambda[j]/S));
  55. Chi_prime_prime[3][3][n] = Chi_prime_prime[3][3][n] + ((2/S)*(1/S));
  56. }
  57. chi[n] = Chi;
  58. cout << "chi squared = " << chi[n] << " ";
  59. prt4 << "chi squared = " << chi[n] << " ";
  60. jacobian << "n = " << n << endl;
  61. jacobian << Chi_prime_prime[1][1][n] << " ";
  62. jacobian << Chi_prime_prime[1][2][n] << " ";
  63. jacobian << Chi_prime_prime[1][3][n] << endl;
  64. jacobian << Chi_prime_prime[2][1][n] << " ";
  65. jacobian << Chi_prime_prime[2][2][n] << " ";
  66. jacobian << Chi_prime_prime[2][3][n] << endl;
  67. jacobian << Chi_prime_prime[3][1][n] << " ";
  68. jacobian << Chi_prime_prime[3][2][n] << " ";
  69. jacobian << Chi_prime_prime[3][3][n] << endl << endl;
  70.  
  71. long double det_J = Chi_prime_prime[1][1][n]*(Chi_prime_prime[2][2][n]*Chi_prime_prime[3][3][n] 
  72. - (Chi_prime_prime[2][3][n]*Chi_prime_prime[3][2][n])) - Chi_prime_prime[1][2][n]
  73. *(Chi_prime_prime[2][1][n]*Chi_prime_prime[3][3][n] - (Chi_prime_prime[2][3][n]*
  74. Chi_prime_prime[3][1][n])) + Chi_prime_prime[1][3][n]*(Chi_prime_prime[2][1][n]*
  75. Chi_prime_prime[3][2][n] - (Chi_prime_prime[2][2][n]*Chi_prime_prime[3][1][n]));
  76.  
  77. long double inv_J[3][3][n];
  78. inv_J[1][1][n] = (1/det_J)*Chi_prime_prime[1][1][n];
  79. inv_J[1][2][n] = (1/det_J)*Chi_prime_prime[1][2][n];
  80. inv_J[1][3][n] = (1/det_J)*Chi_prime_prime[1][3][n];
  81. inv_J[2][1][n] = (1/det_J)*Chi_prime_prime[2][1][n];
  82. inv_J[2][2][n] = (1/det_J)*Chi_prime_prime[2][2][n];
  83. inv_J[2][3][n] = (1/det_J)*Chi_prime_prime[2][3][n]; 
  84. inv_J[3][1][n] = (1/det_J)*Chi_prime_prime[3][1][n];
  85. inv_J[3][2][n] = (1/det_J)*Chi_prime_prime[3][2][n];
  86. inv_J[3][3][n] = (1/det_J)*Chi_prime_prime[3][3][n];
  87. inverse << "n = " << n << endl;
  88. inverse << inv_J[1][1][n] << " ";
  89. inverse << inv_J[1][2][n] << " ";
  90. inverse << inv_J[1][3][n] << endl;
  91. inverse << inv_J[2][1][n] << " ";
  92. inverse << inv_J[2][2][n] << " ";
  93. inverse << inv_J[2][3][n] << endl;
  94. inverse << inv_J[3][1][n] << " ";
  95. inverse << inv_J[3][2][n] << " ";
  96. inverse << inv_J[3][3][n] << endl << endl;
  97.  
  98. long double D[3][1];
  99. for (0 < x; x <= 3; x++)
  100. D[x][1] = 0;
  101. for (0 < x; x <= 3; x++)
  102. D[1][1] = D[1][1] + inv_J[1][x][n]*Chi_prime[1][1];
  103. D[2][1] = D[2][1] + inv_J[2][x][n]*Chi_prime[2][1]; 
  104. D[3][1] = D[3][1] + inv_J[3][x][n]*Chi_prime[3][1];
  105. prt4 << "D[1][1] = " << D[1][1] << ", D[2][1] = " << D[2][1];
  106. prt4 << ", D[3][1] = " << D[3][1] << endl << endl;
  107. cout << "D[1][1] = " << D[1][1] << ", D[2][1] = " << D[2][1];
  108. cout << ", D[3][1] = " << D[3][1] << endl;
  109.  
  110. K_new[1][1] = K_old[1][1] - D[1][1];
  111. K_new[2][1] = K_old[2][1] - D[2][1];
  112. K_new[3][1] = K_old[3][1] - D[3][1];
  113. prt4 << "K[1][1] = " << K_new[1][1] << ", K[2][1] = " << K_new[2][1];
  114. prt4 << ", K[3][1] = " << K_new[3][1] << endl << endl;
  115. cout << "K[1][1] = " << K_new[1][1] << ", K[2][1] = " << K_new[2][1];
  116. cout << ", K[3][1] = " << K_new[3][1] << endl;
  117.  
  118. if (n > 0)
  119. if (Chi_prime_prime[1][1][n] > 0 && Chi_prime_prime[1][1][n - 1] > 0)//concavity stays const.
  120. if (Chi_prime_prime[2][2][n] > 0 && Chi_prime_prime[2][2][n - 1] > 0)
  121. if (Chi_prime_prime[3][3][n] > 0 && Chi_prime_prime[3][3][n - 1] > 0)
  122. {
  123. if (chi_min - chi [n] > 0)
  124. {
  125. chi_min = chi[n];
  126. K_min[1][1] = K_new[1][1];
  127. K_min[2][1] = K_new[2][1];
  128. K_min[3][1] = K_new[3][1];
  129. }
  130. else {
  131. chi_min = chi_min;
  132. K_min[1][1] = K_min[1][1];
  133. K_min[2][1] = K_min[2][1];
  134. K_min[3][1] = K_min[3][1];
  135. }
  136. }
  137. else {
  138. prt4 << "point of inflection passed or in vicinity of local max" << endl << endl;
  139. cout << "minimised chi = " << chi_min << " ";
  140. cout << "corresponding K = " << K_min[1][1];
  141. cout << ", " << K_min[2][1];
  142. cout << ", " << K_min[3][1] << endl << endl;
  143. prt << "minimised chi = " << chi_min;
  144. prt << " corresponding K = " << K_min[1][1];
  145. prt << ", " << K_min[2][1];
  146. prt << ", " << K_min[3][1] << endl << endl;
  147. for (j = 0; j < m; j++)
  148. {
  149. //cout << "i[" << j + 1 << "] = " << i[j] << endl;
  150. //intensity << "i[" << j + 1 << "] = " << i[j] << endl;
  151. intensity << i[j] << endl;
  152. }
  153. system ("pause");
  154. exit(1);
  155. }
  156.  
  157. if (n > 2) //mimimisation criteria
  158. if ((chi[n] = chi_min) && (fabs(chi[n - 1] - chi_min)) < 0.1 && (fabs(chi[n - 2] - chi_min)) < 0.1)
  159. {
  160. cout << "minimised chi = " << chi_min << " " << "corresponding K = ";
  161. cout << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl << endl; 
  162. prt << "minimised chi = " << chi_min << " " << "corresponding K = ";
  163. prt << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl << endl;
  164. for (j = 0; j < m; j++)
  165. {
  166. //cout << "i[" << j + 1 << "] = " << i[j] << endl;
  167. //intensity << "i[" << j + 1 << "] = " << i[j] << endl;
  168. intensity << i[j] << endl;
  169. }
  170. system ("pause");
  171. exit(1);
  172. }
  173.  
  174. K_old[1][1] = K_new[1][1];
  175. K_old[2][1] = K_new[2][1]; 
  176. K_old[3][1] = K_new[3][1];
  177. }
  178. cout << "minimised chi = " << chi_min << " " << "corresponding K = ";
  179. cout << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl; 
  180. prt << "minimised chi = " << chi_min << " " << "corresponding K = ";
  181. prt << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl;
  182. prt.close();
  183. prt4.close();
  184. system ("pause");
  185. }
  186.  
Sep 4 '06 #1
3 1533
D_C
293 100+
The nested loop, which terminates when j < m, you don't ever define m. Perhaps it is 0 by default, or it could be some garbage value.
Sep 4 '06 #2
m is defined in an earlier part of the code as 2048, I couldn't include it all because that would have made the post too large.
Pixie

The nested loop, which terminates when j < m, you don't ever define m. Perhaps it is 0 by default, or it could be some garbage value.
Sep 4 '06 #3
I've got it so the loop ceases to terminate prematurely, those who care to look will find many simple errors which needed to be corrected. The problem now is that the value 'det_J' (the determinate of the jacobian matrix) isn't always well defined. Here are some of the results I get;
--------------------------------------------------------------------------------------------------------------------
n = 0, det = -5.37054e+007
-8.77019e-011 -8.77019e-011 -8.77019e-011
-8.77019e-011 -8.77019e-011 -8.77019e-011
-8.77019e-011 -8.77019e-011 -8.77019e-011

n = 2, det = -5.38309e+007
0.00106051 -0.000932995 -1.26731e-007
-0.000932995 -0.00466498 -6.33653e-007
-1.26731e-007 -2.07463e-007 -8.74974e-011

n = 3, det = 1.#INF
0 0 0
0 0 0
0 4.36697e-007 0

n = 4, det = -5.39565e+007
0.00105804 -0.000930824 -1.26436e-007
-0.000930824 -0.00465412 -6.32178e-007
-1.#IND -2.05964e-007 -8.72938e-011

n = 5, det = 1.#INF
0 0 0
0 0 0
0 4.36697e-007 0

n = 6, det = -5.40821e+007
0.00105559 -0.000928662 -1.26142e-007
-0.000928662 -0.00464331 -6.30709e-007
-1.26142e-007 -2.04471e-007 -8.7091e-011

n = 7, det = -1.#IND
-1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND

n = 8, det = -5.42079e+007
0.00105314 -0.000926508 -1.25849e-007
-0.000926508 -0.00463254 -6.29247e-007
-1.25849e-007 -2.02984e-007 -8.6889e-011

n = 9, det = 1.#INF
0 0 0
0 0 0
-1.#IND 4.36697e-007 0
--------------------------------------------------------------------------------------------------------------------
Pixie
Sep 5 '06 #4

Sign in to post your reply or Sign up for a free account.

Similar topics

1
by: Bob Helber | last post by:
I've got a perl script with a system command, within a for loop, that calls a FORTRAN program. For one of the iterations of the for loop the FORTRAN program quits with an error due to bad input. ...
3
by: zeroDoNotYeSpamtype | last post by:
(I tried to post this earlier, but it seems without success) A similar question to this has been answered many times, to whit the question applying when the index variable of a for loop is...
1
by: Anthony Knittel | last post by:
i've got an expensive loop that keeps running for a long time, whats the best way to put breaks in the program to return a bit of control to the user instead of locking up the system? its just a...
5
by: Blankdraw | last post by:
I can't get this nested loop to break the outer loop at the 5th data value so control can proceed to the next array col and continue pigeon-holing the next 5 in its own column. Why can I not get...
3
by: picknicker187 | last post by:
hi, the following while function(it´s supposed to find word synonyms out of a list) always quits after running through the list in the inner loop once. it´s like the aktuell = aktuell->next on...
1
by: Alvin A. Delagon | last post by:
Greetings, I have to write a python script that would continously monitor and process a queue database. Once the script sees an unprocessed record it will create a thread to process it otherwise...
4
by: Sebastian Loncar | last post by:
Hallo NG, i've big online game with a lots of requests per second and it uses a lot of memory. the game runs as an ASP.NET application and holts all the game data in the memory. The game has 1600...
13
by: Sarath | last post by:
What's the advantage of using for_each than the normal iteration using for loop? Is there any 'tweak' to use stream objects as 'Function' parameter of for_each loop. Seems copy function can do...
7
by: gordon | last post by:
is it possible to send a message to the gui instance while the Tk event loop is running?I mean after i create a gui object like root=Tk() mygui=SomeUI(root) and call root.mainloop() can i...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
0
by: ArrayDB | last post by:
The error message I've encountered is; ERROR:root:Error generating model response: exception: access violation writing 0x0000000000005140, which seems to be indicative of an access violation...
1
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
1
by: Shællîpôpï 09 | last post by:
If u are using a keypad phone, how do u turn on JavaScript, to access features like WhatsApp, Facebook, Instagram....
0
by: af34tf | last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...
0
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 3 Apr 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome former...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.