Python scipy.optimize.nnls() Examples

def _mge_1d(r_array, flux_r, N=20, linspace=False):
    """

    :param r_array:
    :param flux_r:
    :param N:
    :return:
    """
    if linspace is True:
        sigmas = np.linspace(r_array[0], r_array[-1] / 2, N + 2)[1:-1]
    else:
        sigmas = np.logspace(np.log10(r_array[0]), np.log10((r_array[-1] + 0.0000001) / 2.), N + 2)[1:-1]
    # sigmas = np.linspace(r_array[0], r_array[-1]/2, N + 2)[1:-1]

    A = np.zeros((len(flux_r), N))
    for j in np.arange(A.shape[1]):
        A[:, j] = gaussian(r_array, sigmas[j], 1.)
    amplitudes, norm = nnls(A, flux_r)
    return amplitudes, sigmas, norm 

def iter_solver(self, A, W, H, k, it):
        if not sps.issparse(A):
            for j in iter(range(0, H.shape[0])):
                res = opt.nnls(W, A[:, j])
                H[j, :] = res[0]
        else:
            for j in iter(range(0, H.shape[0])):
                res = opt.nnls(W, A[:, j].toarray()[:, 0])
                H[j, :] = res[0]

        if not sps.issparse(A):
            for j in iter(range(0, W.shape[0])):
                res = opt.nnls(H, A[j, :])
                W[j, :] = res[0]
        else:
            for j in iter(range(0, W.shape[0])):
                res = opt.nnls(H, A[j, :].toarray()[0,:])
                W[j, :] = res[0]
        return (W, H) 

def fit(self):
    print "Fitting a model with ", len(self.training_data), " points"
    labels = np.array([row[2] for row in self.training_data])
    data_points = np.array([self._get_features(row) for row in self.training_data])
    self.model = nnls(data_points, labels)
    # TODO: Add a debug logging mode ?
    # print "Residual norm ", self.model[1]
    # print "Model ", self.model[0]
    # Calculate training error
    training_errors = []
    for p in self.training_data:
      predicted = self.predict(p[0], p[1])
      training_errors.append(predicted / p[2])

    print "Average training error %f%%" % ((np.mean(training_errors) - 1.0)*100.0 )
    return self.model[0] 

def optimize(self):
    try:
      prev_cost = self.error()
      prev_w = self.w.copy()
      nz_idcs = self.w > 0
      res = nnls(self.A[:, nz_idcs], self.b, maxiter=100*self.A.shape[1])
      self.w[nz_idcs] = res[0]
      new_cost = self.error()
      if new_cost > prev_cost*(1.+util.TOL):
        raise NumericalPrecisionError('self.optimize() returned a solution with increasing error. Numeric limit possibly reached: preverr = ' + str(prev_cost) + ' err = ' + str(new_cost) + '.\n \
                                        If the two errors are very close, try running bc.util.tolerance(tol) with tol > current tol = ' + str(util.TOL) + ' before running')
    except NumericalPrecisionError as e:
      self.log.warning(e)
      self.w = prev_w
      self.reached_numeric_limit = True
      return 

def fit_nnls(self):
        from scipy.optimize import nnls
        return nnls(self.design_matrix, self.titer_dist)[0] 

def __init__(self, x, y, n_blocks=None, nn=False, separators=None):
        Jackknife.__init__(self, x, y, n_blocks, separators)
        if nn:  # non-negative least squares
            func = lambda x, y: np.atleast_2d(nnls(x, np.array(y).T[0])[0])
        else:
            func = lambda x, y: np.atleast_2d(
                np.linalg.lstsq(x, np.array(y).T[0])[0])

        self.est = func(x, y)
        self.delete_values = self.delete_values(x, y, func, self.separators)
        self.pseudovalues = self.delete_values_to_pseudovalues(
            self.delete_values, self.est)
        (self.jknife_est, self.jknife_var, self.jknife_se, self.jknife_cov) =\
            self.jknife(self.pseudovalues) 

def _support_least_square(X, uv, z, debug=False):
    """WIP, not fonctional!"""
    n_trials, n_channels, n_times = X.shape
    n_atoms, _, n_times_valid = z.shape
    n_times_atom = n_times - n_times_valid + 1

    # Compute DtD
    DtD = compute_DtD(uv, n_channels)
    t0 = n_times_atom - 1
    z_hat = np.zeros(z.shape)

    for idx in range(n_trials):
        Xi = X[idx]
        support_i = z[:, idx].nonzero()
        n_support = len(support_i[0])
        if n_support == 0:
            continue
        rhs = np.zeros((n_support, n_support))
        lhs = np.zeros(n_support)

        for i, (k_i, t_i) in enumerate(zip(*support_i)):
            for j, (k_j, t_j) in enumerate(zip(*support_i)):
                dt = t_i - t_j
                if abs(dt) < n_times_atom:
                    rhs[i, j] = DtD[k_i, k_j, t0 + dt]
            aux_i = np.dot(uv[k_i, :n_channels], Xi[:, t_i:t_i + n_times_atom])
            lhs[i] = np.dot(uv[k_i, n_channels:], aux_i)

        # Solve the non-negative least-square with nnls
        z_star, a = optimize.nnls(rhs, lhs)
        for i, (k_i, t_i) in enumerate(zip(*support_i)):
            z_hat[k_i, idx, t_i] = z_star[i]

    return z_hat 

def calculate(self, x_fit: np.ndarray) -> np.ndarray:
        if x_fit.ndim == 1:
            x_fit = x_fit.reshape(x_fit.shape[0], 1)
        M = self.X.collapse()

        N, p1 = M.shape
        q, p2 = x_fit.T.shape

        self.map = np.zeros((N, q), dtype=np.float32)
        MtM = np.dot(x_fit.T, x_fit)
        for n1 in range(N):
            self.map[n1] = opt.nnls(MtM, np.dot(x_fit.T, M[n1]))[0]
        self.map = self.map.reshape(self.X.shape[:-1] + (x_fit.shape[-1],))

        return self.map 

def __init__(self, X: 'hp.hparray'):
        self.ucls = ucls(X)
        self.nnls = nnls(X)
        self.fcls = fcls(X) 

def nn_omp(X, D, n_nonzero_coefs=None, tol=None):
    """ The Non Negative OMP algorithm of
        'On the Uniqueness of Nonnegative Sparse Solutions to Underdetermined Systems of Equations'"""

    n_samples = X.shape[1]
    n_atoms = D.shape[1]
    Z = np.zeros((n_atoms, n_samples))
    _norm = np.sum(D ** 2, axis=0)
    for i in range(n_samples):

        x = X[:, i]
        r = x
        z = np.zeros(n_atoms)
        Dx = np.array([]).astype(int)
        j = 0
        if n_nonzero_coefs is not None:
            tol = 1e-20

            def cont_criterion():
                not_reached_sparsity = j < n_nonzero_coefs
                return (not_reached_sparsity and norm(r) > tol)
        else:
            cont_criterion = lambda: norm(r) > tol

        while (cont_criterion()):
            a = np.dot(D.T, r)
            a[a < 0] = 0
            e = (norm(r) ** 2) - (a ** 2) / _norm
            k = np.argmin(e)
            Dx = np.append(Dx, k)

            z_est = nnls(D[:, Dx], x)[0]
            r = x - np.dot(D[:, Dx], z_est)
            j += 1

        if j != 0:
            z[Dx] = z_est
        Z[:, i] = z
    return Z 

def fit(self, X, y, sample_weight=None):
        """Fit non-negative linear model.

        Parameters
        ----------
        X : numpy array or sparse matrix of shape [n_samples, n_features]
            Training data
        y : numpy array of shape [n_samples,]
            Target values
        sample_weight : numpy array of shape [n_samples]
            Individual weights for each sample

        Returns
        -------
        self : returns an instance of self.

        """
        X, y = check_X_y(X, y, y_numeric=True, multi_output=False)

        if sample_weight is not None and np.atleast_1d(sample_weight).ndim > 1:
            raise ValueError("Sample weights must be 1D array or scalar")

        X, y, X_offset, y_offset, X_scale = self._preprocess_data(
            X, y, fit_intercept=self.fit_intercept, normalize=self.normalize,
            copy=self.copy_X, sample_weight=sample_weight)

        if sample_weight is not None:
            # Sample weight can be implemented via a simple rescaling.
            X, y = _rescale_data(X, y, sample_weight)

        self.coef_, result = nnls(X, y.squeeze())

        if np.all(self.coef_ == 0):
            raise ConvergenceWarning("All coefficients estimated to be zero in"
                                     " the non-negative least squares fit.")

        self._set_intercept(X_offset, y_offset, X_scale)
        self.opt_result_ = OptimizeResult(success=True, status=0, x=self.coef_,
                                          fun=result)
        return self 

def _train(self, method='nnl1reg',  lam_drop=1.0, lam_pot = 0.5, lam_avi = 3.0, **kwargs):
        '''
        determine the model parameters -- lam_drop, lam_pot, lam_avi are
        the regularization parameters.

        Parameters
        ----------
        method : str, optional
            Description
        lam_drop : float, optional
            Description
        lam_pot : float, optional
            Description
        lam_avi : float, optional
            Description
        **kwargs
            Description
        '''
        self.lam_pot = lam_pot
        self.lam_avi = lam_avi
        self.lam_drop = lam_drop
        if len(self.train_titers)==0:
            raise InsufficientDataException("Error: No titers in training set.")
        else:
            if method=='l1reg':  # l1 regularized fit, no constraint on sign of effect
                self.model_params = self.fit_l1reg()
            elif method=='nnls':  # non-negative least square, not regularized
                self.model_params = self.fit_nnls()
            elif method=='nnl2reg': # non-negative L2 norm regularized fit
                self.model_params = self.fit_nnl2reg()
            elif method=='nnl1reg':  # non-negative fit, branch terms L1 regularized, avidity terms L2 regularized
                self.model_params = self.fit_nnl1reg()

            print('rms deviation on training set=',np.sqrt(self.fit_func()))

        # extract and save the potencies and virus effects. The genetic parameters
        # are subclass specific and need to be process by the subclass
        self.serum_potency = {serum:self.model_params[self.genetic_params+ii]
                              for ii, serum in enumerate(self.sera)}
        self.virus_effect = {strain:self.model_params[self.genetic_params+len(self.sera)+ii]
                             for ii, strain in enumerate(self.test_strains)} 

def nnls(A, B):
    def func(b):
        return optimize.nnls(A, b)[0]
    return np.array(map(func, B)) 

def select(self):
    #do least squares on active set
    res = nnls(self.A[:, self.active_idcs].T, self.b, maxiter=100*self.A.shape[1])

    x_opt = self.A.dot(res[0])
    w_opt = np.zeros(self.w.shape[0])
    w_opt[self.active_idcs] = res[0]
    xw = self.A.dot(self.w)
    sdir = x_opt - xw
    sdir /= np.sqrt((sdir**2).sum())

    #do line search towards x_opt

    #find earliest gamma for which a variable joins the active set
    #anywhere gamma_denom = 0 or gamma < 0, the variable never enters the active set 
    gamma_nums = (sdir - self.x).dot(self.b - xw)
    gamma_denoms = (sdir - self.x).dot(x_opt - xw)
    good_idcs = np.logical_not(np.logical_or(gamma_denoms == 0, gamma_nums*gamma_denoms < 0))
    gammas = np.inf*np.ones(gamma_nums.shape[0])
    gammas[good_idcs] = gamma_nums[good_idcs]/gamma_denoms[good_idcs]
    f_least_angle = gammas.argmin()
    gamma_least_angle = gammas[f_least_angle]

    #find earliest gamma for which a variable leaves the active set
    f_leave_active = -1
    gamma_leave_active = np.inf
    gammas[:] = np.inf
    leave_idcs = w_opt < 0
    gammas[leave_idcs] = self.wts[leave_idcs]/(self.wts[leave_idcs] - w_opt[leave_idcs])
    f_leave_active = gammas.argmin()
    gamma_leave_active = gammas[f_leave_active] 

def _reweight(self, f):
    self.w[f] = 1.
    nz_idcs = self.w > 0
    res = nnls(self.A[:, nz_idcs], self.b, maxiter=100*self.A.shape[1])
    self.w[nz_idcs] = res[0]
    return 

def test_maxiter(self):
        # test that maxiter argument does stop iterations
        # NB: did not manage to find a test case where the default value
        # of maxiter is not sufficient, so use a too-small value
        rndm = np.random.RandomState(1234)
        a = rndm.uniform(size=(100, 100))
        b = rndm.uniform(size=100)
        with assert_raises(RuntimeError):
            nnls(a, b, maxiter=1) 

def test_nnls(self):
        a = arange(25.0).reshape(-1,5)
        x = arange(5.0)
        y = dot(a,x)
        x, res = nnls(a,y)
        assert_(res < 1e-7)
        assert_(norm(dot(a,x)-y) < 1e-7) 

def __init__(self, x, y, n_blocks=None, nn=False, separators=None):
        Jackknife.__init__(self, x, y, n_blocks, separators)
        if nn:  # non-negative least squares
            func = lambda x, y: np.atleast_2d(nnls(x, np.array(y).T[0])[0])
        else:
            func = lambda x, y: np.atleast_2d(
                np.linalg.lstsq(x, np.array(y).T[0])[0])

        self.est = func(x, y)
        self.delete_values = self.delete_values(x, y, func, self.separators)
        self.pseudovalues = self.delete_values_to_pseudovalues(
            self.delete_values, self.est)
        (self.jknife_est, self.jknife_var, self.jknife_se, self.jknife_cov) =\
            self.jknife(self.pseudovalues) 

def test_nnls(self):
        a = arange(25.0).reshape(-1,5)
        x = arange(5.0)
        y = dot(a,x)
        x, res = nnls(a,y)
        assert_(res < 1e-7)
        assert_(norm(dot(a,x)-y) < 1e-7) 

def optimize_weights(num_bones, Tr, C, verts0, K, nnls_soft_constr_weight=1.e+4):
    """ optimize for blending weights according to section 4.4 """
    # reshape Tr according to section 4.3 into 
    # a matrix of d * P of 4-vectors
    # TODO: not sure if this is correct already
    #       during iteration this step sometimes increases the residual :-(
    num_verts = verts0.shape[0]
    Tr1 = Tr.reshape((-1, num_bones, 4))
    new_W = np.zeros((num_bones, num_verts))
    for i in xrange(num_verts):
        # form complete system matrix Sigma * w = y
        Sigma = np.inner(Tr1, V.hom4(verts0[i]))
        #import ipdb; ipdb.set_trace()
        y = C[:,i]
        # choose K columns that individually best explain the residual
        error = ( (Sigma - y[:,np.newaxis]) ** 2 ).sum(axis=0)
        k_best = np.argsort(error)[:K]
        # solve nonlinear least squares problem
        # with additional convexity constraint sum(k_weighs) = 1
        Sigma = Sigma[:,k_best]
        #k_weights0 = W[k_best, i]
        k_weights, _ = nnls(
            np.vstack((Sigma, nnls_soft_constr_weight * np.ones(K))),
            np.append(y, nnls_soft_constr_weight))
        new_W[k_best, i] = k_weights
    return new_W 

def solution_bounds(endmember_occupancies):
    """
    Parameters
    ----------
    endmember_occupancies : 2d array of floats
        A 1D array for each endmember in the solid solution,
        containing the number of atoms of each element on each site.

    Returns
    -------
    solution_bounds : 2d array of floats
        An abbreviated version of endmember_occupancies, 
        where the columns represent the independent compositional 
        bounds on the solution
    """
    # Find bounds for the solution
    i_sorted =zip(*sorted([(i,
                            sum([1 for val in endmember_occupancies.T[i]
                                 if val>1.e-10]))
                           for i in range(len(endmember_occupancies.T))
                                          if np.any(endmember_occupancies.T[i] > 1.e-10)],
                          key=lambda x: x[1]))[0]

    solution_bounds = endmember_occupancies[:,i_sorted[0],np.newaxis]
    for i in i_sorted[1:]:
        if np.abs(nnls(solution_bounds, endmember_occupancies.T[i])[1]) > 1.e-10:
            solution_bounds = np.concatenate((solution_bounds,
                                              endmember_occupancies[:,i,np.newaxis]),
                                             axis=1)
    return solution_bounds 

def change_component_set(self, new_component_list):
        """
        Change the set of basis components without 
        changing the bulk composition. 

        Will raise an exception if the new component set is 
        invalid for the given composition.

        Parameters
        ----------
        new_component_list : list of strings
            New set of basis components.
        """
        composition = np.array([self.atomic_composition[element]
                                for element in self.element_list])
        component_matrix = np.zeros((len(new_component_list), len(self.element_list)))
        
        for i, component in enumerate(new_component_list):
            formula = dictionarize_formula(component)
            for element, n_atoms in formula.items():
                component_matrix[i][self.element_list.index(element)] = n_atoms

        sol = nnls(component_matrix.T, composition)
        if sol[1] < 1.e-12:
            component_amounts = sol[0]
        else:
            raise Exception('Failed to change component set. '
                            'Could not find a non-negative least squares solution. '
                            'Can the bulk composition be described with this set of components?')

        composition = OrderedCounter(dict(zip(new_component_list, component_amounts)))
        self.__init__(composition, 'molar') 

def _test_nnlsm():
    print ('\nTesting nnls routines ...\n')
    m = 100
    n = 10
    k = 200
    rep = 5

    for r in iter(range(0, rep)):
        A = np.random.rand(m, n)
        X_org = np.random.rand(n, k)
        X_org[np.random.rand(n, k) < 0.5] = 0
        B = A.dot(X_org)
        # B = np.random.rand(m,k)
        # A = np.random.rand(m,n/2)
        # A = np.concatenate((A,A),axis=1)
        # A = A + np.random.rand(m,n)*0.01
        # B = np.random.rand(m,k)

        import time
        start = time.time()
        C1, info = nnlsm_blockpivot(A, B)
        elapsed2 = time.time() - start
        rel_norm2 = nla.norm(C1 - X_org) / nla.norm(X_org)
        print ('nnlsm_blockpivot:    ', 'OK  ' if info[0] else 'Fail',\
            'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed2, rel_norm2))

        start = time.time()
        C2, info = nnlsm_activeset(A, B)
        num_backup = 0
        elapsed1 = time.time() - start
        rel_norm1 = nla.norm(C2 - X_org) / nla.norm(X_org)
        print ('nnlsm_activeset:     ', 'OK  ' if info[0] else 'Fail',\
            'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed1, rel_norm1))

        import scipy.optimize as opt
        start = time.time()
        C3 = np.zeros([n, k])
        for i in iter(range(0, k)):
            res = opt.nnls(A, B[:, i])
            C3[:, i] = res[0]
        elapsed3 = time.time() - start
        rel_norm3 = nla.norm(C3 - X_org) / nla.norm(X_org)
        print ('scipy.optimize.nnls: ', 'OK  ',\
            'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed3, rel_norm3))

        if num_backup > 0:
            break
        if rel_norm1 > 10e-5 or rel_norm2 > 10e-5 or rel_norm3 > 10e-5:
            break
        print ('') 

def fit_signatures(W, genome, metric="l2"):
    genome = np.array(genome)
    maxmutation = round(np.sum(genome))
     
    
    
    #initialize the guess
    x0 = np.random.rand(W.shape[1], 1)*maxmutation
    x0= x0/np.sum(x0)*maxmutation
    
    #the optimization step
    #set the bounds and constraints
    bnds = create_bounds([], genome, W.shape[1]) 
    cons1 ={'type': 'eq', 'fun': constraints1, 'args':[genome]}
    
    
    #sol = scipy.optimize.minimize(parameterized_objective2_custom, x0, args=(W, genome),  bounds=bnds, constraints =cons1, tol=1e-30)
    #solution = sol.x
    
    
    ### using NNLS algorithm 
    reg = nnls(W,genome)
    weights = reg[0]
    est_genome = np.dot(W, weights)
    normalised_weights = weights/sum(weights)
    solution = normalised_weights*sum(genome)
    
    
    
    #print(W1)
    #convert the newExposure vector into list type structure
    newExposure = list(solution)
    
    
    # get the maximum value of the new Exposure
    maxcoef = max(newExposure)
    idxmaxcoef = newExposure.index(maxcoef)
    
    newExposure = np.round(newExposure)
    
    # We may need to tweak the maximum value of the new exposure to keep the total number of mutation equal to the original mutations in a genome
    if np.sum(newExposure)!=maxmutation:
        newExposure[idxmaxcoef] = round(newExposure[idxmaxcoef])+maxmutation-sum(newExposure)
     
    
    
    if metric=="cosine":
        newSimilarity = cos_sim(genome, est_genome)
    else:
        newSimilarity = np.linalg.norm(genome-est_genome , ord=2)/np.linalg.norm(genome, ord=2)
    
    return (newExposure, newSimilarity)

# to do the calculation, we need: W, samples(one column), H for the specific sample, tolerence
# Fit signatures with multiprocessing 

def fit_signatures_pool(total_genome, W, index):
    genome = total_genome[:,index]
    genome = np.array(genome)
    maxmutation = round(np.sum(genome))
     
    
    
    #initialize the guess
    x0 = np.random.rand(W.shape[1], 1)*maxmutation
    x0= x0/np.sum(x0)*maxmutation
    
    #the optimization step
    #set the bounds and constraints
    bnds = create_bounds([], genome, W.shape[1]) 
    cons1 ={'type': 'eq', 'fun': constraints1, 'args':[genome]}
    
    
    #sol = scipy.optimize.minimize(parameterized_objective2_custom, x0, args=(W, genome),  bounds=bnds, constraints =cons1, tol=1e-30)
    #solution = sol.x
    
    
    ### using NNLS algorithm 
    reg = nnls(W,genome)
    weights = reg[0]
    est_genome = np.dot(W, weights)
    normalised_weights = weights/sum(weights)
    solution = normalised_weights*sum(genome)
    
    #print(W1)
    #convert the newExposure vector into list type structure
    newExposure = list(solution)
    
    
    # get the maximum value of the new Exposure
    maxcoef = max(newExposure)
    idxmaxcoef = newExposure.index(maxcoef)
    
    newExposure = np.round(newExposure)
    
    # We may need to tweak the maximum value of the new exposure to keep the total number of mutation equal to the original mutations in a genome
    if np.sum(newExposure)!=maxmutation:
        newExposure[idxmaxcoef] = round(newExposure[idxmaxcoef])+maxmutation-sum(newExposure)
     
    
    
    newSimilarity = cos_sim(genome, est_genome)
    
    return (newExposure, newSimilarity) 

def fit_signatures(W, genome, metric="l2"):
    genome = np.array(genome)
    maxmutation = round(np.sum(genome))
     
    
    
    #initialize the guess
    x0 = np.random.rand(W.shape[1], 1)*maxmutation
    x0= x0/np.sum(x0)*maxmutation
    
    #the optimization step
    #set the bounds and constraints
    bnds = create_bounds([], genome, W.shape[1]) 
    cons1 ={'type': 'eq', 'fun': constraints1, 'args':[genome]}
    
    
    #sol = scipy.optimize.minimize(parameterized_objective2_custom, x0, args=(W, genome),  bounds=bnds, constraints =cons1, tol=1e-30)
    #solution = sol.x
    
    
    ### using NNLS algorithm 
    reg = nnls(W,genome)
    weights = reg[0]
    est_genome = np.dot(W, weights)
    normalised_weights = weights/sum(weights)
    solution = normalised_weights*sum(genome)
    
    
    
    #print(W1)
    #convert the newExposure vector into list type structure
    newExposure = list(solution)
    
    
    # get the maximum value of the new Exposure
    maxcoef = max(newExposure)
    idxmaxcoef = newExposure.index(maxcoef)
    
    newExposure = np.round(newExposure)
    
    # We may need to tweak the maximum value of the new exposure to keep the total number of mutation equal to the original mutations in a genome
    if np.sum(newExposure)!=maxmutation:
        newExposure[idxmaxcoef] = round(newExposure[idxmaxcoef])+maxmutation-sum(newExposure)
     
    
    
    if metric=="cosine":
        newSimilarity = cos_sim(genome, est_genome)
    else:
        newSimilarity = np.linalg.norm(genome-est_genome , ord=2)/np.linalg.norm(genome, ord=2)
    
    return (newExposure, newSimilarity)

# to do the calculation, we need: W, samples(one column), H for the specific sample, tolerence
# Fit signatures with multiprocessing 

def fit_signatures_pool(total_genome, W, index):
    genome = total_genome[:,index]
    genome = np.array(genome)
    maxmutation = round(np.sum(genome))
     
    
    
    #initialize the guess
    x0 = np.random.rand(W.shape[1], 1)*maxmutation
    x0= x0/np.sum(x0)*maxmutation
    
    #the optimization step
    #set the bounds and constraints
    bnds = create_bounds([], genome, W.shape[1]) 
    cons1 ={'type': 'eq', 'fun': constraints1, 'args':[genome]}
    
    
    #sol = scipy.optimize.minimize(parameterized_objective2_custom, x0, args=(W, genome),  bounds=bnds, constraints =cons1, tol=1e-30)
    #solution = sol.x
    
    
    ### using NNLS algorithm 
    reg = nnls(W,genome)
    weights = reg[0]
    est_genome = np.dot(W, weights)
    normalised_weights = weights/sum(weights)
    solution = normalised_weights*sum(genome)
    
    #print(W1)
    #convert the newExposure vector into list type structure
    newExposure = list(solution)
    
    
    # get the maximum value of the new Exposure
    maxcoef = max(newExposure)
    idxmaxcoef = newExposure.index(maxcoef)
    
    newExposure = np.round(newExposure)
    
    # We may need to tweak the maximum value of the new exposure to keep the total number of mutation equal to the original mutations in a genome
    if np.sum(newExposure)!=maxmutation:
        newExposure[idxmaxcoef] = round(newExposure[idxmaxcoef])+maxmutation-sum(newExposure)
     
    
    
    newSimilarity = cos_sim(genome, est_genome)
    
    return (newExposure, newSimilarity) 

def _fitting_nnls(data, ref_spectra, *, maxiter=100):
    r"""
    Fitting of multiple spectra using NNLS method.

    Parameters
    ----------

    data : ndarray(float), 2D
        array holding multiple observed spectra, shape (K, N), where K is the number of energy points,
        and N is the number of spectra

    absorption_refs : ndarray(float), 2D
        array of references, shape (K, Q), where Q is the number of references.

    maxiter : int
        maximum number of iterations. Optimization may stop prematurely if convergence criteria are met.

    Returns
    -------

    map_data_fitted : ndarray(float), 2D
        fitting results, shape (Q, N), where Q is the number of references and N is the number of spectra.

    map_rfactor : ndarray(float), 2D
        map that represents R-factor for the fitting, shape (M,N).

    map_residual : ndarray(float), 2D
        residual returned by NNLS algorithm
    """
    assert data.ndim == 2, "Data array 'data' must have 2 dimensions"
    assert ref_spectra.ndim == 2, "Data array 'ref_spectra' must have 2 dimensions"

    n_pts = data.shape[0]
    n_pixels = data.shape[1]
    n_pts_2 = ref_spectra.shape[0]
    n_refs = ref_spectra.shape[1]

    assert n_pts == n_pts_2, f"The number of spectrum points in data ({n_pts}) "\
                             f"and references ({n_pts_2}) do not match."

    assert maxiter > 0, f"The parameter 'maxiter' is zero or negative ({maxiter})"

    map_data_fitted = np.zeros(shape=[n_refs, n_pixels])
    map_rfactor = np.zeros(shape=[n_pixels])
    map_residual = np.zeros(shape=[n_pixels])
    for n in range(n_pixels):
        map_sel = data[:, n]
        result, residual = nnls(ref_spectra, map_sel, maxiter=maxiter)

        rfactor = rfactor_compute(map_sel, result, ref_spectra)

        map_data_fitted[:, n] = result
        map_rfactor[n] = rfactor
        map_residual[n] = residual

    return map_data_fitted, map_rfactor, map_residual