diff --git a/HARK/parallel.py b/HARK/parallel.py index c45db667b..97c8b856f 100644 --- a/HARK/parallel.py +++ b/HARK/parallel.py @@ -129,35 +129,20 @@ def runCommands(agent, command_list): return agent -"""============================================================= - ========= Define a parallel Nelder-Mead algorithm ========== - =============================================================""" - - -def parallelNelderMead( - objFunc, - guess, - perturb=None, - P=1, - ftol=0.000001, - xtol=0.00000001, - maxiter=np.inf, - maxeval=np.inf, - r_param=1.0, - e_param=1.0, - c_param=0.5, - s_param=0.5, - maxcores=None, - name=None, - resume=False, - savefreq=None, - verbose=1, -): +#============================================================= +# ======== Define a parallel Nelder-Mead algorithm ========== +#============================================================= + +def parallelNelderMead(objFunc,guess,perturb=None,P=1,ftol=0.000001,xtol=0.00000001, + maxiter=np.inf,maxeval=np.inf,r_param=1.0,e_param=1.0, + c_param=0.5,s_param=0.5,maxthreads=None,name=None,resume=False, + savefreq=None,verbose=1): """ + A parallel implementation of the Nelder-Mead minimization algorithm, as described in Lee and Wiswall. For long optimization procedures, it can save progress between iterations and resume later. - + Parameters ---------- objFunc : function @@ -169,7 +154,7 @@ def parallelNelderMead( objFunc. If perturb[j] is non-zero, a simplex point will be created that perturbs the j-th element of guess by perturb[j]; if it is zero, then the j-th parameter of objFunc will not be optimized over. By - default, guess=None, indicating that all parameters should be optimized, + default, perturb=None, indicating that all parameters should be optimized, with an initial perturbation of 0.1*guess. P : int Degree of parallelization: the number of vertices of the simplex to try @@ -196,7 +181,7 @@ def parallelNelderMead( Parameter indicating magnitude of the contraction point calculation. s_param: float Parameter indicating magnitude of the shrink calculation. - maxcores : int + maxthreads : int The maximum number of CPU cores that the optimization should use, regardless of the size of the problem. name : string @@ -212,7 +197,7 @@ def parallelNelderMead( verbose : int Indicator for the verbosity of the optimization routine. Higher values generate more text output; verbose=0 produces no text output. - + Returns ------- min_point : np.array @@ -224,31 +209,35 @@ def parallelNelderMead( if resume: simplex, fvals, iters, evals = loadNelderMeadData(name) dim_count = fvals.size - 1 - N = dim_count + 1 # Number of points in simplex - K = simplex.shape[1] # Total number of parameters - # Otherwise, construct the initial simplex and array of function values - else: - if perturb is None: # Default: perturb each parameter by 10% - perturb = 0.1 * guess + N = dim_count+1 # Number of points in simplex + K = simplex.shape[1] # Total number of parameters + + # Otherwise, construct the initial simplex and array of function values + else: + if perturb is None: # Default: perturb each parameter by 10% + perturb = 0.1*guess guess[guess == 0] = 0.1 - params_to_opt = np.where(perturb != 0)[ - 0 - ] # Indices of which parameters to optimize - dim_count = params_to_opt.size # Number of parameters to search over - N = dim_count + 1 # Number of points in simplex - K = guess.size # Total number of parameters - simplex = np.tile(guess, (N, 1)) - for j in range( - dim_count - ): # Perturb each parameter to optimize by the specified distance - simplex[j + 1, params_to_opt[j]] = ( - simplex[j + 1, params_to_opt[j]] + perturb[params_to_opt[j]] - ) - # Initialize a few - fvals = np.zeros(dim_count + 1) + np.nan + + params_to_opt = np.where(perturb != 0)[0] # Indices of which parameters to optimize + dim_count = params_to_opt.size # Number of parameters to search over + N = dim_count+1 # Number of points in simplex + K = guess.size # Total number of parameters + simplex = np.tile(guess,(N,1)) + for j in range(dim_count): # Perturb each parameter to optimize by the specified distance + simplex[j+1,params_to_opt[j]] = simplex[j+1,params_to_opt[j]] + perturb[params_to_opt[j]] + + # Initialize iteration and evaluation counts, plus a 1D array of function values + fvals = np.zeros(dim_count+1) + np.nan + iters = 0 evals = 0 - + + # Make sure degree of parallelization is not illegal + if P > N - 1: + print('Requested degree of simplex parallelization is ' + str(P) + ', but dimension of optimization problem is only ' + str(N-1) + '.') + print('Degree of parallelization has been reduced to ' + str(N-1) + '.') + P = N-1 + # Create the pool of worker processes cpu_cores = multiprocessing.cpu_count() # Total number of available CPU cores cores_to_use = min(cpu_cores, dim_count) @@ -267,48 +256,32 @@ def parallelNelderMead( simplex = simplex[order, :] fmin = fvals[0] f_dist = np.abs(fmin - fvals[-1]) - x_dist = np.max( - np.sqrt( - np.sum(simplex ** 2.0 - np.tile(simplex[0, :], (N, 1)) ** 2.0, axis=1) - ) - ) + x_dist = np.max(np.sqrt(np.sum((simplex - np.tile(simplex[0,:],(N,1)))**2.0,axis=1))) if verbose > 0: - print( - "Evaluated the initial simplex: fmin=" - + str(fmin) - + ", f_dist=" - + str(f_dist) - + ", x_dist=" - + str(x_dist) - ) - else: # Resume an existing search that was cut short + print('Evaluated the initial simplex: fmin=' + str(fmin) + ', f_dist=' + str(f_dist) + ', x_dist=' + str(x_dist)) + if savefreq is not None: + saveNelderMeadData(name, simplex, fvals, iters, evals) + if verbose > 0: + print('Saved search progress in ' + name + '.txt') + else: # Resume an existing search that was cut short if verbose > 0: - print( - "Resuming search after " - + str(iters) - + " iterations and " - + str(evals) - + " function evaluations." - ) - + print('Resuming search after ' + str(iters) + ' iterations and ' + str(evals) + ' function evaluations.') + # Initialize some inputs for the multithreader - j_list = list(range(N - P, N)) - opt_params = [r_param, c_param, e_param] - - # Run the Nelder-Mead algorithm until a terminal condition is met + j_list = range(N-P,N) + opt_params= [r_param,c_param,e_param] + + # Run the Nelder-Mead algorithm until a terminal condition is met go = True while go: t_start = clock() iters += 1 if verbose > 0: - print("Beginning iteration #" + str(iters) + " now.") - - # Update the P worst points of the simplex - output = parallel( - delayed(parallelNelderMeadWorker)(objFunc, simplex, fvals, j, P, opt_params) - for j in j_list - ) - new_subsimplex = np.zeros((P, K)) + np.nan + print('Beginning iteration #' + str(iters) + ' now.') + + # Update the P worst points of the simplex + output = parallel(delayed(parallelNelderMeadWorker)(objFunc,simplex,fvals,j,P,opt_params) for j in j_list) + new_subsimplex = np.zeros((P,K)) + np.nan new_vals = np.zeros(P) + np.nan new_evals = 0 for i in range(P): @@ -316,7 +289,7 @@ def parallelNelderMead( new_vals[i] = output[i][1] new_evals += output[i][2] evals += new_evals - + # Check whether any updates actually happened old_subsimplex = simplex[(N - P) : N, :] if np.max(np.abs(new_subsimplex - old_subsimplex)) == 0: @@ -337,42 +310,22 @@ def parallelNelderMead( if verbose > 0: print("Updated the simplex successfully.") # Otherwise, update the simplex with the new results - simplex[(N - P) : N, :] = new_subsimplex - fvals[(N - P) : N] = new_vals - + simplex[(N-P):N,:] = new_subsimplex + fvals[(N-P):N] = new_vals + # Reorder the simplex from best to worst order = np.argsort(fvals) fvals = fvals[order] simplex = simplex[order, :] fmin = fvals[0] f_dist = np.abs(fmin - fvals[-1]) - x_dist = np.max( - np.sqrt( - np.sum(simplex ** 2.0 - np.tile(simplex[0, :], (N, 1)) ** 2.0, axis=1) - ) - ) + x_dist = np.max(np.sqrt(np.sum((simplex - np.tile(simplex[0,:],(N,1)))**2.0,axis=1))) t_end = clock() if verbose > 0: t_iter = t_end - t_start - print( - "Finished iteration #" - + str(iters) - + " with " - + str(new_evals) - + " evaluations (" - + str(evals) - + " cumulative) in " - + str(t_iter) - + " seconds." - ) - print( - "Simplex status: fmin=" - + str(fmin) - + ", f_dist=" - + str(f_dist) - + ", x_dist=" - + str(x_dist) - ) + print('Finished iteration #' + str(iters) +' with ' + str(new_evals) + ' evaluations (' + str(evals) + ' cumulative) in ' + str(t_iter) + ' seconds.') + print('Simplex status: fmin=' + str(fmin) + ', f_dist=' + str(f_dist) + ', x_dist=' + str(x_dist)) + # Check for terminal conditions if iters >= maxiter: @@ -386,25 +339,25 @@ def parallelNelderMead( print("Function tolerance reached, terminating successfully.") if x_dist < xtol: go = False - print("Parameter tolerance reached, terminating successfully.") - + print('Parameter tolerance reached, terminating successfully.') + # Save the progress of the estimation if desired if savefreq is not None: if (iters % savefreq) == 0: - saveNelderMeadData(name, simplex, fvals, iters, evals) - if verbose > 0: - print("Saved search progress in " + name + ".txt") - + saveNelderMeadData(name, simplex, fvals, iters, evals) + if verbose > 0: + print('Saved search progress in ' + name + '.txt') + # Return the results xopt = simplex[0, :] return xopt, fmin - - + + def saveNelderMeadData(name, simplex, fvals, iters, evals): """ Stores the progress of a parallel Nelder-Mead search in a text file so that it can be resumed later (after manual termination or a crash). - + Parameters ---------- name : string @@ -417,29 +370,33 @@ def saveNelderMeadData(name, simplex, fvals, iters, evals): The number of completed Nelder-Mead iterations. evals : int The cumulative number of function evaluations in the search process. - + Returns ------- - none + None """ - with open(name + ".txt", "w") as f: - my_writer = csv.writer(f, delimiter=" ") + N = simplex.shape[0] # Number of points in simplex + K = simplex.shape[1] # Total number of parameters + + with open(name + '.txt','wb') as f: + my_writer = csv.writer(f, delimiter=',') my_writer.writerow(simplex.shape) my_writer.writerow([iters, evals]) - my_writer.writerow(simplex.flatten()) + for n in range(N): + my_writer.writerow(simplex[n,:]) my_writer.writerow(fvals) - - + + def loadNelderMeadData(name): """ Reads the progress of a parallel Nelder-Mead search from a text file, as created by saveNelderMeadData(). - + Parameters ---------- name : string Name of the txt file from which to read search progress. - + Returns ------- simplex : np.array @@ -451,27 +408,37 @@ def loadNelderMeadData(name): evals : int The cumulative number of function evaluations in the search process. """ - with open(name + ".txt", "rb") as f: - my_reader = csv.reader(f, delimiter=" ") - my_shape_txt = next(my_reader) - shape0 = int(my_shape_txt[0]) - shape1 = int(my_shape_txt[1]) - my_nums_txt = next(my_reader) - iters = int(my_nums_txt[0]) - evals = int(my_nums_txt[1]) - simplex_flat = np.array(next(my_reader), dtype=float) - simplex = np.reshape(simplex_flat, (shape0, shape1)) - fvals = np.array(next(my_reader), dtype=float) + # Open the Nelder-Mead progress file + with open(name + '.txt','rb') as f: + my_reader = csv.reader(f, delimiter=',') - return simplex, fvals, iters, evals + # Get the shape of the simplex and initialize it + my_shape_txt = my_reader.next() + N = int(my_shape_txt[0]) + K = int(my_shape_txt[1]) + simplex = np.zeros((N, K)) + np.nan + + # Get number of iterations and cumulative evaluations from the next line + my_nums_txt = my_reader.next() + iters = int(my_nums_txt[0]) + evals = int(my_nums_txt[1]) + + # Read one line per point of the simplex + for n in range(N): + simplex[n, :] = np.array(my_reader.next(), dtype=float) + # Read the final line to get function values + fvals = np.array(my_reader.next(), dtype=float) -def parallelNelderMeadWorker(objFunc, simplex, f_vals, j, P, opt_params): + return simplex, fvals, iters, evals + + +def parallelNelderMeadWorker(objFunc,simplex,f_vals,j,P,opt_params): """ A worker process for the parallel Nelder-Mead algorithm. Updates one point in the simplex, returning its function value as well. Should basically never be called directly, only by parallelNelderMead(). - + Parameters ---------- objFunc : function @@ -487,7 +454,7 @@ def parallelNelderMeadWorker(objFunc, simplex, f_vals, j, P, opt_params): Degree of parallelization of the algorithm. opt_params : numpy.array Three element array with parameters for reflection, contraction, expansion. - + Returns ------- new_point : numpy.array @@ -506,20 +473,20 @@ def parallelNelderMeadWorker(objFunc, simplex, f_vals, j, P, opt_params): best_val = f_vals[0] # best value in the vertex next_val = f_vals[j - 1] # next best point in the simplex evals = 0 - + # Calculate the centroid of the "good" simplex points - N = simplex.shape[0] # number of points in simplex - centroid = np.mean(simplex[0 : (N - P), :], axis=0) - + N = simplex.shape[0] # number of points in simplex + centroid = np.mean(simplex[0:(N-P),:],axis=0) + # Calculate the reflection point and its function value r_point = centroid + alpha * (centroid - my_point) r_val = objFunc(r_point) evals += 1 - + # Case 1: the reflection point is better than best point - if r_val < best_val: - e_point = r_point + gamma * (r_point - centroid) - e_val = objFunc(e_point) # Calculate expansion point + if r_val < best_val: + e_point = r_point + gamma*(r_point - centroid) + e_val = objFunc(e_point) # Calculate expansion point evals += 1 if e_val < r_val: new_point = e_point @@ -539,8 +506,8 @@ def parallelNelderMeadWorker(objFunc, simplex, f_vals, j, P, opt_params): else: temp_point = my_point temp_val = my_val - c_point = beta * (centroid + temp_point) - c_val = objFunc(c_point) # Calculate contraction point + c_point = temp_point + beta*(centroid - temp_point) + c_val = objFunc(c_point) # Calculate contraction point evals += 1 if c_val < temp_val: new_point = c_point @@ -548,35 +515,14 @@ def parallelNelderMeadWorker(objFunc, simplex, f_vals, j, P, opt_params): else: new_point = temp_point new_val = temp_val - + # Return the outputs return new_point, new_val, evals def main(): print("Sorry, HARKparallel doesn't actually do much on its own.") - print("To see an example of multithreading in HARK, see /Testing/MultithreadDemo.") - print('To ensure full compatibility "out of the box", multithreading is not') - print("used in our models and applications; users can turn it on by modifying") - print("the source code slightly.") - - K = 36 - P = 24 - my_guess = np.random.rand(K) - 0.5 - - def testFunc1(x): - return np.sum(x ** 2.0) / x.size - - xopt, fmin = parallelNelderMead( - testFunc1, - my_guess, - P=P, - maxiter=300, - savefreq=100, - name="testfile", - resume=False, - ) - xopt2, fmin2 = parallelNelderMead(testFunc1, xopt, P=P) + print("You must import functions from it into your workspace.") if __name__ == "__main__":