r"""
        Builds the A and b matrices, and calls the solver specified in the
        ``settings`` attribute.

        Parameters
        ----------
        x : ND-array
            Initial guess of unknown variable
970994
def run(self, x=None):
        

        """
        logger.info('Running ReactiveTransport')
        # Create S1 & S1 for 1st Picard's iteration
        if x is None:
            x = np.zeros(shape=[self.Np, ], dtype=float)
        self[self.settings['quantity']] = x
        self._update_physics()

        x = self._run_reactive(x=x)
        self[self.settings['quantity']] = x
647825
def optimize_NxN(x, y, fit_offset=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = get_weight(x, y, perc).astype(bool)
        if issparse(weights): weights = weights.A
    else:
        weights = None

    x, y = x.astype(np.float64), y.astype(np.float64)

    n_vars = x.shape[1]
    offset, gamma = np.zeros(n_vars), np.zeros(n_vars)

    for i in range(n_vars):
        xi = x[:, i] if weights is None else x[:, i][weights[:, i]]
        yi = y[:, i] if weights is None else y[:, i][weights[:, i]]

        if fit_offset:
            offset[i], gamma[i] = minimize(lambda m: np.sum((-yi + xi * m[1] + m[0])**2), method="L-BFGS-B",
                                           x0=(0, 0.1), bounds=[(0, None), (None, None)]).x
        else:
            gamma[i] = minimize(lambda m: np.sum((-yi + xi * m) ** 2), x0=0.1, method="L-BFGS-B").x
    offset[np.isnan(offset)], gamma[np.isnan(gamma)] = 0, 0
    return offset, gamma
468556
def infer_x(self, y_desired, **kwargs):
        
        
        sigma = kwargs.get('sigma', self.sigma)
        k     = kwargs.get('k', self.k)
         
        xq = self._guess_x(y_desired, k = k, sigma = sigma, **kwargs)[0]    
         
        dists, index = self.fmodel.dataset.nn_x(xq, k = k)

        w = self.fmodel._weights(dists, sigma*sigma)
        
        X   = np.array([self.fmodel.dataset.get_x_padded(i) for i in index])
        Y    = np.array([self.fmodel.dataset.get_y(i) for i in index])
         
        W   = np.diag(w)
        WX  = np.dot(W, X)
        WXT = WX.T
             
        B   = np.dot(np.linalg.pinv(np.dot(WXT, WX)),WXT)
        
        M = np.dot(B, np.dot(W, Y))
        
        _, idx = self.fmodel.dataset.nn_y(y_desired, k=1)        
        ynn = self.fmodel.dataset.get_y(idx[0])
        
        eps = 0.00001
        
        J = np.zeros((len(xq),len(y_desired)), dtype = np.float64)
        for i in range(len(xq)):
            xi = np.array(xq, dtype = np.float64)
            xi[i] = xi[i] + eps
            yi = np.dot(np.append(1.0, xi), M).ravel()
            J[i] = (yi - ynn) / eps
        J = np.transpose(J)
        Jinv = np.linalg.pinv(J)
        x = np.dot(Jinv, y_desired - ynn)
        
        return [x]
1329533
def fit(self, X, y):
    
    X, labels = self._prepare_inputs(X, y, ensure_min_samples=2)
    n, d = X.shape
    num_dims = self.num_dims
    if num_dims is None:
        num_dims = d

    # Measure the total training time
    train_time = time.time()

    # Initialize A to a scaling matrix
    A = np.zeros((num_dims, d))
    np.fill_diagonal(A, 1./(np.maximum(X.max(axis=0)-X.min(axis=0), EPS)))

    # Run NCA
    mask = labels[:, np.newaxis] == labels[np.newaxis, :]
    optimizer_params = {'method': 'L-BFGS-B',
                        'fun': self._loss_grad_lbfgs,
                        'args': (X, mask, -1.0),
                        'jac': True,
                        'x0': A.ravel(),
                        'options': dict(maxiter=self.max_iter),
                        'tol': self.tol
                        }

    # Call the optimizer
    self.n_iter_ = 0
    opt_result = minimize(**optimizer_params)

    self.transformer_ = opt_result.x.reshape(-1, X.shape[1])
    self.n_iter_ = opt_result.nit

    # Stop timer
    train_time = time.time() - train_time
    if self.verbose:
      cls_name = self.__class__.__name__

      # Warn the user if the algorithm did not converge
      if not opt_result.success:
        warnings.warn('[{}] NCA did not converge: {}'.format(
            cls_name, opt_result.message), ConvergenceWarning)

      print('[{}] Training took {:8.2f}s.'.format(cls_name, train_time))

    return self
1325707
def infer_x(self, y):
        
        OptimizedInverseModel.infer_x(self, y)
        if self.fmodel.size() == 0:
            return self._random_x()

        x_guesses = [self._guess_x_simple(y)[0]]
        result = []
        for i, xg in enumerate(x_guesses):
            res = scipy.optimize.minimize(self._error, xg,
                                          args        = (),
                                          method      = self.algo,
                                          bounds      = self.constraints,
                                          options     = self.conf
                                         )



            d = self._error(res.x)
            result.append((d, i, res.x))
        return [self._enforce_bounds(xi) for fi, i, xi in sorted(result)]
194842
def predict_withGradients(self, X):
        
        if X.ndim==1: X = X[None,:]
        ps = self.model.param_array.copy()
        means = []
        stds = []
        dmdxs = []
        dsdxs = []
        for s in self.hmc_samples:
            if self.model._fixes_ is None:
                self.model[:] = s
            else:
                self.model[self.model._fixes_] = s
            self.model._trigger_params_changed()
            m, v = self.model.predict(X)
            std = np.sqrt(np.clip(v, 1e-10, np.inf))
            dmdx, dvdx = self.model.predictive_gradients(X)
            dmdx = dmdx[:,:,0]
            dsdx = dvdx / (2*std)
            means.append(m)
            stds.append(std)
            dmdxs.append(dmdx)
            dsdxs.append(dsdx)
        self.model.param_array[:] = ps
        self.model._trigger_params_changed()
        return means, stds, dmdxs, dsdxs
1361852
def solve_and_plot_series(self, x0, params, varied_data, varied_idx, solver=None, plot_kwargs=None,
                              plot_residuals_kwargs=None, **kwargs):
        
        sol, nfo = self.solve_series(
            x0, params, varied_data, varied_idx, solver=solver, **kwargs)
        ax_sol = self.plot_series(sol, varied_data, varied_idx, info=nfo,
                                  **(plot_kwargs or {}))

        extra = dict(ax_sol=ax_sol, info=nfo)
        if plot_residuals_kwargs:
            extra['ax_resid'] = self.plot_series_residuals_internal(
                varied_data, varied_idx, info=nfo,
                **(plot_residuals_kwargs or {})
            )
        return sol, extra
647199
def init_nn_params(scale, layer_sizes, rs=npr.RandomState(0)):
    """"""
    return [(rs.randn(insize, outsize) * scale,   # weight matrix
             rs.randn(outsize) * scale)           # bias vector
            for insize, outsize in zip(layer_sizes[:-1], layer_sizes[1:])]
1025858
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        X.reshape(blockA.numCols,B.numCols);
        blockB.reshape(B.numRows,B.numCols,false);
        blockX.reshape(X.numRows,X.numCols,false);
        MatrixOps_DDRB.convert(B,blockB);

        alg.solve(blockB,blockX);

        MatrixOps_DDRB.convert(blockX,X);
    }
947061
public static void solveL(DMatrixSparseCSC L , double []x )
    {
        final int N = L.numCols;

        int idx0 = L.col_idx[0];
        for (int col = 0; col < N; col++) {

            int idx1 = L.col_idx[col+1];
            double x_j = x[col] /= L.nz_values[idx0];

            for (int i = idx0+1; i < idx1; i++) {
                int row = L.nz_rows[i];
                x[row] -=  L.nz_values[i]*x_j;
            }

            idx0 = idx1;
        }
    }
966146
def _apply_adj(self, a):
        # is projection the zero operator?
        if self.V.shape[1] == 0:
            return numpy.zeros(a.shape)
        ''''''
        c = inner(self.V, a, ip_B=self.ip_B)
        if self.Q is not None and self.R is not None:
            c = self.Q.dot(scipy.linalg.solve_triangular(self.R.T.conj(), c,
                                                         lower=True))
        return self.W.dot(c)
1369553
def forward(self, x):
        
        fwd = np.zeros((self.num_matrices, self.N, self.N), settings.float_type)
        indices = np.tril_indices(self.N, 0)
        z = np.zeros(len(indices[0])).astype(int)
        for i in range(self.num_matrices):
            fwd[(z + i,) + indices] = x[i, :]
        return fwd.squeeze(axis=0) if self.squeeze else fwd
886370
private double[][] solveL(double[][] X) {
    final int n = L.length, m = X[0].length;
    { // First iteration, simplified.
      final double[] X0 = X[0], L0 = L[0];
      for(int j = 0; j < m; j++) {
        X0[j] /= L0[0];
      }
    }
    for(int k = 1; k < n; k++) {
      final double[] Xk = X[k], Lk = L[k];
      final double iLkk = 1. / Lk[k];
      for(int j = 0; j < m; j++) {
        double Xkj = Xk[j];
        for(int i = 0; i < k; i++) {
          Xkj -= X[i][j] * Lk[i];
        }
        Xk[j] = Xkj * iLkk;
      }
    }
    return X;
  }
212684
def implicit_step(self):
        
        config = self.config
        system = self.system
        dae = self.system.dae

        # constant short names
        In = spdiag([1] * dae.n)
        h = self.h

        while self.err > config.tol and self.niter < config.maxit:
            if self.t - self.t_jac >= 5:
                dae.rebuild = True
                self.t_jac = self.t
            elif self.niter > 4:
                dae.rebuild = True
            elif dae.factorize:
                dae.rebuild = True

            # rebuild Jacobian
            if dae.rebuild:
                exec(system.call.int)
                dae.rebuild = False
            else:
                exec(system.call.int_fg)

            # complete Jacobian matrix dae.Ac
            if config.method == 'euler':
                dae.Ac = sparse(
                    [[In - h * dae.Fx, dae.Gx], [-h * dae.Fy, dae.Gy]],
                    'd')
                dae.q = dae.x - self.x0 - h * dae.f

            elif config.method == 'trapezoidal':
                dae.Ac = sparse([[In - h * 0.5 * dae.Fx, dae.Gx],
                                 [-h * 0.5 * dae.Fy, dae.Gy]], 'd')
                dae.q = dae.x - self.x0 - h * 0.5 * (dae.f + self.f0)

            # windup limiters
            dae.reset_Ac()

            if dae.factorize:
                self.F = self.solver.symbolic(dae.Ac)
                dae.factorize = False
            self.inc = -matrix([dae.q, dae.g])

            try:
                N = self.solver.numeric(dae.Ac, self.F)
                self.solver.solve(dae.Ac, self.F, N, self.inc)
            except ArithmeticError:
                logger.error('Singular matrix')
                dae.check_diag(dae.Gy, 'unamey')
                dae.check_diag(dae.Fx, 'unamex')
                # force quit
                self.niter = config.maxit + 1
                break
            except ValueError:
                logger.warning('Unexpected symbolic factorization')
                dae.factorize = True
                continue
            else:
                inc_x = self.inc[:dae.n]
                inc_y = self.inc[dae.n:dae.m + dae.n]
                dae.x += inc_x
                dae.y += inc_y

            self.err = max(abs(self.inc))
            if np.isnan(self.inc).any():
                logger.error('Iteration error: NaN detected.')
                self.niter = config.maxit + 1
                break

            self.niter += 1

        if self.niter <= config.maxit:
            self.convergence = True
1287251
def _run_opf(self, P, q, AA, ll, uu, xmin, xmax, x0, opt):
        
        AAcvx = tocvx(AA)
        nx = x0.shape[0] # number of variables
        # add var limits to linear constraints
        eyex = spmatrix(1.0, range(nx), range(nx))
        A = sparse([eyex, AAcvx])
        l = matrix([-xmin, ll])
        u = matrix([ xmax, uu])

        Ae, be, Ai, bi = split_linear_constraints(A, l, u)


        if len(P.V) > 0:
            cvx_sol = solvers.qp(P, q, Ai, bi, Ae, be, self.solver, {"x": x0})
        else:
            cvx_sol = solvers.lp(q, Ai, bi, Ae, be, self.solver, {"x": x0})

        return cvx_sol
212602
def xstep(self):
        
        """

        self.YU[:] = self.Y - self.U

        b = self.DSf + self.rho*sl.rfftn(self.YU, None, self.cri.axisN)
        if self.cri.Cd == 1:
            self.Xf[:] = sl.solvedbi_sm(self.Df, self.mu + self.rho,
                                        b, self.c, self.cri.axisM)
        else:
            self.Xf[:] = sl.solvemdbi_ism(self.Df, self.mu + self.rho, b,
                                          self.cri.axisM, self.cri.axisC)

        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            Dop = lambda x: sl.inner(self.Df, x, axis=self.cri.axisM)
            if self.cri.Cd == 1:
                DHop = lambda x: np.conj(self.Df) * x
            else:
                DHop = lambda x: sl.inner(np.conj(self.Df), x,
                                          axis=self.cri.axisC)
            ax = DHop(Dop(self.Xf)) + (self.mu + self.rho)*self.Xf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
437481
@Override
    public boolean setA(DMatrixRMaj A) {
        if( A.numRows < A.numCols )
            throw new IllegalArgumentException("Can't solve for wide systems.  More variables than equations.");
        if( A.numRows > maxRows || A.numCols > maxCols )
            setMaxSize(A.numRows,A.numCols);

        R.reshape(A.numCols,A.numCols);
        a.reshape(A.numRows,1);
        temp.reshape(A.numRows,1);

        _setA(A);
        if( !decomposer.decompose(A) )
            return false;

        gammas = decomposer.getGammas();
        QR = decomposer.getQR();
        decomposer.getR(R,true);
        return true;
    }