Compute the least-squares fit between two arrays.
        A Python translation of 'FITLIN' from 'drutil.f' (Drizzle V2.9).
188936
def fitlin(imgarr,refarr):
    
    # Initialize variables
    _mat = np.zeros((3,3),dtype=np.float64)
    _xorg = imgarr[0][0]
    _yorg = imgarr[0][1]
    _xoorg = refarr[0][0]
    _yoorg = refarr[0][1]
    _sigxox = 0.
    _sigxoy = 0.
    _sigxo = 0.
    _sigyox = 0.
    _sigyoy = 0.
    _sigyo = 0.

    _npos = len(imgarr)
    # Populate matrices
    for i in range(_npos):
        _mat[0][0] += np.power((imgarr[i][0] - _xorg),2)
        _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg)
        _mat[0][2] += (imgarr[i][0] - _xorg)
        _mat[1][1] += np.power((imgarr[i][1] - _yorg),2)
        _mat[1][2] += imgarr[i][1] - _yorg

        _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg)
        _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg)
        _sigxo += refarr[i][0] - _xoorg
        _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg)
        _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg)
        _sigyo += refarr[i][1] - _yoorg

    _mat[2][2] = _npos
    _mat[1][0] = _mat[0][1]
    _mat[2][0] = _mat[0][2]
    _mat[2][1] = _mat[1][2]

    # Now invert this matrix
    _mat = linalg.inv(_mat)

    _a  = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2]
    _b  = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2])
    #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2]

    _c  = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2]
    _d  = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2]
    #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2]

    _xt = _xoorg - _a*_xorg+_b*_yorg
    _yt = _yoorg - _d*_xorg-_c*_yorg

    return [_a,_b,_xt],[_c,_d,_yt]
950122
def fit_null(Y,S_XX,U_XX,nfile,F):
    
    mtSet = limix.MTSet(Y, S_R=S_XX, U_R=U_XX, F=F)

    RV = mtSet.fitNull(cache=False)
    params = np.array([RV['params0_g'],RV['params0_n']])
    np.savetxt(nfile+'.p0',params)
    np.savetxt(nfile+'.nll0',RV['NLL0'])
    np.savetxt(nfile+'.cg0',RV['Cg'])
    np.savetxt(nfile+'.cn0',RV['Cn'])
855283
def fit_texture(layer):
    
    x, y = layer
    x = (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x))
    y = (y - np.nanmin(y)) / (np.nanmax(y) - np.nanmin(y))
    return x, y
965000
def _fit_s(D, a0, logp, tol=1e-7, maxiter=1000):
    ''''''
    N, K = D.shape
    s1 = a0.sum()
    m = a0 / s1
    mlogp = (m*logp).sum()
    for i in xrange(maxiter):
        s0 = s1
        g = psi(s1) - (m*psi(s1*m)).sum() + mlogp
        h = _trigamma(s1) - ((m**2)*_trigamma(s1*m)).sum()

        if g + s1 * h < 0:
            s1 = 1/(1/s0 + g/h/(s0**2))
        if s1 <= 0:
            s1 = s0 * exp(-g/(s0*h + g)) # Newton on log s
        if s1 <= 0:
            s1 = 1/(1/s0 + g/((s0**2)*h + 2*s0*g)) # Newton on 1/s
        if s1 <= 0:
            s1 = s0 - g/h # Newton
        if s1 <= 0:
            raise Exception('Unable to update s from {}'.format(s0))

        a = s1 * m
        if abs(s1 - s0) < tol:
            return a

    raise Exception('Failed to converge after {} iterations, s is {}'
            .format(maxiter, s1))
1167708
def fit(self, X, y=None):
        
        self.unhasher.fit(self._get_terms_iter(X))
        return self
47670
def _evaluate_fits(self, x, fit_key):
        
        fit = self.fits[fit_key]
        fit_params = self._get_fit_params(np.copy(x), fit_key)
        if type(fit) == list:
            res = []
            for i in range(len(fit)):
                res.append(fit[i](fit_params))
            return np.array(res)
        else:
            return fit(fit_params)
758850
def fit(self, x, y):
        
        verbose = self.verbose
        # Create temporary files
        data_file = tempfile.NamedTemporaryFile("w+b", delete=False)
        label_file = tempfile.NamedTemporaryFile("w+b", delete=False)

        start = time.time()
        raw_rules = categorical2pysbrl_data(x, y, data_file.name, label_file.name, supp=self.min_support,
                                            zmin=self.min_rule_len, zmax=self.max_rule_len, method=self.fim_method)
        if verbose > 1:
            print("Info: sbrl data files saved to %s and %s temporarily" % (data_file.name, label_file.name))
        data_file.close()
        label_file.close()
        cat_time = time.time() - start
        if verbose:
            print("Info: time for rule mining: %.4fs" % cat_time)
        n_labels = int(np.max(y)) + 1
        start = time.time()
        _model = train_sbrl(data_file.name, label_file.name, self.lambda_, eta=self.eta,
                            max_iters=self.iters, n_chains=self.n_chains, seed=self.seed,
                            alpha=self.alpha, verbose=verbose)
        train_time = time.time() - start
        if verbose:
            print("Info: training time: %.4fs" % train_time)

        # update model parameters
        self._n_classes = n_labels
        self._n_features = x.shape[1]

        # convert the raw parameters to rules
        self.from_raw(_model[0], _model[1], raw_rules)
        self._supports = self.compute_support(x, y)

        # Close the temp files
        os.unlink(data_file.name)
        os.unlink(label_file.name)
975472
def map_fit(interface, state, label, inp):
    
    import numpy as np
    out = interface.output(0)

    H, J, grad = 0, 0, 0
    for row in inp:
        row = row.strip().split(state["delimiter"])  # split row
        if len(row) > 1:  # check if row is empty
            # add intercept term to every sample
            x = np.array([1] + [(0 if v in state["missing_vals"] else float(v)) for i, v in enumerate(row) if
                                i in state["X_indices"]])
            # map label value to 0 or 1. If label does not match set error
            y = 0 if state["y_map"][0] == row[state["y_index"]] else 1 if state["y_map"][1] == row[
                state["y_index"]] else "Error"

            g = 1. / (1 + np.exp(-np.dot(x, state["thetas"])))  # sigmoid function
            grad += x * (g - y)  # gradient
            H += np.multiply(np.outer(x, x), g * (1 - g))  # Hessian matrix
            J -= np.log(g) if y == 1 else np.log(1 - g)  # J cost function
    out.add("grad", grad)
    out.add("J", J)
    for i, row in enumerate(H):
        out.add(i, row)
1035235
def fit(self, X, y, n_iter=None):
        """"""
        self.n_iter = self.n_iter if n_iter is None else n_iter
        X = getattr(X, 'values', X).reshape(len(X), 1)
        X_1 = self.homogenize(X)
        for i in range(self.n_iter):
            for i in range(0, len(X), 10):  # minibatch learning for numerical stability
                batch = slice(i, min(i + 10, len(X)))
                Xbatch, ybatch = X[batch, :], y[batch]
                X_1_batch = X_1[batch, :]
                self.W += (self.alpha / len(X) ** 1.5) * (
                    self.delta(Xbatch, ybatch).reshape((len(Xbatch), 1)).T.dot(X_1_batch))
        return self
1113095
def _fit_null(self, Y, X_base, scan_onsets=None):
        
        n_V = np.size(Y, axis=1)
        n_T = np.size(Y, axis=0)
        t_start = time.time()

        D, F, run_TRs, n_run = self._prepare_DF(
            n_T, scan_onsets=scan_onsets)
        YTY_diag = np.sum(Y * Y, axis=0)
        YTDY_diag = np.sum(Y * np.dot(D, Y), axis=0)
        YTFY_diag = np.sum(Y * np.dot(F, Y), axis=0)

        tol = self.tol
        n_iter = self.n_iter
        logger.info('Fitting null model'
                    ' for maximum {} times'.format(n_iter))

        # Add DC components capturing run-specific baselines.
        X_DC = self._gen_X_DC(run_TRs)
        X_DC, X_base, idx_DC = self._merge_DC_to_base(
            X_DC, X_base, no_DC=False)
        X_res = None
        param0 = np.zeros(n_V)
        for it in range(0, n_iter):
            if X_res is None:
                X0 = X_base
            else:
                X0 = np.concatenate((X_base, X_res), axis=1)
            n_X0 = X0.shape[1]
            X0TX0, X0TDX0, X0TFX0 = self._make_templates(D, F, X0, X0)
            X0TY, X0TDY, X0TFY = self._make_templates(D, F, X0, Y)
            res_null = scipy.optimize.minimize(
                self._loglike_AR1_null, param0, args=(
                    YTY_diag, YTDY_diag, YTFY_diag,
                    X0TX0, X0TDX0, X0TFX0, X0TY, X0TDY, X0TFY,
                    n_T, n_V, n_run, n_X0),
                method=self.optimizer, jac=True, tol=tol,
                options=self.minimize_options)
            param_change = res_null.x - param0
            param0 = res_null.x.copy()
            est_rho1_AR1_null = 2.0 / np.pi * np.arctan(param0)
            if self.auto_nuisance:
                X0TAX0 = X0TX0[None, :, :] \
                    - est_rho1_AR1_null[:, None, None] \
                    * X0TDX0[None, :, :] \
                    + est_rho1_AR1_null[:, None, None]**2 \
                    * X0TFX0[None, :, :]
                # dimension: space*#baseline*#baseline
                X0TAY = self._make_ar1_quad_form(X0TY, X0TDY, X0TFY,
                                                 est_rho1_AR1_null)
                # dimension: #baseline*space
                beta0s = np.linalg.solve(X0TAX0, X0TAY.T).T
                residuals = Y - np.dot(X_base, beta0s[:np.shape(X_base)[1], :])
                X_res = self.nureg_method(self.n_nureg_).fit_transform(
                    self.preprocess_residual(residuals))
            if np.max(np.abs(param_change)) < self.tol:
                logger.info('The change of parameters is smaller than '
                            'the tolerance value {}. Fitting is finished '
                            'after {} iterations'.format(self.tol, it + 1))
                break
        X0TAX0 = X0TX0[None, :, :] \
            - est_rho1_AR1_null[:, None, None] \
            * X0TDX0[None, :, :] \
            + est_rho1_AR1_null[:, None, None]**2 \
            * X0TFX0[None, :, :]
        # dimension: space*#baseline*#baseline
        X0TAY = self._make_ar1_quad_form(X0TY, X0TDY, X0TFY,
                                         est_rho1_AR1_null)
        # dimension: #baseline*space
        est_beta0_AR1_null = np.linalg.solve(X0TAX0, X0TAY.T).T

        YTAY = self._make_ar1_quad_form(YTY_diag, YTDY_diag, YTFY_diag,
                                        est_rho1_AR1_null)
        # dimension: space,
        YTAcorrY = YTAY - np.sum(X0TAY * est_beta0_AR1_null, axis=0)
        # dimension: space,
        est_sigma_AR1_null = (YTAcorrY / (n_T - n_X0)) ** 0.5
        if idx_DC.size > 1:
            collapsed_DC = np.sum(X0[:, idx_DC], axis=1)
            X0 = np.insert(np.delete(X0, idx_DC, axis=1), 0,
                           collapsed_DC, axis=1)
            collapsed_beta0 = np.mean(est_beta0_AR1_null[idx_DC, :], axis=0)
            est_beta0_AR1_null = np.insert(
                np.delete(est_beta0_AR1_null, idx_DC, axis=0),
                0, collapsed_beta0, axis=0)
        t_finish = time.time()
        logger.info(
            'total time of fitting: {} seconds'.format(t_finish - t_start))
        return est_beta0_AR1_null, est_sigma_AR1_null, est_rho1_AR1_null, X0
53289
def min_func(fun, bounds_min, bounds_max, d=None, rmax=10,
             n0=64, nswp=10, verb=True, smooth_fun=None):
    """"""
    if d is None:
        d = len(bounds_min)
        a = np.asanyarray(bounds_min).copy()
        b = np.asanyarray(bounds_max).copy()
    else:
        a = np.ones(d) * bounds_min
        b = np.ones(d) * bounds_max

    if smooth_fun is None:
        smooth_fun = lambda p, lam: (math.pi / 2 - np.arctan(p - lam))
    #smooth_fun = lambda p, lam: np.exp(-10*(p - lam))

    # We do not need to store the cores, only the interfaces!
    Rx = [[]] * (d + 1)  # Python list for the interfaces
    Rx[0] = np.ones((1, 1))
    Rx[d] = np.ones((1, 1))
    Jy = [np.empty(0, dtype=np.int)] * (d + 1)
    ry = rmax * np.ones(d + 1, dtype=np.int)
    ry[0] = 1
    ry[d] = 1
    n = n0 * np.ones(d, dtype=np.int)
    fun_evals = 0

    grid = [np.reshape(np.linspace(a[i], b[i], n[i]), (n[i], 1))
            for i in xrange(d)]
    for i in xrange(d - 1):
        #cr1 = y[i]
        ry[i + 1] = min(ry[i + 1], n[i] * ry[i])
        cr1 = np.random.randn(ry[i], n[i], ry[i + 1])
        cr1 = reshape(cr1, (ry[i] * n[i], ry[i + 1]))
        q, r = np.linalg.qr(cr1)
        ind = maxvol(q)
        w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
        w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
        Jy[i + 1] = np.hstack((w1, w2))
        Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
        Jy[i + 1] = Jy[i + 1][ind, :]

        # Jy{i+1} = [kron(ones(n(i),1), Jy{i}), kron((1:n(i))', ones(ry(i),1))];
        # Jy{i+1} = Jy{i+1}(ind,:);
    swp = 0
    dirn = -1
    i = d - 1
    lm = float('Inf')
    while swp < nswp:
        # Right-to-left sweep
        # The idea: compute the current core; compute the function of it;
        # Shift locally or globally? Local shift would be the first try
        # Compute the current core

        if np.size(Jy[i]) == 0:
            w1 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w1 = mkron(np.ones((n[i] * ry[i + 1], 1), dtype=np.int), Jy[i])
        w2 = mkron(mkron(np.ones((ry[i + 1], 1), dtype=np.int),
                         grid[i]), np.ones((ry[i], 1), dtype=np.int))
        if np.size(Jy[i + 1]) == 0:
            w3 = np.zeros((ry[i] * n[i] * ry[i + 1], 0), dtype=np.int)
        else:
            w3 = mkron(Jy[i + 1], np.ones((ry[i] * n[i], 1), dtype=np.int))

        J = np.hstack((w1, w2, w3))
        # Just add some random indices to J, which is rnr x d, need to make rn (r + r0) x add,
        # i.e., just generate random r, random n and random multiindex

        cry = fun(J)
        fun_evals += cry.size
        cry = reshape(cry, (ry[i], n[i], ry[i + 1]))
        min_cur = np.min(cry.flatten("F"))
        ind_cur = np.argmin(cry.flatten("F"))
        if lm > min_cur:
            lm = min_cur
            x_full = J[ind_cur, :]
            val = fun(x_full)
            if verb:
                print('New record:', val, 'Point:', x_full, 'fevals:', fun_evals)
        cry = smooth_fun(cry, lm)
        if (dirn < 0 and i > 0):
            cry = reshape(cry, (ry[i], n[i] * ry[i + 1]))
            cry = cry.T
            #q, r = np.linalg.qr(cry)
            u, s, v = mysvd(cry, full_matrices=False)
            ry[i] = min(ry[i], rmax)
            q = u[:, :ry[i]]
            ind = rect_maxvol(q)[0]  # maxvol(q)
            ry[i] = ind.size
            w1 = mkron(np.ones((ry[i + 1], 1), dtype=np.int), grid[i])
            if np.size(Jy[i + 1]) == 0:
                w2 = np.zeros((n[i] * ry[i + 1], 0), dtype=np.int)
            else:
                w2 = mkron(Jy[i + 1], np.ones((n[i], 1), dtype=np.int))
            Jy[i] = np.hstack((w1, w2))
            Jy[i] = reshape(Jy[i], (n[i] * ry[i + 1], -1))
            Jy[i] = Jy[i][ind, :]

        if (dirn > 0 and i < d - 1):
            cry = reshape(cry, (ry[i] * n[i], ry[i + 1]))
            q, r = np.linalg.qr(cry)
            #ind = maxvol(q)
            ind = rect_maxvol(q)[0]
            ry[i + 1] = ind.size
            w1 = mkron(np.ones((n[i], 1), dtype=np.int), Jy[i])
            w2 = mkron(grid[i], np.ones((ry[i], 1), dtype=np.int))
            Jy[i + 1] = np.hstack((w1, w2))
            Jy[i + 1] = reshape(Jy[i + 1], (ry[i] * n[i], -1))
            Jy[i + 1] = Jy[i + 1][ind, :]

        i += dirn
        if i == d or i == -1:
            dirn = -dirn
            i += dirn
            swp = swp + 1
    return val, x_full
1032425
def fit_transform(self, X, y, step_size=0.1, init_weights=None, warm_start=False):
        """"""
        self.fit(X=X, y=y, step_size=step_size, init_weights=init_weights, warm_start=warm_start)

        return self.transform(X=X)
1223021
def fit(self, X, y=None, **fit_params):
        
        if self.transform_type not in self.valid_transforms:
                warnings.warn("Invalid transform type.", stacklevel=2)
        return self
492323
def fit_predict(self, X, y=None, **kwargs):
        
        return self.fit(X, **kwargs).predict(X, **kwargs)
580011
def uflat(self):
          # noqa
        aligned_arrays = [self[key].values[[(slice(None)
                                             if c in coords
                                             else None)
                                            for c in self._coords]].T
                          for key, coords in self.dependent_variables_info]
        return np.vstack(aligned_arrays).flatten("F")
918935
def fit(self, x, y, **kwargs):
        
        self.mean = numpy.mean(y)
        return {}
308328
def compute_stats2(arrayNR, stats, weights):
    
    newshape = list(arrayNR.shape)
    if newshape[1] != len(weights):
        raise ValueError('Got %d weights but %d values!' %
                         (len(weights), newshape[1]))
    newshape[1] = len(stats)  # number of statistical outputs
    newarray = numpy.zeros(newshape, arrayNR.dtype)
    data = [arrayNR[:, i] for i in range(len(weights))]
    for i, func in enumerate(stats):
        newarray[:, i] = apply_stat(func, data, weights)
    return newarray