update

abess-team · May 17, 2021 · 7a685a2 · 7a685a2
1 parent cdca89b
commit 7a685a2
Show file tree

Hide file tree

Showing 5 changed files with 653 additions and 131 deletions.
diff --git a/.coverage b/.coverage
diff --git a/python/pytest/test_all.py b/python/pytest/test_all.py
@@ -29,8 +29,7 @@ def test_gaussian(self):
                         K_max=10, epsilon=10, powell_path=2, s_min=1, s_max=p, lambda_min=0.01, lambda_max=100, is_cv=True, K=5,
                         exchange_num=2, tau=0.1 * np.log(n*p) / n,
                         primary_model_fit_max_iter=10, primary_model_fit_epsilon=1e-6, early_stop=False, approximate_Newton=True, ic_coef=1., thread=5, covariance_update=False)
-        group = np.linspace(1, p, p)
-        model.fit(data.x, data.y, group=group)
+        model.fit(data.x, data.y)
 
         nonzero_true = np.nonzero(data.beta)[0]
         nonzero_fit = np.nonzero(model.beta)[0]
@@ -220,6 +219,13 @@ def test_mulnomial(self):
         group = np.linspace(1, p, p)
         model.fit(data.x, data.y, group=group)
 
+        model2 = abessMultinomial(path_type="seq", sequence=sequence, ic_type='ebic', is_screening=True, screening_size=100,
+                                  K_max=10, epsilon=10, powell_path=2, s_min=1, s_max=p, lambda_min=0.01, lambda_max=100, is_cv=False, K=5,
+                                  exchange_num=2, tau=0.1 * np.log(n*p) / n,
+                                  primary_model_fit_max_iter=30, primary_model_fit_epsilon=1e-6, early_stop=False, approximate_Newton=False, ic_coef=1., thread=5)
+        group = np.linspace(1, p, p)
+        model2.fit(data.x, data.y, group=group)
+
         nonzero_true = np.nonzero(data.beta)[0]
         nonzero_fit = np.nonzero(model.beta)[0]
         print(nonzero_true)

diff --git a/python/src/model_fit.cpp b/python/src/model_fit.cpp
@@ -350,104 +350,6 @@ void multigaussian_fit(Eigen::MatrixXd &x, Eigen::MatrixXd &y, Eigen::VectorXd &
     beta = cg.solveWithGuess(x.adjoint() * y, beta);
 }
 
-Eigen::VectorXd logit_fit(Eigen::MatrixXd &x, Eigen::VectorXd &y, int n, int p, Eigen::VectorXd weights)
-{
-    if (n <= p)
-    {
-        Eigen::MatrixXd X = Eigen::MatrixXd::Ones(n, n);
-        Eigen::VectorXd beta0 = Eigen::VectorXd::Zero(n);
-        Eigen::VectorXd beta1 = Eigen::VectorXd::Zero(n);
-        Eigen::VectorXd one = Eigen::VectorXd::Ones(n);
-        X.rightCols(n - 1) = x.leftCols(n - 1);
-        Eigen::MatrixXd X_new = Eigen::MatrixXd::Zero(n, n);
-        Eigen::VectorXd Pi = pi(X, y, beta0);
-        Eigen::VectorXd log_Pi = Pi.array().log();
-        Eigen::VectorXd log_1_Pi = (one - Pi).array().log();
-        double loglik0 = (y.cwiseProduct(log_Pi) + (one - y).cwiseProduct(log_1_Pi)).dot(weights);
-        Eigen::VectorXd W = Pi.cwiseProduct(one - Pi);
-        Eigen::VectorXd Z = X * beta0 + (y - Pi).cwiseQuotient(W);
-        W = W.cwiseProduct(weights);
-        for (int i = 0; i < n; i++)
-        {
-            X_new.row(i) = X.row(i) * W(i);
-        }
-        beta1 = (X_new.transpose() * X).ldlt().solve(X_new.transpose() * Z);
-
-        double loglik1;
-
-        int j;
-        for (j = 0; j < 30; j++)
-        {
-            Pi = pi(X, y, beta1);
-            log_Pi = Pi.array().log();
-            log_1_Pi = (one - Pi).array().log();
-            loglik1 = (y.cwiseProduct(log_Pi) + (one - y).cwiseProduct(log_1_Pi)).dot(weights);
-            if (abs(loglik0 - loglik1) / (0.1 + abs(loglik1)) < 1e-6)
-            {
-                break;
-            }
-            beta0 = beta1;
-            loglik0 = loglik1;
-            W = Pi.cwiseProduct(one - Pi);
-            Z = X * beta0 + (y - Pi).cwiseQuotient(W);
-            W = W.cwiseProduct(weights);
-            for (int i = 0; i < n; i++)
-            {
-                X_new.row(i) = X.row(i) * W(i);
-            }
-            beta1 = (X_new.transpose() * X).ldlt().solve(X_new.transpose() * Z);
-        }
-        return beta0;
-    }
-
-    else
-    {
-        Eigen::MatrixXd X = Eigen::MatrixXd::Ones(n, p + 1);
-        X.rightCols(p) = x;
-        Eigen::MatrixXd X_new = Eigen::MatrixXd::Zero(n, p + 1);
-        Eigen::VectorXd beta0 = Eigen::VectorXd::Zero(p + 1);
-        Eigen::VectorXd beta1 = Eigen::VectorXd::Zero(p + 1);
-        Eigen::VectorXd one = Eigen::VectorXd::Ones(n);
-        Eigen::VectorXd Pi = pi(X, y, beta0);
-        Eigen::VectorXd log_Pi = Pi.array().log();
-        Eigen::VectorXd log_1_Pi = (one - Pi).array().log();
-        double loglik0 = (y.cwiseProduct(log_Pi) + (one - y).cwiseProduct(log_1_Pi)).dot(weights);
-        Eigen::VectorXd W = Pi.cwiseProduct(one - Pi);
-        Eigen::VectorXd Z = X * beta0 + (y - Pi).cwiseQuotient(W);
-        W = W.cwiseProduct(weights);
-        for (int i = 0; i < n; i++)
-        {
-            X_new.row(i) = X.row(i) * W(i);
-        }
-        beta1 = (X_new.transpose() * X).ldlt().solve(X_new.transpose() * Z);
-        double loglik1;
-
-        int j;
-        for (j = 0; j < 30; j++)
-        {
-            Pi = pi(X, y, beta1);
-            log_Pi = Pi.array().log();
-            log_1_Pi = (one - Pi).array().log();
-            loglik1 = (y.cwiseProduct(log_Pi) + (one - y).cwiseProduct(log_1_Pi)).dot(weights);
-            if (abs(loglik0 - loglik1) / (0.1 + abs(loglik1)) < 1e-6)
-            {
-                break;
-            }
-            beta0 = beta1;
-            loglik0 = loglik1;
-            W = Pi.cwiseProduct(one - Pi);
-            Z = X * beta0 + (y - Pi).cwiseQuotient(W);
-            W = W.cwiseProduct(weights);
-            for (int i = 0; i < n; i++)
-            {
-                X_new.row(i) = X.row(i) * W(i);
-            }
-            beta1 = (X_new.transpose() * X).ldlt().solve(X_new.transpose() * Z);
-        }
-        return beta0;
-    }
-}
-
 void lm_fit(Eigen::MatrixXd &x, Eigen::VectorXd &y, Eigen::VectorXd &weights, Eigen::VectorXd &beta, double &coef0, double loss0, bool approximate_Newton, int primary_model_fit_max_iter, double primary_model_fit_epsilon, double tau, double lambda)
 {
     // beta = (X.adjoint() * X + this->lambda_level * Eigen::MatrixXd::Identity(X.cols(), X.cols())).colPivHouseholderQr().solve(X.adjoint() * y);