main.cpp
// Copyright (c) 2012 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com, Dmitri Maslov, Michele Mosca
//
// This file is part of SQCT.
//
// SQCT is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// SQCT is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with SQCT. If not, see <http://www.gnu.org/licenses/>.
//
#include "sk.h"
#include "exactdecomposer.h"
#include "gcommdecomposer.h"
#include "theoremverification.h"
#include "toptimalitytest.h"
#include "hoptimalitytest.h"
#include "netgenerator.h"
#include "output.h"
#include <fstream>
#include <sstream>
#include <chrono>
#include <boost/program_options.hpp>
#include <boost/timer/timer.hpp>
#include <cmath>
// command line parcing code based on http://www.boost.org/doc/libs/1_49_0/doc/html/program_options/tutorial.html#id2499896
// and BOOST_ROOT/libs/program_options/example/first.cpp
namespace po = boost::program_options;
namespace btm = boost::timer;
using namespace std;
//////////////////////////////////////////////////////////////////////////////
/// \brief Option of SK execution
struct SKOptions
{
SKOptions() :
math_fr( false ), write_dot_qc(false), angle(false) {};
string in_filename; ///< input filename
string out_filename; ///< output statistics filename
bool math_fr; ///< should output be mathematica friendly
int iterations_default; ///< number of SK iterations to perform by default
bool write_dot_qc;///< should we output dot qc or not
string out_dir;///< path where we should output dotqc files
bool angle;///< how we interpret first element \b val of each input line ( Pi/\b val or just \b val)
int max_sde;///< maximal value of sde to use during basic approximation step
};
//////////////////////////////////////////////////////////////////////////////
/// \brief Result of application of the SK algorithm
struct ApplicationResult
{
circuit c; ///< Circuit
int hc; ///< Hadamard counts
int tc; ///< T and inverse of T gates counts
int pc; ///< Phase and inverse of Phase gate counts
int plc; ///< Pauli matrices counts
int total; ///< Total cost using cost function defined by vector gateLibrary::cost
double dst; ///< Trace distance to approximation
double tappr; ///< Time required to do an approximation
double tdecomp; ///< Time required to do a decomposition
int nr; ///< Number of iterations performed
int denom_reduction; ///< Difference between \f$ sde(|\cdot|^2) \f$
/// before and after conversion to canonical form
int denom ; ///< \f$ sde(|\cdot|^2) \f$ of resulting unitary
matrix2x2<mpz_class> exact_uni; /// exact unitary corresponding to the circuit c
/// \brief Summarize gate statistics for Clifford + T library
void updateForCliffordT()
{
typedef gateLibrary gl;
total = c.cost();
auto counts = gl::toCliffordT( c.count() );
hc = counts[ gl::H ];
tc = counts[ gl::T ];
pc = counts[ gl::P ];
plc = counts[ gl::Z ] + counts[ gl::X ] + counts[ gl::Y ];
}
/// \brief Outputs results in CSV format
void toCSVStream( std::ostream& out ) const
{
out.precision(5);
out.setf( ios_base::scientific );
out << nr << ","
<< tc << ","
<< hc << ","
<< pc << ","
<< plc << ","
<< dst << ","
<< tappr << ","
<< tdecomp << "," << denom_reduction << "," << denom ;
}
/// \brief Outputs results in Mathematica friendly format
void toMathematicaStream( std::ostream& out ) const
{
out.precision(10);
out.setf( ios_base::fixed ); // use fixed as mathematica does not understand scientific c++ format
out << "AppResult["
<< nr << ","
<< tc << ","
<< hc << ","
<< pc << ","
<< plc << ",";
out.precision(60);
out << dst << "," ;
out.precision(10);
out << tappr << ","
<< tdecomp << "," << denom_reduction << "," << denom << "]";
}
/// \brief Outputs dot qc file with circuit and adds information about genration process into comments
void generateDotQc( std::string filename, std::string circuit_name, std::string symbolic_form )
{
ofstream ofs(filename.c_str());
if( !ofs )
{
cout << "Unable to open file:" << filename << " for output" <<endl;
return;
}
ofs << "# Computed by SQCT, based on arXiv:1206.5236" << endl;
ofs << "# Published at: http://code.google.com/p/sqct/" << endl;
ofs << "# Symbolic form of unitary to approximate:" << symbolic_form << endl;
ofs << "# Distance between unitary and approximation:" << dst << endl;
ofs << "# (Formula (1) in http://arxiv.org/abs/quant-ph/0411206)" << endl;
ofs << "# Total number of gates:" << tc + hc + pc + plc << endl;
ofs << "# T and T^{Dagger} gates:" << tc << endl;
ofs << "# Hadamard gates:" << hc << endl;
ofs << "# Phase and Phase^{Dagger} gates:" << pc << endl;
ofs << "# Pauli gates:" << plc << endl;
ofs << "# Total cost: " << total << endl;
ofs << "# Number of iterations:" << nr << endl;
ofs << "# Approximation time(seconds):" << tappr << endl;
ofs << "# Decomposition time(seconds):" << tdecomp << endl;
ofs << "# Reduction:" << denom_reduction << endl;
ofs << "# sde(abs(z)^2):" << denom << endl;
ofs << "# Exact unitary:" << endl << exact_uni << endl;
ofs << ".v 1" << endl;
ofs << ".i 1" << endl;
ofs << ".o 1" << endl;
ofs << "BEGIN " << circuit_name <<"(1)" << endl;
c.toStream(ofs);
ofs << "END " << circuit_name << endl;
ofs << "BEGIN" << endl;
ofs << circuit_name << " 1" << endl;
ofs << "END" << endl;
ofs.close();
}
};
///////////////////////////////////////////////////////////////////
/// \brief Process input file with inputs for SK in batch mode, outputs statistics etc
class SKApplication
{
public:
SKApplication( SKOptions& skopt ) :
opt( skopt ), skd(skopt.max_sde + 1)
{}
/// \brief Finds circuit and all supplementary information for unitary
void calculate( const matrix2x2hpr& matrix, int recursion_level, ApplicationResult& ar )
{
boost::timer::cpu_times ct;
boost::timer::cpu_times ct2;
sk::Me res;
{
boost::timer::auto_cpu_timer t;
skd.decompose(matrix,res,recursion_level);
ct = t.elapsed();
}
int gde_before = res.min_gde_abs2();
res.reduce();
int gde_after = res.min_gde_abs2();
ar.denom_reduction = gde_before - gde_after;
ar.denom = res.max_sde_abs2();
ar.exact_uni = res;
{
boost::timer::auto_cpu_timer t;
ed.decompose(res,ar.c);
ct2 = t.elapsed();
}
sk::Ma conv(res);
ar.dst = trace_dist(matrix,conv);
ar.tappr = (double) ct.wall * 1e-9;
ar.tdecomp = (double) ct2.wall * 1e-9;
ar.nr = recursion_level;
ar.updateForCliffordT();
}
/// \brief Process input file defined by opt taking into account specified options
void process()
{
static const double TwoPi = 2. * hprHelpers::toMachine( hprHelpers::pi() );
ifstream ifs( opt.in_filename.c_str() );
ofstream ofs( opt.out_filename.c_str() );
Rotation r;
int iter = 0;
ApplicationResult ar;
int line_number = 0;
while( true )
{
string line;
getline(ifs,line);
iter = opt.iterations_default;
istringstream ss(line,istringstream::in);
if( ! ifs ) break;
double val;
ss >> val >> r.nx >> r.ny >> r.nz >> iter;
if( ! opt.angle )
{
r.num = 1.0;
r.den = val;
}
else
{
r.num = val;
r.den = TwoPi;
}
calculate( r.matrix() , iter, ar );
stringstream fname;
fname << opt.out_dir << "/";
if( r.isSpecial() != 'N' )
fname << r.name() << "-" << iter;
else
fname << opt.out_filename << "." << line_number;
fname << ".qc" ;
if( opt.write_dot_qc )
ar.generateDotQc( fname.str() ,
r.name(),
r.symbolic() );
if( opt.math_fr )
{
ofs << "{" << r.Mathematica() << ",";
ar.toMathematicaStream(ofs);
ofs << "}" << endl;
}
else
{
ofs << r.CSV() << ",";
ar.toCSVStream(ofs);
ofs << endl;
}
line_number++;
}
}
private:
const SKOptions& opt;
sk skd;
exactDecomposer ed;
};
///////////////////////////////////////////////////////////////////
bool print_help( const string& topic )
{
map<string,string> hi;
hi["in"] = "Each line of input file must have the followoing form:\n"
"val,nx,ny,nz,[iter]\n"
"val -- interpreted depending on the value of option --angle\n"
" if --angle is specified then val interpreted as rotation\n"
" angle phi. Otherwise rotation angle phi = 2 Pi / val.\n"
"nx -- projection of rotation axis on, x\n"
"ny -- projection of rotation axis on, y\n"
"nz -- projection of rotation axis on, z\n"
"iter -- optional parameter defining number of iterations used for approximation.\n"
"The matrix that will be approximated is I*cos(phi/2) - \n"
"i sin(phi/2) (nx * X + ny * Y + nz * Z). If the number of iterations is not specified,\n"
"the default values is used. It can be set using --iterations option.\n"
"\n"
"Default output format is CSV and each line of it contains the following:\n"
"num,den,nx,ny,nz,N,Tc,Hc,Pc,Plc,dst,tappr,tdecomp,denom_red,denom \n"
"[num,den,nx,ny,nz] determines input rotation with phi = num * 2Pi / den\n"
"N -- number of iterations of SK algorithm used"
"Tc-- number of T gates in the resulting circuit\n"
"Hc-- number of Hadamard gates in the circuit\n"
"Pc-- number of Phase gates in the circuit\n"
"Plc-- number of Pauli gates in the circuit\n"
"dst-- distance to approximation, see \n"
"(Formula (1) in http://arxiv.org/abs/quant-ph/0411206)\n"
"tappr-- time required for approximation (seconds) \n"
"tdecomp-- time required for exact decomposition (seconds)\n"
"tdecomp-- time required for exact decomposition (seconds)\n"
"denom_red -- change of the power of sqrt(2) in the denominator\n"
"after conversion to canonical form.\n"
"denom -- power of sqrt(2) in the denominator of the exact unitary in\n"
"the canonical form. In the case of Mathematica output numbers are \n"
"decorated with Rot and AppResult list headers.\n";
if( hi.count(topic) )
{
cout << hi[topic] << endl;
return true;
}
return false;
}
/////////////////////////////////////////////////////////////////////
struct ResynthOptions
{
ResynthOptions() :
reverse(false)
{}
vector<string> circuit_files;
bool reverse;
bool math_fr;
};
////////////////////////////////////////////////////////////////////
/// \brief Resyntheisizes circuits using exact decomposition algoritm
struct ResynthApplication
{
ResynthApplication( const ResynthOptions& opt ) :
ro(opt)
{}
void process()
{
for( const auto& fn : ro.circuit_files )
process(fn);
}
void process( const string& filename )
{
const gateLibrary& gl = gateLibrary::instance();
ifstream ifs(filename.c_str());
ofstream ofs( (filename + ".out").c_str() );
// read circuit first
circuit in, out;
in.fromStream(ifs,ro.reverse);
ed.decompose(in,out);
auto in_cost = in.count();
int in_total = in.cost();
in_cost = gl.toCliffordT(in_cost);
auto out_cost = out.count();
int out_total = out.cost();
out_cost = gl.toCliffordT(out_cost);
string cb = "# ";
string ce = "";
if( ro.math_fr )
cb = "(* ",ce = " *)";
ofs << cb << "Computed by SQCT, based on arXiv:1206.5236" << ce << endl;
ofs << cb << "Software published at: http://code.google.com/p/sqct/"<< ce << endl;
ofs << cb << "Gate counts for input circuit [T+Td,H,P+Pd,X,Z,Y] : ["
<< in_cost[gl.T] << "," << in_cost[gl.H] << ","
<< in_cost[gl.P] << "," << in_cost[gl.X] << ","
<< in_cost[gl.Z] << "," << in_cost[gl.Y] << "]" << ce << endl;
ofs << cb << "Total cost of input:" << in_total << ce << endl;
ofs << cb << "Gate counts for output circuit [T+Td,H,P+Pd,X,Z,Y] : ["
<< out_cost[gl.T] << "," << out_cost[gl.H] << ","
<< out_cost[gl.P] << "," << out_cost[gl.X] << ","
<< out_cost[gl.Z] << "," << out_cost[gl.Y] << "]" << ce << endl;
ofs << cb << "Total cost of output:" << out_total << ce << endl;
if( ro.math_fr )
out.toMathStream(ofs);
else
out.toStreamSym(ofs);
ofs.close();
ifs.close();
}
exactDecomposer ed;
const ResynthOptions& ro;
};
////////////////////////////////////////////////////////////////////
/// \brief Options for epsilon net generation
struct enetOptions
{
/// \brief If one element specified -- upper bound for sde of epsilon net to be
/// generated. If two elements specified -- interval of sde to be generated.
vector<int> epsilon_net_layers;
};
////////////////////////////////////////////////////////////////////
struct enetApplication
{
enetApplication( const enetOptions& options ) :
m_options(options)
{}
/// \brief Checks if all layers generated on initial state are available
void check_initial()
{
initial_ok = true;
if( m_layers[0] == 0 )
initial_ok = false;
for( int i = 2; (i < initial_end) && initial_ok ; ++i )
if( m_layers[i] == 0 ) initial_ok = false;
}
/// \brief Collects information about parts of epsilon net available
void check_files()
{
m_layers.resize(100,0);
// collect information about available layers
for( int i = 0; i < 100; ++i )
{
string name = netGenerator::fileName(i);
ifstream ifs(name);
if( !ifs )
m_layers[i] = 0;
else
m_layers[i] = 1;
}
}
/// \brief Generates all layers with \f$ sde(|\cdot|^2) \f$ in the interval [st,end)
void generate( int st, int end )
{
if( ! initial_ok && st < initial_end )
m_ng.generateInitial();
for( int i = max(initial_end,st) ; i < end; ++i )
{
if( m_layers[i] == 0 )
{
epsilonnet base_net;
base_net.loadFromFile( netGenerator::fileName(i-1).c_str() );
unique_ptr<epsilonnet> res( netGenerator::generate( base_net ) );
res->saveToFile( netGenerator::fileName(i).c_str() );
}
}
}
/// \brief Do the work
void process()
{
check_files();
check_initial();
int sz = m_options.epsilon_net_layers.size();
int b = m_options.epsilon_net_layers[0];
switch( sz )
{
case 1:
generate(0,b + 1);
break;
case 2:
generate(b,m_options.epsilon_net_layers[1] + 1);
break;
default:
std::cerr << "Wrong number of parameters. See help." << endl;
}
}
void print_error_message()
{
cout << "Epsilon net files a not available, use --epsilon-net option" << endl
<< "to generate them." << endl;
}
bool check()
{
int st = 0;
int end = m_options.epsilon_net_layers[0] + 1;
check_files();
check_initial();
if( ! initial_ok && st < initial_end )
{
print_error_message();
return false;
}
for( int i = max(initial_end,st) ; i < end; ++i )
{
if( m_layers[i] == 0 )
{
print_error_message();
return false;
}
}
return true;
}
const enetOptions& m_options; ///< Application options
std::vector<int> m_layers; ///< Ones for available layers, zeros for not availible layers
static const int initial_end; ///< Maximal layer generated on initial state
bool initial_ok; ///< If all initial layers are availible
netGenerator m_ng; ///< Epsilon net generator
};
const int enetApplication::initial_end = 21;
////////////////////////////////////////////////////////////////////
struct RandomRotationsOpts
{
string filename;
int samples;
};
////////////////////////////////////////////////////////////////////
struct RotationGeneratorApp
{
RotationGeneratorApp( RandomRotationsOpts& opts) :
m_opt(opts)
{}
/// \brief Generates file with rotations by random angle around random axis
void process()
{
ofstream ofs( m_opt.filename.c_str() );
unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
std::default_random_engine generator (seed);
static const double PI = 3.1415926535897932384626433832795028841971693993751;
normal_distribution<double> nd(0.0,1.0);
uniform_real_distribution<double> urd( 0.0, 2*PI );
for( int i = 0; i < m_opt.samples; ++i )
{
double x = nd(generator);
double y = nd(generator);
double z = nd(generator);
double norm = sqrt(x*x + y*y + z*z);
x /= norm;
y /= norm;
z /= norm;
double phi = urd(generator);
ofs << phi << " " << x << " " << y << " " << z << endl;
}
}
const RandomRotationsOpts& m_opt;
};
////////////////////////////////////////////////////////////////////
struct exactDecompositionOpts
{
std::string filename;
};
////////////////////////////////////////////////////////////////////
struct exactDecompositionApplication
{
exactDecompositionApplication( const exactDecompositionOpts& opts ) :
m_opts(opts)
{
}
void run()
{
ifstream ifs;
ifs.exceptions(std::fstream::badbit | std::fstream::failbit);
ifs.open(m_opts.filename);
matrix2x2<mpz_class> exact_uni;
circuit c;
ifs >> exact_uni;
cout << exact_uni << endl;
ed.decompose(exact_uni,c);
c.toStream(cout);
}
exactDecompositionOpts m_opts;
exactDecomposer ed;
};
////////////////////////////////////////////////////////////////////
void print_about_message()
{
cout << "Copyright (c) 2012 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com" << endl << endl;
cout << "SQCT is free software: you can redistribute it and/or modify" << endl;
cout << "it under the terms of the GNU Lesser General Public License as published by" << endl;
cout << "the Free Software Foundation, either version 3 of the License, or" << endl;
cout << "(at your option) any later version." << endl;
cout << "" << endl;
cout << "SQCT is distributed in the hope that it will be useful," << endl;
cout << "but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl;
cout << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl;
cout << "GNU Lesser General Public License for more details." << endl;
cout << "" << endl;
cout << "You should have received a copy of the GNU Lesser General Public License" << endl;
cout << "along with SQCT. If not, see <http://www.gnu.org/licenses/>." << endl;
cout << "" << endl;
cout << "The program code based on results of http://arxiv.org/abs/1206.5236." << endl;
cout << "It also implements the version of Solovay Kitaev algorithm described" << endl;
cout << "in http://arxiv.org/abs/quant-ph/0505030. " << endl;
cout << "" << endl;
cout << "The list of used libraries, corresponding licenses and authors:" << endl<< endl;
cout << "Boost | Boost Software License" << endl<< endl;
cout << "The GNU Multiple Precision Arithmetic Library | GNU Lesser General Public License v3" << endl;
cout << "(by Free Software Foundation, Inc.)" << endl<< endl;
cout << "The GNU MPFR Library the library | GNU Lesser General Public License v3" << endl;
cout << "(by Free Software Foundation, Inc.," << endl;
cout << " Contributed by the AriC and Caramel projects, INRIA. )" << endl << endl;
cout << "mpfr::real | GNU Lesser General Public License v3" << endl;
cout << "(by Christian Schneider <software(at)chschneider(dot)eu> )" << endl;
cout << "" << endl;
cout << "Source code of this program can be obtained from: " << endl;
}
////////////////////////////////////////////////////////////////////
int main(int ac, char* av[])
{
// help
string help_topic = "";
// sk application parameters
SKOptions sko;
// epsilon net generation
enetOptions eopt;
// resynthesise
ResynthOptions ro;
// randomized rotations
RandomRotationsOpts rro;
// exact decomposition
exactDecompositionOpts edo;
try {
po::options_description desc("Allowed options");
desc.add_options()
("help,H", po::value< string >(&help_topic)->implicit_value(""),
"Produce help message, see help <option name> for more details "
"about specific option.")
("in,I", po::value< string >(&(sko.in_filename)),
"File name with unitaries for approximation.")
("iterations,N", po::value< int >(&(sko.iterations_default))->default_value(3),
"Default number of iteration used by the Solovay Kitaev algorithm, "
"when not specified in input file.")
("angle,A", po::value< bool >(&(sko.angle))->zero_tokens(),
"When true program interprets first entry in each line of input file as "
"as rotation angle, and as denominator of Pi.")
("out,O", po::value< string >(&(sko.out_filename))->default_value("out.csv"),
"File name with summary about "
"unitary approximation results.")
("math-fr,F", po::value< bool >(&(sko.math_fr))->zero_tokens(),
"Should output be Mathematica friendly. "
"Default output format is CSV.")
("dotqc,C", po::value< bool >(&(sko.write_dot_qc))->zero_tokens(),
"Should we output dotqc file for each generated circuit.")
("out-dir,D", po::value< string >(&(sko.out_dir))->default_value("out"),
"Directory to output circuits.")
("max-sde,M", po::value< int >(&(sko.max_sde))->default_value(30),
"Maximal value of sde to use during intial approximation step.")
("epsilon-net,E", po::value< vector<int> >(&(eopt.epsilon_net_layers))->multitoken(),
"Generates epsilon net layers. If one values specified generates all layers "
"with sde(|.|^2) less or equal than this value. If two values specified then sde(|.|^2) "
"of result lies in interval [first,second] ")
("theory-topt",
"Verifies conjecture about T optimality ")
("theory-hopt",
"Performs brute force check necessary for proof of result in Appendix 2 "
"in http://arxiv.org/abs/1206.5236")
("theory-correctness",
"Performs brute force check described by Algoritm 2 "
"in http://arxiv.org/abs/1206.5236 ( proves Lemma 3 ).")
("resynth,S", po::value< vector<string> >(&(ro.circuit_files) )->multitoken(),
"Resynthesize circuit. Output may differ from original by global phase.")
("reverse", po::value< bool >(& (ro.reverse) )->zero_tokens(),
"True if circuit file should be interpreted in matrix multiplication order.")
("random-rot", po::value< int >(&(rro.samples)),
"Generates file with specified number of random rotations.")
("exact", po::value< string >(&(edo.filename)),
"Performs an exact decomposition of a unitary over the ring.")
("about", "Information about the program.")
;
po::variables_map vm;
po::store(po::parse_command_line(ac, av, desc), vm);
po::notify(vm);
if( vm.count("resynth") ) {
ro.math_fr = sko.math_fr;
ResynthApplication ra(ro);
ra.process();
return 0;
}
if(vm.count("epsilon-net"))
{
enetApplication eapp(eopt);
eapp.process();
return 0;
}
if( vm.count("theory-correctness") ) {
cout << is_theorem_true() << endl;
return 0;
}
if( vm.count("theory-topt") ) {
toptimalitytest tt;
return 0;
}
if( vm.count("theory-hopt") ) {
hoptimalitytest ht;
return 0;
}
if ( vm.count("random-rot") ) {
rro.filename = sko.out_filename;
RotationGeneratorApp rga(rro);
rga.process();
return 0;
}
if( vm.count("about") ) {
print_about_message();
return 0;
}
if ( vm.count("in") ) {
eopt.epsilon_net_layers.clear();
eopt.epsilon_net_layers.push_back(sko.max_sde);
enetApplication ea(eopt);
if( ea.check() )
{
SKApplication app(sko);
app.process();
}
return 0;
}
if( vm.count("exact") ) {
exactDecompositionApplication eda(edo);
eda.run();
return 0;
}
if ( true ) {
if( !print_help(help_topic) )
cout << desc << endl;
return 0;
}
}
catch(exception& e) {
cerr << "Error: " << e.what() << endl;
return 1;
}
catch(...) {
cerr << "Exception of unknown type!" << endl;
return 1;
}
}