Skip to content

Commit

Permalink
stp computation revise
Browse files Browse the repository at this point in the history
  • Loading branch information
Zhufei Chu committed Sep 14, 2024
1 parent 73aac51 commit 24a6ac0
Show file tree
Hide file tree
Showing 6 changed files with 117 additions and 158 deletions.
2 changes: 1 addition & 1 deletion docs/conf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

extensions = [
'sphinx.ext.imgmath'
'sphinx.ext.mathjax'
]

templates_path = ['_templates']
Expand Down
40 changes: 36 additions & 4 deletions docs/stp_compute.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,8 @@ According to the definition,

.. math::
A \ltimes B = (A \bigotimes I_{4/4}) \cdot (B \bigotimes I_{4/2}) =
\begin{align}
A \ltimes B = (A \bigotimes I_{4/4}) \cdot (B \bigotimes I_{4/2}) &=
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 1 & 1
Expand All @@ -48,11 +49,13 @@ According to the definition,
0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0
\end{bmatrix} =
\end{bmatrix} \notag \\
&=
\begin{bmatrix}
1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\
0 & 1 & 0 & 1 & 1 & 1 & 0 & 1
\end{bmatrix}.
\end{bmatrix}. \notag
\end{align}
Copy Definition
-----------------
Expand All @@ -63,4 +66,33 @@ compute the STP :math:`A \bigotimes B`, we can examine :math:`B` column by
column from leftside to rightside, if the column in :math:`B` is
:math:`\begin{bmatrix} 1 \\ 0 \end{bmatrix}`, we copy :math:`A_{left}` as a
partial result; otherwise, we copy :math:`A_{right}`. One can verify the
results are exactly the same as the ones computed by native definition.
results are exactly the same as the ones computed by native definition.

Functions
----------------
In header file ``stp/stp_eigen.hpp``, we provide function::

matrix semi_tensor_product( const matrix& A, const matrix& B
const bool verbose = false,
const stp_method method = stp_method::copy_method )

to compute the STP of matrices :math:`A` and :math:`B`, where toggle ``verbose`` is off and toggle ``stp_method``
is used by the copy definition by default.

Example

.. code-block:: c++

matrix A;
matrix B;

//default
auto result = stp::semi_tensor_product( A, B );

//print verbose information
auto result = stp::semi_tensor_product( A, B, true );

//use native definition for STP computation
auto result = stp::semi_tensor_product( A, B, true, stp::stp_method::native_method );

One can find more examples or test cases in ``examples/stp_eigen.cpp`` and ``test/stp_eigen.cpp``.
30 changes: 30 additions & 0 deletions examples/stp_eigen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,4 +27,34 @@ int main()
<< stp::generate_swap_matrix( 2, 4 ) << std::endl;
std::cout << "W[8,2] is \n"
<< stp::generate_swap_matrix( 8, 2 ) << std::endl;

//n % p = 0 type matrix multiplication test
for(int i = 0; i < 10; i++)
{
int m = 200;
int n = 4 << i;
int p = 4;
int q = 100;
matrix A = stp::matrix_random_generation( m, n );
matrix B = stp::matrix_random_generation( p, q );
stp::semi_tensor_product( A, B, true, stp::stp_method::copy_method );
stp::semi_tensor_product( A, B, true, stp::stp_method::native_method );
}

std::cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n";

//p % n = 0 type matrix multiplication test
for(int i = 0; i < 10; i++)
{
int m = 200;
int n = 4;
int p = 4 << i;
int q = 100;
matrix A = stp::matrix_random_generation( m, n );
matrix B = stp::matrix_random_generation( p, q );
stp::semi_tensor_product( A, B, true, stp::stp_method::copy_method );
stp::semi_tensor_product( A, B, true, stp::stp_method::native_method );
}

return 0;
}
47 changes: 0 additions & 47 deletions examples/stp_test.cpp
View file

This file was deleted.

120 changes: 32 additions & 88 deletions include/stp/stp_eigen.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,19 @@ using matrix_chain = std::vector<matrix>; // Defined matrix chain
namespace stp
{

inline matrix matrix_random_generation( int row, int col )
{
matrix result( row, col );
for( int i = 0; i < row; i++ )
{
for( int j = 0; j < col; j++ )
{
result( i, j ) = rand() % 2;
}
}
return result;
}

inline matrix generate_swap_matrix( const int& m, const int& n )
{
matrix swap_matrixXi = matrix::Zero( m * n, m * n );
Expand Down Expand Up @@ -124,18 +137,16 @@ namespace stp
}
}


enum class stp_method : uint8_t
{
self_adaptation = 0,
definition_1 = 1,
definition_2 = 2,
copy_method = 1,
native_method = 2,
};

class semi_tensor_product_impl
{
public:
semi_tensor_product_impl( const matrix& A, const matrix& B, bool verbose, const stp_method method)
semi_tensor_product_impl( const matrix& A, const matrix& B, bool verbose, const stp_method method )
: A( A ), B( B ), verbose( verbose ), method( method )
{
m = A.rows();
Expand All @@ -147,32 +158,14 @@ namespace stp
matrix run()
{
matrix result;
uint64_t complexity_1 = cost1();
uint64_t complexity_2 = cost2();

switch (method)
if( method == stp_method::native_method )
{
case stp_method::definition_1:
result = call_with_stopwatch( time, [&]() { return semi_tensor_product1(); } );
break;

case stp_method::definition_2:
result = call_with_stopwatch( time, [&]() { return semi_tensor_product2(); } );
break;

case stp_method::self_adaptation:
if(complexity_1 < complexity_2)
{
result = call_with_stopwatch( time, [&]() { return semi_tensor_product1(); } );
}
else
{
result = call_with_stopwatch( time, [&]() { return semi_tensor_product2(); } );
}
break;

default:
break;
result = call_with_stopwatch( time, [&]() { return semi_tensor_product_native(); } );
}
else
{
result = call_with_stopwatch( time, [&]() { return semi_tensor_product_copy(); } );
}

if(verbose)
Expand All @@ -189,45 +182,8 @@ namespace stp
return ( m * n ) / std::gcd( m, n );
}

uint64_t cost1()
{
uint64_t t = get_lcm( n, p );
uint64_t kp = 2 * ( ( m * t / n ) * t + t * ( t * q / p ) );
uint64_t total = kp + ( ( 2 + 1 ) * t ) * ( ( m * t / n) * ( t * q / p) );
return total;
}

uint64_t cost2()
{
uint64_t total;
if( n % p == 0 )
{
uint64_t times = n / p;
uint64_t row = m;
uint64_t col = times * q;

total = ( 2 + 1 ) * ( m * times ) * p * q;
}
else if( p % n == 0 )
{
uint64_t times = p / n;
uint64_t row = m * times;
uint64_t col = q;

total = ( 2 + 1 ) * ( times * q ) * m * n;
}
else
{
std::cout << "matrix type error!" << std::endl;
assert(false);
}
return total;
}

matrix semi_tensor_product1()
matrix semi_tensor_product_native()
{
cost = cost1();

unsigned t = get_lcm( n, p );

matrix Ia = matrix::Identity( t / n, t / n );
Expand All @@ -239,7 +195,8 @@ namespace stp
return KPa * KPb;
}

matrix semi_tensor_product2()
//copy method as default
matrix semi_tensor_product_copy()
{
int row, col;
matrix result_matrix;
Expand Down Expand Up @@ -279,39 +236,27 @@ namespace stp
assert(false);
}

cost = cost2();
return result_matrix;
}

void report(const matrix& result)
{
std::cout << "----------------------------------------------------------\n";
std::cout << "--------------------STP Computation----------------------\n";
std::cout << "Dimension A: " << m << " x " << n << "\n";
std::cout << "Dimension B: " << p << " x " << q << "\n";

if( method == stp_method::definition_1 )
if( method == stp_method::native_method )
{
std::cout << "Use the definition 1!\n";
}
else if( method == stp_method::definition_2 )
{
std::cout << "Use the definition 2!\n";
std::cout << "Use the native method to compute STP.\n";
}
else
{
if( cost1() < cost2() )
{
std::cout << "Use the definition 1!\n";
}
else
{
std::cout << "Use the definition 2!\n";
}
std::cout << "Use the copy method to compute STP.\n";
}

std::cout << "Dimensions: " << result.rows() << " x " << result.cols() << "\n";
std::cout << "Result Dimensions: " << result.rows() << " x " << result.cols() << "\n";
std::cout << "Total time: " << to_millisecond( time ) << "ms\n";
std::cout << "----------------------------------------------------------\n";
std::cout << "---------------------------------------------------------\n";
}

private:
Expand All @@ -321,15 +266,14 @@ namespace stp
const stp_method method;
uint32_t m{0u}, n{0u}, p{0u}, q{0u};
stopwatch<>::duration time{0};
uint64_t cost{0u};
};

/*! \brief Compute the semi-tensor product of matrix A and matrix B.
Example:
matrix A, B;
matrix result = semi_tensor_product(A, B);
*/
inline matrix semi_tensor_product( const matrix& A, const matrix& B, const bool verbose = false, const stp_method method = stp_method::self_adaptation)
inline matrix semi_tensor_product( const matrix& A, const matrix& B, const bool verbose = false, const stp_method method = stp_method::copy_method )
{
semi_tensor_product_impl stp_impl( A, B, verbose, method);
return stp_impl.run();
Expand Down
36 changes: 18 additions & 18 deletions test/stp_eigen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,23 +65,23 @@ TEST_CASE( "power reducing matrix", "[eigen]" )
CHECK( stp1 == stp2 );
}

/*
TEST_CASE( "expression to tt", "[eigen]" )
TEST_CASE( "STP calculation by two methods", "[eigen]" )
{
std::string expr1 = "(a & b) | (a & ~c) | (~b & ~c)";
//n % p = 0 type
int m = 200;
int n1 = 4 << 5;
int p1 = 4;
int q = 100;
matrix A1 = stp::matrix_random_generation( m, n1 );
matrix B1 = stp::matrix_random_generation( p1, q );
CHECK( stp::semi_tensor_product( A1, B1, false, stp::stp_method::copy_method ) ==
stp::semi_tensor_product( A1, B1, false, stp::stp_method::native_method ) );

std::vector<std::string> inputs_order1 = { "a", "b", "c" };
MatrixXi mat1 = stp::from_exp_to_nmx( expr1, inputs_order1 );
CHECK( stp::to_binary( mat1 ) == "10001011" );
CHECK( stp::to_hex( mat1 ) == "8B" );
std::vector<std::string> inputs_order2 = { "a", "c", "b" };
MatrixXi mat2 = stp::from_exp_to_nmx( expr1, inputs_order2 );
CHECK( stp::to_binary( mat2 ) == "10001101" );
CHECK( stp::to_hex( mat2 ) == "8D" );
std::vector<std::string> inputs_order3 = { "c", "b", "a" };
MatrixXi mat3 = stp::from_exp_to_nmx( expr1, inputs_order3 );
CHECK( stp::to_binary( mat3 ) == "11010001" );
CHECK( stp::to_hex( mat3 ) == "D1" );
}*/
//p % n = 0 type
int n2 = 4;
int p2 = 4 << 5;
matrix A2 = stp::matrix_random_generation( m, n2 );
matrix B2 = stp::matrix_random_generation( p2, q );
CHECK( stp::semi_tensor_product( A2, B2, false, stp::stp_method::copy_method ) ==
stp::semi_tensor_product( A2, B2, false, stp::stp_method::native_method ) );
}

0 comments on commit 24a6ac0

Please sign in to comment.