Compute solution, multiprocessing implementation.
        Note: estimated speedup approximately by factor of
        number of physical cores. Virtual cores do not improve
        calculation time.
425237
def solve_mp(self, nproc=4, printmat=0, sendback=0, silent=False):
        ''''''
        # Initialize elements
        self.initialize()
        # Compute number of equations
        self.neq = np.sum([e.nunknowns for e in self.elementlist])
        if self.neq == 0: return
        if silent is False:
            print('Number of elements, Number of equations:', len(
                self.elementlist), ',', self.neq)
        if self.neq == 0:
            if silent is False: print('No unknowns. Solution complete')
            return
        mat = np.empty((self.neq, self.neq))
        rhs = np.empty(self.neq)

        # start multiprocessing
        if nproc is None:
            nproc = mp.cpu_count() - 1  # make no. of processes equal to 1 less than no. of cores
        elif nproc > mp.cpu_count():
            print("Given 'nproc' larger than no. of cores on machine. Setting 'nproc' to {}.".format(mp.cpu_count()))
            nproc = mp.cpu_count()

        pool = mp.Pool(processes=nproc)
        results = []
        for e in self.elementlist:
            if e.nunknowns > 0:
                results.append(pool.apply_async(e.equation))
            if silent is False:
                print('.', end='', flush=True)

        pool.close()
        pool.join()

        mat = np.empty((self.neq, self.neq))
        rhs = np.zeros(self.neq)

        ieq = 0

        for p in results:
            imat, irhs = p.get()
            mat[ieq:ieq + imat.shape[0], :] = imat
            rhs[ieq:ieq + irhs.shape[0]] = irhs
            ieq += imat.shape[0]

        # end multiprocessing

        if printmat:
            return mat, rhs
        sol = np.linalg.solve(mat, rhs)
        icount = 0
        for e in self.elementlist:
            if e.nunknowns > 0:
                e.setparams(sol[icount:icount + e.nunknowns])
                icount += e.nunknowns
        if silent is False:
            print()  # needed cause the dots are printed
            print('solution complete')
        elif (silent == 'dot') or (silent == '.'):
            print('.', end='', flush=True)
        if sendback:
            return sol
        return
656651
def infer_process_count(self):
        '''
        '''

        try:
            self.processes = multiprocessing.cpu_count()
        except NotImplementedError:
            self._logger.log(
                'error',
                'Could infer CPU count, setting number of processes back to 4'
            )
            self.processes = 4
215867
def benchmark(N, gates):
    """"""
    qubits = list(range(0, N))
    ket = qf.zero_state(N)

    for n in range(0, N):
        ket = qf.H(n).run(ket)

    for _ in range(0, (gates-N)//3):
        qubit0, qubit1 = random.sample(qubits, 2)
        ket = qf.X(qubit0).run(ket)
        ket = qf.T(qubit1).run(ket)
        ket = qf.CNOT(qubit0, qubit1).run(ket)

    return ket.vec.tensor
524371
protected void optimizeSNE(AffinityMatrix pij, double[][] sol) {
    final int size = pij.size();
    if(size * 3L * dim > 0x7FFF_FFFAL) {
      throw new AbortException("Memory exceeds Java array size limit.");
    }
    // Meta information on each point; joined for memory locality.
    // Gradient, Momentum, and learning rate
    // For performance, we use a flat memory layout!
    double[] meta = new double[size * 3 * dim];
    final int dim3 = dim * 3;
    for(int off = 2 * dim; off < meta.length; off += dim3) {
      Arrays.fill(meta, off, off + dim, 1.); // Initial learning rate
    }
    // Affinity matrix in projected space
    double[][] qij = new double[size][size];

    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Iterative Optimization", iterations, LOG) : null;
    Duration timer = LOG.isStatistics() ? LOG.newDuration(this.getClass().getName() + ".runtime.optimization").begin() : null;
    // Optimize
    for(int it = 0; it < iterations; it++) {
      double qij_sum = computeQij(qij, sol);
      computeGradient(pij, qij, 1. / qij_sum, sol, meta);
      updateSolution(sol, meta, it);
      LOG.incrementProcessed(prog);
    }
    LOG.ensureCompleted(prog);
    if(timer != null) {
      LOG.statistics(timer.end());
    }
  }
163921
def _compute(self):
        
        r = self
        for func, args, kwargs in self._todo:
            r = func(self)(*args, **kwargs)
        return r
819917
def pso(self, n_particles, n_iterations, lower_start=None, upper_start=None, threadCount=1, init_pos=None,
            mpi=False, print_key='PSO'):
        
        if lower_start is None or upper_start is None:
            lower_start, upper_start = np.array(self.lower_limit), np.array(self.upper_limit)
            print("PSO initialises its particles with default values")
        else:
            lower_start = np.maximum(lower_start, self.lower_limit)
            upper_start = np.minimum(upper_start, self.upper_limit)
        if mpi is True:
            pso = MpiParticleSwarmOptimizer(self.chain, lower_start, upper_start, n_particles, threads=1)
            if pso.isMaster():
                print('MPI option chosen')
        else:
            pso = ParticleSwarmOptimizer(self.chain, lower_start, upper_start, n_particles, threads=threadCount)
        if init_pos is None:
            init_pos = (upper_start - lower_start) / 2 + lower_start
        if not init_pos is None:
            pso.gbest.position = init_pos
            pso.gbest.velocity = [0]*len(init_pos)
            pso.gbest.fitness, _ = self.chain.likelihood(init_pos)
        X2_list = []
        vel_list = []
        pos_list = []
        time_start = time.time()
        if pso.isMaster():
            print('Computing the %s ...' % print_key)
        num_iter = 0
        for swarm in pso.sample(n_iterations):
            X2_list.append(pso.gbest.fitness*2)
            vel_list.append(pso.gbest.velocity)
            pos_list.append(pso.gbest.position)
            num_iter += 1
            if pso.isMaster():
                if num_iter % 10 == 0:
                    print(num_iter)
        if not mpi:
            result = pso.gbest.position
        else:
            result = MpiUtil.mpiBCast(pso.gbest.position)

        #if (pso.isMaster() and mpi is True) or self.chain.sampling_option == 'X2_catalogue':
        if mpi is True and not pso.isMaster():
            pass
        else:
            lens_dict, source_dict, lens_light_dict, ps_dict, kwargs_cosmo = self.chain.param.args2kwargs(result)
            print(pso.gbest.fitness * 2 / (self.chain.effectiv_num_data_points(lens_dict, source_dict, lens_light_dict, ps_dict)), 'reduced X^2 of best position')
            print(pso.gbest.fitness, 'logL')
            print(self.chain.effectiv_num_data_points(lens_dict, source_dict, lens_light_dict, ps_dict), 'effective number of data points')
            print(lens_dict, 'lens result')
            print(source_dict, 'source result')
            print(lens_light_dict, 'lens light result')
            print(ps_dict, 'point source result')
            print(kwargs_cosmo, 'cosmo result')
            time_end = time.time()
            print(time_end - time_start, 'time used for PSO', print_key)
            print('===================')
        return result, [X2_list, pos_list, vel_list, []]
299671
private void compute() {

        int n = graph.length;

        // the set 'S' contains the pairs of vertices adjacent to 'y'
        int[] s = new int[n];
        int sizeOfS;

        // order the vertices by degree
        int[] vertices = new int[n];
        for (int v = 0; v < n; v++) {
            vertices[ordering[v]] = v;
        }

        // if the graph is known to be a biconnected component (prepossessing)
        // and there is at least one vertex with a degree > 2 we can skip all
        // vertices of degree 2.
        //
        // otherwise the smallest possible cycle is {0,1,2} (no parallel edges
        // or loops) we can therefore don't need to do the first two shortest
        // paths calculations
        int first = biconnected && nDeg2Vertices < n ? nDeg2Vertices : 2;

        for (int i = first; i < n; i++) {
            final int r = vertices[i];

            ShortestPaths pathsFromR = new ShortestPaths(graph, null, r, limit / 2, ordering);

            // we only check the vertices which belong to the set Vr. this
            // set is vertices reachable from 'r' by only travelling though
            // vertices smaller then r. In the ShortestPaths API this is
            // name 'isPrecedingPathTo'.
            //
            // using Vr allows us to prune the number of vertices to check and
            // discover each cycle exactly once. This is possible as for each
            // simple cycle there is only one vertex with a maximum ordering.
            for (int j = 0; j < i; j++) {
                final int y = vertices[j];
                if (!pathsFromR.isPrecedingPathTo(y)) continue;

                // start refilling set 's' by resetting it's size
                sizeOfS = 0;

                // z is adjacent to y and belong to Vr
                for (final int z : graph[y]) {
                    if (!pathsFromR.isPrecedingPathTo(z)) continue;

                    final int distToZ = pathsFromR.distanceTo(z);
                    final int distToY = pathsFromR.distanceTo(y);

                    // the distance of the path to z is one less then the
                    // path to y. the vertices are adjacent, therefore z must
                    // also belong to the shortest path from r to y.
                    //
                    // we queue up (in 's') all the vertices adjacent to y for
                    // which this holds and then check these once we've processed
                    // all adjacent vertices
                    //
                    //  / ¯ ¯ z1 \          z1 and z2 are added to 's' and
                    // r          y - z3    checked later as p and q (see below)
                    //  \ _ _ z2 /
                    //
                    if (distToZ + 1 == distToY) {
                        s[sizeOfS++] = z;
                    }

                    // if the distances are equal we could have an odd cycle
                    // but we need to check the paths only intersect at 'r'.
                    //
                    // we check the intersect for cases like this, shortest
                    // cycle here is {p .. y, z .. p} not {r .. y, z .. r}
                    //
                    //           / ¯ ¯ y         / ¯ ¯ \ / ¯ ¯ y
                    //  r - - - p      |    or  r       p      |
                    //           \ _ _ z         \ _ _ / \ _ _ z
                    //
                    // if it's the shortest cycle then the intersect is just {r}
                    //
                    //  / ¯ ¯ y
                    // r      |
                    //  \ _ _ z
                    //
                    else if (distToZ == distToY && ordering[z] < ordering[y]) {
                        final int[] pathToY = pathsFromR.pathTo(y);
                        final int[] pathToZ = pathsFromR.pathTo(z);
                        if (singletonIntersect(pathToZ, pathToY)) {
                            Cycle cycle = new OddCycle(pathsFromR, pathToY, pathToZ);
                            add(cycle);
                        }
                    }
                }

                // check each pair vertices adjacent to 'y' for an
                // even cycle, as with the odd cycle we ensure the intersect
                // of the paths {r .. p} and {r .. q} is {r}.
                //
                //  / ¯ ¯ p \
                // r         y
                //  \ _ _ q /
                //
                for (int k = 0; k < sizeOfS; k++) {
                    for (int l = k + 1; l < sizeOfS; l++) {
                        int[] pathToP = pathsFromR.pathTo(s[k]);
                        int[] pathToQ = pathsFromR.pathTo(s[l]);
                        if (singletonIntersect(pathToP, pathToQ)) {
                            Cycle cycle = new EvenCycle(pathsFromR, pathToP, y, pathToQ);
                            add(cycle);
                        }
                    }
                }
            }
        }
    }
719221
def _syncronous_batch_evaluation(self,x):
        
        from multiprocessing import Process, Pipe

        # --- parallel evaluation of the function
        divided_samples = [x[i::self.n_procs] for i in range(self.n_procs)]
        pipe = [Pipe() for i in range(self.n_procs)]
        proc = [Process(target=spawn(self._eval_func),args=(c,k)) for k,(p,c) in zip(divided_samples,pipe)]
        [p.start() for p in proc]
        [p.join() for p in proc]

        # --- time of evaluation is set to constant (=1). This is one of the hypothesis of synchronous batch methods.
        f_evals = np.zeros((x.shape[0],1))
        cost_evals = np.ones((x.shape[0],1))
        i = 0
        for (p,c) in pipe:
            f_evals[i::self.n_procs] = p.recv()[0] # throw away costs
            i += 1
        return f_evals, cost_evals
163144
def _run_parallel_multiprocess(self):
        
        _log.debug("run.parallel.multiprocess.start")
        processes = []
        ProcRunner.instance = self
        for i in range(self._ncores):
            self._status.running(i)
            proc = multiprocessing.Process(target=ProcRunner.run, args=(i,))
            proc.start()
            processes.append(proc)
        for i in range(self._ncores):
            processes[i].join()
            code = processes[i].exitcode
            self._status.success(i) if 0 == code else self._status.fail(i)
        _log.debug("run.parallel.multiprocess.end states={}".format(self._status))
565343
def cpu_count(logical=True):
    
    if logical:
        from multiprocessing import cpu_count
        ncpu=cpu_count()
    else:
        import psutil
        ncpu=psutil.cpu_count(logical=False)
    return ncpu
1050494
def compute(mechanism, subsystem, purviews, cause_purviews,
                effect_purviews):
        
        concept = subsystem.concept(mechanism,
                                    purviews=purviews,
                                    cause_purviews=cause_purviews,
                                    effect_purviews=effect_purviews)
        # Don't serialize the subsystem.
        # This is replaced on the other side of the queue, and ensures
        # that all concepts in the CES reference the same subsystem.
        concept.subsystem = None
        return concept
411954
def parallel(processes, threads):
    
    pool = multithread(threads)
    pool.map(run_process, processes)
    pool.close()
    pool.join()
1281450
function cpuReduce (data, op) {
  var result = op(data[0], data[1])
  for (var i = 2; i < data.length; i++) {
    result = op(result, data[i])
  }
  return result
}
82371
def _process_non_parallel(self, X: Union[pd.DataFrame, np.ndarray], n_refs: int, cluster_array: np.ndarray):
        
        for gap_value, n_clusters in [self._calculate_gap(X, n_refs, n_clusters)
                                      for n_clusters in cluster_array]:
            yield (gap_value, n_clusters)
108949
def mulmod(computation: BaseComputation) -> None:
    
    left, right, mod = computation.stack_pop(num_items=3, type_hint=constants.UINT256)

    if mod == 0:
        result = 0
    else:
        result = (left * right) % mod
    computation.stack_push(result)
918647
def _compute_results(self):
        
        self.Y_best = best_value(self.Y)
        self.x_opt = self.X[np.argmin(self.Y),:]
        self.fx_opt = np.min(self.Y)
384694
private void doPrecalculateProcess() {
		if (discreteProcess != null && discreteProcess.length != 0) {
			return;
		}

		final int numberOfPaths			= this.getNumberOfPaths();
		final int numberOfFactors		= this.getNumberOfFactors();
		final int numberOfComponents	= this.getNumberOfComponents();

		// Allocate Memory
		discreteProcess			= new RandomVariableInterface[getTimeDiscretization().getNumberOfTimeSteps() + 1][getNumberOfComponents()];
		discreteProcessWeights	= new RandomVariableInterface[getTimeDiscretization().getNumberOfTimeSteps() + 1];

		// Set initial Monte-Carlo weights
		discreteProcessWeights[0] = stochasticDriver.getRandomVariableForConstant(1.0 / numberOfPaths);

		// Set initial value
		RandomVariableInterface[] initialState = getInitialState();
		final RandomVariableInterface[] currentState = new RandomVariableInterface[numberOfComponents];
		for (int componentIndex = 0; componentIndex < numberOfComponents; componentIndex++) {
			currentState[componentIndex] = initialState[componentIndex];
			discreteProcess[0][componentIndex] = applyStateSpaceTransform(componentIndex, currentState[componentIndex]);
		}

		/*
		 * Evolve the process using an Euler scheme.
		 * The evolution is performed multi-threadded.
		 * Each component of the vector runs in its own thread.
		 */

		// We do not allocate more threads the twice the number of processors.
		int numberOfThreads = Math.min(Math.max(2 * Runtime.getRuntime().availableProcessors(),1),numberOfComponents);
		executor = Executors.newFixedThreadPool(numberOfThreads);

		// Evolve process
		for (int timeIndex2 = 1; timeIndex2 < getTimeDiscretization().getNumberOfTimeSteps()+1; timeIndex2++) {
			final int timeIndex = timeIndex2;
			// Generate process from timeIndex-1 to timeIndex
			final double deltaT = getTime(timeIndex) - getTime(timeIndex - 1);

			// Fetch drift vector
			RandomVariableInterface[] drift = null;
			try {
				drift = getDrift(timeIndex - 1, discreteProcess[timeIndex - 1], null);
			}
			catch(Exception e) {
				throw new RuntimeException("Drift calculaton failed at time index " + timeIndex + " (time=" + getTime(timeIndex - 1) + ") . See cause of this exception for details.", e);
			}

			// Calculate new realization
			Vector<Future<RandomVariableInterface>> discreteProcessAtCurrentTimeIndex = new Vector<Future<RandomVariableInterface>>(numberOfComponents);
			discreteProcessAtCurrentTimeIndex.setSize(numberOfComponents);
			for (int componentIndex2 = 0; componentIndex2 < numberOfComponents; componentIndex2++) {
				final int componentIndex = componentIndex2;

				final RandomVariableInterface	driftOfComponent	= drift[componentIndex];

				// Check if the component process has stopped to evolve
				if (driftOfComponent == null) {
					continue;
				}

				Callable<RandomVariableInterface> worker = new  Callable<RandomVariableInterface>() {
					public RandomVariableInterface call() {
						if(scheme == Scheme.EULER_FUNCTIONAL) {
							currentState[componentIndex] = applyStateSpaceTransformInverse(componentIndex, discreteProcess[timeIndex - 1][componentIndex]);
						}

						RandomVariableInterface[]	factorLoadings		= getFactorLoading(timeIndex - 1, componentIndex, discreteProcess[timeIndex - 1]);

						// Check if the component process has stopped to evolve
						if (factorLoadings == null) {
							return null;
						}

						// Apply drift
						if(driftOfComponent != null) {
							currentState[componentIndex] = currentState[componentIndex].addProduct(driftOfComponent, deltaT);
						}

						// Apply diffusion
						RandomVariableInterface[] brownianIncrement	= stochasticDriver.getIncrement(timeIndex - 1);
						currentState[componentIndex] = currentState[componentIndex].addSumProduct(Arrays.asList(factorLoadings), Arrays.asList(brownianIncrement));

						// Transform the state space to the value space and return it.
						return applyStateSpaceTransform(componentIndex, currentState[componentIndex]).cache();
					}
				};


				/*
				 * Optional multi-threadding (asyncronous calculation of the components)
				 */
				Future<RandomVariableInterface> result = null;
				try {
					if(isUseMultiThreadding) {
						result = executor.submit(worker);
					} else {
						result = new FutureWrapper<RandomVariableInterface>(worker.call());
					}
				} catch (Exception e) {
					throw new RuntimeException("Euler step failed at time index " + timeIndex + " (time=" + getTime(timeIndex) + "). See cause of this exception for details.", e);
				}

				// The following line will add the result of the calculation to the vector discreteProcessAtCurrentTimeIndex
				discreteProcessAtCurrentTimeIndex.set(componentIndex, result);
			}

			// Fetch results and move to discreteProcess[timeIndex]
			for (int componentIndex = 0; componentIndex < numberOfComponents; componentIndex++) {
				try {
					Future<RandomVariableInterface> discreteProcessAtCurrentTimeIndexAndComponent = discreteProcessAtCurrentTimeIndex.get(componentIndex);
					if(discreteProcessAtCurrentTimeIndexAndComponent != null) {
						discreteProcess[timeIndex][componentIndex] = discreteProcessAtCurrentTimeIndexAndComponent.get().cache();
					} else {
						discreteProcess[timeIndex][componentIndex] = null;
					}
				} catch (InterruptedException e) {
					throw new RuntimeException("Euler step failed at time index " + timeIndex + " (time=" + getTime(timeIndex) + "). See cause of this exception for details.", e);
				} catch (ExecutionException e) {
					throw new RuntimeException("Euler step failed at time index " + timeIndex + " (time=" + getTime(timeIndex) + "). See cause of this exception for details.", e);
				}
			}

			if (scheme == Scheme.PREDICTOR_CORRECTOR) {
				// Apply corrector step to realizations at next time step

				RandomVariableInterface[] driftWithPredictor = getDrift(timeIndex - 1, discreteProcess[timeIndex], null);

				for (int componentIndex = 0; componentIndex < getNumberOfComponents(); componentIndex++) {
					RandomVariableInterface driftWithPredictorOfComponent		= driftWithPredictor[componentIndex];
					RandomVariableInterface driftWithoutPredictorOfComponent	= drift[componentIndex];

					if (driftWithPredictorOfComponent == null || driftWithoutPredictorOfComponent == null) {
						continue;
					}

					// Calculated the predictor corrector drift adjustment
					RandomVariableInterface driftAdjustment = driftWithPredictorOfComponent.sub(driftWithoutPredictorOfComponent).div(2.0).mult(deltaT);

					// Add drift adjustment
					currentState[componentIndex] = currentState[componentIndex].add(driftAdjustment);

					// Re-apply state space transform
					discreteProcess[timeIndex][componentIndex] = applyStateSpaceTransform(componentIndex, currentState[componentIndex]);
				} // End for(componentIndex)
			} // End if(scheme == Scheme.PREDICTOR_CORRECTOR)

			// Set Monte-Carlo weights
			discreteProcessWeights[timeIndex] = discreteProcessWeights[timeIndex - 1];
		} // End for(timeIndex)

		try {
			executor.shutdown();
		}
		catch(SecurityException e) {
			// @TODO Improve exception handling here
		}
	}