Details

    • Type: New Feature New Feature
    • Status: Open
    • Priority: Major Major
    • Resolution: Unresolved
    • Affects Version/s: 3.0
    • Fix Version/s: 4.0
    • Labels:
      None

      Description

      During experiments with space flight trajectory optimizations I recently
      observed, that the direct optimization algorithm BOBYQA
      http://plato.asu.edu/ftp/other_software/bobyqa.zip
      from Mike Powell is significantly better than the simple Powell algorithm
      already in commons.math. It uses significantly lower function calls and is
      more reliable for high dimensional problems. You can replace CMA-ES in many
      more application cases by BOBYQA than by the simple Powell optimizer.
      I would like to contribute a Java port of the algorithm.
      I maintained the structure of the original FORTRAN code, so the
      code is fast but not very nice.

      License status: Michael Powell has sent the agreement via snail mail

      • it hasn't arrived yet.

      Progress: The attached patch relative to the trunk contains both the
      optimizer and the related unit tests - which are all green now.

      Performance:
      Performance difference (number of function evaluations)
      PowellOptimizer / BOBYQA for different test functions (taken from
      the unit test of BOBYQA, dimension=13 for most of the
      tests.

      Rosen = 9350 / 1283
      MinusElli = 118 / 59
      Elli = 223 / 58
      ElliRotated = 8626 / 1379
      Cigar = 353 / 60
      TwoAxes = 223 / 66
      CigTab = 362 / 60
      Sphere = 223 / 58
      Tablet = 223 / 58
      DiffPow = 421 / 928
      SsDiffPow = 614 / 219
      Ackley = 757 / 97
      Rastrigin = 340 / 64

      The number for DiffPow should be discussed with Michael Powell,
      I will send him the details.

      Open Problems

      • Checkstyle violations[1] because of the original Fortran source:
        • Original method comments were copied: Javadoc standard documentation should be added, but the original documentation should stay (as a reference to what the original intended behaviour was) untouched until we are sure that the code behaves as expected.
        • Multiple variable declarations per line.
        • "goto" conversions:
          • "goto"s not convertible in loops were translated into a finite automaton (switch statement)
          • "no default in switch"
          • "fall through from previous case in switch"
      • Unexplored code paths: "throw" statements have been introduced in the code. Each should be triggered by at least one unit test. They are currently commented out in provision of the 3.0 release (cf. MATH-712) but should be re-enabled afterwards.

      [1] Once the violations are solved, the following lines should be removed from the source file:

      • // CHECKSTYLE: stop all
      • //CHECKSTYLE: resume all
      1. bobyqa_convert.pl
        0.8 kB
        Gilles
      2. BOBYQA.math.patch
        158 kB
        Dr. Dietmar Wolz
      3. BOBYQA.v02.math.patch
        158 kB
        Dr. Dietmar Wolz
      4. bobyqa.zip
        39 kB
        Dr. Dietmar Wolz
      5. BOBYQAOptimizer.java.patch
        3 kB
        Dr. Dietmar Wolz
      6. bobyqaoptimizer0.4.zip
        129 kB
        Dr. Dietmar Wolz
      7. BOBYQAOptimizer0.4.zip
        26 kB
        Dr. Dietmar Wolz
      8. bobyqav0.3.zip
        39 kB
        Dr. Dietmar Wolz

        Issue Links

          Activity

          Hide
          Dr. Dietmar Wolz added a comment -

          Initial code contribution (BOBYQAOptimizer.java + Unit tests)

          Show
          Dr. Dietmar Wolz added a comment - Initial code contribution (BOBYQAOptimizer.java + Unit tests)
          Hide
          Gilles added a comment -

          Oops, this is indeed FORTRAN code in Java clothes. There is quite a lot of work to make it look like Java...
          Personally, I'd rather not commit it as is, because it is unmaintainable (it does not match anything in CM in terms of design and programming style).

          Show
          Gilles added a comment - Oops, this is indeed FORTRAN code in Java clothes. There is quite a lot of work to make it look like Java... Personally, I'd rather not commit it as is, because it is unmaintainable (it does not match anything in CM in terms of design and programming style).
          Hide
          Dr. Dietmar Wolz added a comment -

          Tests with the original code confirmed a problem of BOBYQA with
          DiffPow and SsDiffPow. But the tests revealed a problem with the
          Java code: Initialization of the bounds with Double.MIN_VALUE and
          Double.MAX_VALUE (representing no bounds)
          caused trouble, I changed the values to -1e300, 1e300.
          The old performance numbers a not valid, BOBYQA
          changed the initial guess because of the bounds problem.

          Adjusted accuracy level for PowellOptimizer and the new
          BOBYQA code - now initialized correctly, lead to the
          following performance difference (number of function evaluations)
          PowellOptimizer / BOBYQA for different test functions (taken from
          the unit test of BOBYQA, dimension=13 for most of the
          tests.

          Rosen = 9831 / 1404
          MinusElli = 183 / 115
          Elli = 399 / 147
          ElliRotated = 8964 / 7165
          Cigar = 348 / 56
          TwoAxes = 223 / 70
          CigTab = 513 / 61
          Sphere = 223 / 56
          Tablet = 223 / 57
          DiffPow = 308 / 10139
          SsDiffPow = 667 / 4520
          Ackley = 814 / 458
          Rastrigin = 327 / 167

          Essentially the same picture as before.

          Show
          Dr. Dietmar Wolz added a comment - Tests with the original code confirmed a problem of BOBYQA with DiffPow and SsDiffPow. But the tests revealed a problem with the Java code: Initialization of the bounds with Double.MIN_VALUE and Double.MAX_VALUE (representing no bounds) caused trouble, I changed the values to -1e300, 1e300. The old performance numbers a not valid, BOBYQA changed the initial guess because of the bounds problem. Adjusted accuracy level for PowellOptimizer and the new BOBYQA code - now initialized correctly, lead to the following performance difference (number of function evaluations) PowellOptimizer / BOBYQA for different test functions (taken from the unit test of BOBYQA, dimension=13 for most of the tests. Rosen = 9831 / 1404 MinusElli = 183 / 115 Elli = 399 / 147 ElliRotated = 8964 / 7165 Cigar = 348 / 56 TwoAxes = 223 / 70 CigTab = 513 / 61 Sphere = 223 / 56 Tablet = 223 / 57 DiffPow = 308 / 10139 SsDiffPow = 667 / 4520 Ackley = 814 / 458 Rastrigin = 327 / 167 Essentially the same picture as before.
          Hide
          Dr. Dietmar Wolz added a comment -

          "FORTRAN code in Java clothes" was the best I could do so far.
          Perhaps someone else voluteers to transfer the code into a more
          readable version using the commons.math matrix operations?
          The problem is that the original code heavily relies on pointers and gotos.
          It is indeed probably only maintainable by the author.
          The current code could simplify the development of a "clean" version
          by enabling parallel debugging using only a Java development environment.
          I started with embedding a C++ port via JNI treating the optimizer
          as a "black box". Since JNI is not possible in commons.math, I tried
          to convert the code to Java.
          Unfortunately there is nothing publicly available which
          has comparable qualities in terms of usability.
          From the developer perspective we have a maintanance problem, but from
          the user perspective, it is not very satisfactory to abandon BOBYQA.

          Show
          Dr. Dietmar Wolz added a comment - "FORTRAN code in Java clothes" was the best I could do so far. Perhaps someone else voluteers to transfer the code into a more readable version using the commons.math matrix operations? The problem is that the original code heavily relies on pointers and gotos. It is indeed probably only maintainable by the author. The current code could simplify the development of a "clean" version by enabling parallel debugging using only a Java development environment. I started with embedding a C++ port via JNI treating the optimizer as a "black box". Since JNI is not possible in commons.math, I tried to convert the code to Java. Unfortunately there is nothing publicly available which has comparable qualities in terms of usability. From the developer perspective we have a maintanance problem, but from the user perspective, it is not very satisfactory to abandon BOBYQA.
          Hide
          Gilles added a comment -

          No problem, as I said on the ML. Thanks for the contribution.

          [...] the tests revealed a problem with the Java code [...]

          This is another reason why I don't like a diff/patch for the initial version. I've already made small changes so that it doesn't apply. Could you please attach the Java file itself (so that I don't have to work out the diff between diffs...)? Thanks.

          Show
          Gilles added a comment - No problem, as I said on the ML. Thanks for the contribution. [...] the tests revealed a problem with the Java code [...] This is another reason why I don't like a diff/patch for the initial version. I've already made small changes so that it doesn't apply. Could you please attach the Java file itself (so that I don't have to work out the diff between diffs...)? Thanks.
          Hide
          Dr. Dietmar Wolz added a comment -

          As requested the source files of the second patch.

          Show
          Dr. Dietmar Wolz added a comment - As requested the source files of the second patch.
          Hide
          Gilles added a comment -

          It seems that none of the tests exercise the "rescue" method.

          Show
          Gilles added a comment - It seems that none of the tests exercise the "rescue" method.
          Hide
          Dr. Dietmar Wolz added a comment -

          Confirmed - this also happens with the original code, I will ask Mike Powell to provide us with a function exercising rescue. You are probably suspicious whether the Java rescue works, as am I.

          Show
          Dr. Dietmar Wolz added a comment - Confirmed - this also happens with the original code, I will ask Mike Powell to provide us with a function exercising rescue. You are probably suspicious whether the Java rescue works, as am I.
          Hide
          Gilles added a comment -

          In fact I'm more worried about introducing bugs that will go unnoticed.

          Show
          Gilles added a comment - In fact I'm more worried about introducing bugs that will go unnoticed.
          Hide
          Dr. Dietmar Wolz added a comment -

          Now I remember that I debugged rescue (and removed a nasty infinite loop) testing MinusElli with GoalType.MAXIMIZE before I fixed the
          GoalType.MAXIMIZE handling. Testing MinusElli with GoalType.MINIMIZE still calls rescue:

          @Test
          public void testRescue()

          { double[] startPoint = point(DIM,1.0); double[][] boundaries = null; RealPointValuePair expected = new RealPointValuePair(point(DIM,0.0),-7.5072566333E53); doTest(new MinusElli(), startPoint, boundaries, GoalType.MINIMIZE, 1e45, 1e24, 1000, expected); }

          But parallel testing with the original revealed a problem, so please wait a few hours before working on rescue, I will fix it soon.

          Show
          Dr. Dietmar Wolz added a comment - Now I remember that I debugged rescue (and removed a nasty infinite loop) testing MinusElli with GoalType.MAXIMIZE before I fixed the GoalType.MAXIMIZE handling. Testing MinusElli with GoalType.MINIMIZE still calls rescue: @Test public void testRescue() { double[] startPoint = point(DIM,1.0); double[][] boundaries = null; RealPointValuePair expected = new RealPointValuePair(point(DIM,0.0),-7.5072566333E53); doTest(new MinusElli(), startPoint, boundaries, GoalType.MINIMIZE, 1e45, 1e24, 1000, expected); } But parallel testing with the original revealed a problem, so please wait a few hours before working on rescue, I will fix it soon.
          Hide
          Gilles added a comment -

          There's no hurry. I'm not going to stop sleeping until the conversion is done; there is so much to do...
          I've also started reading the reference paper in the hope to figure the algorithm bits out of the Fortran gigantic functions.

          Show
          Gilles added a comment - There's no hurry. I'm not going to stop sleeping until the conversion is done; there is so much to do... I've also started reading the reference paper in the hope to figure the algorithm bits out of the Fortran gigantic functions.
          Hide
          Dr. Dietmar Wolz added a comment -

          rescue was as I feared not correct, the fixed version v0.3 behaves exactly as the original.

          A breakup of the working area in separate double arrays should be straightforward, but still some work because of the nasty loop counting starting with 1.

          Something like
          double[] xbase = new double[n];
          double[] xpt = new double[n * npt];
          double[] fval = new double[npt];
          double[] xopt = new double[n];
          double[] gopt = new double[n];
          double[] hq = new double[n * np / 2];
          double[] pq = new double[npt];
          double[] bmat = new double[ndim * n];
          double[] zmat = new double[npt * (npt - np)];
          double[] sl = new double[n];
          double[] su = new double[n];
          double[] xnew = new double[n];
          double[] xalt = new double[n];
          double[] d__ = new double[n];
          double[] vlag = new double[ndim];
          double[] w = new double[?];

          Class ScopePtr could then be removed and the code should get even faster.

          The "gotos"/ switches are much harder to eliminate.

          Show
          Dr. Dietmar Wolz added a comment - rescue was as I feared not correct, the fixed version v0.3 behaves exactly as the original. A breakup of the working area in separate double arrays should be straightforward, but still some work because of the nasty loop counting starting with 1. Something like double[] xbase = new double [n] ; double[] xpt = new double [n * npt] ; double[] fval = new double [npt] ; double[] xopt = new double [n] ; double[] gopt = new double [n] ; double[] hq = new double [n * np / 2] ; double[] pq = new double [npt] ; double[] bmat = new double [ndim * n] ; double[] zmat = new double [npt * (npt - np)] ; double[] sl = new double [n] ; double[] su = new double [n] ; double[] xnew = new double [n] ; double[] xalt = new double [n] ; double[] d__ = new double [n] ; double[] vlag = new double [ndim] ; double[] w = new double [?] ; Class ScopePtr could then be removed and the code should get even faster. The "gotos"/ switches are much harder to eliminate.
          Hide
          Luc Maisonobe added a comment -

          The Software Grant from prof. M.J.D. Powell has been filed in Apache Software Foundation records.
          So this administrative step is over, we can include the code as soon as we want to.

          Show
          Luc Maisonobe added a comment - The Software Grant from prof. M.J.D. Powell has been filed in Apache Software Foundation records. So this administrative step is over, we can include the code as soon as we want to.
          Hide
          Dr. Dietmar Wolz added a comment -

          The attached new version contains three variants of BOBYQA implementing the
          MultivariateRealOptimizer together with their associated unit tests:

          BOBYQAOptimizerOld - essentially the previous version 0.3 with minor adaptions
          in the unit tests added for comparison.

          BOBYQAOptimizer - the actual version I use myself - all "pointers" are replaced
          by Java arrays. This version is much easier to read and faster than the old
          one. Nevertheless the "finite state machine logic" realized by the original
          code using gotos is still there - in Java realized using switch/case statements.

          BOBYQAOptimizerC - a version showing how MultivariateRealOptimizer can be implemented using the C-port of BOBYQA from the dlib library. The C+-code doing the mapping from Java to C+ to Java via JNI is also included.

          Some remarks:

          • all three versions are semantically equivalent - at least from the point of
            the actual unit test coverage.
          • the new Java version is about 30% faster than the old one - array access is
            cheaper than method calls simulating pointer arithmetic.
          • I observed a significant change of behaviour when statements like
            a = a+(b...) are replaced by a += (b...). Subtle accuracy differences add up
            during the optimization process.
          • The C version was tested both on Mac (Snow Leopard) and Windows7,
            both for 32 and 64-bit Java. For 64bit the Microsoft C++-compiler produced correct
            results only if optimization is switched off completely, on the Mac the
            gnu compiler works correctly using -O3.
          • On the Mac (a Macbook pro with a new sandy bridge processor) the C-version
            performs the unit tests three times faster than the (new) Java version.
            The difference on Win7 is much smaller, since I couldn't use compiler optimization
            there yet. The performance difference is not significant for my own application of the algo (space flight trajectory optimization) because of the high cost for the eval function. But this is perhaps not true for all other applications.

          So if we want to include BOBYQA in commons.math in my opinion we shouldn't
          neglect performance issues generally. Even a direct port of the code has a significant
          performance disadvantage, any refactoring should try to avoid adding too much
          additional performance bottlenecks.

          Show
          Dr. Dietmar Wolz added a comment - The attached new version contains three variants of BOBYQA implementing the MultivariateRealOptimizer together with their associated unit tests: BOBYQAOptimizerOld - essentially the previous version 0.3 with minor adaptions in the unit tests added for comparison. BOBYQAOptimizer - the actual version I use myself - all "pointers" are replaced by Java arrays. This version is much easier to read and faster than the old one. Nevertheless the "finite state machine logic" realized by the original code using gotos is still there - in Java realized using switch/case statements. BOBYQAOptimizerC - a version showing how MultivariateRealOptimizer can be implemented using the C-port of BOBYQA from the dlib library. The C+ -code doing the mapping from Java to C + to Java via JNI is also included. Some remarks: all three versions are semantically equivalent - at least from the point of the actual unit test coverage. the new Java version is about 30% faster than the old one - array access is cheaper than method calls simulating pointer arithmetic. I observed a significant change of behaviour when statements like a = a+(b...) are replaced by a += (b...). Subtle accuracy differences add up during the optimization process. The C version was tested both on Mac (Snow Leopard) and Windows7, both for 32 and 64-bit Java. For 64bit the Microsoft C++-compiler produced correct results only if optimization is switched off completely, on the Mac the gnu compiler works correctly using -O3. On the Mac (a Macbook pro with a new sandy bridge processor) the C-version performs the unit tests three times faster than the (new) Java version. The difference on Win7 is much smaller, since I couldn't use compiler optimization there yet. The performance difference is not significant for my own application of the algo (space flight trajectory optimization) because of the high cost for the eval function. But this is perhaps not true for all other applications. So if we want to include BOBYQA in commons.math in my opinion we shouldn't neglect performance issues generally. Even a direct port of the code has a significant performance disadvantage, any refactoring should try to avoid adding too much additional performance bottlenecks.
          Hide
          Gilles added a comment -

          BOBYQAOptimizer - the actual version I use myself - all "pointers" are replaced
          by Java arrays. This version is much easier to read and faster than the old
          one.

          It's a pity that you didn't mention this version earlier; I already spent quite a few hours replacing the "ScopedPtr" variables. Only a few of them remains in my working version: namely, those that are created as offsets in "w" before calling "trsbox", "rescue", "altmov" and "update".
          Since I also made a few other changes along the way, I don't feel like starting all (almost) over again...
          Hence, I'll continue with my incremental changes; but, at some point, I could use some help to convert the state machine code into proper function calls.

          IMO, we should first arrive at a clearer code before worrying about performance (the more so that, as you pointed out, this algorithm will probably be put to use when function evaluation is expensive, overwhelming the optimizer's code running time).

          The refactoring from efficient Fortran to (very probably less efficient) Java is a big effort but it's indispensable: If it were not to become understandable and maintainable, I don't see the point in including it in CM; you could just provide the straight translation in a JAR file and people would use it as a black box.

          Show
          Gilles added a comment - BOBYQAOptimizer - the actual version I use myself - all "pointers" are replaced by Java arrays. This version is much easier to read and faster than the old one. It's a pity that you didn't mention this version earlier; I already spent quite a few hours replacing the "ScopedPtr" variables. Only a few of them remains in my working version: namely, those that are created as offsets in "w" before calling "trsbox", "rescue", "altmov" and "update". Since I also made a few other changes along the way, I don't feel like starting all (almost) over again... Hence, I'll continue with my incremental changes; but, at some point, I could use some help to convert the state machine code into proper function calls. IMO, we should first arrive at a clearer code before worrying about performance (the more so that, as you pointed out, this algorithm will probably be put to use when function evaluation is expensive, overwhelming the optimizer's code running time). The refactoring from efficient Fortran to (very probably less efficient) Java is a big effort but it's indispensable: If it were not to become understandable and maintainable, I don't see the point in including it in CM; you could just provide the straight translation in a JAR file and people would use it as a black box.
          Hide
          Phil Steitz added a comment -

          I agree strongly with Gilles comment on maintainability. At apache, our goal is to develop software that can be maintained by a group of volunteers that changes over time. In [math], that means algorithms need to be documented and it has to be possible for a Java developer with sufficient mathematical knowledge to read the references and javadoc and make sense of the code.

          Show
          Phil Steitz added a comment - I agree strongly with Gilles comment on maintainability. At apache, our goal is to develop software that can be maintained by a group of volunteers that changes over time. In [math] , that means algorithms need to be documented and it has to be possible for a Java developer with sufficient mathematical knowledge to read the references and javadoc and make sense of the code.
          Hide
          Gilles added a comment -

          Dietmar,

          I think I've come across a bug in function "rescue". At line 2443 (of the Java translation I use for comparison, i.e. the one with "ScopedPtr"):

          i__2 = ndim;
          for (ip = npt; ip <= i__2; ip++) {
            sum += bmat.get(ip + j * bmat_dim1) * w.get(ip);
          }
          

          Whereas, in the original Fortran code (at line 1544):

                DO 220 IP=NPT+1,NDIM
            220 SUM=SUM+BMAT(IP,J)*W(IP)
          

          Can you confirm?

          If indeed, the loop should start at "npt + 1", I've made the change; but tests still all pass! Does this mean that the one exercising rescue is too lenient?
          Trying to remove the "w" ("split" work space) from "rescue" provides many possibilities for my own mistakes; thus I am a little scared that they would also go unnoticed...

          Show
          Gilles added a comment - Dietmar, I think I've come across a bug in function "rescue". At line 2443 (of the Java translation I use for comparison, i.e. the one with "ScopedPtr"): i__2 = ndim; for (ip = npt; ip <= i__2; ip++) { sum += bmat.get(ip + j * bmat_dim1) * w.get(ip); } Whereas, in the original Fortran code (at line 1544): DO 220 IP=NPT+1,NDIM 220 SUM=SUM+BMAT(IP,J)*W(IP) Can you confirm? If indeed, the loop should start at "npt + 1", I've made the change; but tests still all pass! Does this mean that the one exercising rescue is too lenient? Trying to remove the "w" ("split" work space) from "rescue" provides many possibilities for my own mistakes; thus I am a little scared that they would also go unnoticed...
          Hide
          Dr. Dietmar Wolz added a comment -

          Hi Gilles,
          bug is confirmed - and it is also in the new array based version.
          Debugging the unit test reveals sum is 0 here, so for rescue the
          unit tests are not sufficient at the moment.

          Creating the array based version caused a lot of these bugs initially, but I
          would favour this version
          over the old one because it is much cleaner and easier to read and maintain.

          My experience is that if you introduce this kind of bug in any part beside
          rescue, you
          are guaranteed to see a difference in the output I put into the junit tests form
          my newest attachment.

          The number of function calls will differ from the C/Fortran version. But you
          definitinely have
          to compare these function call numbers with the C/Fortran version, to be sure
          everything is

          as intended. Even slight changes like a = a+b replaced by a += b can cause
          trouble
          because of subtle accuracy problems accumulating over time.

          Rescue is a completely different story, therfore I asked Mike Powell whether he
          could
          provide us with meaningful samples/testsm his answer was:

          ---------------
          Excerpt from a mail from Mike Powell:
          Concerning RESCUE, it is present in BOBYQA because, in some tests
          with unattainable accuracy requirements, a little more progress could
          be made by invoking RESCUE. Checking the correctness of RESCUE was
          possible only by modifying the package in a way that forced it to be
          called in situations that were not contaminated severely by loss of
          accuracy. I may decide to delete RESCUE from the Fortran package that
          I send to people who ask for BOBYQA.
          -------------

          To make progress we have a lot of options to decide:

          1) Keep rescue? May be removing it from the initial release in commons.math
          doesn't do too much harm.
          When keeping it we should find useful applications before and add them to the
          unit tests.

          2) Using arrays instead of pointers and don't use the large shared working
          space? I would prefer so, since the
          arrays-version is faster in Java (pointers need to be emulated) and the code
          becomes much clearer.

          3) Complete redesign/refactor the code?
          A difficult issue. I would say that this is very hard to achieve.
          Instead I would try to build something equivalent from scratch.
          Problem is that Mike Powell put 50 years of experience developing optimization
          algos
          into BOBYQA and you can see that comparing the number of cost function calls
          needed with other
          optimization algos. So from the user perspective BOBYQA has a huge value even if
          not completely

          refactored/redesigned and it is not an easy decision to keep BOBYQA out of
          commons.math for
          code-aesthetic reasons.

          I would like to see other opinions on these options.

          Gilles commented on MATH-621:
          -----------------------------

          Dietmar,

          I think I've come across a bug in function "rescue". At line 2443 (of the Java
          translation I use for comparison, i.e. the one with "ScopedPtr"):

          i__2 = ndim;
          for (ip = npt; ip <= i__2; ip++) {
            sum += bmat.get(ip + j * bmat_dim1) * w.get(ip);
          }
          

          Whereas, in the original Fortran code (at line 1544):

                DO 220 IP=NPT+1,NDIM
            220 SUM=SUM+BMAT(IP,J)*W(IP)
          

          Can you confirm?

          If indeed, the loop should start at "npt + 1", I've made the change; but tests
          still all pass! Does this mean that the one exercising rescue is too lenient?
          Trying to remove the "w" ("split" work space) from "rescue" provides many
          possibilities for my own mistakes; thus I am a little scared that they would
          also go unnoticed...


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - Hi Gilles, bug is confirmed - and it is also in the new array based version. Debugging the unit test reveals sum is 0 here, so for rescue the unit tests are not sufficient at the moment. Creating the array based version caused a lot of these bugs initially, but I would favour this version over the old one because it is much cleaner and easier to read and maintain. My experience is that if you introduce this kind of bug in any part beside rescue, you are guaranteed to see a difference in the output I put into the junit tests form my newest attachment. The number of function calls will differ from the C/Fortran version. But you definitinely have to compare these function call numbers with the C/Fortran version, to be sure everything is as intended. Even slight changes like a = a+b replaced by a += b can cause trouble because of subtle accuracy problems accumulating over time. Rescue is a completely different story, therfore I asked Mike Powell whether he could provide us with meaningful samples/testsm his answer was: --------------- Excerpt from a mail from Mike Powell: Concerning RESCUE, it is present in BOBYQA because, in some tests with unattainable accuracy requirements, a little more progress could be made by invoking RESCUE. Checking the correctness of RESCUE was possible only by modifying the package in a way that forced it to be called in situations that were not contaminated severely by loss of accuracy. I may decide to delete RESCUE from the Fortran package that I send to people who ask for BOBYQA. ------------- To make progress we have a lot of options to decide: 1) Keep rescue? May be removing it from the initial release in commons.math doesn't do too much harm. When keeping it we should find useful applications before and add them to the unit tests. 2) Using arrays instead of pointers and don't use the large shared working space? I would prefer so, since the arrays-version is faster in Java (pointers need to be emulated) and the code becomes much clearer. 3) Complete redesign/refactor the code? A difficult issue. I would say that this is very hard to achieve. Instead I would try to build something equivalent from scratch. Problem is that Mike Powell put 50 years of experience developing optimization algos into BOBYQA and you can see that comparing the number of cost function calls needed with other optimization algos. So from the user perspective BOBYQA has a huge value even if not completely refactored/redesigned and it is not an easy decision to keep BOBYQA out of commons.math for code-aesthetic reasons. I would like to see other opinions on these options. Gilles commented on MATH-621 : ----------------------------- Dietmar, I think I've come across a bug in function "rescue". At line 2443 (of the Java translation I use for comparison, i.e. the one with "ScopedPtr"): i__2 = ndim; for (ip = npt; ip <= i__2; ip++) { sum += bmat.get(ip + j * bmat_dim1) * w.get(ip); } Whereas, in the original Fortran code (at line 1544): DO 220 IP=NPT+1,NDIM 220 SUM=SUM+BMAT(IP,J)*W(IP) Can you confirm? If indeed, the loop should start at "npt + 1", I've made the change; but tests still all pass! Does this mean that the one exercising rescue is too lenient? Trying to remove the "w" ("split" work space) from "rescue" provides many possibilities for my own mistakes; thus I am a little scared that they would also go unnoticed... – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Dr. Dietmar Wolz added a comment -

          The new version was not mentioned earlier because it didn't exist. I did
          that on my vacation and uploaded it at the day I got the tests through.

          I suspected you tried something similar but I thought that it would it
          make easier for you to locate bugs in your version if you had something
          working using arrays for comparison.

          Conversion of the state machine is far from trivial, if you want to do
          it properly. Problem is that you either define a lot of global variables
          or have to transfer a huge number of parameters.

          What would the jar-solution mean for the users?
          How will a user find the jar, who will host these kind of black box extensions
          to CM?

          Instead of providing a jar you could also provide a dll/shared libraries
          together with a JNI interface. It is also a black box -
          but up to six times faster according to my tests. But I agree - performance
          is not essentially here because you usually have expensive evaluation functions.

          Having the code published could inspire others to work and improve on it.

          Gilles commented on MATH-621:
          -----------------------------

          It's a pity that you didn't mention this version earlier; I already spent
          quite a few hours replacing the "ScopedPtr" variables. Only a few of them
          remains in my working version: namely, those that are created as offsets in "w"
          before calling "trsbox", "rescue", "altmov" and "update".
          Since I also made a few other changes along the way, I don't feel like starting
          all (almost) over again...
          Hence, I'll continue with my incremental changes; but, at some point, I could
          use some help to convert the state machine code into proper function calls.

          IMO, we should first arrive at a clearer code before worrying about performance
          (the more so that, as you pointed out, this algorithm will probably be put to
          use when function evaluation is expensive, overwhelming the optimizer's code
          running time).

          The refactoring from efficient Fortran to (very probably less efficient) Java is
          a big effort but it's indispensable: If it were not to become understandable and
          maintainable, I don't see the point in including it in CM; you could just
          provide the straight translation in a JAR file and people would use it as a
          black box.


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - The new version was not mentioned earlier because it didn't exist. I did that on my vacation and uploaded it at the day I got the tests through. I suspected you tried something similar but I thought that it would it make easier for you to locate bugs in your version if you had something working using arrays for comparison. Conversion of the state machine is far from trivial, if you want to do it properly. Problem is that you either define a lot of global variables or have to transfer a huge number of parameters. What would the jar-solution mean for the users? How will a user find the jar, who will host these kind of black box extensions to CM? Instead of providing a jar you could also provide a dll/shared libraries together with a JNI interface. It is also a black box - but up to six times faster according to my tests. But I agree - performance is not essentially here because you usually have expensive evaluation functions. Having the code published could inspire others to work and improve on it. Gilles commented on MATH-621 : ----------------------------- It's a pity that you didn't mention this version earlier; I already spent quite a few hours replacing the "ScopedPtr" variables. Only a few of them remains in my working version: namely, those that are created as offsets in "w" before calling "trsbox", "rescue", "altmov" and "update". Since I also made a few other changes along the way, I don't feel like starting all (almost) over again... Hence, I'll continue with my incremental changes; but, at some point, I could use some help to convert the state machine code into proper function calls. IMO, we should first arrive at a clearer code before worrying about performance (the more so that, as you pointed out, this algorithm will probably be put to use when function evaluation is expensive, overwhelming the optimizer's code running time). The refactoring from efficient Fortran to (very probably less efficient) Java is a big effort but it's indispensable: If it were not to become understandable and maintainable, I don't see the point in including it in CM; you could just provide the straight translation in a JAR file and people would use it as a black box. – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Phil Steitz added a comment -

          Dietmar,

          Thanks for working on this. Here are answers to your questions about "black box" solutions.

          What would the jar-solution mean for the users?
          How will a user find the jar, who will host these kind of black box extensions
          to CM?

          One solution is Apache Extras[1]. GitHub or SourceForge are others. We do not host code not maintained by the ASF on ASF infrastructure, so it would have to be external.

          Instead of providing a jar you could also provide a dll/shared libraries
          together with a JNI interface. It is also a black box -
          but up to six times faster according to my tests. But I agree - performance
          is not essentially here because you usually have expensive evaluation functions.

          Again, this is possible, but would have to be hosted externally to Apache.

          Having the code published could inspire others to work and improve on it.

          +1 for this and that we could do here. The key is to get the code to the point where others can follow it so they get "inspired" rather than confused and frustrated.

          [1] http://code.google.com/a/apache-extras.org/hosting/

          Show
          Phil Steitz added a comment - Dietmar, Thanks for working on this. Here are answers to your questions about "black box" solutions. What would the jar-solution mean for the users? How will a user find the jar, who will host these kind of black box extensions to CM? One solution is Apache Extras [1] . GitHub or SourceForge are others. We do not host code not maintained by the ASF on ASF infrastructure, so it would have to be external. Instead of providing a jar you could also provide a dll/shared libraries together with a JNI interface. It is also a black box - but up to six times faster according to my tests. But I agree - performance is not essentially here because you usually have expensive evaluation functions. Again, this is possible, but would have to be hosted externally to Apache. Having the code published could inspire others to work and improve on it. +1 for this and that we could do here. The key is to get the code to the point where others can follow it so they get "inspired" rather than confused and frustrated. [1] http://code.google.com/a/apache-extras.org/hosting/
          Hide
          Dr. Dietmar Wolz added a comment -

          The code is now fully documented - even the internals of the private functions.
          Gilles and me have already removed the main obstacle - the use of pointers und the
          nasty huge working area.
          A further refactoring - splitting the private functions into smaller ones - now is much more work.
          Maybe Gilles can give an estimation of the effort?
          I started and was a bit frustrated soon since I had to introduce so many global variables, I felt the
          code even lost readability. Now the variables are at least local to the (admittedly large) functions.
          I got the feeling that it would be easier to invent something equivalent from scratch describing the
          algorithm mathematically as a set of matrix operations first - similar as it was done for CMA-ES.
          Then a clean implementation (maybe starting with a matlab prototype) could be done.
          I can perform the second step - but for the first one an active researcher in the field would
          be required - a researcher willing to license his work for Apache. As discussed in the mailing
          list there are candidates and Luc offered to approach one of them.

          To summarize the opinions so far:
          Besides Luc and me everyone seems to favor to keep out BOBYQA from CM. It could be
          hosted at Apache Extras, GitHub or SourceForge.
          Gilles continues to try to further enhance the code with the aim to make it "suitable" for
          inclusion into CM.

          Phil Steitz commented on MATH-621:
          ----------------------------------
          +1 for this and that we could do here. The key is to get the code to the point where others can follow it so they get "inspired" rather than confused and frustrated.

          Show
          Dr. Dietmar Wolz added a comment - The code is now fully documented - even the internals of the private functions. Gilles and me have already removed the main obstacle - the use of pointers und the nasty huge working area. A further refactoring - splitting the private functions into smaller ones - now is much more work. Maybe Gilles can give an estimation of the effort? I started and was a bit frustrated soon since I had to introduce so many global variables, I felt the code even lost readability. Now the variables are at least local to the (admittedly large) functions. I got the feeling that it would be easier to invent something equivalent from scratch describing the algorithm mathematically as a set of matrix operations first - similar as it was done for CMA-ES. Then a clean implementation (maybe starting with a matlab prototype) could be done. I can perform the second step - but for the first one an active researcher in the field would be required - a researcher willing to license his work for Apache. As discussed in the mailing list there are candidates and Luc offered to approach one of them. To summarize the opinions so far: Besides Luc and me everyone seems to favor to keep out BOBYQA from CM. It could be hosted at Apache Extras, GitHub or SourceForge. Gilles continues to try to further enhance the code with the aim to make it "suitable" for inclusion into CM. Phil Steitz commented on MATH-621 : ---------------------------------- +1 for this and that we could do here. The key is to get the code to the point where others can follow it so they get "inspired" rather than confused and frustrated.
          Hide
          Gilles added a comment -

          I certainly also would like the code to be in CM. If just for the selfish reason that I've been working on the conversion for many days now . I don't think that anyone wants it out. However, my view of the CM software library is that it should be a repository of best practices for scientific computing in Java, and not "just" a toolbox.
          This implies IMHO that we try our best that the implementations are coded in a proper OO way...

          That said, I'm a bit dismayed that we don't seem to combine our forces to achieve that goal.
          I cannot use your new version, and you prefer it to the old one; that is the situation that I wanted to avoid when I proposed to do a first pass on the code before committing it to the repository!
          Please allow me a few more days to get to a "consistent" state (e.g. removing the big work space). A that point I propose to

          1. commit your original "straight" translation as provided in your first patch
          2. commit my current version of that code

          Then we can use that as a common base to refactor:

          1. Replace the "goto"/state machine constructs
          2. Use matrix operations (from package "linear") whenever possible, in place of the explicit loops

          I don't quite see what seems to bother you with "global" variables versus long argument list. One of the "problems" of the original (Fortran-like) implementation is that many of the function arguments are either

          • input-output (i.e. they are changed in-place), or
          • work space

          The former will be more suitably turned into instance variables, and the latter should become a local variable. I've already started this transformation; it reduces the the length of the argument lists and make the code clearer namely because one will be able to tell which variables are input (function arguments) and which are output (instance variable).

          The clean, from scratch, implementation is also an option. However, a couple of weeks ago, I had the impression that no one stepped forward to work on it. That's also why I took the pedestrian way, having the feeling that we could slowly but surely (thanks to your reference implementation and unit tests) walk towards the goal of having BOBYQA included in CM before the 3.0 release (while waiting for the best person for the task would seem to make it unlikely).

          Once the 2 commits I talked about above are done, the code will still be far from my ideal OO but I hope that it will be in a state that will make it more likely that it will evolve towards that ideal.

          How does that sound?

          Show
          Gilles added a comment - I certainly also would like the code to be in CM. If just for the selfish reason that I've been working on the conversion for many days now . I don't think that anyone wants it out. However, my view of the CM software library is that it should be a repository of best practices for scientific computing in Java, and not "just" a toolbox. This implies IMHO that we try our best that the implementations are coded in a proper OO way... That said, I'm a bit dismayed that we don't seem to combine our forces to achieve that goal. I cannot use your new version, and you prefer it to the old one; that is the situation that I wanted to avoid when I proposed to do a first pass on the code before committing it to the repository! Please allow me a few more days to get to a "consistent" state (e.g. removing the big work space). A that point I propose to commit your original "straight" translation as provided in your first patch commit my current version of that code Then we can use that as a common base to refactor: Replace the "goto"/state machine constructs Use matrix operations (from package "linear") whenever possible, in place of the explicit loops I don't quite see what seems to bother you with "global" variables versus long argument list. One of the "problems" of the original (Fortran-like) implementation is that many of the function arguments are either input-output (i.e. they are changed in-place), or work space The former will be more suitably turned into instance variables, and the latter should become a local variable. I've already started this transformation; it reduces the the length of the argument lists and make the code clearer namely because one will be able to tell which variables are input (function arguments) and which are output (instance variable). The clean, from scratch, implementation is also an option. However, a couple of weeks ago, I had the impression that no one stepped forward to work on it. That's also why I took the pedestrian way, having the feeling that we could slowly but surely (thanks to your reference implementation and unit tests) walk towards the goal of having BOBYQA included in CM before the 3.0 release (while waiting for the best person for the task would seem to make it unlikely). Once the 2 commits I talked about above are done, the code will still be far from my ideal OO but I hope that it will be in a state that will make it more likely that it will evolve towards that ideal. How does that sound?
          Hide
          Gilles added a comment -

          Another problem: The conditionals, at line 567 and at line 570, are never entered:

          if (xnew.get(j) == sl.get(j)) {
            bdtest = w.get(j);
          }
          if (xnew.get(j) == su.get(j)) {
            bdtest = -w.get(j);
          }
          

          In some sense, it would be fine if they could be removed; because otherwise, in my converted code, there would the extreme inconvenience that the content of "w" is accessed before it is initialized...

          Show
          Gilles added a comment - Another problem: The conditionals, at line 567 and at line 570, are never entered: if (xnew.get(j) == sl.get(j)) { bdtest = w.get(j); } if (xnew.get(j) == su.get(j)) { bdtest = -w.get(j); } In some sense, it would be fine if they could be removed; because otherwise, in my converted code, there would the extreme inconvenience that the content of "w" is accessed before it is initialized...
          Hide
          Gilles added a comment -

          Original version (with "ScopedPtr" variables) committed in revision 1154543.

          Show
          Gilles added a comment - Original version (with "ScopedPtr" variables) committed in revision 1154543.
          Hide
          Gilles added a comment -

          My first steps in the (hopefully) good direction were committed in revision 1154550.

          I think that the next step would be to replace all 1-based loops by 0-based loops, together with changing all occurrences of "FortranArray" and "FortranMatrix" by "ArrayRealVector" and "Array2DRealMatrix", respectively. This is a little delicate because all changes must be done at once, and also because some indices are hard-coded and there are some tests on various index bounds (e.g. in "prelim" and "rescue")... However, I hope that Dietmar's new version with arrays (which I haven't looked at yet) will help sort this out.

          Show
          Gilles added a comment - My first steps in the (hopefully) good direction were committed in revision 1154550. I think that the next step would be to replace all 1-based loops by 0-based loops, together with changing all occurrences of "FortranArray" and "FortranMatrix" by "ArrayRealVector" and "Array2DRealMatrix", respectively. This is a little delicate because all changes must be done at once, and also because some indices are hard-coded and there are some tests on various index bounds (e.g. in "prelim" and "rescue")... However, I hope that Dietmar's new version with arrays (which I haven't looked at yet) will help sort this out.
          Hide
          Dr. Dietmar Wolz added a comment -

          Hi Gilles,
          it was not my intention to "separate forces", but the only thing I could do
          during my
          vacation partly without internet access was to focus on one single topic, I
          chose
          the "replace all 1-based loops by 0-based loops" task. Of course I didn't expect
          that you restart from scratch, but that we somehow integrate what we have done
          so far.
          Question is now: should we really go to CM matrices in one step, or use 0-based
          Arrays
          as an intermediate step? I could for instance use your code as a basis and try
          to come
          up with an 0-based equivalent as an intermediate step. What do you think?
          By the way, thanks for your significant efforts.

          Gilles commented on MATH-621:
          That said, I'm a bit dismayed that we don't seem to combine our forces to
          achieve that goal.

          ----- Ursprüngliche Mail ----
          Von: Gilles (JIRA) <jira@apache.org>
          An: drdietmarwolz@yahoo.de
          Gesendet: Samstag, den 6. August 2011, 19:16:27 Uhr
          Betreff: [jira] [Commented] (MATH-621) BOBYQA is missing in optimization

          [
          https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13080431#comment-13080431
          ]

          Gilles commented on MATH-621:
          -----------------------------

          My first steps in the (hopefully) good direction were committed in revision
          1154550.

          I think that the next step would be to replace all 1-based loops by 0-based
          loops, together with changing all occurrences of "FortranArray" and
          "FortranMatrix" by "ArrayRealVector" and "Array2DRealMatrix", respectively. This
          is a little delicate because all changes must be done at once, and also because
          some indices are hard-coded and there are some tests on various index bounds
          (e.g. in "prelim" and "rescue")... However, I hope that Dietmar's new version
          with arrays (which I haven't looked at yet) will help sort this out.

          Show
          Dr. Dietmar Wolz added a comment - Hi Gilles, it was not my intention to "separate forces", but the only thing I could do during my vacation partly without internet access was to focus on one single topic, I chose the "replace all 1-based loops by 0-based loops" task. Of course I didn't expect that you restart from scratch, but that we somehow integrate what we have done so far. Question is now: should we really go to CM matrices in one step, or use 0-based Arrays as an intermediate step? I could for instance use your code as a basis and try to come up with an 0-based equivalent as an intermediate step. What do you think? By the way, thanks for your significant efforts. Gilles commented on MATH-621 : That said, I'm a bit dismayed that we don't seem to combine our forces to achieve that goal. ----- Ursprüngliche Mail ---- Von: Gilles (JIRA) <jira@apache.org> An: drdietmarwolz@yahoo.de Gesendet: Samstag, den 6. August 2011, 19:16:27 Uhr Betreff: [jira] [Commented] ( MATH-621 ) BOBYQA is missing in optimization [ https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13080431#comment-13080431 ] Gilles commented on MATH-621 : ----------------------------- My first steps in the (hopefully) good direction were committed in revision 1154550. I think that the next step would be to replace all 1-based loops by 0-based loops, together with changing all occurrences of "FortranArray" and "FortranMatrix" by "ArrayRealVector" and "Array2DRealMatrix", respectively. This is a little delicate because all changes must be done at once, and also because some indices are hard-coded and there are some tests on various index bounds (e.g. in "prelim" and "rescue")... However, I hope that Dietmar's new version with arrays (which I haven't looked at yet) will help sort this out.
          Hide
          Gilles added a comment -

          Hello Dietmar.

          Of course I didn't expect
          that you restart from scratch, but that we somehow integrate what we have done
          so far.

          But I did not restart from scratch; I used your Java translation and went from there!
          In fact, it looks like I used the original Fortran because I indeed did transform the "ScopedPtr" that were inherently bi-dimensional (i.e. matrices like "bmat" and "zmat") into (new auxiliary) "FortranMatrix" objects.
          It was only when I checked for the bug in "rescue" that I looked into the original Fortran, and discovered that it was using matrices! I still don't understand why they were translated into 1-d arrays in your code...

          Question is now: should we really go to CM matrices in one step, or use 0-based
          Arrays as an intermediate step? I could for instance use your code as a basis and try
          to come up with an 0-based equivalent as an intermediate step. What do you think?

          Going to 0-based loops and replace "FortranArray" and "FortranMatrix" with their CM-equivalent ("ArrayRealVector" and "Array2DRowRealMatrix"), which are 0-based, can only be done in one step, if I'm not mistaken.

          I've started to write a script that would do the translation of everything that can be spotted automatically i.e. construct like

          for (int i = 1; i < n; i++)
          

          and

          zmat.getEntry(2, j)
          

          and

          FortranArray glag = new FortranArray(n);
          

          etc.

          But there are some contructs that must be changed concomitantly (like some test on array bounds in "prelim") and cannot be done automatically. I also suspect that some simplification could be done there because of the use of a variable containing the "number of evaluations" + 1. Unfortunately, my first attempt was not successful

          So, what I suggest is that

          • I finish my (Perl) script; I'll test that it makes all the intended replacements (but obviously the resulting code will not pass the tests anymore until all the non-trivial replacements have been correctly performed).
          • I'll post it here so that you can run it on your working copy
          • you could then attempt to perform the rest of the 1-based to 0-based conversion.

          Once we are left with "ArrayRealVector" and "Array2DRowRealMatrix", we can see whether some of the explicit loops can replaced by matrix operation methods from CM's "RealMatrix" interface.
          After this "low-level clean-up" we can discuss how to introduce the "bounded optimization" concept into the CM API (I've already marked the code with "XXX" to that purpose).
          Concurrently to the preceding point, it would be nice to gradually tackle the "goto" problem.

          What do you think?

          Show
          Gilles added a comment - Hello Dietmar. Of course I didn't expect that you restart from scratch, but that we somehow integrate what we have done so far. But I did not restart from scratch; I used your Java translation and went from there! In fact, it looks like I used the original Fortran because I indeed did transform the "ScopedPtr" that were inherently bi-dimensional (i.e. matrices like "bmat" and "zmat") into (new auxiliary) "FortranMatrix" objects. It was only when I checked for the bug in "rescue" that I looked into the original Fortran, and discovered that it was using matrices! I still don't understand why they were translated into 1-d arrays in your code... Question is now: should we really go to CM matrices in one step, or use 0-based Arrays as an intermediate step? I could for instance use your code as a basis and try to come up with an 0-based equivalent as an intermediate step. What do you think? Going to 0-based loops and replace "FortranArray" and "FortranMatrix" with their CM-equivalent ("ArrayRealVector" and "Array2DRowRealMatrix"), which are 0-based, can only be done in one step, if I'm not mistaken. I've started to write a script that would do the translation of everything that can be spotted automatically i.e. construct like for ( int i = 1; i < n; i++) and zmat.getEntry(2, j) and FortranArray glag = new FortranArray(n); etc. But there are some contructs that must be changed concomitantly (like some test on array bounds in "prelim") and cannot be done automatically. I also suspect that some simplification could be done there because of the use of a variable containing the "number of evaluations" + 1. Unfortunately, my first attempt was not successful So, what I suggest is that I finish my (Perl) script; I'll test that it makes all the intended replacements (but obviously the resulting code will not pass the tests anymore until all the non-trivial replacements have been correctly performed). I'll post it here so that you can run it on your working copy you could then attempt to perform the rest of the 1-based to 0-based conversion. Once we are left with "ArrayRealVector" and "Array2DRowRealMatrix", we can see whether some of the explicit loops can replaced by matrix operation methods from CM's "RealMatrix" interface. After this "low-level clean-up" we can discuss how to introduce the "bounded optimization" concept into the CM API (I've already marked the code with "XXX" to that purpose). Concurrently to the preceding point, it would be nice to gradually tackle the "goto" problem. What do you think?
          Hide
          Dr. Dietmar Wolz added a comment -

          If you already started writing scripts we should see what they generate.

          I needed several days fixing the indices / removing all newly introduced bugs
          when I moved to 0-based Java arrays.

          My hope is that the indexing problems can be solved by looking into my 0-based
          Java array version.

          maybe because I used the C++-port in dlib?
          http://dlib.net/dlib/optimization/optimization_bobyqa.h.html
          there are no matrices.

          But I agree, starting with matrices and automatically transform them into
          a Array2DRowRealMatrix seems
          a good idea.

          ----- Ursprüngliche Mail ----
          Von: Gilles (JIRA) <jira@apache.org>
          An: drdietmarwolz@yahoo.de
          Gesendet: Mittwoch, den 10. August 2011, 16:19:27 Uhr
          Betreff: [jira] [Commented] (MATH-621) BOBYQA is missing in optimization

          [
          https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13082361#comment-13082361
          ]

          Gilles commented on MATH-621:
          -----------------------------

          Hello Dietmar.

          Of course I didn't expect
          that you restart from scratch, but that we somehow integrate what we have done
          so far.

          But I did not restart from scratch; I used your Java translation and went from
          there!
          In fact, it looks like I used the original Fortran because I indeed did
          transform the "ScopedPtr" that were inherently bi-dimensional (i.e. matrices
          like "bmat" and "zmat") into (new auxiliary) "FortranMatrix" objects.
          It was only when I checked for the bug in "rescue" that I looked into the
          original Fortran, and discovered that it was using matrices! I still don't
          understand why they were translated into 1-d arrays in your code...

          Question is now: should we really go to CM matrices in one step, or use 0-based
          Arrays as an intermediate step? I could for instance use your code as a basis
          and try
          to come up with an 0-based equivalent as an intermediate step. What do you
          think?

          Going to 0-based loops and replace "FortranArray" and "FortranMatrix" with their
          CM-equivalent ("ArrayRealVector" and "Array2DRowRealMatrix"), which are 0-based,
          can only be done in one step, if I'm not mistaken.

          I've started to write a script that would do the translation of everything that
          can be spotted automatically i.e. construct like

          for (int i = 1; i < n; i++)
          

          and

          zmat.getEntry(2, j)
          

          and

          FortranArray glag = new FortranArray(n);
          

          etc.

          But there are some contructs that must be changed concomitantly (like some test
          on array bounds in "prelim") and cannot be done automatically. I also suspect
          that some simplification could be done there because of the use of a variable
          containing the "number of evaluations" + 1. Unfortunately, my first attempt was
          not successful

          So, what I suggest is that

          • I finish my (Perl) script; I'll test that it makes all the intended
            replacements (but obviously the resulting code will not pass the tests anymore
            until all the non-trivial replacements have been correctly performed).
          • I'll post it here so that you can run it on your working copy
          • you could then attempt to perform the rest of the 1-based to 0-based
            conversion.

          Once we are left with "ArrayRealVector" and "Array2DRowRealMatrix", we can see
          whether some of the explicit loops can replaced by matrix operation methods from
          CM's "RealMatrix" interface.
          After this "low-level clean-up" we can discuss how to introduce the "bounded
          optimization" concept into the CM API (I've already marked the code with "XXX"
          to that purpose).
          Concurrently to the preceding point, it would be nice to gradually tackle the
          "goto" problem.

          What do you think?


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - If you already started writing scripts we should see what they generate. I needed several days fixing the indices / removing all newly introduced bugs when I moved to 0-based Java arrays. My hope is that the indexing problems can be solved by looking into my 0-based Java array version. maybe because I used the C++-port in dlib? http://dlib.net/dlib/optimization/optimization_bobyqa.h.html there are no matrices. But I agree, starting with matrices and automatically transform them into a Array2DRowRealMatrix seems a good idea. ----- Ursprüngliche Mail ---- Von: Gilles (JIRA) <jira@apache.org> An: drdietmarwolz@yahoo.de Gesendet: Mittwoch, den 10. August 2011, 16:19:27 Uhr Betreff: [jira] [Commented] ( MATH-621 ) BOBYQA is missing in optimization [ https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13082361#comment-13082361 ] Gilles commented on MATH-621 : ----------------------------- Hello Dietmar. Of course I didn't expect that you restart from scratch, but that we somehow integrate what we have done so far. But I did not restart from scratch; I used your Java translation and went from there! In fact, it looks like I used the original Fortran because I indeed did transform the "ScopedPtr" that were inherently bi-dimensional (i.e. matrices like "bmat" and "zmat") into (new auxiliary) "FortranMatrix" objects. It was only when I checked for the bug in "rescue" that I looked into the original Fortran, and discovered that it was using matrices! I still don't understand why they were translated into 1-d arrays in your code... Question is now: should we really go to CM matrices in one step, or use 0-based Arrays as an intermediate step? I could for instance use your code as a basis and try to come up with an 0-based equivalent as an intermediate step. What do you think? Going to 0-based loops and replace "FortranArray" and "FortranMatrix" with their CM-equivalent ("ArrayRealVector" and "Array2DRowRealMatrix"), which are 0-based, can only be done in one step, if I'm not mistaken. I've started to write a script that would do the translation of everything that can be spotted automatically i.e. construct like for ( int i = 1; i < n; i++) and zmat.getEntry(2, j) and FortranArray glag = new FortranArray(n); etc. But there are some contructs that must be changed concomitantly (like some test on array bounds in "prelim") and cannot be done automatically. I also suspect that some simplification could be done there because of the use of a variable containing the "number of evaluations" + 1. Unfortunately, my first attempt was not successful So, what I suggest is that I finish my (Perl) script; I'll test that it makes all the intended replacements (but obviously the resulting code will not pass the tests anymore until all the non-trivial replacements have been correctly performed). I'll post it here so that you can run it on your working copy you could then attempt to perform the rest of the 1-based to 0-based conversion. Once we are left with "ArrayRealVector" and "Array2DRowRealMatrix", we can see whether some of the explicit loops can replaced by matrix operation methods from CM's "RealMatrix" interface. After this "low-level clean-up" we can discuss how to introduce the "bounded optimization" concept into the CM API (I've already marked the code with "XXX" to that purpose). Concurrently to the preceding point, it would be nice to gradually tackle the "goto" problem. What do you think? – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Gilles added a comment -

          Here is a script that converts "for" loops from 1-based to 0-based.

          Show
          Gilles added a comment - Here is a script that converts "for" loops from 1-based to 0-based.
          Hide
          Gilles added a comment -

          Introducing the "INDEX_OFFSET" constant in my last commit made the script simpler; it only converts "for" loops (hopefully all of them). It also changes the value of "INDEX_OFFSET" from 1 to 0, allowing to leave the auxiliary "FortranArray" and "FortranMatrix" in place (they just become no-op wrappers).

          So, do you want to try correcting the remaining indices?
          The first error occurs in "prelim" (at line 1779)...

          A good strategy might be to insert "INDEX_OFFSET" at the appropriate places so that when reverting "INDEX_OFFSET" to 1, we can still test that the behaviour was as before. I thought that it could then help refactoring the "do while" loop into something more understandable.

          Show
          Gilles added a comment - Introducing the "INDEX_OFFSET" constant in my last commit made the script simpler; it only converts "for" loops (hopefully all of them). It also changes the value of "INDEX_OFFSET" from 1 to 0, allowing to leave the auxiliary "FortranArray" and "FortranMatrix" in place (they just become no-op wrappers). So, do you want to try correcting the remaining indices? The first error occurs in "prelim" (at line 1779)... A good strategy might be to insert "INDEX_OFFSET" at the appropriate places so that when reverting "INDEX_OFFSET" to 1, we can still test that the behaviour was as before. I thought that it could then help refactoring the "do while" loop into something more understandable.
          Hide
          Dr. Dietmar Wolz added a comment -

          Will try to fix the offsets, but this will probably take some time.

          Show
          Dr. Dietmar Wolz added a comment - Will try to fix the offsets, but this will probably take some time.
          Hide
          Dr. Dietmar Wolz added a comment -

          No changes from the perl generated code
          beside the ones necessary to get INDEX_OFFSET=0 working. Introduced INDEX_OFFSET where possible but there were
          many other adaptions necessary (just compare the perl generated code with the attachment). Version 0.3 had some useful
          additional minor changes/refactorings missing here (see remarks below),
          but the main work for 0.3 was the index change, and this we have here again. Remarks:

          1) The perl script has damaged the "for loop" intendation

          2) n, npt and nptm should be global variables and not set separately
          in each method

          3) "System generated locals": Declare variables in the scope they are needed and
          not method-globally if not necessary

          4) testDiagonalRosen() is a copy/paste leftover from CMAES, should be removed

          5) We should shink about removing rescue as proposed by Mike Powell.

          Show
          Dr. Dietmar Wolz added a comment - No changes from the perl generated code beside the ones necessary to get INDEX_OFFSET=0 working. Introduced INDEX_OFFSET where possible but there were many other adaptions necessary (just compare the perl generated code with the attachment). Version 0.3 had some useful additional minor changes/refactorings missing here (see remarks below), but the main work for 0.3 was the index change, and this we have here again. Remarks: 1) The perl script has damaged the "for loop" intendation 2) n, npt and nptm should be global variables and not set separately in each method 3) "System generated locals": Declare variables in the scope they are needed and not method-globally if not necessary 4) testDiagonalRosen() is a copy/paste leftover from CMAES, should be removed 5) We should shink about removing rescue as proposed by Mike Powell.
          Hide
          Gilles added a comment -

          Thanks for the work.
          However, if I change the "INDEX_OFFSET" constant (setting it back to "1"), the tests fail.
          I see that you hard-coded the offset in most places instead of using "INDEX_OFFSET". I still think that this place-holder would be useful to keep track of places where the index variables might have been set to fit with the Fortran 1-based counting... Don't you?

          The perl script has damaged the "for loop" intendation

          Sorry, I didn't see that. But that's easy to fix. I'll do it after the issue with INDEX_OFFSET is settled.

          n, npt and nptm should be global variables and not set separately
          in each method

          Yes, I agree. But there are probably many other variables for which this is true ("zmat", "bmat", etc).

          "System generated locals": Declare variables in the scope they are needed [...]

          Agreed, of course. I had started to do that mainly with "d__1"; then there are many cases where the same variable was reused whereas we would prefer to create yet another one with a more explicit name.

          testDiagonalRosen() is a copy/paste leftover from CMAES, should be removed

          OK, I'll do it in the next commit.

          We should shink about removing rescue as proposed by Mike Powell.

          I'm all for anything that leads to removing unnecessary lines of code
          If you are indeed confident that, in most cases, the added complexity is not worth it, I'll just delete it.

          Show
          Gilles added a comment - Thanks for the work. However, if I change the "INDEX_OFFSET" constant (setting it back to "1"), the tests fail. I see that you hard-coded the offset in most places instead of using "INDEX_OFFSET". I still think that this place-holder would be useful to keep track of places where the index variables might have been set to fit with the Fortran 1-based counting... Don't you? The perl script has damaged the "for loop" intendation Sorry, I didn't see that. But that's easy to fix. I'll do it after the issue with INDEX_OFFSET is settled. n, npt and nptm should be global variables and not set separately in each method Yes, I agree. But there are probably many other variables for which this is true ("zmat", "bmat", etc). "System generated locals": Declare variables in the scope they are needed [...] Agreed, of course. I had started to do that mainly with "d__1"; then there are many cases where the same variable was reused whereas we would prefer to create yet another one with a more explicit name. testDiagonalRosen() is a copy/paste leftover from CMAES, should be removed OK, I'll do it in the next commit. We should shink about removing rescue as proposed by Mike Powell. I'm all for anything that leads to removing unnecessary lines of code If you are indeed confident that, in most cases, the added complexity is not worth it, I'll just delete it.
          Hide
          Dr. Dietmar Wolz added a comment -

          I see that you hard-coded the offset in most places instead of using "INDEX_OFFSET". I still think that this place-holder would be useful to keep track of places where the index variables might have been set to fit with the Fortran 1-based counting... Don't you?

          I am not convinced yet. I thought INDEX_OFFSET as a tool to support the conversion. If you don't use
          INDEX_OFFSET in the for loops (for int i = INDEX_OFFSET ...) I don't see why to introduce it artificially
          in other places. The final aim should be to get rid of the Fortran-Arrays/Matrices and have 0-based access. I don't see
          it essential to maintain INDEX_OFFSET as a kind of back reference to the old Fortran code in the future.
          We have the unit tests as regression test.

          Just try to convert one method - lets say prelim - the way you want to have it.
          The working 0-based version 0.4 should make this easy. Then lets have a look at it.
          I suspect it to become rather ugly using INDEX_OFFSET in all places. But then we
          also should convert the for loops as (for int i = INDEX_OFFSET ...) so that the code runs
          again with INDEX_OFFSET=1. If you then really think it is better this way, I will help to
          convert the other methods.

          Show
          Dr. Dietmar Wolz added a comment - I see that you hard-coded the offset in most places instead of using "INDEX_OFFSET". I still think that this place-holder would be useful to keep track of places where the index variables might have been set to fit with the Fortran 1-based counting... Don't you? I am not convinced yet. I thought INDEX_OFFSET as a tool to support the conversion. If you don't use INDEX_OFFSET in the for loops (for int i = INDEX_OFFSET ...) I don't see why to introduce it artificially in other places. The final aim should be to get rid of the Fortran-Arrays/Matrices and have 0-based access. I don't see it essential to maintain INDEX_OFFSET as a kind of back reference to the old Fortran code in the future. We have the unit tests as regression test. Just try to convert one method - lets say prelim - the way you want to have it. The working 0-based version 0.4 should make this easy. Then lets have a look at it. I suspect it to become rather ugly using INDEX_OFFSET in all places. But then we also should convert the for loops as (for int i = INDEX_OFFSET ...) so that the code runs again with INDEX_OFFSET=1. If you then really think it is better this way, I will help to convert the other methods.
          Hide
          Gilles added a comment -

          OK. Keeping INDEX_OFFSET might be more work than really useful. I'll remove it also.

          Show
          Gilles added a comment - OK. Keeping INDEX_OFFSET might be more work than really useful. I'll remove it also.
          Hide
          Gilles added a comment -

          1-based indexing issue solved in revision 1158015.

          Show
          Gilles added a comment - 1-based indexing issue solved in revision 1158015.
          Hide
          Gilles added a comment -

          Removed "testDiagonalRosen" unit test in revision 1158017.

          Show
          Gilles added a comment - Removed "testDiagonalRosen" unit test in revision 1158017.
          Hide
          Gilles added a comment -

          Commenting out "rescue" (line 671) makes the "testRescue" test fail, as expected. So, if I also remove the test, we are fine. However, do you know whether I can also remove the whole "case 190" (lines 667-697) as well as any code that references that "state" (e.g. lines 791-796, 846-851, 2597-2599, etc.)?

          Show
          Gilles added a comment - Commenting out "rescue" (line 671) makes the "testRescue" test fail, as expected. So, if I also remove the test, we are fine. However, do you know whether I can also remove the whole "case 190" (lines 667-697) as well as any code that references that "state" (e.g. lines 791-796, 846-851, 2597-2599, etc.)?
          Hide
          Dr. Dietmar Wolz added a comment -

          "also remove the whole "case 190" (lines 667-697) "

          My gut feeling is "yes". We could just do it and rely on our unit tests.
          Alternatively I could ask Mike Powell. It doesn't seem that case 190
          does anything beside adaptions after rescue.

          ----- Ursprüngliche Mail ----
          Von: Gilles (JIRA) <jira@apache.org>
          An: drdietmarwolz@yahoo.de
          Gesendet: Montag, den 15. August 2011, 23:27:27 Uhr
          Betreff: [jira] [Commented] (MATH-621) BOBYQA is missing in optimization

          [
          https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13085340#comment-13085340
          ]

          Gilles commented on MATH-621:
          -----------------------------

          Commenting out "rescue" (line 671) makes the "testRescue" test fail, as
          expected. So, if I also remove the test, we are fine. However, do you know
          whether I can also remove the whole "case 190" (lines 667-697) as well as any
          code that references that "state" (e.g. lines 791-796, 846-851, 2597-2599,
          etc.)?


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - "also remove the whole "case 190" (lines 667-697) " My gut feeling is "yes". We could just do it and rely on our unit tests. Alternatively I could ask Mike Powell. It doesn't seem that case 190 does anything beside adaptions after rescue. ----- Ursprüngliche Mail ---- Von: Gilles (JIRA) <jira@apache.org> An: drdietmarwolz@yahoo.de Gesendet: Montag, den 15. August 2011, 23:27:27 Uhr Betreff: [jira] [Commented] ( MATH-621 ) BOBYQA is missing in optimization [ https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13085340#comment-13085340 ] Gilles commented on MATH-621 : ----------------------------- Commenting out "rescue" (line 671) makes the "testRescue" test fail, as expected. So, if I also remove the test, we are fine. However, do you know whether I can also remove the whole "case 190" (lines 667-697) as well as any code that references that "state" (e.g. lines 791-796, 846-851, 2597-2599, etc.)? – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Gilles added a comment -

          "rescue" removed in revision 1158448.

          Show
          Gilles added a comment - "rescue" removed in revision 1158448.
          Hide
          Dr. Dietmar Wolz added a comment -

          So finally we use the cm vectors / matrices, fine.
          Whats next? My feeling is that before we look at the case statements we should do some easy stuff:
          1) remove unused variables (eclipse shows some ugly yellow warnings)
          2) reorganize variable use (limit their scope). Should be done with the case statement removal in mind.
          3) would be great if we could use/identify some of the ArrayRealVector / Array2DRowRealMatrix operations (beside getEntry/setEntry)
          We should coordinate 1) + 2), 3) may perhaps be done independently -
          and later merged.

          Show
          Dr. Dietmar Wolz added a comment - So finally we use the cm vectors / matrices, fine. Whats next? My feeling is that before we look at the case statements we should do some easy stuff: 1) remove unused variables (eclipse shows some ugly yellow warnings) 2) reorganize variable use (limit their scope). Should be done with the case statement removal in mind. 3) would be great if we could use/identify some of the ArrayRealVector / Array2DRowRealMatrix operations (beside getEntry/setEntry) We should coordinate 1) + 2), 3) may perhaps be done independently - and later merged.
          Hide
          Gilles added a comment -

          remove unused variables

          Would you take that one?

          reorganize variable use (limit their scope)

          I've done this for "altmov" and "update".
          Let me know which function(s) you'd like to deal with.

          would be great if we could use/identify some of the ArrayRealVector [...]

          This can be left for after the other two points I think.

          When you make changes, please provide a patch for each step separately: i.e. one patch for point 1, one patch for each function, for point 2.

          Show
          Gilles added a comment - remove unused variables Would you take that one? reorganize variable use (limit their scope) I've done this for "altmov" and "update". Let me know which function(s) you'd like to deal with. would be great if we could use/identify some of the ArrayRealVector [...] This can be left for after the other two points I think. When you make changes, please provide a patch for each step separately: i.e. one patch for point 1, one patch for each function, for point 2.
          Hide
          Dr. Dietmar Wolz added a comment -

          >>remove unused variables
          >Would you take that one?

          patch related to "remove unused variables"
          and "remove unused imports"

          Show
          Dr. Dietmar Wolz added a comment - >>remove unused variables >Would you take that one? patch related to "remove unused variables" and "remove unused imports"
          Hide
          Gilles added a comment -

          The patch does not apply cleanly on the last revision. Are you submitted to the "commits" ML?
          If not, we should really coordinate. In my mind, that meant that while waiting for your answer to my questions here, I can continue updating the code...

          "remove unused imports"

          Please don't do that. Or in a separate patch. Or not with Eclipse: Because it reordered the "import" lines, out of 12 changed lines, only 2 were really suppressing unused import statements.

          Show
          Gilles added a comment - The patch does not apply cleanly on the last revision. Are you submitted to the "commits" ML? If not, we should really coordinate . In my mind, that meant that while waiting for your answer to my questions here, I can continue updating the code... "remove unused imports" Please don't do that. Or in a separate patch. Or not with Eclipse: Because it reordered the "import" lines, out of 12 changed lines, only 2 were really suppressing unused import statements.
          Hide
          Dr. Dietmar Wolz added a comment -

          I checked out immediately before producing the patch immediately before
          attaching it here, strange that
          there was a conflict. Is it a good idea to maintain a import order inconsistent
          with "organize imports" in
          eclipse? This makes "organize imports" useless. But ok, so should I again
          produce a patch without
          the imports? But I cannot avoid a conflict if you check in just as I produce the
          patch. By the way: such situations
          are much easier to handle using git, but this is no option here.

          ----- Ursprüngliche Mail ----
          Von: Gilles (JIRA) <jira@apache.org>
          An: drdietmarwolz@yahoo.de
          Gesendet: Donnerstag, den 18. August 2011, 13:42:27 Uhr
          Betreff: [jira] [Commented] (MATH-621) BOBYQA is missing in optimization

          [
          https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13086967#comment-13086967
          ]

          Gilles commented on MATH-621:
          -----------------------------

          The patch does not apply cleanly on the last revision. Are you submitted to the
          "commits" ML?
          If not, we should really coordinate. In my mind, that meant that while waiting
          for your answer to my questions here, I can continue updating the code...

          "remove unused imports"

          Please don't do that. Or in a separate patch. Or not with Eclipse: Because it
          reordered the "import" lines, out of 12 changed lines, only 2 were really
          suppressing unused import statements.


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - I checked out immediately before producing the patch immediately before attaching it here, strange that there was a conflict. Is it a good idea to maintain a import order inconsistent with "organize imports" in eclipse? This makes "organize imports" useless. But ok, so should I again produce a patch without the imports? But I cannot avoid a conflict if you check in just as I produce the patch. By the way: such situations are much easier to handle using git, but this is no option here. ----- Ursprüngliche Mail ---- Von: Gilles (JIRA) <jira@apache.org> An: drdietmarwolz@yahoo.de Gesendet: Donnerstag, den 18. August 2011, 13:42:27 Uhr Betreff: [jira] [Commented] ( MATH-621 ) BOBYQA is missing in optimization [ https://issues.apache.org/jira/browse/MATH-621?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13086967#comment-13086967 ] Gilles commented on MATH-621 : ----------------------------- The patch does not apply cleanly on the last revision. Are you submitted to the "commits" ML? If not, we should really coordinate . In my mind, that meant that while waiting for your answer to my questions here, I can continue updating the code... "remove unused imports" Please don't do that. Or in a separate patch. Or not with Eclipse: Because it reordered the "import" lines, out of 12 changed lines, only 2 were really suppressing unused import statements. – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Gilles added a comment -

          Oops, sorry; I had not updated the to latest revision! Working on two computers, I forgot to coordinate with myself.
          Patch committed in revision 1159203.

          Is it a good idea to maintain a import order inconsistent with "organize imports" in eclipse?

          Not everybody use Eclipse...

          Show
          Gilles added a comment - Oops, sorry; I had not updated the to latest revision! Working on two computers, I forgot to coordinate with myself. Patch committed in revision 1159203. Is it a good idea to maintain a import order inconsistent with "organize imports" in eclipse? Not everybody use Eclipse...
          Hide
          Dr. Dietmar Wolz added a comment -

          So you want to force Eclipse users not to use the "organize import"
          functionality?
          What is the advantage for the non-Eclipse-users? Are there other tools with
          "organize import" functionality which produce a different order of imports?


          This message is automatically generated by JIRA.
          For more information on JIRA, see: http://www.atlassian.com/software/jira

          Show
          Dr. Dietmar Wolz added a comment - So you want to force Eclipse users not to use the "organize import" functionality? What is the advantage for the non-Eclipse-users? Are there other tools with "organize import" functionality which produce a different order of imports? – This message is automatically generated by JIRA. For more information on JIRA, see: http://www.atlassian.com/software/jira
          Hide
          Gilles added a comment -

          So you want to force Eclipse users not to use the "organize import" functionality?

          No; just that, if you want to perform such cleanup, you should provide it in a separate patch. Otherwise I can get complaints: The rule is that a commit should aim at fixing one distinct issue at a time.

          Show
          Gilles added a comment - So you want to force Eclipse users not to use the "organize import" functionality? No; just that, if you want to perform such cleanup, you should provide it in a separate patch. Otherwise I can get complaints: The rule is that a commit should aim at fixing one distinct issue at a time.
          Hide
          Gilles added a comment -

          An unexplored code path: See revision 1159438.

          Show
          Gilles added a comment - An unexplored code path: See revision 1159438.
          Hide
          Gilles added a comment -

          I tried to replace a loop with a matrix operation: The computation of the first "sum" in "case 90" (line 564). Although the new and old computations differ by less than 1e-15, they induce failures:

          Failed tests: 
            testAckley(org.apache.commons.math.optimization.direct.BOBYQAOptimizerTest): expected:<0.0> but was:<1.0731970423449866E-8>
          
          Tests in error: 
            testDiffPow(org.apache.commons.math.optimization.direct.BOBYQAOptimizerTest): illegal state: maximal count (12,000) exceeded: evaluations
          

          The first is possibly acceptable due to numerical errors but the second seems more problematic.

          Anyways, it looks like these kinds of code transformation should not be attempted at this point.

          Show
          Gilles added a comment - I tried to replace a loop with a matrix operation: The computation of the first "sum" in "case 90" (line 564). Although the new and old computations differ by less than 1e-15, they induce failures: Failed tests: testAckley(org.apache.commons.math.optimization.direct.BOBYQAOptimizerTest): expected:<0.0> but was:<1.0731970423449866E-8> Tests in error: testDiffPow(org.apache.commons.math.optimization.direct.BOBYQAOptimizerTest): illegal state: maximal count (12,000) exceeded: evaluations The first is possibly acceptable due to numerical errors but the second seems more problematic. Anyways, it looks like these kinds of code transformation should not be attempted at this point.
          Hide
          Gilles added a comment -

          "prelim" cleaned up in revision 1159792.

          Another occurrence of a supposedly innocuous change (dividing by a number vs multiplying by the inverse of that number) induces the same two unit tests to fail.

          Show
          Gilles added a comment - "prelim" cleaned up in revision 1159792. Another occurrence of a supposedly innocuous change (dividing by a number vs multiplying by the inverse of that number) induces the same two unit tests to fail.
          Hide
          Gilles added a comment -

          Many code branches in "prelim" unexplored by the current set of unit tests.
          The unused code paths are marked by throwing an exception as their last statement.

          Show
          Gilles added a comment - Many code branches in "prelim" unexplored by the current set of unit tests. The unused code paths are marked by throwing an exception as their last statement.
          Hide
          Dr. Dietmar Wolz added a comment -

          Problems with acurracy after refactoring are common with bobyqa. I always monitored also the number of function evaluations for each unit test. It is nontrivial to estimate which change is tolerable and which not. But if we plan further refactorings like a more extensive use of the CM matrix methods, we have to accept that we modify also the behaviour of the algorithm and have to adapt the tests. People who prefer the original can use the attachments here.

          Show
          Dr. Dietmar Wolz added a comment - Problems with acurracy after refactoring are common with bobyqa. I always monitored also the number of function evaluations for each unit test. It is nontrivial to estimate which change is tolerable and which not. But if we plan further refactorings like a more extensive use of the CM matrix methods, we have to accept that we modify also the behaviour of the algorithm and have to adapt the tests. People who prefer the original can use the attachments here.
          Hide
          Gilles added a comment -

          Looking at the unit tests I see that, in the "assertEquals", tolerance for the function values are not very stringent (e.g. 1e-8 instead of, say, 1e-15); so this means that your initial translation was already not expected to reproduce exactly the "reference" results (Fortran or C++ ). Am I correct?

          So, in the "Ackley" case, the tolerance is set to 1e-8, and the alternative computation produces "1.07319e-8" instead of "0.0". Thus the absolute change in tolerance is smaller than 1e-10, which might be put on the account of error propagation .

          However, in "DiffPow", the number of evaluations with the your original code is 6016, while the alternative computation makes it fail even as 12000 are allowed. This looks like a much more drastic change in behaviour!
          Pardon my possible ignorance, but it makes we wonder whether the test is not too fragile or whether there is a bug in the implementation... In other words it would have been a "lucky accident" that the original computation succeeds.

          Show
          Gilles added a comment - Looking at the unit tests I see that, in the "assertEquals", tolerance for the function values are not very stringent (e.g. 1e-8 instead of, say, 1e-15); so this means that your initial translation was already not expected to reproduce exactly the "reference" results (Fortran or C++ ). Am I correct? So, in the "Ackley" case, the tolerance is set to 1e-8, and the alternative computation produces "1.07319e-8" instead of "0.0". Thus the absolute change in tolerance is smaller than 1e-10, which might be put on the account of error propagation . However, in "DiffPow", the number of evaluations with the your original code is 6016, while the alternative computation makes it fail even as 12000 are allowed. This looks like a much more drastic change in behaviour! Pardon my possible ignorance, but it makes we wonder whether the test is not too fragile or whether there is a bug in the implementation... In other words it would have been a "lucky accident" that the original computation succeeds.
          Hide
          Dr. Dietmar Wolz added a comment -

          For me the unit tests initially had a different porpose: to ensure we reproduce the original exactly in java. I compared the eval numbers and ensured they didn't change from the original. Now we are in a different phase. Failing tests should be more seen as a warning and we could adapt the limits. But it is not easy to decide what can betolerated.

          Show
          Dr. Dietmar Wolz added a comment - For me the unit tests initially had a different porpose: to ensure we reproduce the original exactly in java. I compared the eval numbers and ensured they didn't change from the original. Now we are in a different phase. Failing tests should be more seen as a warning and we could adapt the limits. But it is not easy to decide what can betolerated.
          Hide
          Nigel Goodwin added a comment -

          I'm new to this group (but not to optimization) and I have been trying out BOBYQA. Some of the issues which have arisen are:

          • yes, I want the Java source code as part of commons - maths, source code is why I use Apache.
          • there is no trapping of max evaluations exception TooManyEvaluationsException line 833
          • I am struggling with starting points, and this may be an issue with the original algorithm. If a starting poin is given with a value between lower bound and lower bound + trust region, the value is evaluated at lower bound + trust region. But if the starting point is at lower bound, it is evaluated at lower bound. This is important for me because I am using a good starting point from a GA method, so I want that exact starting point to be used in the initial derivative calculations.

          I am using a GA based method, followed by CMAES and/or BOBYQA. I generally find BOBYQA better than CMAES, but not so much if the function does not have smooth derivatives.

          BTW, I find CMAES not good at all with bounds, and it looks like CMAES java code is not using the latest thinking.

          Show
          Nigel Goodwin added a comment - I'm new to this group (but not to optimization) and I have been trying out BOBYQA. Some of the issues which have arisen are: yes, I want the Java source code as part of commons - maths, source code is why I use Apache. there is no trapping of max evaluations exception TooManyEvaluationsException line 833 I am struggling with starting points, and this may be an issue with the original algorithm. If a starting poin is given with a value between lower bound and lower bound + trust region, the value is evaluated at lower bound + trust region. But if the starting point is at lower bound, it is evaluated at lower bound. This is important for me because I am using a good starting point from a GA method, so I want that exact starting point to be used in the initial derivative calculations. I am using a GA based method, followed by CMAES and/or BOBYQA. I generally find BOBYQA better than CMAES, but not so much if the function does not have smooth derivatives. BTW, I find CMAES not good at all with bounds, and it looks like CMAES java code is not using the latest thinking.
          Hide
          Nigel Goodwin added a comment -

          I'm finding my way around this....

          The fortran has:

          W(JSL)=XL(J)-X(J)
          W(JSU)=XU(J)-X(J)
          IF (W(JSL) .GE. -RHOBEG) THEN
          IF (W(JSL) .GE. ZERO) THEN
          X(J)=XL(J)
          W(JSL)=ZERO
          W(JSU)=TEMP
          ELSE
          X(J)=XL(J)+RHOBEG
          W(JSL)=-RHOBEG
          W(JSU)=DMAX1(XU(J)-X(J),RHOBEG)
          END IF
          ELSE IF (W(JSU) .LE. RHOBEG) THEN
          IF (W(JSU) .LE. ZERO) THEN
          X(J)=XU(J)
          W(JSL)=-TEMP
          W(JSU)=ZERO
          ELSE
          X(J)=XU(J)-RHOBEG
          W(JSL)=DMIN1(XL(J)-X(J),-RHOBEG)
          W(JSU)=RHOBEG
          END IF
          END IF

          so the original also has this problem, and resets X(J) in the line

          X(J)=XL(J)+RHOBEG

          It is non-trivial to change this behaviour, i fear.

          Any thoughts? It is quite common to try a local minimization starting from a good points from a stochastic minimization. Or else do some local minimization as a stochastic operator. See for example genoud.

          I will investigate some more.

          Show
          Nigel Goodwin added a comment - I'm finding my way around this.... The fortran has: W(JSL)=XL(J)-X(J) W(JSU)=XU(J)-X(J) IF (W(JSL) .GE. -RHOBEG) THEN IF (W(JSL) .GE. ZERO) THEN X(J)=XL(J) W(JSL)=ZERO W(JSU)=TEMP ELSE X(J)=XL(J)+RHOBEG W(JSL)=-RHOBEG W(JSU)=DMAX1(XU(J)-X(J),RHOBEG) END IF ELSE IF (W(JSU) .LE. RHOBEG) THEN IF (W(JSU) .LE. ZERO) THEN X(J)=XU(J) W(JSL)=-TEMP W(JSU)=ZERO ELSE X(J)=XU(J)-RHOBEG W(JSL)=DMIN1(XL(J)-X(J),-RHOBEG) W(JSU)=RHOBEG END IF END IF so the original also has this problem, and resets X(J) in the line X(J)=XL(J)+RHOBEG It is non-trivial to change this behaviour, i fear. Any thoughts? It is quite common to try a local minimization starting from a good points from a stochastic minimization. Or else do some local minimization as a stochastic operator. See for example genoud. I will investigate some more.
          Hide
          Gilles added a comment -

          there is no trapping of max evaluations exception TooManyEvaluationsException line 833

          Indeed, that's intended behaviour: The user sets the limit, and must catch the exception if he wants to perform something instead of letting the application abort.

          As for changing the "qualitative" behaviour of the algorithm itself, we are not there yet.
          This (Java) code is brand-new and in a state of flux. You are very much welcome to help make it more Java-like: The main task is to figure out how to replace the big "switch" statements (i.e. a set of disguised "goto"s) into a sequence of method calls.

          Show
          Gilles added a comment - there is no trapping of max evaluations exception TooManyEvaluationsException line 833 Indeed, that's intended behaviour: The user sets the limit, and must catch the exception if he wants to perform something instead of letting the application abort. As for changing the "qualitative" behaviour of the algorithm itself, we are not there yet. This (Java) code is brand-new and in a state of flux. You are very much welcome to help make it more Java-like: The main task is to figure out how to replace the big "switch" statements (i.e. a set of disguised "goto"s) into a sequence of method calls.
          Hide
          Nigel Goodwin added a comment -

          I've found what i think is a bug, no doubt arising from the change to zero based indexing.

          Line 541 should read:

          curv = hq.getEntry((j + 1 + (j 1) * (j 1)) / 2 - 1);

          It is clear that, when this loop is encountered with j=0, the original code:

          curv = hq.getEntry((j + j * j) / 2 - 1);

          will cause an array out of bounds exception.

          As far as i can make out, it is indexing the diagonal elements of the second derivative matrix.

          There may be a similar bug in line 1777, but I think I wexercise this code and have not had an exception. Worth checking.

          Show
          Nigel Goodwin added a comment - I've found what i think is a bug, no doubt arising from the change to zero based indexing. Line 541 should read: curv = hq.getEntry((j + 1 + (j 1) * (j 1)) / 2 - 1); It is clear that, when this loop is encountered with j=0, the original code: curv = hq.getEntry((j + j * j) / 2 - 1); will cause an array out of bounds exception. As far as i can make out, it is indexing the diagonal elements of the second derivative matrix. There may be a similar bug in line 1777, but I think I wexercise this code and have not had an exception. Worth checking.
          Hide
          Gilles added a comment -

          You are right that there is a bug, but comparing with the original Fortran, I think that the line should rather read:

          curv = hq.getEntry((j + j * j) / 2);
          

          Unfortunately, this statement is in a code path that is not explored by the current unit tests.

          By the way, when writing code excerpts here, don't forget to enclose them within the appropriate tag:

          {code}
            // Code goes here
          {code}
          

          Otherwise it can get scrambled...

          For the other line, I have no idea.

          Show
          Gilles added a comment - You are right that there is a bug, but comparing with the original Fortran, I think that the line should rather read: curv = hq.getEntry((j + j * j) / 2); Unfortunately, this statement is in a code path that is not explored by the current unit tests. By the way, when writing code excerpts here, don't forget to enclose them within the appropriate tag: {code} // Code goes here {code} Otherwise it can get scrambled... For the other line, I have no idea.
          Hide
          Gilles added a comment -

          Fixed first reported bug in revision 1165656. Thanks.

          Show
          Gilles added a comment - Fixed first reported bug in revision 1165656. Thanks.
          Hide
          Nigel Goodwin added a comment - - edited

          let me think.

          Assume hq is a lower triangular matrix, stored in order of columns.

          So, using a 1-index, we have

          hq(1) is mat(1,1)
          hq(3) is mat(2,2)
          hq(6) is mat(3,3) etc.

          So in 1-index, the diagonal element of the j'th row is

          h((j + j*j)/2)
          

          as in the Fortran

          So in 0-index, we have

          hq(0) is mat(0,0)
          hq(2) is mat(1,1)
          hq(5) is mat(2,2)

          So the diagonal element is given by

          ((j+1) + (j + 1) * ( j+ 1))/2 - 1

          and I think

          (j + j*j)/2

          is incorrect.

          I'm pretty sure about this - I do a lot of work with triangular matrices.

          Try another test - look at the last diagonal element - when j is n - 1, then the index is (n + n*n)/2 -1 which is n*(n+1)/2 - 1, which is compatible with the dimension when declared, n*np/2.

          Show
          Nigel Goodwin added a comment - - edited let me think. Assume hq is a lower triangular matrix, stored in order of columns. So, using a 1-index, we have hq(1) is mat(1,1) hq(3) is mat(2,2) hq(6) is mat(3,3) etc. So in 1-index, the diagonal element of the j'th row is h((j + j*j)/2) as in the Fortran So in 0-index, we have hq(0) is mat(0,0) hq(2) is mat(1,1) hq(5) is mat(2,2) So the diagonal element is given by ((j+1) + (j + 1) * ( j+ 1))/2 - 1 and I think (j + j*j)/2 is incorrect. I'm pretty sure about this - I do a lot of work with triangular matrices. Try another test - look at the last diagonal element - when j is n - 1, then the index is (n + n*n)/2 -1 which is n*(n+1)/2 - 1, which is compatible with the dimension when declared, n*np/2.
          Hide
          Nigel Goodwin added a comment - - edited

          forgot to close the cod tag, can't see where to edit my comment...

          What I meant was

          
          h((j + j*j)/2)
          
          
          Show
          Nigel Goodwin added a comment - - edited forgot to close the cod tag, can't see where to edit my comment... What I meant was h((j + j*j)/2)
          Hide
          Sebb added a comment -

          Hover over the comment you want to edit; there will be an edit icon in the top rhs of the grey box.

          Show
          Sebb added a comment - Hover over the comment you want to edit; there will be an edit icon in the top rhs of the grey box.
          Hide
          Nigel Goodwin added a comment -

          I've checked line 1777 carefully and compared the logic with the Fortran, and it looks OK. I haven't checked line 1812.

          Note - it is well known that correcting a bug is likely to create a new bug, and if there is a bug in one area of code, there is likely to be a bug in similar areas. Also, incorrect indexing of the second derivative matrix may not cause an exception, but may simply make the algorithm perform less well. So it's worth a few more eyes to look at indexing of hq and check it.

          Show
          Nigel Goodwin added a comment - I've checked line 1777 carefully and compared the logic with the Fortran, and it looks OK. I haven't checked line 1812. Note - it is well known that correcting a bug is likely to create a new bug, and if there is a bug in one area of code, there is likely to be a bug in similar areas. Also, incorrect indexing of the second derivative matrix may not cause an exception, but may simply make the algorithm perform less well. So it's worth a few more eyes to look at indexing of hq and check it.
          Hide
          Nigel Goodwin added a comment -

          Let me be clear - revision 1165656 is incorrect, IMHO.

          Show
          Nigel Goodwin added a comment - Let me be clear - revision 1165656 is incorrect, IMHO.
          Hide
          Gilles added a comment -

          OK, I trust you

          There were many other cases of matrices disguised as vectors. This is the last one, I think, which I left for the end because of the tricky indexing.
          So, instead of the current line 236:

          final ArrayRealVector hq = new ArrayRealVector(n * np / 2);
          

          I'd like to have:

          final Array2DRowRealMatrix hq = new Array2DRowRealMatrix(n, n);
          

          Would you like to provide the patch with the necessary changes (1d to 2d indexing)?

          Show
          Gilles added a comment - OK, I trust you There were many other cases of matrices disguised as vectors. This is the last one, I think, which I left for the end because of the tricky indexing. So, instead of the current line 236: final ArrayRealVector hq = new ArrayRealVector(n * np / 2); I'd like to have: final Array2DRowRealMatrix hq = new Array2DRowRealMatrix(n, n); Would you like to provide the patch with the necessary changes (1d to 2d indexing)?
          Hide
          Nigel Goodwin added a comment - - edited

          No I would not!

          I do a lot of work with positive definite symmetric matrices, and I ALWAYS store and manipulate them as upper (or lower) triangular matrices stored as a vector or 1D array. This is the best way to do it, and nothing would motivate me to make a retrograde step and store symmetric matrices as a matrix or 2D array (unless maybe underlying storages was a 1 D array, and even then I would have grave misgivings).

          I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage. Plus i do a bit of loop unrolling etc etc but that is another story.

          I first looked at code provided by Prof Powell back in the early 1980's, his rank one (or was it rank two) cholesky factor update, and no doubt this has influenced my tendency to always store upper (or lower) triangular matrices as 1D arrays.

          I even hate it when there is a special class for symmetric matrices, and it has a last step which puts the lower triangular elements equal to the upper triangular elements. No, don't do it, do it properly!

          In most of my application code, cholesky factorisation and back substitution are the bottlenecks. I do tens of thousands of back substitutions.

          Now, if there was a symmetric matrix class with an underlying storage of a 1D array, and some suitable get and set functions, you might persuade me, so long as the get and set were just utilities to be used in non critical parts of the code.

          Store a symmetric matrix as a 2d array? No, never, no!

          ps soory about the capitals. But you may note I feel strongly. It is why all Java libraries I have looked at (including apache commons maths) suck for symmetric matrices, and i had to write my own.

          It is late, and I may be wrong, but I think LAPACK symmetric code (or is it BLAS?) stores symmetric matrics as 1D triangular arrays.

          Show
          Nigel Goodwin added a comment - - edited No I would not! I do a lot of work with positive definite symmetric matrices, and I ALWAYS store and manipulate them as upper (or lower) triangular matrices stored as a vector or 1D array. This is the best way to do it, and nothing would motivate me to make a retrograde step and store symmetric matrices as a matrix or 2D array (unless maybe underlying storages was a 1 D array, and even then I would have grave misgivings). I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage. Plus i do a bit of loop unrolling etc etc but that is another story. I first looked at code provided by Prof Powell back in the early 1980's, his rank one (or was it rank two) cholesky factor update, and no doubt this has influenced my tendency to always store upper (or lower) triangular matrices as 1D arrays. I even hate it when there is a special class for symmetric matrices, and it has a last step which puts the lower triangular elements equal to the upper triangular elements. No, don't do it, do it properly! In most of my application code, cholesky factorisation and back substitution are the bottlenecks. I do tens of thousands of back substitutions. Now, if there was a symmetric matrix class with an underlying storage of a 1D array, and some suitable get and set functions, you might persuade me, so long as the get and set were just utilities to be used in non critical parts of the code. Store a symmetric matrix as a 2d array? No, never, no! ps soory about the capitals. But you may note I feel strongly. It is why all Java libraries I have looked at (including apache commons maths) suck for symmetric matrices, and i had to write my own. It is late, and I may be wrong, but I think LAPACK symmetric code (or is it BLAS?) stores symmetric matrics as 1D triangular arrays.
          Hide
          Nigel Goodwin added a comment -

          Sorrt for saying the same thing 3 times! I'm still finding my way around here, although I intend to post very rarely - I am a user of Apache rather than a contributor, unless there is something critical to which I can provide some help.

          Show
          Nigel Goodwin added a comment - Sorrt for saying the same thing 3 times! I'm still finding my way around here, although I intend to post very rarely - I am a user of Apache rather than a contributor, unless there is something critical to which I can provide some help.
          Hide
          Gilles added a comment -

          I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage.

          The point is that this issue deals with an optimization algorithm; not how to properly store a symmetric matrix. Algorithms that use matrices should not have to re-implement matrix operations.

          A matrix represented by a "...Matrix" object certainly will help in understanding what the BOBYQA code does and how to further adapt it to the Java programming style.

          Your knowledge might be a useful contribution (on the "dev" ML) to the discussion on how to improve the design and efficiency of the various matrix classes. Bringing some of your concrete arguments could lead to the creation of an efficient "SymmetricMatrix" to be used in this algorithm.

          Show
          Gilles added a comment - I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage. The point is that this issue deals with an optimization algorithm; not how to properly store a symmetric matrix. Algorithms that use matrices should not have to re-implement matrix operations. A matrix represented by a "...Matrix" object certainly will help in understanding what the BOBYQA code does and how to further adapt it to the Java programming style. Your knowledge might be a useful contribution (on the "dev" ML) to the discussion on how to improve the design and efficiency of the various matrix classes. Bringing some of your concrete arguments could lead to the creation of an efficient "SymmetricMatrix" to be used in this algorithm.
          Hide
          Nigel Goodwin added a comment -

          Various responses, then i'm out!

          a) you have no idea how grateful I am to you and others for the whole set of apache commons maths code

          b) I use the code for end user applications

          c) if it ain't broke, don;t fix it - so my main concern at present is finding any bugs

          d) if you want to redesign Powell's code, using other existing apache code, good luck! It will not be easy, and will introduce bugs, and will take some time to settle down. There may be better ways to use your skills and effort.

          e) surely the main priority is to create a java BOBYQA which is robust and is used. Adapting it to conform to some pre-existing or debatable style is surely secondary?

          f) yes, a nice symmetric matrix class with some suitable algorithms would be nice.

          g) I am new to apache, I look around for a lot of mathematical library java code, there is a lot out there, I use for example JMetal for multi objective optimisation. I have looked at jama. I have tried out JCuda. Etc etc. What makes me choose a particular package? Flexibility, simplicity, best of breed algorithms, easy to understand, easy to extend. Now, coming back to BOBYQA, it is complex enough that I will probably never attempt to understand or modify it. So I will just take it as a black box and hope it works. No doubt some other algorithm will appear in a year or two, and some good person like your self will also modify it. BOBYQA is complex, and I leave it to the experts, and will leave well alone unless I have studied it full time for 6 months or more and know every detail. Not going to happen!

          h) I have already moved on and thinking about hybrid quasi newton or radial basis function/genetic algorithm approaches, but will leave further study to next year.

          i) Oh yes, the one thing I could not find any Java code for was a BesselK function of second (or third?) order, used to define a Matern function. That could be a useful extension to apache special functions, more useful work than 'prettifying' a complex algorithm like BOBYQA.

          j) Hope you like robust debates! But I fear I am finished here, got to move on to some more application oriented development.

          k) look me up on LinkedIn

          Show
          Nigel Goodwin added a comment - Various responses, then i'm out! a) you have no idea how grateful I am to you and others for the whole set of apache commons maths code b) I use the code for end user applications c) if it ain't broke, don;t fix it - so my main concern at present is finding any bugs d) if you want to redesign Powell's code, using other existing apache code, good luck! It will not be easy, and will introduce bugs, and will take some time to settle down. There may be better ways to use your skills and effort. e) surely the main priority is to create a java BOBYQA which is robust and is used. Adapting it to conform to some pre-existing or debatable style is surely secondary? f) yes, a nice symmetric matrix class with some suitable algorithms would be nice. g) I am new to apache, I look around for a lot of mathematical library java code, there is a lot out there, I use for example JMetal for multi objective optimisation. I have looked at jama. I have tried out JCuda. Etc etc. What makes me choose a particular package? Flexibility, simplicity, best of breed algorithms, easy to understand, easy to extend. Now, coming back to BOBYQA, it is complex enough that I will probably never attempt to understand or modify it. So I will just take it as a black box and hope it works. No doubt some other algorithm will appear in a year or two, and some good person like your self will also modify it. BOBYQA is complex, and I leave it to the experts, and will leave well alone unless I have studied it full time for 6 months or more and know every detail. Not going to happen! h) I have already moved on and thinking about hybrid quasi newton or radial basis function/genetic algorithm approaches, but will leave further study to next year. i) Oh yes, the one thing I could not find any Java code for was a BesselK function of second (or third?) order, used to define a Matern function. That could be a useful extension to apache special functions, more useful work than 'prettifying' a complex algorithm like BOBYQA. j) Hope you like robust debates! But I fear I am finished here, got to move on to some more application oriented development. k) look me up on LinkedIn
          Hide
          Dr. Dietmar Wolz added a comment -

          I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage.

          A previous version of the code (attached as version 0.4) directly uses 1D arrays instead of Matrices and is definitely faster. But here in CM the focus is to produce
          code which is easily understandable. This way we have a chance to understand and improve the algorithm itself. If performance of the optimization algo is important
          (which in many applications is not the case because "interesting" cost functions are normally expensive to evaluate) then you could call the C-version of
          BOBYQA via JNI which gives you a performance boost of factor 6-10. This is possible still using the CM interfaces and when the cost function
          is implemented in Java.

          Show
          Dr. Dietmar Wolz added a comment - I have written a lot of symmetric matrix code, and it ALWAYS outperforms any Java library I can find, simply because I use the proper 1D array form of storage. A previous version of the code (attached as version 0.4) directly uses 1D arrays instead of Matrices and is definitely faster. But here in CM the focus is to produce code which is easily understandable. This way we have a chance to understand and improve the algorithm itself. If performance of the optimization algo is important (which in many applications is not the case because "interesting" cost functions are normally expensive to evaluate) then you could call the C-version of BOBYQA via JNI which gives you a performance boost of factor 6-10. This is possible still using the CM interfaces and when the cost function is implemented in Java.
          Hide
          Gilles added a comment -

          Hi Dietmar.

          Could you please have a look at revision 1185917?
          I'm worried that this seems like a sensitivity to the limited precision of the floating point representation of numbers.
          I recall that you mentioned you had detected this kind of behaviour.
          I think that, if this is intended behaviour, its rationale should be documented; otherwise, one can only wonder why a logically equivalent statement might break the algorithm.
          However, I'd prefer to think that the unit tests are too stringent, i.e. that the mathematical algorithm is actually robust and that the the tolerance in the failing unit tests should be adapted to take into account the limited precision. Or probably even better: that it would be possible to detect that the numerical procedure is in trouble (due to a too low tolerance) and throw an appropriate exception instead of fail after exhausting the evaluations counter (there was a similar problem in issue MATH-631).

          Equally worrisome is the fact that not all the branches are covered by the "baseline" unit tests, so that my inserting bugs in statements there will go unnoticed!
          Do you have ideas on how to increase the coverage (in a controlled way)?

          Show
          Gilles added a comment - Hi Dietmar. Could you please have a look at revision 1185917? I'm worried that this seems like a sensitivity to the limited precision of the floating point representation of numbers. I recall that you mentioned you had detected this kind of behaviour. I think that, if this is intended behaviour, its rationale should be documented; otherwise, one can only wonder why a logically equivalent statement might break the algorithm. However, I'd prefer to think that the unit tests are too stringent, i.e. that the mathematical algorithm is actually robust and that the the tolerance in the failing unit tests should be adapted to take into account the limited precision. Or probably even better: that it would be possible to detect that the numerical procedure is in trouble (due to a too low tolerance) and throw an appropriate exception instead of fail after exhausting the evaluations counter (there was a similar problem in issue MATH-631 ). Equally worrisome is the fact that not all the branches are covered by the "baseline" unit tests, so that my inserting bugs in statements there will go unnoticed! Do you have ideas on how to increase the coverage (in a controlled way)?
          Hide
          Dr. Dietmar Wolz added a comment -

          Hi Gilles,
          what do you mean with "logically equivalent"? Computers cannot perform exact arithmetic,
          there is always a rounding error, which can accumulate differently when you choose different
          (but mathematically equivalent) algorithms. See for instance http://en.wikipedia.org/wiki/Numerical_stability.
          The existing unit tests were initially designed by me to detect fortran/java conversion errors.
          After the initial conversion phase - when we try to "beautify" the algorithm by using CM Matrix
          functionality we can loosen the accuracy requirements - which seems ok, as long as we maintain
          the applicability to real problems. If we run into real problems, perhaps the underlying CM matrix
          function implementation needs a redesign to increase accuracy, other algos would then also profit from
          an improvement here.

          The unit test problems are quite artificial, which could explain
          their limited coverage. I see two options to improve coverage:
          1) Understand why some parts of the code is not covered and design unit tests using this knowledge.
          This should be done by researchers in the area of optimization algorithms - or by the original
          author M. Powell.
          2) Collect practical applications of the algo by users of CM. Check their coverage and add a unit test
          if a newly covered code section is detected.

          Show
          Dr. Dietmar Wolz added a comment - Hi Gilles, what do you mean with "logically equivalent"? Computers cannot perform exact arithmetic, there is always a rounding error, which can accumulate differently when you choose different (but mathematically equivalent) algorithms. See for instance http://en.wikipedia.org/wiki/Numerical_stability . The existing unit tests were initially designed by me to detect fortran/java conversion errors. After the initial conversion phase - when we try to "beautify" the algorithm by using CM Matrix functionality we can loosen the accuracy requirements - which seems ok, as long as we maintain the applicability to real problems. If we run into real problems, perhaps the underlying CM matrix function implementation needs a redesign to increase accuracy, other algos would then also profit from an improvement here. The unit test problems are quite artificial, which could explain their limited coverage. I see two options to improve coverage: 1) Understand why some parts of the code is not covered and design unit tests using this knowledge. This should be done by researchers in the area of optimization algorithms - or by the original author M. Powell. 2) Collect practical applications of the algo by users of CM. Check their coverage and add a unit test if a newly covered code section is detected.
          Hide
          Gilles added a comment -

          "logically equivalent" == "mathematically equivalent"

          What I'm asking is whether the test failures are meaningful.
          I understand that one cannot expect the numbers to be the same all through the last decimal places when reordering some operations. When such a case arises, we should probably increase the tolerances so that the test passes.
          But I'm wondering whether it is normal that reordering should lead to an increase of the number of evaluations.

          I think that the accuracy thresholds should take rounding into account, in the sense that the results of two logically/mathematically equivalent computations should be considered equal (unless there is an intrinsic feature of the algorithm causing "really" different results, in which case a comment should make it clear).
          In this instance,

           a + 2 * dx
          

          and

           a + dx + dx
          

          give different results.
          One explanation could be that "a + dx" is still "a". But IMO, that means that the algorithm is fragile: An addition was meant but nothing has actually happened. Hence, I'd tend to say that further computations are doubtful...
          That's what I mean by "detect that the numerical procedure is in trouble"; the input data (e.g. tolerance value) leads it to ineffectiveness, which should be detected and reported as such.

          In fact, I thought that the unit tests came from an original test suite used by BOBYQA's author! Does such a suite exist?
          Alternatively, (related to your point 2), we can try and set up our own suite using "standard" problems; there has been some attempt in this sense with the "BatteryNISTTest" class introduced recently. This class already led to exercising a code path not covered by the existing tests; however, I was hoping that someone would be more systematic in the selection of a test suite of "well-known" (by the optimization community) problems.
          Of course, this is going back to the discussion we had a few weeks ago: Do we wait for a hypothetical expert, or do we do something now?

          Show
          Gilles added a comment - "logically equivalent" == "mathematically equivalent" What I'm asking is whether the test failures are meaningful . I understand that one cannot expect the numbers to be the same all through the last decimal places when reordering some operations. When such a case arises, we should probably increase the tolerances so that the test passes. But I'm wondering whether it is normal that reordering should lead to an increase of the number of evaluations. I think that the accuracy thresholds should take rounding into account, in the sense that the results of two logically/mathematically equivalent computations should be considered equal (unless there is an intrinsic feature of the algorithm causing "really" different results, in which case a comment should make it clear). In this instance, a + 2 * dx and a + dx + dx give different results. One explanation could be that "a + dx" is still "a". But IMO, that means that the algorithm is fragile: An addition was meant but nothing has actually happened. Hence, I'd tend to say that further computations are doubtful... That's what I mean by "detect that the numerical procedure is in trouble"; the input data (e.g. tolerance value) leads it to ineffectiveness, which should be detected and reported as such. In fact, I thought that the unit tests came from an original test suite used by BOBYQA's author! Does such a suite exist? Alternatively, (related to your point 2), we can try and set up our own suite using "standard" problems; there has been some attempt in this sense with the "BatteryNISTTest" class introduced recently. This class already led to exercising a code path not covered by the existing tests; however, I was hoping that someone would be more systematic in the selection of a test suite of "well-known" (by the optimization community) problems. Of course, this is going back to the discussion we had a few weeks ago: Do we wait for a hypothetical expert, or do we do something now?
          Hide
          Dr. Dietmar Wolz added a comment -

          > from an original test suite used by BOBYQA's author!
          No, they are derived from the CMA-ES testsuite which itself was derived from samples from Nikolaus Hansen (the author of CMA-ES). Nikolaus Hansen did quite extensive testing because of the stochastic nature of CMA-ES. I have no objections to loosen the accuracy thresholds and to add "systematic" tests. Regarding the
          "hypothetical expert" we perhaps have to actively search for one - as far as I remember there were some guys in south america extending BOBYQA at the FORTRAN level. The problem is to find someone familiar with unit tests and the concepts behind BOBYQA. Perhaps we can raise interest at some university to work with us?

          Show
          Dr. Dietmar Wolz added a comment - > from an original test suite used by BOBYQA's author! No, they are derived from the CMA-ES testsuite which itself was derived from samples from Nikolaus Hansen (the author of CMA-ES). Nikolaus Hansen did quite extensive testing because of the stochastic nature of CMA-ES. I have no objections to loosen the accuracy thresholds and to add "systematic" tests. Regarding the "hypothetical expert" we perhaps have to actively search for one - as far as I remember there were some guys in south america extending BOBYQA at the FORTRAN level. The problem is to find someone familiar with unit tests and the concepts behind BOBYQA. Perhaps we can raise interest at some university to work with us?
          Hide
          Gilles added a comment -

          Anybody willing to help is welcome!
          Maybe Mike Powell could provide some unit tests .
          Can you contact the people in South America?

          Do you have an idea on how to actively search for an expert?
          I assume that telling that we try to rewrite an algorithm in Java will not look very exciting...

          Show
          Gilles added a comment - Anybody willing to help is welcome! Maybe Mike Powell could provide some unit tests . Can you contact the people in South America? Do you have an idea on how to actively search for an expert? I assume that telling that we try to rewrite an algorithm in Java will not look very exciting...
          Hide
          Dr. Dietmar Wolz added a comment -

          Wrote a letter to Prof. Pilotta, one of the authors of http://www.scielo.br/pdf/cam/v30n1/09.pdf.
          Commented on http://forum.openopt.org/viewtopic.php?pid=1422#p1422 referring to MATH-621.

          Show
          Dr. Dietmar Wolz added a comment - Wrote a letter to Prof. Pilotta, one of the authors of http://www.scielo.br/pdf/cam/v30n1/09.pdf . Commented on http://forum.openopt.org/viewtopic.php?pid=1422#p1422 referring to MATH-621 .
          Hide
          Gilles added a comment -

          No further work will be achieved before release 3.0: Postponing.

          Show
          Gilles added a comment - No further work will be achieved before release 3.0: Postponing.
          Hide
          Thomas Neidhart added a comment -

          Findbugs reports a missing break statement in the BOBYQAOptimizer at line 1880. In case it is intended, it should be documents in the code.

          Show
          Thomas Neidhart added a comment - Findbugs reports a missing break statement in the BOBYQAOptimizer at line 1880. In case it is intended, it should be documents in the code.
          Hide
          Gilles added a comment - - edited

          As you can figure out from reading this issue, there is a lot of work to make this code readable and maintainable.
          The "missing break" is most probably intended, but the code should not be written that way nevertheless; so I would refrain adding comments on constructs that should anyways be refactored. I much prefer that FindBugs keeps reminding us that this class is not representative of Commons Math coding standards.

          Show
          Gilles added a comment - - edited As you can figure out from reading this issue, there is a lot of work to make this code readable and maintainable. The "missing break" is most probably intended, but the code should not be written that way nevertheless; so I would refrain adding comments on constructs that should anyways be refactored. I much prefer that FindBugs keeps reminding us that this class is not representative of Commons Math coding standards.
          Hide
          Thomas Neidhart added a comment -

          I am perfectly aware of this, I just wanted to put a pointer here so it's not missed.

          Show
          Thomas Neidhart added a comment - I am perfectly aware of this, I just wanted to put a pointer here so it's not missed.

            People

            • Assignee:
              Unassigned
              Reporter:
              Dr. Dietmar Wolz
            • Votes:
              0 Vote for this issue
              Watchers:
              2 Start watching this issue

              Dates

              • Created:
                Updated:

                Development