Add C++ of balanced binary tree butterfly factorization + C++ tests.

misc/test/CMakeLists.txt

@@ -199,7 +199,7 @@ if(MATIO_LIB_FILE AND MATIO_INC_DIR AND BUILD_READ_MAT_FILE AND NOT NOCPPTESTS)
 	# faust_multiplication : time comparison between Faust-vector product and Dense matrix-vector product
 
 
-	list(APPEND tests hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication  faust_matdense_conjugate GivensFGFT GivensFGFTSparse GivensFGFTParallel GivensFGFTParallelSparse test_MatDiag faust_matsparse_mul faust_matsparse_index_op GivensFGFTComplex GivensFGFTComplexSparse GivensFGFTParallelComplex faust_toeplitz faust_circ faust_hankel faust_sparse_prox_sp palm4msa_2020 hierarchical2020 hierarchical2020Hadamard hierarchicalFactorizationHadamard hierarchicalFactorizationButterfly test_MatBSR)
+	list(APPEND tests hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication  faust_matdense_conjugate GivensFGFT GivensFGFTSparse GivensFGFTParallel GivensFGFTParallelSparse test_MatDiag faust_matsparse_mul faust_matsparse_index_op GivensFGFTComplex GivensFGFTComplexSparse GivensFGFTParallelComplex faust_toeplitz faust_circ faust_hankel faust_sparse_prox_sp palm4msa_2020 hierarchical2020 hierarchical2020Hadamard hierarchicalFactorizationHadamard hierarchicalFactorizationButterfly hierarchicalFactorizationButterflyBalanced test_MatBSR)
 
 	if(FAUST_TORCH)
 		list(APPEND tests faust_torch)

misc/test/src/C++/hierarchicalFactorizationButterfly.cpp.in

@@ -13,8 +13,8 @@ using namespace std;
 int main(int argc, char** argv)
 {
 	int log2dim;
-	if(argc > 2)
-		log2dim = std::atoi(argv[2]);
+	if(argc >= 2)
+		log2dim = std::atoi(argv[1]);
 	else
 		log2dim = 10;
 	auto wht = TransformHelper<FPP,Cpu>::hadamardFaust(log2dim, false);
@@ -100,6 +100,14 @@ int main(int argc, char** argv)
 	std::cout << "Factorization dir: " << dir << std::endl;
 	auto th = butterfly_hierarchical(A, dir);
 	th->display();
+//	auto support = support_DFT<FPP>((int) log2dim);
+//	for(int i=0;i<support.size();i++)
+//	{
+//		auto fac = th->get_fact(i);
+//		fac.setNZtoOne();
+//		fac -= Faust::MatDense<FPP,Cpu>(*support[i]);
+//		std::cout << " fac " << i << " err:" << fac.norm()/ th->get_fact(i).norm() << std::endl;
+//	}
 	auto M = th->get_product();
 	M -= A;
 	auto err = M.norm()/A.norm();

misc/test/src/C++/hierarchicalFactorizationButterflyBalanced.cpp.in

+#include "faust_init_from_matio_params.h"
+#include "faust_init_from_matio_core.h"
+#include "faust_MatDense.h"
+#include "faust_TransformHelper.h"
+#include "faust_butterfly.h"
+
+typedef @TEST_FPP@ FPP;
+
+using namespace Faust;
+using namespace std;
+
+
+int main(int argc, char** argv)
+{
+	int log2dim;
+	if(argc >= 2)
+		log2dim = std::atoi(argv[1]);
+	else
+		log2dim = 10;
+	auto wht = TransformHelper<FPP,Cpu>::hadamardFaust(log2dim, false);
+	Faust::MatDense<FPP,Cpu> A;
+	/********** NOTE: a lot of garbage is left here, it was useful during the first writing of butterfly module (it could be again) */
+	//	Faust::MatDense<FPP,Cpu> support1, support2, A;
+	//	Faust::MatDense<FPP,Cpu> submatrixA, bestx, besty, X, Y, refX, refY;
+	//	init_faust_mat_from_matio(support1, "supports.mat", "support1");
+	//	init_faust_mat_from_matio(support2, "supports.mat", "support2");
+	//	init_faust_mat_from_matio(submatrixA, "svd.mat", "submatrixA");
+	//	init_faust_mat_from_matio(bestx, "svd.mat", "bestx");
+	//	init_faust_mat_from_matio(besty, "svd.mat", "besty");
+	//	init_faust_mat_from_matio(A, "A.mat", "A");
+	//	init_faust_mat_from_matio(refX, "XY.mat", "X");
+	//	init_faust_mat_from_matio(refY, "XY.mat", "Y");
+
+	//	support1.Display();
+	//	support2.Display();
+	cout << "matrix to factorize:" << std::endl;
+	A = wht->get_product();
+	A.Display();
+
+
+	//	vector<vector<faust_unsigned_int>*> cec, noncec;
+	////	vector<vector<faust_unsigned_int>> cec, noncec;
+	//	retrieveCEC(support1, support2, cec, noncec);
+	//	std::cout << "cec:" << std::endl;
+	//	print_classes(cec);
+	//	std::cout << "noncec:" << std::endl;
+	//	print_classes(noncec);
+	//
+	//	Faust::MatDense<FPP,Cpu> _bestx(bestx.getNbRow(), bestx.getNbCol());
+	//	Faust::MatDense<FPP,Cpu> _besty(besty.getNbRow(), besty.getNbCol());
+	//	int r = bestx.getNbCol();
+	//	submatrixA.best_low_rank(r, _bestx, _besty);
+	//	std::cout << bestx.norm() << " " << _bestx.norm() << std::endl;
+	//	std::cout << besty.norm() << " " << _besty.norm() << std::endl;
+	//	bestx.Display();
+	//	_bestx.Display();
+	//	besty.Display();
+	//	_besty.Display();
+	//
+	//
+	//	r = 1;
+	//	auto RP = support1.col_nonzero_inds(r);
+	//	for(int i=0;i<RP.size();i++)
+	//		std::cout << RP[i] << " ";
+	//	std::cout << std::endl;
+	//	auto CP = support2.row_nonzero_inds(r);
+	//	for(int i=0;i<CP.size();i++)
+	//		std::cout << CP[i] << " ";
+	//	std::cout << std::endl;
+	//
+	//	Faust::MatDense<FPP, Cpu> submat;
+	//	A.submatrix(RP, CP, submat);
+	//	submat.Display();
+	//
+	////    colx, rowx = np.meshgrid(ce, RP)
+	////    print("ce=", ce, "RP=", RP, "colx=", colx, "rowx=", rowx)
+	////    coly, rowy = np.meshgrid(CP, ce)
+	////    print("ce=", ce, "CP=", CP, "coly=", coly, "rowy=", rowy)
+	////    bestx, besty = best_low_rank(submatrixA, len(ce))
+	////    if not exists("svd.mat"):
+	////        savemat("svd.mat", {'submatrixA': submatrixA, 'bestx': bestx, 'besty':
+	////                            besty})
+	////    X[rowx, colx] = bestx
+	////    Y[rowy, coly] = besty
+	//
+	//	solveDTO(A, support1, support2, X, Y);
+	//	std::cout << X.norm() << " " << refX.norm() << std::endl;
+	//	std::cout << Y.norm() << " " << refY.norm() << std::endl;
+	//
+	//	Faust::MatDense<FPP, Cpu> ones(2, 2);
+	//	Faust::MatDense<FPP, Cpu> id(5, 5);
+	//	ones.setOnes();
+	//	id.setEyes();
+	//
+	//	Faust::MatDense<FPP, Cpu> out;
+	//	Faust::kron(ones, id, out);
+	//	std::cout << out.norm() << std::endl;
+
+	auto th = butterfly_hierarchical(A, BALANCED);
+	th->display();
+//	auto support = support_DFT<FPP>((int) log2dim);
+//	for(int i=0;i<support.size();i++)
+//	{
+//		auto fac = th->get_fact(i);
+//		fac.setNZtoOne();
+//		fac -= Faust::MatDense<FPP,Cpu>(*support[i]);
+//		std::cout << " fac " << i << " err:" << fac.norm()/ th->get_fact(i).norm() << std::endl;
+//	}
+	auto M = th->get_product();
+	M -= A;
+	auto err = M.norm()/A.norm();
+	std::cout << "err=" << err << std::endl;
+	if(err > 1e-6)
+		exit(1);
+}

src/algorithm/factorization/faust_butterfly.h

@@ -6,6 +6,7 @@ namespace Faust
 	{
 		LEFT,
 		RIGHT,
+		BALANCED
 	};
 
 	template<typename FPP>
@@ -26,6 +27,9 @@ namespace Faust
 	template<typename FPP>
 		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical(const Faust::MatDense<FPP, Cpu>& A, const std::vector<Faust::MatSparse<FPP, Cpu>*> &supports, const ButterflyFactDir& dir=RIGHT);
 
+	template<typename FPP>
+		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical_balanced(const Faust::MatDense<FPP, Cpu>& A, const std::vector<Faust::MatSparse<FPP, Cpu>*> &supports);
+
 	template<typename FPP>
 		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical(const Faust::MatDense<FPP, Cpu>& A, const ButterflyFactDir &dir=RIGHT);
 

src/algorithm/factorization/faust_butterfly.hpp

@@ -289,6 +289,89 @@ namespace Faust
 			return th;
 		}
 
+	template<typename FPP>
+		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical_balanced(const Faust::MatDense<FPP, Cpu>& A, const std::vector<Faust::MatSparse<FPP, Cpu>*> &supports)
+		{
+			auto th = new TransformHelper<FPP, Cpu>();
+			double l2_nfacts = std::log2((double)supports.size());
+			size_t d = ceil(l2_nfacts); // depth of factorization tree
+			/** initialize the support tree structure (the size of each vector -- one per level of tree) **/
+			std::vector<std::vector<MatDense<FPP, Cpu>>> tree_supports(d);
+			// set supports for each tree level
+			size_t i=0;
+			while(i<d-1)
+			{
+				tree_supports[i] = std::vector<MatDense<FPP, Cpu>>(1 << (i+1));
+//				std::cout << "tree_supports["<<i<<"].size():" << tree_supports[i].size() << std::endl;
+				i++;
+			}
+			// the number of leafs of tree is not necessarily a power of two (remaining factorization(s))
+			size_t n = supports.size();
+			if(d > 1)
+				tree_supports[d-1] = std::vector<MatDense<FPP, Cpu>>(2*(n - tree_supports[d-2].size()));
+			else // factorization in two factors only
+				tree_supports[d-1] = std::vector<MatDense<FPP, Cpu>>(2);
+//			std::cout << "tree_supports[" << d-1 << "].size():" << tree_supports[d-1].size() << std::endl;
+			/** initialize the supports in the tree starting from the leafs -- that are the supports arguments **/
+			for(int i=d-1;i>=0;i--)
+			{
+				#pragma omp parallel for schedule(dynamic)
+				for(int j = 0;j < tree_supports[i].size(); j++)
+				{
+					if(i == d - 1)
+					{
+						// tree_supports[i][j] is a leaf on the last level
+//						std::cout << "tree_supports[" << i << "," << j << "] received support[" << j << "]" << std::endl;
+						tree_supports[i][j] = *supports[j];
+					}
+					else if (2*(j+1) > tree_supports[i+1].size())
+					{
+						// tree_supports[i][j] is a leaf
+//						std::cout << "tree_supports[" << i << "," << j << "] received support[" << j+tree_supports[i+1].size()/2 << "]" << std::endl;
+						tree_supports[i][j] = *supports[j+tree_supports[i+1].size()/2];
+					}
+					else
+					{
+						// tree_supports[i][j] is not a leaf so it's calculated using child nodes
+						// tree_supports[i][j] = tree_supports[i+1][2*j] * tree_supports[i+1][2*j+]
+//						std::cout << "tree_supports[" << i << "," << j << "] received tree_supports[" << i+1 << "," << 2*j << "]*tree_supports[" << i+1 << "," << 2*j+1 << "]" << std::endl;
+						gemm(tree_supports[i+1][2*j], tree_supports[i+1][2*j+1], tree_supports[i][j], FPP(1.0), FPP(0), 'N', 'N');
+					}
+				}
+			}
+			// the supports are now set for all factorization tree nodes
+			// factorize according to the tree of supports, level by level
+			std::vector<MatDense<FPP, Cpu>> next_lvl_facs;
+			std::vector<MatDense<FPP, Cpu>> input_facs;
+			// first factorization
+			next_lvl_facs.resize(2);
+			solveDTO(A, tree_supports[0][0], tree_supports[0][1], next_lvl_facs[0], next_lvl_facs[1]);
+			input_facs = next_lvl_facs;
+			for(size_t i=0;i<d-1;i++)
+			{
+				next_lvl_facs.resize(tree_supports[i+1].size());
+				auto j_bound = std::min(tree_supports[i].size(), tree_supports[i+1].size()/2);
+				#pragma omp parallel for schedule(static)
+				for(size_t j=0;j < j_bound;j++)
+					solveDTO(input_facs[j], tree_supports[i+1][2*j], tree_supports[i+1][2*j+1], next_lvl_facs[j*2], next_lvl_facs[j*2+1]);
+				if(i < d - 2)
+				{
+					input_facs = next_lvl_facs;
+				}
+			}
+			// build the resulting Faust
+			for(auto f: next_lvl_facs)
+			{
+				th->push_back(new MatSparse<FPP, Cpu>(f), false, false);
+			}
+			if(next_lvl_facs.size() < supports.size())
+				for(int i=next_lvl_facs.size()/2;i<input_facs.size();i++)
+				{
+					th->push_back(new MatSparse<FPP, Cpu>(input_facs[i]), false, false);
+				}
+			return th;
+		}
+
 	template<typename FPP>
 		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical(const Faust::MatDense<FPP, Cpu>& A, const ButterflyFactDir &dir/*=RIGHT*/)
 		{
@@ -301,7 +384,11 @@ namespace Faust
 			//	std::cout << "support norms" << std::endl;
 			//	for(auto s: support)
 			//		std::cout << s->norm() << std::endl;
-			auto th = butterfly_hierarchical(A, support, dir);
+			TransformHelper<FPP, Cpu>* th = nullptr;
+			if(dir == BALANCED)
+				th = butterfly_hierarchical_balanced(A, support);
+			else
+				th = butterfly_hierarchical(A, support, dir);
 			return th;
 		}