home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mesaaiok.zip / source / solver.cpp < prev    next >
C/C++ Source or Header  |  1995-12-01  |  8KB  |  299 lines
  1.  
  2. #define INCL_DOSPROCESS
  3. #define INCL_DOSERRORS
  4. #include <os2.h>
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <math.h>
  8. #include "range.h"
  9. #include "model.h"
  10. #include "solver.h"
  11.  
  12. double waitForNewValue(Solver &sv, Model &model, Address &ad, double oldval, int &maxtime)
  13. {
  14.     int totaltime = 1;
  15.     double newValue;
  16.     model.getCellValue(ad,newValue);
  17.  
  18.     while (( totaltime < maxtime)  && (newValue == oldval)) {
  19.          DosSleep(50);
  20.          totaltime += 50;
  21.          model.getCellValue(ad,newValue);
  22.          if (sv.shouldStop()) {
  23.              DosExit(0,0);  //kill the thread
  24.          } // endif
  25.     }
  26.  
  27.     if (totaltime < maxtime) {
  28.        maxtime = 0;
  29.     }
  30.  
  31.     return newValue;
  32. }
  33.  
  34. double func(Solver & sv,Model &model, Address & x, Address & y, double newX)
  35. {
  36.     double oldX;
  37.     double oldY;
  38.     double newY;
  39.     int maxtime = 30000;  // 30 second wait time
  40.  
  41.     model.getCellValue(x,oldX);
  42.     model.getCellValue(y,oldY);
  43.  
  44.     if (oldX == newX) return oldY;   // the new X value is the same as the old one
  45.  
  46.     model.setCellValue(x,newX);
  47.     newY = waitForNewValue(sv,model,y,oldY,maxtime);
  48.     return newY;
  49. }
  50.  
  51. int bracketRoot(Solver & sv, Model &model, Address & x, Address & y, double &lft, double &rght, double offset)
  52. {
  53.     int j;
  54.     double f1,f2;
  55.     double gap;
  56.  
  57.     if (lft > rght) {
  58.        f1 = lft;
  59.        lft = rght;
  60.        rght = f1;
  61.     }
  62.     if (lft == rght) rght += 0.1;
  63.     gap = rght - lft;
  64.     f1 = func(sv,model,x,y,lft) - offset;
  65.     f2 = func(sv,model,x,y,rght) - offset;
  66.     for (j = 0 ; j < 100 ; j++) {
  67.          if (f1 * f2 < 0.0) return 1;
  68.          if (fabs(f1) < fabs(f2)) {
  69.              gap *= 2;
  70.              rght = lft;
  71.              lft -= gap;
  72.              f2 = f1;
  73.              f1 = func(sv,model,x,y,lft) - offset;
  74.          } else {
  75.              gap *=2;
  76.              lft = rght;
  77.              rght += gap;
  78.              f1 = f2;
  79.              f2 = func(sv,model,x,y,rght) - offset;
  80.          }
  81.     }
  82.     return 0;
  83. }
  84.  
  85. int bisectToRoot(Solver & sv,Model &model,Address & x, Address & y, double &lft, double &rght, double offset, double acc)
  86. {
  87.    int j;
  88.    double dx,f,fmid,xmid,rtb;
  89.  
  90.     f = func(sv,model,x,y,lft) - offset;
  91.     dx = rght - lft;
  92.     xmid = (rght+lft)/2;   // find midpoint;
  93.     for ( j = 0 ; j < 40 ; j ++ ) {
  94.          fmid=func(sv,model,x,y,xmid) - offset;
  95.          if (((fmid < 0.0)  && (f < 0.0)) || ((fmid > 0.0)  && (f > 0.0))){
  96.             lft = xmid;
  97.             f = fmid;
  98.          } else {
  99.             rght = xmid;
  100.          }
  101.          dx = rght - lft;
  102.          if (fabs(dx) < acc || fmid == 0) {
  103.              lft=rght=xmid;
  104.              return 1;
  105.          }
  106.          xmid = (rght+lft)/2;   // find midpoint;
  107.     }
  108.     return 0;
  109. }
  110.  
  111. void _System solveThread(unsigned long p)
  112. {
  113.     Solver * sp = (Solver *)p;
  114.     double left,right;
  115.  
  116.     sp->model().getCellValue(sp->changeCell(),left);
  117.     right = left;
  118.  
  119.     if (!bracketRoot(*sp,sp->model(),sp->changeCell(), sp->finalCell(), left, right, sp->finalValue())) {
  120.        sp->setError("Could not bracket value!  Try a better guess.");
  121.        return;
  122.     }
  123.     if (!bisectToRoot(*sp,sp->model(),sp->changeCell(),sp->finalCell(),
  124.         left,right,sp->finalValue(),sp->precision())) {
  125.         sp->setError("Could not find value! Try a better guess or lowering the precision.");
  126.         return;
  127.     }
  128.     sp->setPercentComplete(100);
  129. }
  130.  
  131.  
  132.  
  133. int bracketMinMax(Solver &sv,Model &model, Address & x, Address & y,
  134.                double &lft, double &fl,double &rght, double &fr, int min)
  135. {
  136.     int j;
  137.     double fmid, mid;
  138.     double gap;
  139.  
  140.     if (lft > rght) {
  141.        fl = lft;
  142.        lft = rght;
  143.        rght = fl;
  144.     }
  145.     if (lft == rght) rght += 0.1;
  146.     fl = func(sv,model,x,y,lft);
  147.     fr = func(sv,model,x,y,rght);
  148.     mid = (rght + lft) / 2;
  149.     gap = (rght - lft) / 2;
  150.     fmid = func(sv,model,x,y,mid);
  151.  
  152.     for (j = 0 ; j < 100 ; j++) {
  153.          if (( fl > fmid ) && ( fr > fmid ) && min) return 1;
  154.          if (( fl < fmid ) && ( fr < fmid ) && !min) return 1;
  155.          if (((fl < fr) && min) || ((fl > fr) && !min)) {
  156.              gap *= 2;
  157.              rght = mid;
  158.              mid = lft;
  159.              lft -= gap;
  160.              fr = fmid;
  161.              fmid = fl;
  162.              fl = func(sv,model,x,y,lft);
  163.          } else {
  164.              gap *=2;
  165.              lft = mid;
  166.              mid = rght;
  167.              rght += gap;
  168.              fl = fmid;
  169.              fmid = fr;
  170.              fr = func(sv,model,x,y,rght);
  171.          }
  172.     }
  173.     return 0;
  174. }
  175.  
  176.  
  177.  
  178. int bisectToMinMax(Solver & sv,Model &model,Address & x, Address & y,
  179.                 double &lft,double & fl,double &rght,double &fr, double acc, int min)
  180. {
  181.     double gap = (fabs(rght - lft)) / 7;
  182.     double dx = lft;
  183.     double fx;
  184.     double oldfx, oldoldfx;
  185.  
  186.     oldfx = fx = oldoldfx = fl;
  187.  
  188.     for (int j = 0; j < 40; j++ ) {
  189.         while ((dx < rght) && (((oldfx >= fx) && min) ||((oldfx <= fx) && !min)))  {
  190.             oldoldfx = oldfx;
  191.             oldfx = fx;
  192.             dx += gap;
  193.             fx = func(sv,model,x,y,dx);
  194.         } /* endwhile */
  195.  
  196.         rght = dx;
  197.         fr = fx;
  198.  
  199.         dx -= 2*gap;
  200.  
  201.         if (dx < lft) {
  202.            dx = lft;
  203.         } /* endif */
  204.  
  205.         lft = dx;
  206.  
  207.         if (fabs(lft-rght) < acc) {
  208.            lft = rght;
  209.             return 1;
  210.         } /* endif */
  211.  
  212.         fl = oldfx = fx = oldoldfx;
  213.         gap = (fabs(rght - lft)) / 7;
  214.     } /* endfor */
  215.     return 0;
  216. }
  217.  
  218.  
  219. void _System minimizeThread(unsigned long p)
  220. {
  221.     Solver * sp = (Solver *)p;
  222.     double left,right, fl,fr;
  223.  
  224.     sp->model().getCellValue(sp->changeCell(),left);
  225.     right = left;
  226.  
  227.     if (!bracketMinMax(*sp,sp->model(),sp->changeCell(), sp->finalCell(), left, fl,right,fr, 1)) {
  228.        sp->setError("Could not bracket minimum!  Try a better guess.");
  229.        return;
  230.     }
  231.     if (!bisectToMinMax(*sp,sp->model(),sp->changeCell(),sp->finalCell(),
  232.         left,fl,right,fr,sp->precision(),1)) {
  233.         sp->setError("Could not find minimum! Try a better guess or lowering the precision.");
  234.         return;
  235.     }
  236.     sp->setPercentComplete(100);
  237.  
  238. }
  239.  
  240. void _System maximizeThread(unsigned long p)
  241. {
  242.     Solver * sp = (Solver *)p;
  243.     double left,right, fl,fr;
  244.  
  245.     sp->model().getCellValue(sp->changeCell(),left);
  246.     right = left;
  247.  
  248.     if (!bracketMinMax(*sp,sp->model(),sp->changeCell(), sp->finalCell(), left, fl ,right, fr, 0)) {
  249.        sp->setError("Could not bracket maximum!  Try a better guess.");
  250.        return;
  251.     }
  252.     if (!bisectToMinMax(*sp,sp->model(),sp->changeCell(),sp->finalCell(),
  253.         left,fl,right,fr,sp->precision(),0)) {
  254.         sp->setError("Could not find maximum! Try a better guess or lowering the precision.");
  255.         return;
  256.     }
  257.     sp->setPercentComplete(100);
  258.  
  259. }
  260.  
  261.  
  262. Solver :: Solver(Model & theModel, Address & change, Address &final,
  263.                double precision, int type, double value) :
  264.     _model(theModel),chng(change),finl(final),prec(precision),
  265.     tp(type),fval(value)
  266. {
  267.     error[0] = 0;
  268.     percentComp = 0;
  269.     dostop = 0;
  270.     switch (tp) {
  271.     case 0:
  272.        DosCreateThread(&tid,&solveThread,(ULONG)this,0L,16384);
  273.        break;
  274.     case 1:
  275.        DosCreateThread(&tid,&minimizeThread,(ULONG)this,0L,16384);
  276.        break;
  277.     case 2:
  278.        DosCreateThread(&tid,&maximizeThread,(ULONG)this,0L,16384);
  279.        break;
  280.     } /* endswitch */
  281. }
  282.  
  283.     // stop the solving process, blocks until the solving thread is done
  284. Solver & Solver :: stop()
  285. {
  286. //    int rc = DosKillThread(tid);
  287.       // don't use DosKillThread.  it can leave MesaCore in an unknown state
  288.       // and possibly cause crashes
  289.     dostop = 1;
  290.     setError("Stopped before completion!");
  291.     return *this;
  292. }
  293.  
  294. Solver :: ~Solver()
  295. {
  296.     dostop = 1;
  297.     DosWaitThread(&tid,DCWW_WAIT);
  298. }
  299.