grortir.externals.pyswarm._is_feasible_wrapper() - Code Metrics - Inspection of "Merge pull request #30 from qbahn/qbahn-patch-1" - qbahn/grortir - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — master ( 8050ab...8ebc53 )

by Wojtek

created 2016-04-05 19:53 UTC

grortir.externals.pyswarm._is_feasible_wrapper() A

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Metric	Value
cc	1
dl	0
loc	2
rs	10

# pylint: skip-file
from functools import partial
import numpy as np

def _obj_wrapper(func, args, kwargs, x):
    return func(x, *args, **kwargs)

def _is_feasible_wrapper(func, x):
    return np.all(func(x)>=0)

def _cons_none_wrapper(x):
    return np.array([0])

def _cons_ieqcons_wrapper(ieqcons, args, kwargs, x):
    return np.array([y(x, *args, **kwargs) for y in ieqcons])

def _cons_f_ieqcons_wrapper(f_ieqcons, args, kwargs, x):
    return np.array(f_ieqcons(x, *args, **kwargs))

def pso(func, lb, ub, ieqcons=[], f_ieqcons=None, args=(), kwargs={},
        swarmsize=100, omega=0.5, phip=0.5, phig=0.5, maxiter=100,
        minstep=1e-8, minfunc=1e-8, debug=False, processes=1,
        particle_output=False):
    """
    Perform a particle swarm optimization (PSO)

    Parameters
    ==========
    func : function
        The function to be minimized
    lb : array
        The lower bounds of the design variable(s)
    ub : array
        The upper bounds of the design variable(s)

    Optional
    ========
    ieqcons : list
        A list of functions of length n such that ieqcons[j](x,*args) >= 0.0 in
        a successfully optimized problem (Default: [])
    f_ieqcons : function
        Returns a 1-D array in which each element must be greater or equal
        to 0.0 in a successfully optimized problem. If f_ieqcons is specified,
        ieqcons is ignored (Default: None)
    args : tuple
        Additional arguments passed to objective and constraint functions
        (Default: empty tuple)
    kwargs : dict
        Additional keyword arguments passed to objective and constraint
        functions (Default: empty dict)
    swarmsize : int
        The number of particles in the swarm (Default: 100)
    omega : scalar
        Particle velocity scaling factor (Default: 0.5)
    phip : scalar
        Scaling factor to search away from the particle's best known position
        (Default: 0.5)
    phig : scalar
        Scaling factor to search away from the swarm's best known position
        (Default: 0.5)
    maxiter : int
        The maximum number of iterations for the swarm to search (Default: 100)
    minstep : scalar
        The minimum stepsize of swarm's best position before the search
        terminates (Default: 1e-8)
    minfunc : scalar
        The minimum change of swarm's best objective value before the search
        terminates (Default: 1e-8)
    debug : boolean
        If True, progress statements will be displayed every iteration
        (Default: False)
    processes : int
        The number of processes to use to evaluate objective function and
        constraints (default: 1)
    particle_output : boolean
        Whether to include the best per-particle position and the objective
        values at those.

    Returns
    =======
    g : array
        The swarm's best known position (optimal design)
    f : scalar
        The objective value at ``g``
    p : array
        The best known position per particle
    pf: arrray
        The objective values at each position in p

    """

    assert len(lb)==len(ub), 'Lower- and upper-bounds must be the same length'
    assert hasattr(func, '__call__'), 'Invalid function handle'
    lb = np.array(lb)
    ub = np.array(ub)
    assert np.all(ub>lb), 'All upper-bound values must be greater than lower-bound values'

    vhigh = np.abs(ub - lb)
    vlow = -vhigh

    # Initialize objective function
    obj = partial(_obj_wrapper, func, args, kwargs)

    # Check for constraint function(s) #########################################
    if f_ieqcons is None:
        if not len(ieqcons):
            if debug:
                print('No constraints given.')
            cons = _cons_none_wrapper
        else:
            if debug:
                print('Converting ieqcons to a single constraint function')
            cons = partial(_cons_ieqcons_wrapper, ieqcons, args, kwargs)
    else:
        if debug:
            print('Single constraint function given in f_ieqcons')
        cons = partial(_cons_f_ieqcons_wrapper, f_ieqcons, args, kwargs)
    is_feasible = partial(_is_feasible_wrapper, cons)

    # Initialize the multiprocessing module if necessary
    if processes > 1:
        import multiprocessing
        mp_pool = multiprocessing.Pool(processes)

    # Initialize the particle swarm ############################################
    S = swarmsize
    D = len(lb)  # the number of dimensions each particle has
    x = np.random.rand(S, D)  # particle positions
    v = np.zeros_like(x)  # particle velocities
    p = np.zeros_like(x)  # best particle positions
    fx = np.zeros(S)  # current particle function values
    fs = np.zeros(S, dtype=bool)  # feasibility of each particle
    fp = np.ones(S)*np.inf  # best particle function values
    g = []  # best swarm position
    fg = np.inf  # best swarm position starting value

    # Initialize the particle's position
    x = lb + x*(ub - lb)

    # Calculate objective and constraints for each particle
    if processes > 1:
        fx = np.array(mp_pool.map(obj, x))
        fs = np.array(mp_pool.map(is_feasible, x))
    else:
        for i in range(S):
            fx[i] = obj(x[i, :])
            fs[i] = is_feasible(x[i, :])

    # Store particle's best position (if constraints are satisfied)
    i_update = np.logical_and((fx < fp), fs)
    p[i_update, :] = x[i_update, :].copy()
    fp[i_update] = fx[i_update]

    # Update swarm's best position
    i_min = np.argmin(fp)
    if fp[i_min] < fg:
        fg = fp[i_min]
        g = p[i_min, :].copy()
    else:
        # At the start, there may not be any feasible starting point, so just
        # give it a temporary "best" point since it's likely to change
        g = x[0, :].copy()

    # Initialize the particle's velocity
    v = vlow + np.random.rand(S, D)*(vhigh - vlow)

    # Iterate until termination criterion met ##################################
    it = 1
    while it <= maxiter:
        rp = np.random.uniform(size=(S, D))
        rg = np.random.uniform(size=(S, D))

        # Update the particles velocities
        v = omega*v + phip*rp*(p - x) + phig*rg*(g - x)
        # Update the particles' positions
        x = x + v
        # Correct for bound violations
        maskl = x < lb
        masku = x > ub
        x = x*(~np.logical_or(maskl, masku)) + lb*maskl + ub*masku

        # Update objectives and constraints
        if processes > 1:
            fx = np.array(mp_pool.map(obj, x))
            fs = np.array(mp_pool.map(is_feasible, x))
        else:
            for i in range(S):
                fx[i] = obj(x[i, :])
                fs[i] = is_feasible(x[i, :])

        # Store particle's best position (if constraints are satisfied)
        i_update = np.logical_and((fx < fp), fs)
        p[i_update, :] = x[i_update, :].copy()
        fp[i_update] = fx[i_update]

        # Compare swarm's best position with global best position
        i_min = np.argmin(fp)
        if fp[i_min] < fg:
            if debug:
                print('New best for swarm at iteration {:}: {:} {:}'\
                    .format(it, p[i_min, :], fp[i_min]))

            p_min = p[i_min, :].copy()
            stepsize = np.sqrt(np.sum((g - p_min)**2))

            if np.abs(fg - fp[i_min]) <= minfunc:
                print('Stopping search: Swarm best objective change less than {:}'\
                    .format(minfunc))
                if particle_output:
                    return p_min, fp[i_min], p, fp
                else:
                    return p_min, fp[i_min]
            elif stepsize <= minstep:
                print('Stopping search: Swarm best position change less than {:}'\
                    .format(minstep))
                if particle_output:
                    return p_min, fp[i_min], p, fp
                else:
                    return p_min, fp[i_min]
            else:
                g = p_min.copy()
                fg = fp[i_min]

        if debug:
            print('Best after iteration {:}: {:} {:}'.format(it, g, fg))
        it += 1

    print('Stopping search: maximum iterations reached --> {:}'.format(maxiter))

    if not is_feasible(g):
        print("However, the optimization couldn't find a feasible design. Sorry")
    if particle_output:
        return g, fg, p, fp
    else:
        return g, fg

1			# pylint: skip-file
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:
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 = omegav + phiprp(p - x) + phigrg*(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)) + lbmaskl + ub*masku
181
182			# Update objectives and constraints
183			if processes > 1:
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

qbahn / grortir

Push — master ( 8050ab...8ebc53 )

grortir.externals.pyswarm._is_feasible_wrapper() A

Complexity

Size

Duplication

Duplication Side-by-Side

Filter issues like