Solves optimization problem.

        Parameters
        ----------
        problem : Object
762458
def solve(self,problem):
        
    
        # Local vars
        norm2 = self.norm2
        norminf = self.norminf
        parameters = self.parameters
        
        # Parameters
        tol = parameters['tol']
        maxiter = parameters['maxiter']
        quiet = parameters['quiet']
        sigma = parameters['sigma']
        eps = parameters['eps']
        eps_cold = parameters['eps_cold']

        # Problem
        if not isinstance(problem,QuadProblem):
            problem = cast_problem(problem)
            quad_problem = QuadProblem(None,None,None,None,None,None,problem=problem)
        else:
            quad_problem = problem
        self.problem = problem
        self.quad_problem = quad_problem

        # Linsolver
        self.linsolver = new_linsolver(parameters['linsolver'],'symmetric')

        # Reset
        self.reset()

        # Checks
        if not np.all(problem.l <= problem.u):
            raise OptSolverError_NoInterior(self)
    
        # Data
        self.H = quad_problem.H
        self.g = quad_problem.g
        self.A = quad_problem.A
        self.AT = quad_problem.A.T
        self.b = quad_problem.b
        self.l = quad_problem.l-tol/10.
        self.u = quad_problem.u+tol/10.
        self.n = quad_problem.H.shape[0]
        self.m = quad_problem.A.shape[0]
        self.e = np.ones(self.n)
        self.I = eye(self.n,format='coo')
        self.Onm = coo_matrix((self.n,self.m))
        self.Omm = coo_matrix((self.m,self.m))

        # Initial primal
        if quad_problem.x is None:
            self.x = (self.u + self.l)/2.
        else:
            self.x = np.maximum(np.minimum(quad_problem.x,problem.u),problem.l)

        # Initial duals
        if quad_problem.lam is None:
            self.lam = np.zeros(self.m)
        else:
            self.lam = quad_problem.lam.copy()
        if quad_problem.mu is None:
            self.mu = np.ones(self.x.size)*eps_cold
        else:
            self.mu = np.maximum(quad_problem.mu,eps)
        if quad_problem.pi is None:
            self.pi = np.ones(self.x.size)*eps_cold
        else:
            self.pi = np.maximum(quad_problem.pi,eps)

        # Check interior
        try:
            assert(np.all(self.l < self.x)) 
            assert(np.all(self.x < self.u))
            assert(np.all(self.mu > 0))
            assert(np.all(self.pi > 0))
        except AssertionError:
            raise OptSolverError_Infeasibility(self)

        # Init vector
        self.y = np.hstack((self.x,self.lam,self.mu,self.pi))

        # Complementarity measures
        self.eta_mu = np.dot(self.mu,self.u-self.x)/self.x.size
        self.eta_pi = np.dot(self.pi,self.x-self.l)/self.x.size

        # Objective scaling
        fdata = self.func(self.y)
        self.obj_sca = np.maximum(norminf(self.g+self.H*self.x)/10.,1.)
        self.H = self.H/self.obj_sca
        self.g = self.g/self.obj_sca
        fdata = self.func(self.y)

        # Header
        if not quiet:
            print('\nSolver: IQP')
            print('-----------')
                                   
        # Outer
        s = 0.
        self.k = 0
        while True:

            # Complementarity measures
            self.eta_mu = np.dot(self.mu,self.u-self.x)/self.x.size
            self.eta_pi = np.dot(self.pi,self.x-self.l)/self.x.size
            
            # Init eval
            fdata = self.func(self.y)
            fmax = norminf(fdata.f)
            gmax = norminf(fdata.GradF)
            
            # Done
            if fmax < tol and sigma*np.maximum(self.eta_mu,self.eta_pi) < tol:
                self.set_status(self.STATUS_SOLVED)
                self.set_error_msg('')
                return

            # Target
            tau = sigma*norminf(fdata.GradF)
           
            # Header
            if not quiet:
                if self.k > 0:
                    print('')
                print('{0:^3s}'.format('iter'), end=' ')
                print('{0:^9s}'.format('phi'), end=' ')
                print('{0:^9s}'.format('fmax'), end=' ')
                print('{0:^9s}'.format('gmax'), end=' ')
                print('{0:^8s}'.format('cu'), end=' ')
                print('{0:^8s}'.format('cl'), end=' ')
                print('{0:^8s}'.format('s'))
 
            # Inner
            while True:
                
                # Eval
                fdata = self.func(self.y)
                fmax = norminf(fdata.f)
                gmax = norminf(fdata.GradF)
                compu = norminf(self.mu*(self.u-self.x))
                compl = norminf(self.pi*(self.x-self.l))
                phi = (0.5*np.dot(self.x,self.H*self.x)+np.dot(self.g,self.x))*self.obj_sca
                
                # Show progress
                if not quiet:
                    print('{0:^3d}'.format(self.k), end=' ')
                    print('{0:^9.2e}'.format(phi), end=' ')
                    print('{0:^9.2e}'.format(fmax), end=' ')
                    print('{0:^9.2e}'.format(gmax), end=' ')
                    print('{0:^8.1e}'.format(compu), end=' ')
                    print('{0:^8.1e}'.format(compl), end=' ')
                    print('{0:^8.1e}'.format(s))
                
                # Done
                if gmax < tau:
                    break

                # Done
                if fmax < tol and np.maximum(compu,compl) < tol:
                    break

                # Maxiters
                if self.k >= maxiter:
                    raise OptSolverError_MaxIters(self)
                    
                # Search direction
                ux = self.u-self.x
                xl = self.x-self.l
                D1 = spdiags(self.mu/ux,0,self.n,self.n,format='coo')
                D2 = spdiags(self.pi/xl,0,self.n,self.n,format='coo')
                fbar = np.hstack((-fdata.rd+fdata.ru/ux-fdata.rl/xl,fdata.rp))
                if self.A.shape[0] > 0:
                    Jbar = bmat([[tril(self.H)+D1+D2,None],
                                 [-self.A,self.Omm]],format='coo')
                else:
                    Jbar = bmat([[tril(self.H)+D1+D2]],
                                format='coo')
                try:
                    if not self.linsolver.is_analyzed():
                        self.linsolver.analyze(Jbar)
                    pbar = self.linsolver.factorize_and_solve(Jbar,fbar)
                except RuntimeError:
                    raise OptSolverError_BadLinSystem(self)
                px = pbar[:self.n]
                pmu = (-fdata.ru + self.mu*px)/ux
                ppi = (-fdata.rl - self.pi*px)/xl
                p = np.hstack((pbar,pmu,ppi))

                # Steplength bounds
                indices = px > 0
                s1 = np.min(np.hstack(((1.-eps)*(self.u-self.x)[indices]/px[indices],np.inf)))
                indices = px < 0
                s2 = np.min(np.hstack(((eps-1.)*(self.x-self.l)[indices]/px[indices],np.inf)))
                indices = pmu < 0
                s3 = np.min(np.hstack(((eps-1.)*self.mu[indices]/pmu[indices],np.inf)))
                indices = ppi < 0
                s4 = np.min(np.hstack(((eps-1.)*self.pi[indices]/ppi[indices],np.inf)))
                smax = np.min([s1,s2,s3,s4])
                
                # Line search
                s,fdata = self.line_search(self.y,p,fdata.F,fdata.GradF,self.func,smax)

                # Update x
                self.y += s*p
                self.k += 1
                self.x,self.lam,self.mu,self.pi = self.extract_components(self.y)

                # Check
                try:
                    assert(np.all(self.x < self.u))
                    assert(np.all(self.x > self.l))
                    assert(np.all(self.mu > 0))
                    assert(np.all(self.pi > 0))
                except AssertionError:
                    raise OptSolverError_Infeasibility(self)
1249857
def _solve(self, sense=None):
        
        # Remove the constraints from the last run
        while len(self._remove_constr) > 0:
            self._remove_constr.pop().delete()

        try:
            return self._prob.solve(sense=sense)
        except lp.SolverError as e:
            raise_from(MOMAError(text_type(e)), e)
        finally:
            self._remove_constr = []
1244593
def _optimize(self, period, tmpid, use_gradient=True):
        """"""
        theta_0 = [self.y.min(), self.y.max() - self.y.min(), 0]
        result = minimize(self._chi2, theta_0, jac=bool(use_gradient),
                          bounds=[(None, None), (0, None), (None, None)],
                          args=(period, tmpid, use_gradient))
        return result.x
523958
def optimize(self, piter=3, pmaxf=300, ppert=0.1):
        '''
        

        '''

        self._save_npz()
        optimized = pPLD(self.ID, piter=piter, pmaxf=pmaxf,
                         ppert=ppert, debug=True, clobber=True)
        optimized.publish()
        self.reset()
1195695
def _optimize_A(self, A):
        
        right_eigenvectors = self.right_eigenvectors_[:, :self.n_macrostates]
        flat_map, square_map = get_maps(A)
        alpha = to_flat(1.0 * A, flat_map)

        def obj(x):
            return -1 * self._objective_function(
                    x, self.transmat_, right_eigenvectors, square_map,
                    self.populations_
            )

        alpha = scipy.optimize.basinhopping(
                obj, alpha, niter_success=1000,
        )['x']

        alpha = scipy.optimize.fmin(
                obj, alpha, full_output=True, xtol=1E-4, ftol=1E-4,
                maxfun=5000, maxiter=100000
        )[0]

        if np.isneginf(obj(alpha)):
            raise ValueError(
                    "Error: minimization has not located a feasible point.")

        A = to_square(alpha, square_map)
        return A
1003680
def _solve_location_param(self):
        
        params = copy.copy(self.params)
        del params['loc']

        f = lambda location: self.distr_f.ppf(0.5, loc=location, **params) - 1
        return optimize.brentq(f, -10, 10)
131525
def _create_optimizer(self, X, y, status):
        """"""
        if self.optimizer is None:
            self.optimizer = 'avltree'

        times, ranks = y

        if self.optimizer == 'simple':
            optimizer = SimpleOptimizer(X, status, self.alpha, self.rank_ratio, timeit=self.timeit)
        elif self.optimizer == 'PRSVM':
            optimizer = PRSVMOptimizer(X, status, self.alpha, self.rank_ratio, timeit=self.timeit)
        elif self.optimizer == 'direct-count':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            SurvivalCounter(X, ranks, status, len(ranks), times), timeit=self.timeit)
        elif self.optimizer == 'rbtree':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            OrderStatisticTreeSurvivalCounter(X, ranks, status, RBTree, times),
                                            timeit=self.timeit)
        elif self.optimizer == 'avltree':
            optimizer = LargeScaleOptimizer(self.alpha, self.rank_ratio, self.fit_intercept,
                                            OrderStatisticTreeSurvivalCounter(X, ranks, status, AVLTree, times),
                                            timeit=self.timeit)
        else:
            raise ValueError('unknown optimizer: {0}'.format(self.optimizer))

        return optimizer
818780
def prepare_for_optimization(self):
    """"""
    if self.params['eig_type'] == 'TF':
      self.eig_vec_estimate = self.get_min_eig_vec_proxy()
    elif self.params['eig_type'] == 'LZS':
      self.eig_vec_estimate = self.dual_object.m_min_vec
    else:
      self.eig_vec_estimate = tf.placeholder(tf.float32, shape=(self.dual_object.matrix_m_dimension, 1))
    self.stopped_eig_vec_estimate = tf.stop_gradient(self.eig_vec_estimate)
    # Eig value is v^\top M v, where v is eigen vector
    self.eig_val_estimate = tf.matmul(
        tf.transpose(self.stopped_eig_vec_estimate),
        self.dual_object.get_psd_product(self.stopped_eig_vec_estimate))
    # Penalizing negative of min eigen value because we want min eig value
    # to be positive
    self.total_objective = (
        self.dual_object.unconstrained_objective
        + 0.5 * tf.square(
            tf.maximum(-self.penalty_placeholder * self.eig_val_estimate, 0)))
    global_step = tf.Variable(0, trainable=False)
    # Set up learning rate as a placeholder
    self.learning_rate = tf.placeholder(tf.float32, shape=[])

    # Set up the optimizer
    if self.params['optimizer'] == 'adam':
      self.optimizer = tf.train.AdamOptimizer(learning_rate=self.learning_rate)
    elif self.params['optimizer'] == 'adagrad':
      self.optimizer = tf.train.AdagradOptimizer(learning_rate=self.learning_rate)
    elif self.params['optimizer'] == 'momentum':
      self.optimizer = tf.train.MomentumOptimizer(
          learning_rate=self.learning_rate,
          momentum=self.params['momentum_parameter'],
          use_nesterov=True)
    else:
      self.optimizer = tf.train.GradientDescentOptimizer(
          learning_rate=self.learning_rate)

    # Write out the projection step
    self.train_step = self.optimizer.minimize(
        self.total_objective, global_step=global_step)

    self.sess.run(tf.global_variables_initializer())

    # Projecting the dual variables
    proj_ops = []
    for i in range(self.dual_object.nn_params.num_hidden_layers + 1):
      # Lambda_pos is non negative for switch indices,
      # Unconstrained for positive indices
      # Zero for negative indices
      proj_ops.append(self.dual_object.lambda_pos[i].assign(
          tf.multiply(self.dual_object.positive_indices[i],
                      self.dual_object.lambda_pos[i])+
          tf.multiply(self.dual_object.switch_indices[i],
                      tf.nn.relu(self.dual_object.lambda_pos[i]))))
      proj_ops.append(self.dual_object.lambda_neg[i].assign(
          tf.multiply(self.dual_object.negative_indices[i],
                      self.dual_object.lambda_neg[i])+
          tf.multiply(self.dual_object.switch_indices[i],
                      tf.nn.relu(self.dual_object.lambda_neg[i]))))
      # Lambda_quad is only non zero and positive for switch
      proj_ops.append(self.dual_object.lambda_quad[i].assign(
          tf.multiply(self.dual_object.switch_indices[i],
                      tf.nn.relu(self.dual_object.lambda_quad[i]))))
      # Lambda_lu is always non negative
      proj_ops.append(self.dual_object.lambda_lu[i].assign(
          tf.nn.relu(self.dual_object.lambda_lu[i])))

    self.proj_step = tf.group(proj_ops)

    # Create folder for saving stats if the folder is not None
    if (self.params.get('stats_folder') and
        not tf.gfile.IsDirectory(self.params['stats_folder'])):
      tf.gfile.MkDir(self.params['stats_folder'])
41533
def solve_fba(self, objective):
        
        self._prob.set_objective(self._v_wt[objective])
        return self._solve(lp.ObjectiveSense.Maximize)
1243584
def solve(self):
        
        if self.world_coords is None:
            log.warning("solving for Assembly without world coordinates set: %r", self)

        for (component, world_coords) in solver(self.constraints, self.world_coords):
            component.world_coords = world_coords
404965
def optimize(self, angles0, target):
        """"""
        def new_objective(angles):
            a = angles - angles0
            if isinstance(self.smooth_factor, (np.ndarray, list)):
                if len(a) == len(self.smooth_factor):
                    return (self.f(angles, target) +
                            np.sum(self.smooth_factor * np.power(a, 2)))
                else:
                    raise ValueError('len(smooth_factor) != number of joints')
            else:
                return (self.f(angles, target) +
                        self.smooth_factor * np.sum(np.power(a, 2)))

        return scipy.optimize.minimize(
            new_objective,
            angles0,
            **self.optimizer_opt).x
600200
def solve(self):
        '''
        
        '''
        self.defUtilityFuncs()
        self.DiscFacEff = self.DiscFac*self.LivPrb
        self.makePFcFunc()
        self.defValueFuncs()
        solution = ConsumerSolution(cFunc=self.cFunc, vFunc=self.vFunc, vPfunc=self.vPfunc,
                                    mNrmMin=self.mNrmMin, hNrm=self.hNrmNow,
                                    MPCmin=self.MPC, MPCmax=self.MPC)
        #solution = self.addSSmNrm(solution)
        return solution
113349
def optimise(self, data, initial_theta=None, full_output=False, **kwargs):
        

        data = self._format_data(data)

        if initial_theta is None:
            initial_theta = self.estimate(data)

        # Which parameters will be optimised, and which will be fixed?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)

        # For the Cannon model, the normalised fluxes are a multiplicative
        # operation, so we don't actually need to optimise the astrophysical
        # parameters, we only need to optimise the other ones. At each set of
        # non-astrophysical parameters we can determine the optimal
        # astrophysical parameters by a matrix inversion.
        parameters = [p for p in self.parameters if p not in ignore_parameters \
            and p not in self.grid_points.dtype.names]

        # What model wavelength ranges will be required?
        wavelengths_required = []
        for channel, spectrum in zip(matched_channels, data):
            if channel is None: continue
            z = initial_theta.get("z", initial_theta.get("z_{}".format(channel),
                0))
            wavelengths_required.append(
                [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

        # (If necessary,..) train the Cannon model around the closest point.
        subset_bounds = self._initialise_approximator(closest_point=\
            [initial_theta[p] for p in self.grid_points.dtype.names], 
            wavelengths_required=wavelengths_required, **kwargs)
        
        # Get the optimisation keyword arguments.
        op_kwargs = self._configuration.get("optimise", {}).copy()
        op_kwargs.update(kwargs)

        # Get fixed keywords.
        fixed = op_kwargs.pop("fixed", {})
        if fixed is None: fixed = {}
        if fixed:
            # Remove non-parameters from the 'fixed' keywords.
            keys = set(fixed).intersection(self.parameters)
            ignored_keys = set(fixed).difference(parameters)
            if ignored_keys:
                logger.warn("Ignoring optimise.fixed keys because they are not"
                    " model parameters: {0}".format(", ".join(ignored_keys)))
            # If the 'fixed' value is provided, use that. Otherwise if it is
            # None then use the initial_theta value.
            fixed = dict(zip(keys, 
                [(fixed[k], initial_theta.get(k, None))[fixed[k] is None] \
                    for k in keys]))

            logger.info("Fixing keyword arguments (these will not be optimised)"
                ": {}".format(fixed))

        # Remove fixed parameters from the parameters to be optimised
        #parameters = list(set(parameters).difference(fixed))
        parameters = [p for p in parameters if p not in fixed]

        # Translate input bounds.
        nbs = (None, None) # No boundaries.
        input_bounds = op_kwargs.pop("bounds", {})
        # Because astrophysical parameters are not being 'optimised', we may not
        # actually need the subset_bounds.get call, but there may be fringe
        # cases that I cannot immediately think of.
        op_kwargs["bounds"] = [input_bounds.get(p, subset_bounds.get(p, nbs)) \
            for p in parameters]

        # Apply data masks now so we don't have to do it on the fly.
        masked_data, pixels_affected = self._apply_data_mask(data)

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters)

        logger.info("Optimising parameters: {0}".format(", ".join(parameters)))
        logger.info("Optimisation keywords: {0}".format(op_kwargs))

        # Create the objective function. The objective function needs to
        # continuum-normalise the observed spectra, put it at rest, then solve
        # for the best astrophysical parameters, then *generate* that spectrum.
        debug = kwargs.get("debug", False)
        def nlp(t, return_labels=False):
            # Apply fixed keywords.
            theta = dict(zip(parameters, t))
            theta.update(fixed)

            observed_variances = []
            observed_intensities = []
            for channel, spectrum in zip(matched_channels, data):

                # For each channel:
                # 1) Put the observed spectrum at rest.
                z = theta.get("z", theta.get("z_{}".format(channel), 0))
                rest_observed_disp = spectrum.disp * (1 - z)

                # 2) Correct for continuum.
                j, coeffs = 0, []
                while theta.get("continuum_{0}_{1}".format(channel, j), None):
                    coeffs.append(theta["continuum_{0}_{1}".format(channel, j)])
                    j += 1

                # Remember: model.__call__ calculates continuum based on the
                # *observed* wavelength points, so here we do the same (e.g.,
                # not those that have potentially been corrected for redshift)
                if not coeffs: continuum = 1.0
                else:
                    continuum = np.polyval(coeffs[::-1], spectrum.disp) 
                
                # 3) Put the observed data onto the self.wavelengths scale.
                # [TODO] do the resampling correctly.
                rebinned_observed_intensities = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.flux / continuum,
                    left=np.nan, right=np.nan)

                rebinned_observed_variances = np.interp(self.wavelengths,
                    rest_observed_disp, spectrum.variance / continuum,
                    left=np.nan, right=np.nan)

                observed_variances.append(rebinned_observed_variances)
                observed_intensities.append(rebinned_observed_intensities)

            # [TODO] This may be the wrong thing to do.
            observed_variances \
                = np.nanmean(np.vstack(observed_variances), axis=0)
            observed_intensities \
                = np.nanmean(np.vstack(observed_intensities), axis=0)

            # Solve for the astrophysical parameters.
            try:
                labels = self._solve_labels(observed_intensities,
                    observed_variances)

            except:
                logger.exception("Could not determine labels:")
                if debug: raise
                return np.inf

            if return_labels:
                return labels

            # Put the optimal astrophysical parameters into our theta.
            theta.update(labels)

            # Make the log-probability call.
            # Note: The inference._ln_probability call takes a dictionary.
            return -inference._ln_probability(theta, self, data, debug,
                matched_channels=matched_channels)

        # Do the optimisation.
        p0 = np.array([initial_theta[p] for p in parameters])
        x_opt = optimise.minimise(nlp, p0, **op_kwargs)

        # Put the result into a usable form.
        result = dict(zip(parameters, x_opt))
        result.update(nlp(x_opt, return_labels=True))
        result.update(fixed)
        
        result = OrderedDict([(k, result[k]) \
            for k in self.parameters if k in result])

        if full_output:
            # Create model fluxes and calculate some metric.
            chi_sq, dof, model_fluxes = self._chi_sq(result, data)
            result = (result, chi_sq, dof, model_fluxes)

        # Remove any prepared convolution functions.
        self._destroy_convolution_functions()

        return result
64910
def _update_solution_data(self, s, HH, CC, C0):
        
        x = s["x"]
        Va_v = self.om.get_var("Va")
        Pg_v = self.om.get_var("Pg")

        Va = x[Va_v.i1:Va_v.iN + 1]
        Pg = x[Pg_v.i1:Pg_v.iN + 1]
#        f = 0.5 * dot(x.T * HH, x) + dot(CC.T, x)

        s["f"] = s["f"] + C0

        # Put the objective function value in the solution.
#        solution["f"] = f

        return Va, Pg
95195
def solve(self, solver_klass=None):
        
        # Start the clock.
        t0 = time()

        # Build an OPF model with variables and constraints.
        om = self._construct_opf_model(self.case)
        if om is None:
            return {"converged": False, "output": {"message": "No Ref Bus."}}

        # Call the specific solver.
#        if self.opt["verbose"]:
#            print '\nPYLON Version %s, %s', "0.4.2", "April 2010"
        if solver_klass is not None:
            result = solver_klass(om, opt=self.opt).solve()
        elif self.dc:
#            if self.opt["verbose"]:
#                print ' -- DC Optimal Power Flow\n'
            result = DCOPFSolver(om, opt=self.opt).solve()
        else:
#            if self.opt["verbose"]:
#                print ' -- AC Optimal Power Flow\n'
            result = PIPSSolver(om, opt=self.opt).solve()

        result["elapsed"] = time() - t0

        if self.opt.has_key("verbose"):
            if self.opt["verbose"]:
                logger.info("OPF completed in %.3fs." % result["elapsed"])

        return result
1266432
def solve_with_cvxopt(sdp, solverparameters=None):
    
    P = convert_to_picos(sdp)
    P.set_option("solver", "cvxopt")
    P.set_option("verbose", sdp.verbose)
    if solverparameters is not None:
        for key, value in solverparameters.items():
            P.set_option(key, value)
    solution = P.solve()
    x_mat = [np.array(P.get_valued_variable('X'))]
    y_mat = [np.array(P.get_constraint(i).dual)
             for i in range(len(P.constraints))]
    return -solution["cvxopt_sol"]["primal objective"] + \
        sdp.constant_term, \
        -solution["cvxopt_sol"]["dual objective"] + \
        sdp.constant_term, \
        x_mat, y_mat, solution["status"]
404947
def _solve_tsp(dist, niter):
    

    # number of nodes
    N = dist.shape[0]

    # tsp path for quick calculation of cost
    ii = np.arange(N)
    jj = np.hstack((np.arange(1, N), 0))

    # for each node, cache a sorted list of all other nodes in order of
    # increasing distance.
    dsort = [np.argsort(d) for d in dist]
    dsort = [d[d != i] for i, d in enumerate(dsort)]

    # randomly initialize path through graph
    path = np.random.permutation(N)
    idx = np.argsort(path)
    cost = np.sum(dist[path[ii], path[jj]])

    # keep track of objective function over time
    cost_hist = [cost]

    # optimization loop
    node = 0
    itercount = 0
    n = 0

    while n < N and itercount < niter:

        # count iterations
        itercount += 1

        # we'll try breaking the connection i -> j
        i = path[node]
        j = path[(node+1) % N]

        # We are breaking i -> j so we can remove the cost of that connection.
        c = cost - dist[i, j]

        # Search over nodes k that are closer to j than i.
        for k in dsort[j]:

            # Can safely continue if dist[i,j] < dist[k,j] for the remaining k.
            if k == i:
                n += 1
                break

            # Break connection k -> p.
            # Add connection j -> p.
            # Add connection i -> k.
            p = path[(idx[k]+1) % N]
            new_cost = c - dist[k, p] + dist[j, p] + dist[i, k]

            # If this swap improves the cost, implement it and move to next i.
            if new_cost < cost:
                path = reverse_segment(path, idx[j], idx[k])
                idx = np.argsort(path)
                cost = new_cost
                # Restart from the begining of the graph.
                cost_hist.append(cost)
                n = 0
                break

        # move to next node
        node = (node + 1) % N

    return path, cost_hist