factoring for multivariate
422610
function(symbol, factors) {  
                if(symbol.group === FN) { 
                    if(symbol.fname === 'sqrt') {
                        var factors2 = new Factors(),
                            arg = __.Factor.common(symbol.args[0].clone(), factors2);
                        arg = __.Factor.coeffFactor(arg, factors2);
                        symbol = _.symfunction('sqrt', [arg]);
                        factors2.each(function(x) {
                            symbol = _.multiply(symbol, _.parse(core.Utils.format('sqrt({0})', x)));
                        });
                    }
                    else
                        factors.add(symbol);
                }
                else { 
                    //symbol = __.Factor.common(symbol, factors);
                    symbol = __.Factor.mSqfrFactor(symbol, factors);
                    
                    var vars = variables(symbol),
                        symbols = symbol.collectSymbols().map(function(x) {
                            return Symbol.unwrapSQRT(x);
                        }),
                        sorted = {},
                        maxes = {},
                        l = vars.length, n = symbols.length;
                    //take all the variables in the symbol and organize by variable name
                    //e.g. a^2+a^2+b*a -> {a: {a^3, a^2, b*a}, b: {b*a}}
                    
                    for(var i=0; i<l; i++) {
                        var v = vars[i];
                        sorted[v] = new Symbol(0);
                        for(var j=0; j<n; j++) {
                            var s = symbols[j];
                            if(s.contains(v)) {
                                var p = s.value === v ? s.power.toDecimal() : s.symbols[v].power.toDecimal();
                                if(!maxes[v] || p < maxes[v]) maxes[v] = p;
                                sorted[v] = _.add(sorted[v], s.clone());
                            }
                        }
                    }
                    
                    for(var x in sorted) {
                        var r = _.parse(x+'^'+maxes[x]); 
                        var div = _.divide(sorted[x], r);
                        var new_factor = _.expand(div); 
                        
                        if(new_factor.equals(1))
                            break; //why divide by one. Just move 
                        var divided = __.div(symbol.clone(), new_factor); 
                        if(divided[0].equals(0)) { 
                            //cant factor anymore
                            break;
                        }
                        
                        var neg_numeric_factor = isInt(new_factor) && new_factor.lessThan(0);
                        
                        if(divided[1].equals(0) && !neg_numeric_factor) { //we found at least one factor
                            //factors.add(new_factor);
                            var d = __.div(symbol.clone(), divided[0].clone());
                            var r = d[0];
                            symbol = d[1];
                            //we don't want to just flip the sign. If the remainder is -1 then we accomplished nothing
                            //and we just return the symbol;
                            //If r equals zero then there's nothing left to do so we're done
                            if(r.equals(-1))
                                return symbol;

                            var factor = divided[0]; 
                            if(symbol.equals(factor)) {
                                var rem = __.Factor.reduce(factor, factors);
                                
                                if(!symbol.equals(rem)) 
                                    return __.Factor.mfactor(rem, factors);
                            }
                            else {
                                factors.add(factor); 
                                //if the remainder of the symbol is zero then we're done. TODO: Rethink this logic a bit.
                                if(symbol.equals(0))
                                    return r;
                            }
                            
                            if(r.isConstant('all')) { 
                                factors.add(r);
                                return r;
                            }
                            
                            return __.Factor.mfactor(r, factors);
                        }
                    }
                    
                }
                
                //difference of squares factorization
                symbol = __.Factor.sqdiff(symbol, factors);
                //factors by fishing for zeroes
                symbol = __.Factor.zeroes(symbol, factors);
                
                return symbol;
            }
1348579
public Set<Integer> getAdjacentFactors(int factorNum) {
    Set<Integer> adjacentFactors = Sets.newHashSet();
    for (Integer variableNum : factorVariableMap.get(factorNum)) {
      adjacentFactors.addAll(variableFactorMap.get(variableNum));
    }
    return adjacentFactors;
  }
759096
private List<UnitFactor> list(UnitFactor factor) {
    List<UnitFactor> list = new ArrayList<>();
    list.add(factor);
    return list;
  }
529179
def factor_list_to_map(factors):
        
        ret = {}
        for k, g in itertools.groupby(factors):
            ret[k] = len(list(g))
        return ret
1184742
def get_all_monomials(variables, extramonomials, substitutions, degree,
                      removesubstitutions=True):
    
    monomials = get_monomials(variables, degree)
    if extramonomials is not None:
        monomials.extend(extramonomials)
    if removesubstitutions and substitutions is not None:
        monomials = [monomial for monomial in monomials if monomial not
                     in substitutions]
        monomials = [remove_scalar_factor(apply_substitutions(monomial,
                                                              substitutions))
                     for monomial in monomials]
    monomials = unique(monomials)
    return monomials
674639
public void computeBasis( int numComponents ) {
        if( numComponents > A.getNumCols() )
            throw new IllegalArgumentException("More components requested that the data's length.");
        if( sampleIndex != A.getNumRows() )
            throw new IllegalArgumentException("Not all the data has been added");
        if( numComponents > sampleIndex )
            throw new IllegalArgumentException("More data needed to compute the desired number of components");

        this.numComponents = numComponents;

        // compute the mean of all the samples
        for( int i = 0; i < A.getNumRows(); i++ ) {
            for( int j = 0; j < mean.length; j++ ) {
                mean[j] += A.get(i,j);
            }
        }
        for( int j = 0; j < mean.length; j++ ) {
            mean[j] /= A.getNumRows();
        }

        // subtract the mean from the original data
        for( int i = 0; i < A.getNumRows(); i++ ) {
            for( int j = 0; j < mean.length; j++ ) {
                A.set(i,j,A.get(i,j)-mean[j]);
            }
        }

        // Compute SVD and save time by not computing U
        SingularValueDecomposition<DMatrixRMaj> svd =
                DecompositionFactory_DDRM.svd(A.numRows, A.numCols, false, true, false);
        if( !svd.decompose(A) )
            throw new RuntimeException("SVD failed");

        V_t = svd.getV(null,true);
        DMatrixRMaj W = svd.getW(null);

        // Singular values are in an arbitrary order initially
        SingularOps_DDRM.descendingOrder(null,false,W,V_t,true);

        // strip off unneeded components and find the basis
        V_t.reshape(numComponents,mean.length,true);
    }
342838
def upsample(self, factor):
        """"""
        self.command.append("upsample")
        self.command.append(factor)
        return self
23054
def alternating_least_squares(Ciu, factors, **kwargs):
    
    log.warning("This method is deprecated. Please use the AlternatingLeastSquares"
                " class instead")

    model = AlternatingLeastSquares(factors=factors, **kwargs)
    model.fit(Ciu)
    return model.item_factors, model.user_factors
305165
def get_support(variables, polynomial):
    
    support = []
    if is_number_type(polynomial):
        support.append([0] * len(variables))
        return support
    for monomial in polynomial.expand().as_coefficients_dict():
        tmp_support = [0] * len(variables)
        mon, _ = __separate_scalar_factor(monomial)
        symbolic_support = flatten(split_commutative_parts(mon))
        for s in symbolic_support:
            if isinstance(s, Pow):
                base = s.base
                if is_adjoint(base):
                    base = base.adjoint()
                tmp_support[variables.index(base)] = s.exp
            elif is_adjoint(s):
                tmp_support[variables.index(s.adjoint())] = 1
            elif isinstance(s, (Operator, Symbol)):
                tmp_support[variables.index(s)] = 1
        support.append(tmp_support)
    return support
347096
function(x, y, z, w) {
      var perm = this.perm;
      var grad4 = this.grad4;

      var n0, n1, n2, n3, n4; // Noise contributions from the five corners
      // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
      var s = (x + y + z + w) * F4; // Factor for 4D skewing
      var i = Math.floor(x + s);
      var j = Math.floor(y + s);
      var k = Math.floor(z + s);
      var l = Math.floor(w + s);
      var t = (i + j + k + l) * G4; // Factor for 4D unskewing
      var X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
      var Y0 = j - t;
      var Z0 = k - t;
      var W0 = l - t;
      var x0 = x - X0; // The x,y,z,w distances from the cell origin
      var y0 = y - Y0;
      var z0 = z - Z0;
      var w0 = w - W0;
      // For the 4D case, the simplex is a 4D shape I won't even try to describe.
      // To find out which of the 24 possible simplices we're in, we need to
      // determine the magnitude ordering of x0, y0, z0 and w0.
      // Six pair-wise comparisons are performed between each possible pair
      // of the four coordinates, and the results are used to rank the numbers.
      var rankx = 0;
      var ranky = 0;
      var rankz = 0;
      var rankw = 0;
      if (x0 > y0) rankx++;
      else ranky++;
      if (x0 > z0) rankx++;
      else rankz++;
      if (x0 > w0) rankx++;
      else rankw++;
      if (y0 > z0) ranky++;
      else rankz++;
      if (y0 > w0) ranky++;
      else rankw++;
      if (z0 > w0) rankz++;
      else rankw++;
      var i1, j1, k1, l1; // The integer offsets for the second simplex corner
      var i2, j2, k2, l2; // The integer offsets for the third simplex corner
      var i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
      // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
      // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
      // impossible. Only the 24 indices which have non-zero entries make any sense.
      // We use a thresholding to set the coordinates in turn from the largest magnitude.
      // Rank 3 denotes the largest coordinate.
      i1 = rankx >= 3 ? 1 : 0;
      j1 = ranky >= 3 ? 1 : 0;
      k1 = rankz >= 3 ? 1 : 0;
      l1 = rankw >= 3 ? 1 : 0;
      // Rank 2 denotes the second largest coordinate.
      i2 = rankx >= 2 ? 1 : 0;
      j2 = ranky >= 2 ? 1 : 0;
      k2 = rankz >= 2 ? 1 : 0;
      l2 = rankw >= 2 ? 1 : 0;
      // Rank 1 denotes the second smallest coordinate.
      i3 = rankx >= 1 ? 1 : 0;
      j3 = ranky >= 1 ? 1 : 0;
      k3 = rankz >= 1 ? 1 : 0;
      l3 = rankw >= 1 ? 1 : 0;
      // The fifth corner has all coordinate offsets = 1, so no need to compute that.
      var x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords
      var y1 = y0 - j1 + G4;
      var z1 = z0 - k1 + G4;
      var w1 = w0 - l1 + G4;
      var x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in (x,y,z,w) coords
      var y2 = y0 - j2 + 2.0 * G4;
      var z2 = z0 - k2 + 2.0 * G4;
      var w2 = w0 - l2 + 2.0 * G4;
      var x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in (x,y,z,w) coords
      var y3 = y0 - j3 + 3.0 * G4;
      var z3 = z0 - k3 + 3.0 * G4;
      var w3 = w0 - l3 + 3.0 * G4;
      var x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in (x,y,z,w) coords
      var y4 = y0 - 1.0 + 4.0 * G4;
      var z4 = z0 - 1.0 + 4.0 * G4;
      var w4 = w0 - 1.0 + 4.0 * G4;
      // Work out the hashed gradient indices of the five simplex corners
      var ii = i & 255;
      var jj = j & 255;
      var kk = k & 255;
      var ll = l & 255;
      // Calculate the contribution from the five corners
      var t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
      if (t0 < 0) n0 = 0.0;
      else {
        var gi0 = (perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32) * 4;
        t0 *= t0;
        n0 = t0 * t0 * (grad4[gi0] * x0 + grad4[gi0 + 1] * y0 + grad4[gi0 + 2] * z0 + grad4[gi0 + 3] * w0);
      }
      var t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
      if (t1 < 0) n1 = 0.0;
      else {
        var gi1 = (perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32) * 4;
        t1 *= t1;
        n1 = t1 * t1 * (grad4[gi1] * x1 + grad4[gi1 + 1] * y1 + grad4[gi1 + 2] * z1 + grad4[gi1 + 3] * w1);
      }
      var t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
      if (t2 < 0) n2 = 0.0;
      else {
        var gi2 = (perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32) * 4;
        t2 *= t2;
        n2 = t2 * t2 * (grad4[gi2] * x2 + grad4[gi2 + 1] * y2 + grad4[gi2 + 2] * z2 + grad4[gi2 + 3] * w2);
      }
      var t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
      if (t3 < 0) n3 = 0.0;
      else {
        var gi3 = (perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32) * 4;
        t3 *= t3;
        n3 = t3 * t3 * (grad4[gi3] * x3 + grad4[gi3 + 1] * y3 + grad4[gi3 + 2] * z3 + grad4[gi3 + 3] * w3);
      }
      var t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
      if (t4 < 0) n4 = 0.0;
      else {
        var gi4 = (perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32) * 4;
        t4 *= t4;
        n4 = t4 * t4 * (grad4[gi4] * x4 + grad4[gi4 + 1] * y4 + grad4[gi4 + 2] * z4 + grad4[gi4 + 3] * w4);
      }
      // Sum up and scale the result to cover the range [-1,1]
      return 27.0 * (n0 + n1 + n2 + n3 + n4);
    }
86206
def factorization_machine_model(factor_size, num_features,
                                lr_mult_config, wd_mult_config, init_config):
    
    x = mx.symbol.Variable("data", stype='csr')
    # factor, linear and bias terms
    v = mx.symbol.Variable("v", shape=(num_features, factor_size), stype='row_sparse',
                           init=init_config['v'], lr_mult=lr_mult_config['v'],
                           wd_mult=wd_mult_config['v'])
    w = mx.symbol.Variable('w', shape=(num_features, 1), stype='row_sparse',
                      init=init_config['w'], lr_mult=lr_mult_config['w'],
                      wd_mult=wd_mult_config['w'])
    w0 = mx.symbol.Variable('w0', shape=(1,), init=init_config['w0'],
                       lr_mult=lr_mult_config['w0'], wd_mult=wd_mult_config['w0'])
    w1 = mx.symbol.broadcast_add(mx.symbol.dot(x, w), w0)

    # squared terms for subtracting self interactions
    v_s = mx.symbol._internal._square_sum(data=v, axis=1, keepdims=True)
    x_s = x.square()
    bd_sum = mx.sym.dot(x_s, v_s)

    # interactions
    w2 = mx.symbol.dot(x, v)
    w2_squared = 0.5 * mx.symbol.square(data=w2)

    # putting everything together
    w_all = mx.symbol.Concat(w1, w2_squared, dim=1)
    sum1 = w_all.sum(axis=1, keepdims=True)
    sum2 = -0.5 * bd_sum
    model = sum1 + sum2

    y = mx.symbol.Variable("softmax_label")
    model = mx.symbol.LogisticRegressionOutput(data=model, label=y)
    return model
203750
def __push_facvar_sparse(self, polynomial, block_index, row_offset, i, j):
        
        width = self.block_struct[block_index - 1]
        # Preprocess the polynomial for uniform handling later
        # DO NOT EXPAND THE POLYNOMIAL HERE!!!!!!!!!!!!!!!!!!!
        # The simplify_polynomial bypasses the problem.
        # Simplifying here will trigger a bug in SymPy related to
        # the powers of daggered variables.
        # polynomial = polynomial.expand()
        if is_number_type(polynomial) or polynomial.is_Mul:
            elements = [polynomial]
        else:
            elements = polynomial.as_coeff_mul()[1][0].as_coeff_add()[1]
        # Identify its constituent monomials
        for element in elements:
            results = self._get_index_of_monomial(element)
            # k identifies the mapped value of a word (monomial) w
            for (k, coeff) in results:
                if k > -1 and coeff != 0:
                    self.F[row_offset + i * width + j, k] += coeff
724666
private static int[] getFactors(int _value, int _max) {
      final int factors[] = new int[MAX_GROUP_SIZE];
      int factorIdx = 0;

      for (int possibleFactor = 1; possibleFactor <= _max; possibleFactor++) {
         if ((_value % possibleFactor) == 0) {
            factors[factorIdx++] = possibleFactor;
         }
      }

      return (Arrays.copyOf(factors, factorIdx));
   }
530827
def _process_monomial(self, monomial, n_vars):
        
        coeff, monomial = monomial.as_coeff_Mul()
        k = 0
        # Have we seen this monomial before?
        conjugate = False
        try:
            # If yes, then we improve sparsity by reusing the
            # previous variable to denote this entry in the matrix
            k = self.monomial_index[monomial]
        except KeyError:
            # An extra round of substitutions is granted on the conjugate of
            # the monomial if all the variables are Hermitian
            daggered_monomial = \
                apply_substitutions(Dagger(monomial), self.substitutions,
                                    self.pure_substitution_rules)
            try:
                k = self.monomial_index[daggered_monomial]
                conjugate = True
            except KeyError:
                # Otherwise we define a new entry in the associated
                # array recording the monomials, and add an entry in
                # the moment matrix
                k = n_vars + 1
                self.monomial_index[monomial] = k
        if conjugate:
            k = -k
        return k, coeff
474741
def cofactor_n(x, i, j):
    '''
    '''
    m = x.shape[1]
    mr = np.arange(m)
    i_idx = mr[mr != i]
    j_idx = mr[mr != j]
    return det_n(x[:, i_idx[:, np.newaxis],
                      j_idx[np.newaxis, :]])
176879
def factors(self, rank):
        
        self._check_rank(rank)
        return [result.factors for result in self.results[rank]]
1152866
private Pair<Factor, Assignment> getRegressableFactorForFamily(BoostingFactorFamily family, VariableMatch match,
      List<Integer> conditionedVariables, MarginalSet marginals) {
    VariableNumMap marginalVars = match.getMatchedVariables().removeAll(conditionedVariables);
    VariableNumMap fixedVars = match.getMatchedVariables().intersection(conditionedVariables);
    Factor factorMarginal = marginals.getMarginal(marginalVars.getVariableNums())
        .relabelVariables(match.getMappingToTemplate());

    Assignment assignment = marginals.getConditionedValues().intersection(fixedVars)
        .mapVariables(match.getMappingToTemplate().getVariableIndexReplacementMap());

    VariableNumMap outerProductVars = family.getUnconditionalVariables().removeAll(factorMarginal.getVars());
    if (outerProductVars.size() > 0) {
      factorMarginal = factorMarginal.outerProduct(TableFactor.pointDistribution(outerProductVars,
          assignment.intersection(outerProductVars)));
      assignment = assignment.removeAll(outerProductVars.getVariableNumsArray());
    }

    return Pair.of(factorMarginal, assignment);
  }