grortir.externals.pyswarm._cons_f_ieqcons_wrapper() - Code Metrics - Inspection of "Add pySwarm code." - qbahn/grortir - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#29)

by Wojtek

created 2016-04-04 22:31 UTC

grortir.externals.pyswarm._cons_f_ieqcons_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
			0 ignored issues – show introduced 2016-04-04 22:33 UTC by Report Bug Copy Issue Report 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 2016-04-04 22:33 UTC by Report Bug Copy Issue Report 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 = 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:
			0 ignored issues – show Duplication introduced 2016-04-04 22:33 UTC by Report Bug Copy Issue Report 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

qbahn / grortir

Pull Request — master (#29)

grortir.externals.pyswarm._cons_f_ieqcons_wrapper() A

Complexity

Size

Duplication

Duplication Side-by-Side

Filter issues like