Completed
Push — master ( 1d0a34...8050ab )
by Wojtek
8s
created

grortir.externals.pyswarm._cons_none_wrapper()   A

Complexity

Conditions 1

Size

Total Lines 2

Duplication

Lines 0
Ratio 0 %
Metric Value
cc 1
dl 0
loc 2
rs 10
1
# pylint: skip-file
0 ignored issues
show
introduced by
Ignoring entire file
Loading history...
2
from functools import partial
3
import numpy as np
4
5
def _obj_wrapper(func, args, kwargs, x):
6
    return func(x, *args, **kwargs)
7
8
def _is_feasible_wrapper(func, x):
9
    return np.all(func(x)>=0)
10
11
def _cons_none_wrapper(x):
12
    return np.array([0])
13
14
def _cons_ieqcons_wrapper(ieqcons, args, kwargs, x):
15
    return np.array([y(x, *args, **kwargs) for y in ieqcons])
16
17
def _cons_f_ieqcons_wrapper(f_ieqcons, args, kwargs, x):
18
    return np.array(f_ieqcons(x, *args, **kwargs))
19
20
def pso(func, lb, ub, ieqcons=[], f_ieqcons=None, args=(), kwargs={},
21
        swarmsize=100, omega=0.5, phip=0.5, phig=0.5, maxiter=100,
22
        minstep=1e-8, minfunc=1e-8, debug=False, processes=1,
23
        particle_output=False):
24
    """
25
    Perform a particle swarm optimization (PSO)
26
27
    Parameters
28
    ==========
29
    func : function
30
        The function to be minimized
31
    lb : array
32
        The lower bounds of the design variable(s)
33
    ub : array
34
        The upper bounds of the design variable(s)
35
36
    Optional
37
    ========
38
    ieqcons : list
39
        A list of functions of length n such that ieqcons[j](x,*args) >= 0.0 in
40
        a successfully optimized problem (Default: [])
41
    f_ieqcons : function
42
        Returns a 1-D array in which each element must be greater or equal
43
        to 0.0 in a successfully optimized problem. If f_ieqcons is specified,
44
        ieqcons is ignored (Default: None)
45
    args : tuple
46
        Additional arguments passed to objective and constraint functions
47
        (Default: empty tuple)
48
    kwargs : dict
49
        Additional keyword arguments passed to objective and constraint
50
        functions (Default: empty dict)
51
    swarmsize : int
52
        The number of particles in the swarm (Default: 100)
53
    omega : scalar
54
        Particle velocity scaling factor (Default: 0.5)
55
    phip : scalar
56
        Scaling factor to search away from the particle's best known position
57
        (Default: 0.5)
58
    phig : scalar
59
        Scaling factor to search away from the swarm's best known position
60
        (Default: 0.5)
61
    maxiter : int
62
        The maximum number of iterations for the swarm to search (Default: 100)
63
    minstep : scalar
64
        The minimum stepsize of swarm's best position before the search
65
        terminates (Default: 1e-8)
66
    minfunc : scalar
67
        The minimum change of swarm's best objective value before the search
68
        terminates (Default: 1e-8)
69
    debug : boolean
70
        If True, progress statements will be displayed every iteration
71
        (Default: False)
72
    processes : int
73
        The number of processes to use to evaluate objective function and
74
        constraints (default: 1)
75
    particle_output : boolean
76
        Whether to include the best per-particle position and the objective
77
        values at those.
78
79
    Returns
80
    =======
81
    g : array
82
        The swarm's best known position (optimal design)
83
    f : scalar
84
        The objective value at ``g``
85
    p : array
86
        The best known position per particle
87
    pf: arrray
88
        The objective values at each position in p
89
90
    """
91
92
    assert len(lb)==len(ub), 'Lower- and upper-bounds must be the same length'
93
    assert hasattr(func, '__call__'), 'Invalid function handle'
94
    lb = np.array(lb)
95
    ub = np.array(ub)
96
    assert np.all(ub>lb), 'All upper-bound values must be greater than lower-bound values'
97
98
    vhigh = np.abs(ub - lb)
99
    vlow = -vhigh
100
101
    # Initialize objective function
102
    obj = partial(_obj_wrapper, func, args, kwargs)
103
104
    # Check for constraint function(s) #########################################
105
    if f_ieqcons is None:
106
        if not len(ieqcons):
107
            if debug:
108
                print('No constraints given.')
109
            cons = _cons_none_wrapper
110
        else:
111
            if debug:
112
                print('Converting ieqcons to a single constraint function')
113
            cons = partial(_cons_ieqcons_wrapper, ieqcons, args, kwargs)
114
    else:
115
        if debug:
116
            print('Single constraint function given in f_ieqcons')
117
        cons = partial(_cons_f_ieqcons_wrapper, f_ieqcons, args, kwargs)
118
    is_feasible = partial(_is_feasible_wrapper, cons)
119
120
    # Initialize the multiprocessing module if necessary
121
    if processes > 1:
122
        import multiprocessing
123
        mp_pool = multiprocessing.Pool(processes)
124
125
    # Initialize the particle swarm ############################################
126
    S = swarmsize
127
    D = len(lb)  # the number of dimensions each particle has
128
    x = np.random.rand(S, D)  # particle positions
129
    v = np.zeros_like(x)  # particle velocities
130
    p = np.zeros_like(x)  # best particle positions
131
    fx = np.zeros(S)  # current particle function values
132
    fs = np.zeros(S, dtype=bool)  # feasibility of each particle
133
    fp = np.ones(S)*np.inf  # best particle function values
134
    g = []  # best swarm position
135
    fg = np.inf  # best swarm position starting value
136
137
    # Initialize the particle's position
138
    x = lb + x*(ub - lb)
139
140
    # Calculate objective and constraints for each particle
141
    if processes > 1:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
142
        fx = np.array(mp_pool.map(obj, x))
143
        fs = np.array(mp_pool.map(is_feasible, x))
144
    else:
145
        for i in range(S):
146
            fx[i] = obj(x[i, :])
147
            fs[i] = is_feasible(x[i, :])
148
149
    # Store particle's best position (if constraints are satisfied)
150
    i_update = np.logical_and((fx < fp), fs)
151
    p[i_update, :] = x[i_update, :].copy()
152
    fp[i_update] = fx[i_update]
153
154
    # Update swarm's best position
155
    i_min = np.argmin(fp)
156
    if fp[i_min] < fg:
157
        fg = fp[i_min]
158
        g = p[i_min, :].copy()
159
    else:
160
        # At the start, there may not be any feasible starting point, so just
161
        # give it a temporary "best" point since it's likely to change
162
        g = x[0, :].copy()
163
164
    # Initialize the particle's velocity
165
    v = vlow + np.random.rand(S, D)*(vhigh - vlow)
166
167
    # Iterate until termination criterion met ##################################
168
    it = 1
169
    while it <= maxiter:
170
        rp = np.random.uniform(size=(S, D))
171
        rg = np.random.uniform(size=(S, D))
172
173
        # Update the particles velocities
174
        v = omega*v + phip*rp*(p - x) + phig*rg*(g - x)
175
        # Update the particles' positions
176
        x = x + v
177
        # Correct for bound violations
178
        maskl = x < lb
179
        masku = x > ub
180
        x = x*(~np.logical_or(maskl, masku)) + lb*maskl + ub*masku
181
182
        # Update objectives and constraints
183
        if processes > 1:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
184
            fx = np.array(mp_pool.map(obj, x))
185
            fs = np.array(mp_pool.map(is_feasible, x))
186
        else:
187
            for i in range(S):
188
                fx[i] = obj(x[i, :])
189
                fs[i] = is_feasible(x[i, :])
190
191
        # Store particle's best position (if constraints are satisfied)
192
        i_update = np.logical_and((fx < fp), fs)
193
        p[i_update, :] = x[i_update, :].copy()
194
        fp[i_update] = fx[i_update]
195
196
        # Compare swarm's best position with global best position
197
        i_min = np.argmin(fp)
198
        if fp[i_min] < fg:
199
            if debug:
200
                print('New best for swarm at iteration {:}: {:} {:}'\
201
                    .format(it, p[i_min, :], fp[i_min]))
202
203
            p_min = p[i_min, :].copy()
204
            stepsize = np.sqrt(np.sum((g - p_min)**2))
205
206
            if np.abs(fg - fp[i_min]) <= minfunc:
207
                print('Stopping search: Swarm best objective change less than {:}'\
208
                    .format(minfunc))
209
                if particle_output:
210
                    return p_min, fp[i_min], p, fp
211
                else:
212
                    return p_min, fp[i_min]
213
            elif stepsize <= minstep:
214
                print('Stopping search: Swarm best position change less than {:}'\
215
                    .format(minstep))
216
                if particle_output:
217
                    return p_min, fp[i_min], p, fp
218
                else:
219
                    return p_min, fp[i_min]
220
            else:
221
                g = p_min.copy()
222
                fg = fp[i_min]
223
224
        if debug:
225
            print('Best after iteration {:}: {:} {:}'.format(it, g, fg))
226
        it += 1
227
228
    print('Stopping search: maximum iterations reached --> {:}'.format(maxiter))
229
230
    if not is_feasible(g):
231
        print("However, the optimization couldn't find a feasible design. Sorry")
232
    if particle_output:
233
        return g, fg, p, fp
234
    else:
235
        return g, fg