{"action":"created","issue":{"url":"https://api.github.com/repos/scipy/scipy/issues/5101","labels_url":"https://api.github.com/repos/scipy/scipy/issues/5101/labels{/name}","comments_url":"https://api.github.com/repos/scipy/scipy/issues/5101/comments","events_url":"https://api.github.com/repos/scipy/scipy/issues/5101/events","html_url":"https://github.com/scipy/scipy/issues/5101","id":98384142,"number":5101,"title":"Bug in Differential Evolution (Array-length with floating point rounding error)","user":{"login":"Hvass-Labs","id":13588114,"avatar_url":"https://avatars.githubusercontent.com/u/13588114?v=3","gravatar_id":"","url":"https://api.github.com/users/Hvass-Labs","html_url":"https://github.com/Hvass-Labs","followers_url":"https://api.github.com/users/Hvass-Labs/followers","following_url":"https://api.github.com/users/Hvass-Labs/following{/other_user}","gists_url":"https://api.github.com/users/Hvass-Labs/gists{/gist_id}","starred_url":"https://api.github.com/users/Hvass-Labs/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/Hvass-Labs/subscriptions","organizations_url":"https://api.github.com/users/Hvass-Labs/orgs","repos_url":"https://api.github.com/users/Hvass-Labs/repos","events_url":"https://api.github.com/users/Hvass-Labs/events{/privacy}","received_events_url":"https://api.github.com/users/Hvass-Labs/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":7,"created_at":"2015-07-31T12:47:59Z","updated_at":"2015-08-01T01:44:21Z","closed_at":null,"body":"Hello,\r\n\r\nThanks for making scipy!\r\n\r\nThere is a bug in scipy 0.15.1 in the file _differentialevolution.py\r\n\r\nThe bug is in the function init_population_lhs(self) in this line:\r\n\r\n    rdrange += np.atleast_2d(np.arange(0., 1., segsize)).T\r\n\r\nand specifically in this part:\r\n\r\n    np.arange(0., 1., segsize)\r\n\r\nThe problem is that floating-point rounding errors can sometimes create arrays of the wrong length. For example, if the dimensionality of the search-space is 29 and the population size is 35 then the bug is triggered.\r\n\r\nI have fixed the source-code and attached it below. I have also revised and cleaned up the code a bit. The commenting was also confusing and misleading. I'm not sure mine is much better because I haven't studied exactly what the function is supposed to do (again, the commenting was lacking).\r\n\r\nThere are still some things in the function that should be improved, but this would require changes to e.g. the class init as well. I have mentioned these issues in the comments.\r\n\r\nI have also noticed other strange things in other parts of scipy's DE source-code so I can imagine there are more bugs.\r\n\r\nDE is a very useful optimizer and many people probably appreciate that it is now included in scipy, but the implementation needs to be improved. I would suggest fresh and critical eyes from the scipy team go over the source-code to simplify, test and document it properly. I would have liked to help you with this but I'm new to Python.\r\n\r\nI have attached my own DE implementation for Matlab as a reference so you can see how simple it can be done.\r\n\r\nFurthermore, please use the typical DE terminology. You have renamed the DE control parameters to 'mutation' and 'recombination'. Please use their typical names in the research literature 'CR' and 'F' so a practitioner can easily see which is which, or use the full names such as crossover_probability (instead of CR) and differential_weight (instead of F).\r\n\r\nNot everything is bad with scipy's DE implementation and there are some pretty nice features! I particularly like the convergence detection system and that L-BFGS-B optimization can be continued from DE's best found solution - assuming it isn't buggy! :-)\r\n\r\nThanks again for making scipy!\r\n\r\nCheers, Magnus.\r\n\r\n\r\n###########################\r\n# _differentialevolution.py\r\n# OLD Version:\r\n\r\n    def init_population_lhs(self):\r\n        \"\"\"\r\n        Initializes the population with Latin Hypercube Sampling\r\n        Latin Hypercube Sampling ensures that the sampling of parameter space\r\n        is maximised.\r\n        \"\"\"\r\n        samples = np.size(self.population, 0)\r\n        N = np.size(self.population, 1)\r\n        rng = self.random_number_generator\r\n\r\n        # Generate the intervals\r\n        segsize = 1.0 / samples\r\n\r\n        # Fill points uniformly in each interval\r\n        rdrange = rng.rand(samples, N) * segsize\r\n        rdrange += np.atleast_2d(np.arange(0., 1., segsize)).T\r\n\r\n        # Make the random pairings\r\n        self.population = np.zeros_like(rdrange)\r\n\r\n        for j in range(N):\r\n            order = rng.permutation(range(samples))\r\n            self.population[:, j] = rdrange[order, j]\r\n\r\n\r\n###########################\r\n# _differentialevolution.py\r\n# NEW version:\r\n\r\n    def init_population_lhs(self):\r\n        \"\"\"\r\n        Initializes the population with Latin Hypercube Sampling\r\n        Latin Hypercube Sampling ensures that the sampling of parameter space\r\n        is maximised.\r\n\r\n        NOTE: This should describe with a math formula what is being calculated!\r\n        \"\"\"\r\n\r\n        # Get number of samples for each candidate solution.\r\n        # This is messy and should be self.num_samples instead.\r\n        num_samples = np.size(self.population, 0)\r\n\r\n        # Get search-space dimensionality.\r\n        # This is messy and should be self.dim instead.\r\n        dim = np.size(self.population, 1)\r\n\r\n        # Get RNG.\r\n        rng = self.random_number_generator\r\n\r\n        # Generate random samples.\r\n        samples = rng.rand(num_samples, dim)\r\n\r\n        # Scale used for random samples.\r\n        scale = 1.0 / num_samples\r\n\r\n        # Scale the samples.\r\n        samples *= scale\r\n\r\n        # Offset each individual sample so they are in different ranges.\r\n        # This spreads out the samples evenly.\r\n        offset = scale * np.arange(1, num_samples + 1)\r\n        samples += np.atleast_2d(offset).T\r\n\r\n        # Allocate array for population of candidate solutions and initialize to zero.\r\n        self.population = np.zeros_like(samples)\r\n\r\n        # Initialize population of candidate solutions by permutation of the random samples.\r\n        for j in range(dim):\r\n            index = rng.permutation(range(num_samples))\r\n            self.population[:, j] = samples[index, j]\r\n\r\n###########################\r\n\r\n% ------------------------------------------------------\r\n% SwarmOps - Heuristic optimization for Matlab\r\n% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.\r\n% Please see the file license.txt for license details.\r\n% SwarmOps on the internet: http://www.Hvass-Labs.org/\r\n% ------------------------------------------------------\r\n\r\n% Differential Evolution (DE) by Storner and Price (1).\r\n% DE is an optimization method that does not use the\r\n% gradient of the problem. This is the DE/rand/1/bin variant.\r\n% Literature references:\r\n% (1) R. Storn and K. Price. Differential evolution - a simple\r\n%     and efficient heuristic for global optimization over\r\n%     continuous spaces. Journal of Global Optimization,\r\n%     11:341-359, 1997.\r\n% Parameters:\r\n%     problem; name or handle of optimization problem, e.g. @myproblem.\r\n%     data; data-struct, see e.g. the file myproblemdata.m\r\n%     parameters; behavioural parameters for optimizer,\r\n%                 see file deparameters.m\r\n% Returns:\r\n%     bestX; best found position in the search-space.\r\n%     bestFitness; fitness of bestX.\r\n%     evaluations; number of fitness evaluations performed.\r\n\r\nfunction [bestX, bestFitness, evaluations] = de(problem, data, parameters)\r\n\r\n    % Copy data contents to local variables for convenience.\r\n    n = data.Dim;\r\n    acceptableFitness = data.AcceptableFitness;\r\n    maxEvaluations = data.MaxEvaluations;\r\n    lowerBound = data.LowerBound;\r\n    upperBound = data.UpperBound;\r\n\r\n    % Behavioural parameters for this optimizer.\r\n    np = parameters(1);       % Population size\r\n    cr = parameters(2);       % Crossover probability.\r\n    f = parameters(3);        % Differential weight.\r\n\r\n    % Initialize population.\r\n    x = initpopulation(np, n, data.LowerInit, data.UpperInit);\r\n\r\n    % Compute fitness for all agents in population.\r\n    fitness = zeros(1, np); % Pre-allocate for efficiency.\r\n    for i=1:np\r\n        fitness(i) = feval(problem, x(i,:), data);\r\n    end\r\n\r\n    % Determine fitness and index of best agent.\r\n    [bestFitness, bestIndex] = min(fitness);\r\n\r\n    % Perform optimization iterations until acceptable fitness\r\n    % is achieved or the maximum number of fitness evaluations\r\n    % has been performed.\r\n    evaluations = np; % Fitness evaluations above count as iterations.\r\n    while (evaluations < maxEvaluations) && (bestFitness > acceptableFitness)\r\n\r\n        % Pick random and distinct agents from population.\r\n        indices = randperm(np);\r\n        i = indices(1);        % Agent to update.\r\n        a = indices(2);        % Other agent a.\r\n        b = indices(3);        % Other agent b.\r\n        c = indices(4);        % Other agent c.\r\n\r\n        % Pick random index.\r\n        R = randindex(n);\r\n\r\n        % Copy old agent as a basis for the new agent.\r\n        y = x(i,:);\r\n\r\n        % Compute new agent using DE formula.\r\n        for j=1:n\r\n            if (j==R) || (rand(1,1) < cr)\r\n                y(j) = x(a,j) + f * (x(b,j) - x(c,j));\r\n            end\r\n        end\r\n\r\n        % Bound position to search-space.\r\n        y = bound(y, lowerBound, upperBound);\r\n\r\n        % Compute fitness.\r\n        newFitness = feval(problem, y, data);\r\n\r\n        % Update agent if fitness improvement.\r\n        if newFitness < fitness(i)\r\n            % Update fitness.\r\n            fitness(i) = newFitness;\r\n\r\n            % Update agent position.\r\n            x(i,:) = y;\r\n\r\n            % Update population's best, if improvement.\r\n            if newFitness < bestFitness\r\n                bestFitness = newFitness;\r\n                bestIndex = i;\r\n            end\r\n        end\r\n\r\n        % Increment counter.\r\n        evaluations = evaluations + 1;\r\n    end\r\n\r\n    % Set best found solution (bestFitness already set).\r\n    bestX = x(bestIndex,:);\r\nend\r\n\r\n% ------------------------------------------------------\r\n"},"comment":{"url":"https://api.github.com/repos/scipy/scipy/issues/comments/126848469","html_url":"https://github.com/scipy/scipy/issues/5101#issuecomment-126848469","issue_url":"https://api.github.com/repos/scipy/scipy/issues/5101","id":126848469,"user":{"login":"andyfaff","id":702934,"avatar_url":"https://avatars.githubusercontent.com/u/702934?v=3","gravatar_id":"","url":"https://api.github.com/users/andyfaff","html_url":"https://github.com/andyfaff","followers_url":"https://api.github.com/users/andyfaff/followers","following_url":"https://api.github.com/users/andyfaff/following{/other_user}","gists_url":"https://api.github.com/users/andyfaff/gists{/gist_id}","starred_url":"https://api.github.com/users/andyfaff/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/andyfaff/subscriptions","organizations_url":"https://api.github.com/users/andyfaff/orgs","repos_url":"https://api.github.com/users/andyfaff/repos","events_url":"https://api.github.com/users/andyfaff/events{/privacy}","received_events_url":"https://api.github.com/users/andyfaff/received_events","type":"User","site_admin":false},"created_at":"2015-08-01T01:44:21Z","updated_at":"2015-08-01T01:44:21Z","body":"I had a look at your matlab code and the general algorithm in the two is very similar for rand/1/bin.  However, I did have one comment.  Normally the entire population is iterated through in sequence. The original Storn and Price paper uses sequential iteration of the population, as does the Wormington paper I referred to above. The *algorithm* section of the Wikipedia page does as well.  However, I think the *sample code* section visits the population in a random manner.\r\n\r\n I'm not sure that your code is performing sequential iteration of the population. You have:\r\n\r\n    while (evaluations < maxEvaluations) && (bestFitness > acceptableFitness)\r\n        % Pick random and distinct agents from population.\r\n        indices = randperm(np);\r\n        i = indices(1);        % Agent to update.\r\n        a = indices(2);        % Other agent a.\r\n        b = indices(3);        % Other agent b.\r\n        c = indices(4);        % Other agent c.\r\n\r\n`i` should be calculated as something like `i = evaluations % np`, which would correspond to the entire population being iterated in sequence.  The `randperm` function could well be doing this, but it seems like  `i` is also a random number. If it is then random entries of the population are being mutated on each evaluation, which I don't think is correct.\r\n"}}