docs: private <A> double[] MStepLevenbergMarquardt(double a, double b, long[] t, A array, NumberArrayAdapter<?, A> adapter) {
      final int size = adapter.size(array);
      final int prior1 = BitsUtil.cardinality(t);
      final int prior0 = size - prior1;
  
      final int maxiter = 10;
      final double minstep = 1e-8, sigma = 1e-12;
      // target value for "set" objects
      final double loTarget = (prior1 + 1.0) / (prior1 + 2.0);
      // target value for "unset" objects
      final double hiTarget = 1.0 / (prior0 + 2.0);
      // t[i] := t.get(i) ? hiTarget : loTarget.
  
      // Reset, or continue with previous values?
      // a = 0.0;
      // b = FastMath.log((prior0 + 1.0) / (prior1 + 1.0));
      double fval = 0.0;
      for(int i = 0; i < size; i++) {
        final double val = adapter.getDouble(array, i);
        final double fApB = val * a + b;
        final double ti = BitsUtil.get(t, i) ? hiTarget : loTarget;
        if(fApB >= 0) {
          fval += ti * fApB + FastMath.log(1 + FastMath.exp(-fApB));
        }
        else {
          fval += (ti - 1) * fApB + FastMath.log(1 + FastMath.exp(fApB));
        }
      }
      for(int it = 0; it < maxiter; it++) {
        // Update Gradient and Hessian (use H’ = H + sigma I)
        double h11 = sigma, h22 = sigma, h21 = 0.0;
        double g1 = 0.0, g2 = 0.0;
        for(int i = 0; i < size; i++) {
          final double val = adapter.getDouble(array, i);
          final double fApB = val * a + b;
          final double p, q;
          if(fApB >= 0) {
            p = FastMath.exp(-fApB) / (1.0 + FastMath.exp(-fApB));
            q = 1.0 / (1.0 + FastMath.exp(-fApB));
          }
          else {
            p = 1.0 / (1.0 + FastMath.exp(fApB));
            q = FastMath.exp(fApB) / (1.0 + FastMath.exp(fApB));
          }
          final double d2 = p * q;
          h11 += val * val * d2;
          h22 += d2;
          h21 += val * d2;
          final double d1 = (BitsUtil.get(t, i) ? hiTarget : loTarget) - p;
          g1 += val * d1;
          g2 += d1;
        }
        // Stop condition
        if(Math.abs(g1) < 1e-5 && Math.abs(g2) < 1e-5) {
          break;
        }
        // Compute modified Newton directions
        final double det = h11 * h22 - h21 * h21;
        final double dA = -(h22 * g1 - h21 * g2) / det;
        final double dB = -(-h21 * g1 + h11 * g2) / det;
        final double gd = g1 * dA + g2 * dB;
        double stepsize = 1.0;
        while(stepsize >= minstep) { // Line search
          final double newA = a + stepsize * dA;
          final double newB = b + stepsize * dB;
          double newf = 0.0;
          for(int i = 0; i < size; i++) {
            final double val = adapter.getDouble(array, i);
            final double fApB = val * newA + newB;
            final double ti = BitsUtil.get(t, i) ? hiTarget : loTarget;
            if(fApB >= 0) {
              newf += ti * fApB + FastMath.log(1 + FastMath.exp(-fApB));
            }
            else {
              newf += (ti - 1) * fApB + FastMath.log(1 + FastMath.exp(fApB));
            }
          }
          if(newf < fval + 0.0001 * stepsize * gd) {
            a = newA;
            b = newB;
            fval = newf;
            break; // Sufficient decrease satisfied
          }
          else {
            stepsize /= 2.0;
          }
          if(stepsize < minstep) {
            LOG.debug("Minstep hit.");
            break;
          }
        }
      }
      return new double[] { a, b };
    }