sort zeta precompute optimizations

with reference implementation for comparison

src/algorithms/path_sgd.cpp

@@ -121,12 +121,18 @@ namespace odgi {
                                                                     iter_with_max_learning_rate,
                                                                     eps);
 
-                // cache zipf zetas for our full path space (heavy, but one-off)
+                // cache zipf zetas for our full path space
                 if (progress) {
                     std::cerr << "[odgi::path_linear_sgd] calculating zetas for " << (space <= space_max ? space : space_max + (space - space_max) / space_quantization_step + 1) << " zipf distributions" << std::endl;
                 }
-
-                std::vector<double> zetas((space <= space_max ? space : space_max + (space - space_max) / space_quantization_step + 1)+1);
+                std::cout << "space: " << space << std::endl;
+                std::cout << "space_max: " << space_max << std::endl;
+                std::cout << "space_quantization_step: " << space_quantization_step << std::endl;
+
+                auto start_zeta = std::chrono::high_resolution_clock::now();
+                /*
+                // reference zeta computation
+                std::vector<double> zetas_ref((space <= space_max ? space : space_max + (space - space_max) / space_quantization_step + 1)+1);
                 uint64_t last_quantized_i = 0;
 #pragma omp parallel for schedule(static,1)
                 for (uint64_t i = 1; i < space+1; ++i) {
@@ -139,11 +145,40 @@ namespace odgi {
 
                     if (quantized_i != last_quantized_i){
                         dirtyzipf::dirty_zipfian_int_distribution<uint64_t>::param_type z_p(1, compressed_space, theta);
-                        zetas[quantized_i] = z_p.zeta();
+                        zetas_ref[quantized_i] = z_p.zeta();
 
                         last_quantized_i = quantized_i;
                     }
                 }
+                */
+                // fast zeta computation
+                std::vector<double> zetas((space <= space_max ? space : space_max + (space - space_max) / space_quantization_step + 1)+1);
+                double zeta_tmp = 0.0;
+                for (uint64_t i = 1; i < space + 1; i++) {
+                    zeta_tmp += dirtyzipf::fast_precise_pow(1.0 / i, theta);
+                    if (i <= space_max) {
+                        zetas[i] = zeta_tmp;
+                    }
+                    if (i >= space_max && (i - space_max) % space_quantization_step == 0) {
+                        zetas[space_max + 1 + (i - space_max) / space_quantization_step] = zeta_tmp;
+                    }
+                }
+                auto end_zeta = std::chrono::high_resolution_clock::now();
+                uint32_t duration_zeta_ms = std::chrono::duration_cast<std::chrono::milliseconds>(end_zeta - start_zeta).count();
+                std::cout << "zeta precomputation took " << duration_zeta_ms << "ms" << std::endl;
+
+
+                /*
+                for (int i = 1; i < zetas.size(); i++) {
+                    if (zetas[i] != zetas_ref[i]) {
+                        std::cout << "WARNING[" << i << "]: " << zetas[i] << " vs " << zetas_ref[i] << std::endl;
+                    }
+                }
+
+                std::cout << zetas[zetas.size() - 1] << std::endl;
+                std::vector<double> X_tmp(X.size());
+                return X_tmp;
+                */
 
                 // how many term updates we make
                 std::atomic<uint64_t> term_updates;