.. _program_listing_file_blackbird_cpp_BlackbirdProgram.h:

Program Listing for File BlackbirdProgram.h
===========================================

|exhale_lsh| :ref:`Return to documentation for file <file_blackbird_cpp_BlackbirdProgram.h>` (``blackbird_cpp/BlackbirdProgram.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   // Copyright 2019 Xanadu Quantum Technologies Inc.
   
   // Licensed under the Apache License, Version 2.0 (the "License");
   // you may not use this file except in compliance with the License.
   // You may obtain a copy of the License at
   
   //     http://www.apache.org/licenses/LICENSE-2.0
   
   // Unless required by applicable law or agreed to in writing, software
   // distributed under the License is distributed on an "AS IS" BASIS,
   // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   // See the License for the specific language governing permissions and
   // limitations under the License.
   #include <unordered_map>
   
   typedef std::vector<std::vector<std::complex<double>>> complexmat;
   typedef std::vector<std::complex<double>> complexvec;
   
   typedef std::vector<std::vector<double>> floatmat;
   typedef std::vector<double> floatvec;
   
   typedef std::vector<std::vector<int>> intmat;
   typedef std::vector<int> intvec;
   
   namespace blackbird {
   
       enum class Device {
           Chip0,
           Gaussian,
           Fock
       };
   
       enum class Gate {
           // State preparations
           Vacuum, Coherent, Squeezed, Thermal, Fock, Catstate,
           // one mode gates
           Rgate, Dgate, Xgate, Zgate, Sgate, Pgate, Vgate,
           // two mode gates
           BSgate, S2gate, CXgate, CZgate, CKgate,
           // channels
           LossChannel, ThermalLossChannel,
           // decompositions
           Interferometer, GaussianTransform, Gaussian,
           // measurements
           MeasureFock, MeasureHomodyne, MeasureHeterodyne, MeasureIntensity};
   
       enum class ParDomain {
           None,
           Int,
           Float,
           Complex,
           String,
           Bool,
           ArrayFloat,
           ArrayComplex
       };
   
       //================================
       // Operation definitions
       //================================
   
       class Operation {
           public:
               Gate gate;
               intvec modes;
   
               std::string name;
               int num_params;
               ParDomain domain;
               int num_modes;
   
               int i1;
   
               double f1;
               double f2;
               double f3;
   
               std::complex<double> c1;
   
               complexmat U1;
   
               floatmat S1;
               floatmat S2;
   
               void print_op() {
                   std::cout << name << "() | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               template <typename T>
               void print_op(T p1) {
                   std::cout << name << "(" << p1 << ") | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               template <typename T, typename S>
               void print_op(T p1, S p2) {
                   std::cout << name << "(" << p1 << ", " << p2 << ") | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               template <typename T, typename S, typename O>
               void print_op(T p1, S p2, O p3) {
                   std::cout << name << "(" << p1 << ", " << p2 << ", " << p3 << ") | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               template <typename T>
               void print_array_op(T p1) {
                   std::cout << name << "(" << std::endl;
                   for (auto j : p1) {
                       for (auto k : j) {
                           std::cout << k << ',' << ' ';
                       }
                       std::cout << std::endl;
                   }
                   std::cout << ") | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               template <typename T, typename S>
               void print_array_op(T p1, S p2) {
                   std::cout << name << "(" << std::endl;
                   for (auto j : p1) {
                       for (auto k : j) {
                           std::cout << k << ',' << ' ';
                       }
                       std::cout << std::endl;
                   }
                   std::cout << ", " << std::endl;
                   for (auto j : p2) {
                       for (auto k : j) {
                           std::cout << k << ',' << ' ';
                       }
                       std::cout << std::endl;
                   }
                   std::cout << ") | ";
                   for (auto j : modes) {
                       std::cout << j << ',' << ' ';
                   }
                   std::cout << std::endl;
               };
   
               void check_num_args(int args_size) {
                   if (args_size != num_params) {
                       std::string np = std::to_string(num_params);
                       std::string as = std::to_string(args_size);
                       throw std::invalid_argument(name+" is a "+np+" parameter operation, but "
                           +as+" parameters were given.");
                   }
               };
   
               void check_num_modes(int mode_size) {
                   if (mode_size != num_modes) {
                       std::string nm = std::to_string(num_modes);
                       std::string ms = std::to_string(mode_size);
                       throw std::invalid_argument(name+" is a "+nm+" mode operation, but "
                           +ms+" modes were given.");
                   }
               };
       };
   
   
       //================================
       // State preparation
       //================================
   
       class Vacuum : public Operation {
           public:
               Vacuum(intvec m) {
                   gate = Gate::Vacuum;
                   name = "Vacuum";
                   num_params = 0;
                   num_modes = 1;
                   domain = ParDomain::None;
   
                   // put input validation here
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       class Coherent : public Operation {
           public:
               Coherent(floatvec displacement_phase, intvec m) {
                   gate = Gate::Coherent;
                   name = "Coherent";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = displacement_phase[0];
                   f2 = displacement_phase[1];
   
                   // put input validation here
                   int args_size = static_cast<int>(displacement_phase.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Coherent(complexvec alpha, intvec m) {
                   gate = Gate::Coherent;
                   name = "Coherent";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = std::abs(alpha[0]);
                   f2 = std::arg(alpha[0]);
   
                   // put input validation here
                   int args_size = static_cast<int>(alpha.size())+1;
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Coherent(intvec alpha, intvec m) {
                   throw std::invalid_argument("Coherent must have real arguments");
               };
       };
   
   
       class Squeezed : public Operation {
           public:
               Squeezed(floatvec r_phase, intvec m) {
                   gate = Gate::Squeezed;
                   name = "Squeezed";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   if (r_phase.size() == 1) {
                       r_phase.push_back(0);
                   }
                   f1 = r_phase[0];
                   f2 = r_phase[1];
   
                   // put input validation here
                   int args_size = static_cast<int>(r_phase.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Squeezed(complexvec alpha, intvec m) {
                   throw std::invalid_argument("Squeezed must have real arguments");
               };
               Squeezed(intvec alpha, intvec m) {
                   throw std::invalid_argument("Squeezed must have real arguments");
               };
       };
   
   
       class Thermal : public Operation {
           public:
               Thermal(floatvec nbar, intvec m) {
                   gate = Gate::Thermal;
                   name = "Thermal";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = nbar[0];
   
                   // put input validation here
                   if (f1 < 0) {
                       throw std::invalid_argument("Thermal state: cannot have negative thermal population.");
                   }
                   int args_size = static_cast<int>(nbar.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Thermal(complexvec nbar, intvec m) {
                   throw std::invalid_argument("Thermal must have real arguments");
               };
               Thermal(intvec nbar, intvec m) {
                   throw std::invalid_argument("Thermal must have real arguments");
               };
       };
   
   
       class Fock : public Operation {
           public:
               Fock(intvec n, intvec m) {
                   gate = Gate::Fock;
                   name = "Fock";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Int;
   
                   i1 = n[0];
   
                   // put input validation here
                   if (i1 < 0) {
                       throw std::invalid_argument("Fock state: cannot have negative Fock state.");
                   }
                   int args_size = static_cast<int>(n.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Fock(floatvec n, intvec m) {
                   throw std::invalid_argument("Fock must have integer arguments");
               };
   
               Fock(complexvec n, intvec m) {
                   throw std::invalid_argument("Fock must have integer arguments");
               };
       };
   
       class Catstate : public Operation {
           // Cat state
           public:
               Catstate(complexvec alpha, intvec m, float parity=0.) {
                   gate = Gate::Catstate;
                   name = "Catstate";
                   num_params = 3;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = std::abs(alpha[0]);
                   f2 = std::arg(alpha[0]);
                   f3 = parity;
   
                   // put input validation here
                   int args_size = static_cast<int>(alpha.size())+2;
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Catstate(floatvec alpha, intvec m) {
                   throw std::invalid_argument("Catstate alpha parameter must be complex");
               };
               Catstate(intvec alpha, intvec m) {
                   throw std::invalid_argument("Catstate alpha parameter must be complex");
               };
       };
   
   
       //================================
       // Single mode gates
       //================================
   
       class Rgate : public Operation {
           // Rotation gate
           public:
               Rgate(complexvec phi, intvec m) {
                   throw std::invalid_argument("Rgate must have real arguments");
               };
               Rgate(intvec phi, intvec m) {
                   throw std::invalid_argument("Rgate must have real arguments");
               };
               Rgate(floatvec phi, intvec m) {
                   gate = Gate::Rgate;
                   name = "Rgate";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = phi[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(phi.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class Dgate : public Operation {
           public:
               Dgate(floatvec displacement_phase, intvec m) {
                   gate = Gate::Dgate;
                   name = "Dgate";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = displacement_phase[0];
                   f2 = displacement_phase[1];
   
                   // put input validation here
                   int args_size = static_cast<int>(displacement_phase.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Dgate(complexvec alpha, intvec m) {
                   gate = Gate::Dgate;
                   name = "Dgate";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = std::abs(alpha[0]);
                   f2 = std::arg(alpha[0]);
   
                   // put input validation here
                   int args_size = static_cast<int>(alpha.size())+1;
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
               Dgate(intvec x, intvec m) {
                   throw std::invalid_argument("Dgate must have real arguments");
               };
       };
   
       class Xgate : public Operation {
           public:
               Xgate(complexvec x, intvec m) {
                   throw std::invalid_argument("Xgate must have real arguments");
               };
               Xgate(intvec x, intvec m) {
                   throw std::invalid_argument("Xgate must have real arguments");
               };
               Xgate(floatvec x, intvec m) {
                   gate = Gate::Xgate;
                   name = "Xgate";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = x[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(x.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class Zgate : public Operation {
           public:
               Zgate(complexvec z, intvec m) {
                   throw std::invalid_argument("Zgate must have real arguments");
               };
               Zgate(intvec z, intvec m) {
                   throw std::invalid_argument("Zgate must have real arguments");
               };
               Zgate(floatvec z, intvec m) {
                   gate = Gate::Zgate;
                   name = "Zgate";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = z[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(z.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class Sgate : public Operation {
           public:
               Sgate(floatvec r_phase, intvec m) {
                   gate = Gate::Sgate;
                   name = "Sgate";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   if (r_phase.size() == 1) {
                       r_phase.push_back(0);
                   }
                   f1 = r_phase[0];
                   f2 = r_phase[1];
   
                   // put input validation here
                   int args_size = static_cast<int>(r_phase.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               Sgate(complexvec alpha, intvec m) {
                   throw std::invalid_argument("Sgate must have real arguments");
               };
               Sgate(intvec alpha, intvec m) {
                   throw std::invalid_argument("Sgate must have real arguments");
               };
       };
   
       class Pgate : public Operation {
           public:
               Pgate(complexvec s, intvec m) {
                   throw std::invalid_argument("Pgate must have real arguments");
               };
               Pgate(intvec s, intvec m) {
                   throw std::invalid_argument("Pgate must have real arguments");
               };
               Pgate(floatvec s, intvec m) {
                   gate = Gate::Pgate;
                   name = "Pgate";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = s[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(s.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class Vgate : public Operation {
           public:
               Vgate(complexvec s, intvec m) {
                   throw std::invalid_argument("Vgate must have real arguments");
               };
               Vgate(intvec s, intvec m) {
                   throw std::invalid_argument("Vgate must have real arguments");
               };
               Vgate(floatvec s, intvec m) {
                   gate = Gate::Vgate;
                   name = "Vgate";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = s[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(s.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       // ==================================
       // Multimode gates
       // ==================================
   
       class BSgate : public Operation {
           public:
               BSgate(complexvec theta_phi, intvec m) {
                   throw std::invalid_argument("BSgate must have real arguments");
               };
               BSgate(intvec theta_phi, intvec m) {
                   throw std::invalid_argument("BSgate must have real arguments");
               };
               BSgate(floatvec theta_phi, intvec m) {
                   gate = Gate::BSgate;
                   name = "BSgate";
                   num_params = 2;
                   num_modes = 2;
                   domain = ParDomain::Float;
   
                   f1 = theta_phi[0];
                   f2 = theta_phi[1];
   
                   // put input validation here
                   int args_size = static_cast<int>(theta_phi.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       class S2gate : public Operation {
           public:
               S2gate(complexvec r_phi, intvec m) {
                   throw std::invalid_argument("S2gate must have real arguments");
               };
               S2gate(intvec r_phi, intvec m) {
                   throw std::invalid_argument("S2gate must have real arguments");
               };
               S2gate(floatvec r_phi, intvec m) {
                   gate = Gate::S2gate;
                   name = "S2gate";
                   num_params = 2;
                   num_modes = 2;
                   domain = ParDomain::Float;
   
                   int args_size = static_cast<int>(r_phi.size());
   
                   if (r_phi.size() == 1) {
                       r_phi.push_back(0);
                       args_size = 2;
                   }
   
                   f1 = r_phi[0];
                   f2 = r_phi[1];
   
                   // put input validation here
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class CXgate : public Operation {
           public:
               CXgate(complexvec s, intvec m) {
                   throw std::invalid_argument("CXgate must have real arguments");
               };
               CXgate(intvec s, intvec m) {
                   throw std::invalid_argument("CXgate must have real arguments");
               };
               CXgate(floatvec s, intvec m) {
                   gate = Gate::CXgate;
                   name = "CXgate";
                   num_params = 1;
                   num_modes = 2;
                   domain = ParDomain::Float;
   
                   f1 = s[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(s.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       class CZgate : public Operation {
           // Beamsplitter
           public:
               CZgate(complexvec s, intvec m) {
                   throw std::invalid_argument("CZgate must have real arguments");
               };
               CZgate(intvec s, intvec m) {
                   throw std::invalid_argument("CZgate must have real arguments");
               };
               CZgate(floatvec s, intvec m) {
                   gate = Gate::CZgate;
                   name = "CZgate";
                   num_params = 1;
                   num_modes = 2;
                   domain = ParDomain::Float;
   
                   f1 = s[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(s.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       class CKgate : public Operation {
           // Beamsplitter
           public:
               CKgate(complexvec s, intvec m) {
                   throw std::invalid_argument("CKgate must have real arguments");
               };
               CKgate(intvec s, intvec m) {
                   throw std::invalid_argument("CKgate must have real arguments");
               };
               CKgate(floatvec s, intvec m) {
                   gate = Gate::CKgate;
                   name = "CKgate";
                   num_params = 1;
                   num_modes = 2;
                   domain = ParDomain::Float;
   
                   f1 = s[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(s.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       // ==================================
       // Loss channels
       // ==================================
   
       class LossChannel : public Operation {
           public:
               LossChannel(complexvec phi, intvec m) {
                   throw std::invalid_argument("LossChannel must have real arguments");
               };
               LossChannel(intvec phi, intvec m) {
                   throw std::invalid_argument("LossChannel must have real arguments");
               };
               LossChannel(floatvec T, intvec m) {
                   gate = Gate::LossChannel;
                   name = "LossChannel";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = T[0];
   
                   // put input validation here
                   if (0 <= f1 <= 1) {
                       throw std::invalid_argument("LossChannel: loss parameter must be between 0 and 1");
                   }
                   int args_size = static_cast<int>(T.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       class ThermalLossChannel : public Operation {
           public:
               ThermalLossChannel(complexvec phi, intvec m) {
                   throw std::invalid_argument("ThermalLossChannel must have real arguments");
               };
               ThermalLossChannel(intvec phi, intvec m) {
                   throw std::invalid_argument("ThermalLossChannel must have real arguments");
               };
               ThermalLossChannel(floatvec T, intvec m) {
                   gate = Gate::ThermalLossChannel;
                   name = "ThermalLossChannel";
                   num_params = 2;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = T[0];
                   f2 = T[1];
   
                   // put input validation here
                   if (0 <= f1 <= 1) {
                       throw std::invalid_argument("ThermalLossChannel: loss parameter must be between 0 and 1");
                   }
                   if (f2 < 0) {
                       throw std::invalid_argument("ThermalLossChannel: thermal population must be positive");
                   }
                   int args_size = static_cast<int>(T.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       // ==================================
       // Decomposition gates
       // ==================================
   
       class Interferometer : public Operation {
           public:
               Interferometer(complexmat unitary, intvec m) {
                   gate = Gate::Interferometer;
                   name = "Interferometer";
                   num_params = 1;
                   num_modes = 0;
                   domain = ParDomain::ArrayComplex;
   
                   U1 = unitary;
   
                   // put input validation here
                   int U_rows = static_cast<int>(U1.size());
                   int U_cols = static_cast<int>(U1[0].size());
                   int mode_size = static_cast<int>(m.size());
   
                   if (U_rows != U_cols) {
                       throw std::invalid_argument("Interferometer unitary matrix should be square");
                   }
                   else if (U_rows != mode_size) {
                       throw std::invalid_argument("Interferometer unitary matrix must have the same size as the number of modes");
                   }
                   modes = m;
               };
       };
   
   
       class GaussianTransform : public Operation {
           public:
               GaussianTransform(floatmat symplectic, intvec m) {
                   gate = Gate::GaussianTransform;
                   name = "GaussianTransform";
                   num_params = 1;
                   num_modes = 0;
                   domain = ParDomain::ArrayFloat;
   
                   S1 = symplectic;
   
                   // put input validation here
                   int rows = static_cast<int>(S1.size());
                   int cols = static_cast<int>(S1[0].size());
                   int mode_size = static_cast<int>(m.size());
   
                   if (rows != cols) {
                       throw std::invalid_argument("GaussianTransform symplectic matrix should be square");
                   }
                   else if (rows != 2*mode_size) {
                       throw std::invalid_argument("GaussianTransform symplectic matrix should have size double the number of modes");
                   }
                   modes = m;
               };
       };
   
   
       class Gaussian : public Operation {
           public:
               Gaussian(floatmat cov, floatmat means, intvec m) {
                   gate = Gate::Gaussian;
                   name = "Gaussian";
                   num_params = 2;
                   num_modes = 0;
                   domain = ParDomain::ArrayFloat;
   
                   S1 = cov;
                   S2 = means;
   
                   // put input validation here
                   int rows = static_cast<int>(S1.size());
                   int cols = static_cast<int>(S1[0].size());
                   int mode_size = static_cast<int>(m.size());
   
                   if (rows != cols) {
                       throw std::invalid_argument("Gaussian: covariance matrix should be square");
                   }
                   else if (rows != 2*mode_size) {
                       throw std::invalid_argument("Gaussian: covariance matrix should have size double the number of modes");
                   }
   
                   rows = static_cast<int>(S2.size());
                   cols = static_cast<int>(S2[0].size());
   
                   if (rows != 1) {
                       throw std::invalid_argument("Gaussian: means vector should only have 1 row");
                   }
                   else if (cols != 2*mode_size) {
                       throw std::invalid_argument("Gaussian: means vector should have size double the number of modes");
                   }
   
                   modes = m;
               };
   
               Gaussian(floatmat cov, intvec m) {
                   gate = Gate::Gaussian;
                   name = "Gaussian";
                   num_params = 2;
                   num_modes = 0;
                   domain = ParDomain::ArrayFloat;
   
                   S1 = cov;
   
                   // put input validation here
                   int rows = static_cast<int>(S1.size());
                   int cols = static_cast<int>(S1[0].size());
                   int mode_size = static_cast<int>(m.size());
   
                   floatmat means_vector(1, floatvec(2*mode_size, 0.0));
                   S2 = means_vector;
   
                   if (rows != cols) {
                       throw std::invalid_argument("Gaussian: covariance matrix should be square");
                   }
                   else if (rows != 2*mode_size) {
                       throw std::invalid_argument("Gaussian: covariance matrix should have size double the number of modes");
                   }
   
                   modes = m;
               };
       };
   
       // ==================================
       // Measurements
       // ==================================
   
       class MeasureFock : public Operation {
           public:
               MeasureFock(intvec m) {
                   gate = Gate::MeasureFock;
                   name = "MeasureFock";
                   num_params = 0;
                   num_modes = 1;
                   domain = ParDomain::None;
   
                   // put input validation here
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class MeasureHomodyne : public Operation {
           public:
               MeasureHomodyne(intvec m) {
                   gate = Gate::MeasureHomodyne;
                   name = "MeasureHomodyne";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = 0.0;
   
                   // put input validation here
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               MeasureHomodyne(floatvec phi, intvec m) {
                   gate = Gate::MeasureHomodyne;
                   name = "MeasureHomodyne";
                   num_params = 1;
                   num_modes = 1;
                   domain = ParDomain::Float;
   
                   f1 = phi[0];
   
                   // put input validation here
                   int args_size = static_cast<int>(phi.size());
                   check_num_args(args_size);
   
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
   
               MeasureHomodyne(complexvec phi, intvec m) {
                   throw std::invalid_argument("MeasureHomodyne must have real arguments");
               };
               MeasureHomodyne(intvec phi, intvec m) {
                   throw std::invalid_argument("MeasureHomodyne must have real arguments");
               };
       };
   
       class MeasureHeterodyne : public Operation {
           public:
               MeasureHeterodyne(intvec m) {
                   gate = Gate::MeasureHeterodyne;
                   name = "MeasureHeterodyne";
                   num_params = 0;
                   num_modes = 1;
                   domain = ParDomain::None;
   
                   // put input validation here
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
       class MeasureIntensity : public Operation {
           public:
               MeasureIntensity(intvec m) {
                   gate = Gate::MeasureIntensity;
                   name = "MeasureIntensity";
                   num_params = 0;
                   num_modes = 1;
                   domain = ParDomain::None;
   
                   // put input validation here
                   int mode_size = static_cast<int>(m.size());
                   check_num_modes(mode_size);
                   modes = m;
               };
       };
   
   
       //================================
       // Device definitions
       //================================
   
       class Program {
           public:
               Device name;
               std::vector<Operation*> operations;
   
               virtual void print_device_info() {};
   
               void print_operations() {
                   for (auto i : operations) {
                       if (i->num_params == 0) {
                           i->print_op();
                       }
                       else if (i->num_params == 1) {
                           if (i->domain == ParDomain::Float) {
                               i->print_op(i->f1);
                           }
                           else if (i->domain == ParDomain::Int) {
                               i->print_op(i->i1);
                           }
                           else if (i->domain == ParDomain::Complex) {
                               i->print_op(i->c1);
                           }
                           else if (i->domain == ParDomain::ArrayComplex) {
                               i->print_array_op(i->U1);
                           }
                           else if (i->domain == ParDomain::ArrayFloat) {
                               i->print_array_op(i->S1);
                           }
                       }
                       else if (i->num_params == 2) {
                           if (i->domain == ParDomain::Float) {
                               i->print_op(i->f1, i->f2);
                           }
                           else if (i->domain == ParDomain::ArrayFloat) {
                               i->print_array_op(i->S1, i->S2);
                           }
                       }
                       else if (i->num_params == 3) {
                           i->print_op(i->f1, i->f2, i->f3);
                       }
                   }
               };
       };
   
   
       class Chip0 : public Program {
           public:
               int shots = 1;
   
               Chip0(Program dev) {
                   name = Device::Chip0;
               };
               Chip0(int s=1) {
                   name = Device::Chip0;
                   shots = s;
               };
               void print_device_info() {
                   std::cout << "Device: Chip0" << std::endl;
                   std::cout << "Shots: " << shots << std::endl;
               };
       };
   
   
       class GaussianSimulator : public Program {
           public:
               int ns;
               int shots = 1;
               double hb = 2.0;
   
               GaussianSimulator() {
                   name = Device::Gaussian;
               };
               GaussianSimulator(Program dev) {
                   name = Device::Gaussian;
               };
               GaussianSimulator(int num_subsystems, int s=0, double hbar=2) {
                   name = Device::Gaussian;
                   ns = num_subsystems;
                   shots = s;
                   hb = hbar;
               };
   
               void print_device_info() {
                   std::cout << "Device: Gaussian (simulator)" << std::endl;
                   std::cout << "Number of subsystems: " << ns << std::endl;
                   std::cout << "Shots: " << shots << std::endl;
                   std::cout << "hbar: " << hb << std::endl;
               };
       };
   
   
       class FockSimulator : public Program {
           public:
               int ns;
               int cutoff;
               int shots = 1;
               double hb = 2.0;
   
               FockSimulator() {
                   name = Device::Fock;
               };
               FockSimulator(Program dev) {
                   name = Device::Fock;
               };
               FockSimulator(int num_subsystems, int cutoff_dim, int s=0, double hbar=2) {
                   name = Device::Fock;
                   cutoff = cutoff_dim;
                   ns = num_subsystems;
                   shots = s;
                   hb = hbar;
               };
   
               void print_device_info() {
                   std::cout << "Device: Fock (simulator)" << std::endl;
                   std::cout << "Number of subsystems: " << ns << std::endl;
                   std::cout << "Cutoff dimension: " << cutoff << std::endl;
                   std::cout << "Shots: " << shots << std::endl;
                   std::cout << "hbar: " << hb << std::endl;
               };
       };
   }