Compute the Deterministic Shared Response Model

        Parameters
        ----------
        X : list of 2D arrays, element i has shape=[voxels_i, samples]
            Each element in the list contains the fMRI data of one subject.

        y : not used
770751
def fit(self, X, y=None):
        
        logger.info('Starting Deterministic SRM')

        # Check the number of subjects
        if len(X) <= 1:
            raise ValueError("There are not enough subjects "
                             "({0:d}) to train the model.".format(len(X)))

        # Check for input data sizes
        if X[0].shape[1] < self.features:
            raise ValueError(
                "There are not enough samples to train the model with "
                "{0:d} features.".format(self.features))

        # Check if all subjects have same number of TRs
        number_trs = X[0].shape[1]
        number_subjects = len(X)
        for subject in range(number_subjects):
            assert_all_finite(X[subject])
            if X[subject].shape[1] != number_trs:
                raise ValueError("Different number of samples between subjects"
                                 ".")

        # Run SRM
        self.w_, self.s_ = self._srm(X)

        return self
669111
def MeanVarEstEQ(y, x, covariates, tol=1e-8):
    

    pcount = covariates.shape[0] + 2
    N = y.shape[0]
    beta_count = pcount * 2
    X = numpy.ones((N, pcount))
    X[:, 1] = x

    aprime = [1, 0]
    idx = 2
    for cov in covariates:
        aprime.append(-numpy.mean(cov)/numpy.std(cov))
        X[:, idx] = cov
        idx += 1
    aprime = numpy.matrix(aprime)

    # http://stackoverflow.com/questions/7267226/range-for-floats
    def frange(x, y, jump):
        while x < y:
            yield x
            x += jump
        yield 0

    class PhiReturn(object):
        def __init__(self, phi, dtheta):
            self.phi = phi
            self.dtheta = dtheta

    def dot_diag(x, y):
        """This is a conveniance function to perform the following \
        in a more efficient manner:

        :param x: arr1
        :param y: arr2

        x.dot(numpy.diag(y))

        y must be a single dimensioned array

        """
        if len(y.shape) != 1:
            print >> sys.stderr, "You can't pass an array of shape %s to dot_diag" % (y.shape)
            sys.exit(1)

        result = numpy.empty(x.shape)
        for i in range(0, x.shape[0]):
            result[i] = x[i] * y

        return result

    def diag_dot(x, y):
        """diagonal multiply
        """
        result = numpy.empty(y.shape)
        for i in range(0, y.shape[0]):
            result[i] = y[i] * x[i]
        return result

    def Phi(theta):
        MM = y - numpy.dot(X, theta[0:pcount])
        SS = numpy.exp(numpy.dot(-X, theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        Phi = numpy.hstack((DD1.dot(X), (DD2-0.5).dot(X)))

        tX = X.transpose()

        t1 = numpy.empty(tX.shape)
        t2 = numpy.empty(tX.shape)
        t3 = numpy.empty(tX.shape)
        for i in range(0,tX.shape[0]):
            t1[i] = tX[i]*-SS
            t2[i] = tX[i]*-DD1
            t3[i] = tX[i]*-DD2
        t1 = t1.dot(X)
        t2 = t2.dot(X)
        t3 = t3.dot(X)

        dtheta = numpy.hstack((numpy.vstack((t1, t2)), numpy.vstack((t2, t3))))
        return PhiReturn(Phi, dtheta)



    class MvSolveReturn(object):
        def __init__(self, theta, dtheta):
            self.theta = theta
            self.dtheta = dtheta

    def MVsolve(theta_new):
        solution_found = False

        mvsolve_iterations = 0

        while not solution_found:
            theta_old = theta_new.copy()
            tmp = Phi(theta_old)
            #print "ITR", mvsolve_iterations, numpy.mean(theta_new),numpy.sum(numpy.absolute(theta_new - theta_old))
            theta_new = theta_old - scipy.linalg.solve(tmp.dtheta, tmp.phi)
            mvsolve_iterations += 1

            if (numpy.sum(numpy.absolute(theta_new - theta_old)) < tol):
                tmp = Phi(theta_new)
                solution_found = True
            else:
                if mvsolve_iterations > 25000:
                    #print >> sys.stderr, mvsolve_iterations, "failures"
                    raise UnsolvedLocus("")

        return MvSolveReturn(theta_new, tmp.dtheta), mvsolve_iterations

    def MVcalcB(theta):
        MM = y - X.dot(theta[0:pcount])
        SS = numpy.exp(-X.dot(theta[pcount:]))
        DD1 = MM * SS
        DD2 = 0.5 * MM**2 * SS
        AA = numpy.hstack((diag_dot(DD1, X), diag_dot(DD2-0.5,X)))

        return numpy.transpose(AA).dot(AA)/ N


    mod = None
    itr = 0
    total_iterations = 0
    for i in frange(0.00, 1.0, 0.05):
        theta=numpy.empty((beta_count))
        theta[:] = i

        try:
            mod, iterations = MVsolve(theta)
            total_iterations += iterations
            if i > 0.05:
                print >> sys.stderr, "Completed: ", total_iterations, itr
            break
        except exceptions.ValueError as e:
            pass
        except numpy.linalg.linalg.LinAlgError as e:
            pass
        except Exception as inst:
            #print type(inst)
            pass
        itr += 1

    if not mod:
        raise UnsolvedLocus("")

    try:
        ainv = scipy.linalg.inv(mod.dtheta) * N
    except:
        raise ValueError("Singular Matrix Encountered")
    B = MVcalcB(mod.theta)
    V = ainv.dot(B).dot(ainv.transpose())



    # Focus on the two parameters of interest
    theta2 = numpy.array([mod.theta[1], mod.theta[pcount+1]])
    V2 = V[1:beta_count:pcount,1:beta_count:pcount]
    pvalt = 1 - scipy.stats.chi2.cdf(theta2.dot(scipy.linalg.inv(V2)).dot(theta2) * N, 2)


    ## From Chun's updated code:
    theta= mod.theta
    se = numpy.sqrt(numpy.diag(V)/N)
    pval = 2*scipy.stats.norm.cdf(-numpy.absolute(mod.theta/se))
    return pvalt, theta, pval, se, V/N
1239871
def linear_model(x, y):
        
        x1, x2 = x
        y1, y2 = y
        m = (y2 - y1) / (x2 - x1)
        b = y1 - (m * x1)
        return lambda x: m * x + b
109773
def _construct_sparse_features(self, x):
        
        I, J, K = x.shape
        new_array_height = (x != 0).sum(axis=2).max()
        index_array = -np.ones((I, J, new_array_height), dtype='int64')
        value_array = -np.ones((I, J, new_array_height), dtype='float64')
        populate_sparse_features(x, index_array, value_array, I, J, K)
        return index_array, value_array
470185
@Override
    public int predict(double[] x) {
        double minDist = Double.MAX_VALUE;
        int bestCluster = 0;

        int i = 0;
        for (Node neuron : nodes) {
            double dist = Math.squaredDistance(x, neuron.w);
            if (dist < minDist) {
                minDist = dist;
                bestCluster = i;
            }
            i++;
        }

        if (y == null || y.length != nodes.size()) {
            return bestCluster;
        } else {

            return y[bestCluster];
        }
    }
57512
def K(self, X, X2=None):
                """"""        
                X,slices = X[:,:-1],index_to_slices(X[:,-1])
                if X2 is None:
                        X2,slices2 = X,slices
                        K = np.zeros((X.shape[0], X.shape[0]))
                else:
                        X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
                        K = np.zeros((X.shape[0], X2.shape[0]))

                tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
                ttdist = (X[:,0][:,None] - X2[:,0][None,:])
                
                vyt = self.variance_Yt
                
                lyt=1/(2*self.lengthscale_Yt)

                a = -self.a
                c = self.c

                kyy = lambda tdist: np.exp(-lyt*(tdist))

                k1 = lambda tdist: (2*lyt - 4*lyt**2 *(tdist) )

                k4 = lambda tdist: 2*lyt*(tdist)

                for i, s1 in enumerate(slices):
                        for j, s2 in enumerate(slices2):
                                for ss1 in s1:
                                    for ss2 in s2:
                                        if i==0 and j==0:
                                            K[ss1,ss2] = vyt*kyy(tdist[ss1,ss2])
                                        elif i==0 and j==1:
                                            K[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                            #K[ss1,ss2] = (2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                        elif i==1 and j==1:
                                            K[ss1,ss2] = ( k1(tdist[ss1,ss2]) + 1. )*vyt* kyy(tdist[ss1,ss2])+self.ubias
                                        else:
                                            K[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                            #K[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                #stop
                return K
1094541
def compute_samples(self):
        
        shape = [self.batch_size, self.eq_samples, self.iw_samples, self.num_latent]
        mu = tf.nn.sigmoid(self.mu)
        z = tf.select(tf.random_uniform(shape) - mu > 0, tf.ones(shape), tf.zeros(shape))
        return z
1302048
def set_data_matrix(self, X):
        
        X = check_array(X, accept_sparse=sparse_formats)
        self.X = X
557439
def _check_forward_mode_input_array(self, X: np.ndarray) -> int:
        
        # Find the length of each variable to infer T
        if not isinstance(X, np.ndarray):
            raise ValueError('X must be a numpy array, dict, or scalar')
        # Get the shape and tensor rank
        shape = X.shape
        tensor_rank = len(shape)
        T = 0
        # Only 1D and 2D arrays are supported
        if tensor_rank not in (0, 1, 2):
            raise ValueError(f'Shape of X = {X.shape}. Numpy array must be a 1D vector or 2D matrix')        
        if tensor_rank == 0:
            T = 1
        # If the input was a 1D vector, its length must EITHER (1) be T, or (2) m, with T == 1
        if tensor_rank == 1 and (shape[0] != self.m) and self.m != 1:
            raise ValueError(f'Error: X has shape {X.shape}, incompatible with m = {self.m} on fluxion.')
        # Return the value of T in this situation
        if tensor_rank == 1 and shape[0] == self.m:
            T = 1
        if tensor_rank == 1 and self.m == 1:
            T = shape[0]
        # If the input was a 2D vector, it must be of shape Txn
        if tensor_rank == 2 and (shape[1] != self.m):
            raise ValueError(f'Error: X has shape {X.shape}, incompatible with m = {self.m} on fluxion.')
        if tensor_rank == 2:
            T = shape[0]
        return T
216480
def infer_x(self, y):
        
        assert len(y) == self.fmodel.dim_y, "Wrong dimension for y. Expected %i, got %i" % (self.fmodel.dim_y, len(y))
        if len(self.fmodel.dataset) == 0:
            return [[0.0]*self.dim_x]
        else:
            _, index = self.fmodel.dataset.nn_y(y, k=1)
            return [self.fmodel.dataset.get_x(index[0])]
135053
def fit(self, X):
        
        D = self._initialize(X)
        for i in range(self.max_iter):
            gamma = self._transform(D, X)
            e = np.linalg.norm(X - gamma.dot(D))
            if e < self.tol:
                break
            D, gamma = self._update_dict(X, D, gamma)

        self.components_ = D
        return self
895183
def compute(self, x):
        '''
        
        '''
        q_learning = copy(self.__greedy_q_learning)
        q_learning.epsilon_greedy_rate = x[0]
        q_learning.alpha_value = x[1]
        q_learning.gamma_value = x[2]
        if self.__init_state_key is not None:
            q_learning.learn(state_key=self.__init_state_key, limit=int(x[3]))
        else:
            q_learning.learn(limit=x[3])
        q_sum = q_learning.q_df.q_value.sum()
        if q_sum != 0:
            cost = q_learning.q_df.shape[0] / q_sum
        else:
            cost = q_learning.q_df.shape[0] / 1e-4

        return cost
1116330
def _transform(self, X):
        
        X = self._matrix_adjust(X)

        X = check_array(X, accept_sparse='csc', force_all_finite=False,
                        dtype=int)
        if X.min() < 0:
            raise ValueError("X needs to contain only non-negative integers.")
        n_samples, n_features = X.shape

        indices = self.feature_indices_
        if n_features != indices.shape[0] - 1:
            raise ValueError("X has different shape than during fitting."
                             " Expected %d, got %d."
                             % (indices.shape[0] - 1, n_features))

        # Replace all indicators which were below `minimum_fraction` in the
        # training set by 'other'
        if self.minimum_fraction is not None:
            for column in range(X.shape[1]):
                if sparse.issparse(X):
                    indptr_start = X.indptr[column]
                    indptr_end = X.indptr[column + 1]
                    unique = np.unique(X.data[indptr_start:indptr_end])
                else:
                    unique = np.unique(X[:, column])

                for unique_value in unique:
                    if unique_value not in self.do_not_replace_by_other_[column]:
                        if sparse.issparse(X):
                            indptr_start = X.indptr[column]
                            indptr_end = X.indptr[column + 1]
                            X.data[indptr_start:indptr_end][
                                X.data[indptr_start:indptr_end] ==
                                unique_value] = SPARSE_ENCODINGS['OTHER']
                        else:
                            X[:, column][X[:, column] == unique_value] = SPARSE_ENCODINGS['OTHER']

        if sparse.issparse(X):
            n_values_check = X.max(axis=0).toarray().flatten() + 1
        else:
            n_values_check = np.max(X, axis=0) + 1

        # Replace all indicators which are out of bounds by 'other' (index 0)
        if (n_values_check > self.n_values_).any():
            # raise ValueError("Feature out of bounds. Try setting n_values.")
            for i, n_value_check in enumerate(n_values_check):
                if (n_value_check - 1) >= self.n_values_[i]:
                    if sparse.issparse(X):
                        indptr_start = X.indptr[i]
                        indptr_end = X.indptr[i+1]
                        X.data[indptr_start:indptr_end][X.data[indptr_start:indptr_end] >= self.n_values_[i]] = 0
                    else:
                        X[:, i][X[:, i] >= self.n_values_[i]] = 0

        if sparse.issparse(X):
            row_indices = X.indices
            column_indices = []
            for i in range(len(X.indptr) - 1):
                nbr = X.indptr[i + 1] - X.indptr[i]
                column_indices_ = [indices[i]] * nbr
                column_indices_ += X.data[X.indptr[i]:X.indptr[i + 1]]
                column_indices.extend(column_indices_)
            data = np.ones(X.data.size)
        else:
            column_indices = (X + indices[:-1]).ravel()
            row_indices = np.repeat(np.arange(n_samples, dtype=np.int32),
                                    n_features)
            data = np.ones(n_samples * n_features)
        out = sparse.coo_matrix((data, (row_indices, column_indices)),
                                shape=(n_samples, indices[-1]),
                                dtype=self.dtype).tocsc()

        out = out[:, self.active_features_]
        return out.tocsr() if self.sparse else out.toarray()
990446
def XanyKXany(self):
        
        result = np.empty((self.P,self.F_any.shape[1],self.F_any.shape[1]), order='C')
        for p in range(self.P):
            X1D = self.Fstar_any * self.D[:,p:p+1]
            X1X2 = X1D.T.dot(self.Fstar_any)
            result[p] = X1X2
        return result
115856
private Map<Integer, Double> predict(final double[] x) {
        Map<Integer, Double> result = new HashMap<>();
        for (int i = 0; i < model.weights.length; i++) {
            double y = VectorUtils.dotProduct(x, model.weights[i]);
            y += model.bias[i];
            result.put(i, y);
        }

        return result;
    }
956023
def _call(self, x):
        """"""

        # Convert to array with most
        arr = self._asarray(x)
        svd_diag = np.linalg.svd(arr, compute_uv=False)

        # Rotate the axes so the svd-direction is first
        s_reordered = moveaxis(svd_diag, -1, 0)

        # Return nuclear norm
        return self.outernorm(self.pwisenorm(s_reordered))
1367421
def solve
      set_parameters
      # Vectors holding the values used for calculating the weights:
      a = NArray.float(@n)
      b = NArray.float(@n)
      # Set an initial estimate for the probabilities of true segmentation:
      @n.times do |i|
        @weights_current[i] = @decisions[i, true].mean
      end
      # Proceed iteratively until we have converged to a local optimum:
      k = 0
      while k < max_iterations do
        # Copy weights:
        @weights_previous = @weights_current.dup
        # E-step: Estimation of the conditional expectation of the complete data log likelihood function.
        # Deriving the estimator for the unobserved true segmentation (T).
        @n.times do |i|
          voxel_decisions = @decisions[i, true]
          # Find the rater-indices for this voxel where the raters' decisions equals 1 and 0:
          positive_indices, negative_indices = (voxel_decisions.eq 1).where2
          # Determine ai:
          # Multiply by corresponding sensitivity (or 1 - sensitivity):
          a_decision1_factor = (positive_indices.length == 0 ? 1 : @p[positive_indices].prod)
          a_decision0_factor = (negative_indices.length == 0 ? 1 : (1 - @p[negative_indices]).prod)
          a[i] = @weights_previous[i] * a_decision1_factor * a_decision0_factor
          # Determine bi:
          # Multiply by corresponding specificity (or 1 - specificity):
          b_decision0_factor = (negative_indices.length == 0 ? 1 : @q[negative_indices].prod)
          b_decision1_factor = (positive_indices.length == 0 ? 1 : (1 - @q[positive_indices]).prod)
          b[i] = @weights_previous[i] * b_decision0_factor * b_decision1_factor
          # Determine Wi: (take care not to divide by zero)
          if a[i] > 0 or b[i] > 0
            @weights_current[i] = a[i] / (a[i] + b[i])
          else
            @weights_current[i] = 0
          end
        end
        # M-step: Estimation of the performance parameters by maximization.
        # Finding the values of the expert performance level parameters that maximize the conditional expectation
        # of the complete data log likelihood function (phi - p,q).
        @r.times do |j|
          voxel_decisions = @decisions[true, j]
          # Find the voxel-indices for this rater where the rater's decisions equals 1 and 0:
          positive_indices, negative_indices = (voxel_decisions.eq 1).where2
          # Determine sensitivity:
          # Sum the weights for the indices where the rater's decision equals 1:
          sum_positive = (positive_indices.length == 0 ? 0 : @weights_current[positive_indices].sum)
          @p[j] = sum_positive / @weights_current.sum
          # Determine specificity:
          # Sum the weights for the indices where the rater's decision equals 0:
          sum_negative = (negative_indices.length == 0 ? 0 : (1 - @weights_current[negative_indices]).sum)
          @q[j] = sum_negative / (1 - @weights_current).sum
        end
        # Bump our iteration index:
        k += 1
        # Abort if we have reached the local optimum: (there is no change in the sum of weights)
        if @weights_current.sum - @weights_previous.sum == 0
#puts "Iteration aborted as optimum solution was found!" if @verbose
#logger.info("Iteration aborted as optimum solution was found!")
          break
        end
      end
      # Set the true segmentation:
      @true_segmentation_vector = @weights_current.round
      # Set the weights attribute:
      @weights = @weights_current
      # As this vector doesn't make much sense to the user, it must be converted to a volume. If volume reduction has
      # previously been performed, this must be taken into account when transforming it to a volume:
      construct_segmentation_volume
      # Construct a BinVolume instance for the true segmentation and add it as a master volume to the BinMatcher instance.
      update_bin_matcher
      # Set the phi variable:
      @phi[0, true] = @p
      @phi[1, true] = @q
    end