Skip to content

Commit

Permalink
Renames resolution to cutoff
Browse files Browse the repository at this point in the history
  • Loading branch information
@nquesada
nquesada committed Jan 17, 2020
1 parent 0a5406a commit b8c9059
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 65 deletions.
118 changes: 59 additions & 59 deletions include/hermite_multidimensional.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,17 +29,17 @@ typedef unsigned long long int ullint;
* Returns the index of the one dimensional flattened vector corresponding to the multidimensional tensor
*
* @param pos
* @param resolution
* @param cutoff
*
* @return index on flattened vector
*/
ullint vec2index(std::vector<int> &pos, int resolution) {
ullint vec2index(std::vector<int> &pos, int cutoff) {
int dim = pos.size();
ullint nextCoordinate = 0;

nextCoordinate = pos[0]-1;
for(int ii = 0; ii < dim-1; ii++) {
nextCoordinate = nextCoordinate*resolution + (pos[ii+1]-1);
nextCoordinate = nextCoordinate*cutoff + (pos[ii+1]-1);
}

return nextCoordinate;
Expand All @@ -52,12 +52,12 @@ ullint vec2index(std::vector<int> &pos, int resolution) {
* @param nextPos a vector of integers
* @param jumpFrom a vector of integers
* @param jump integer specifying whether to jump to the next index
* @param resolution integer specifying the cuotff
* @param cutoff integer specifying the cuotff
* @dim dimension of the R matrix
*
* @k index necessary for knowing which elements are needed from the input vector y and matrix R
*/
int update_iterator(std::vector<int> &nextPos, std::vector<int> &jumpFrom, int &jump, const int &resolution, const int &dim) {
int update_iterator(std::vector<int> &nextPos, std::vector<int> &jumpFrom, int &jump, const int &cutoff, const int &dim) {
if (jump > 0) {
jumpFrom[jump] += 1;
jump = 0;
Expand All @@ -66,7 +66,7 @@ int update_iterator(std::vector<int> &nextPos, std::vector<int> &jumpFrom, int &
std::vector<int> forwardStep(dim, 0);
forwardStep[ii] = 1;

if ( forwardStep[ii] + nextPos[ii] > resolution) {
if ( forwardStep[ii] + nextPos[ii] > cutoff) {
nextPos[ii] = 1;
jumpFrom[ii] = 1;
jump = ii+1;
Expand Down Expand Up @@ -95,14 +95,14 @@ namespace libwalrus {
* @param R a flattened vector of size \f$n^2\f$, representing a
* \f$n\times n\f$ symmetric matrix.
* @param y a flattened vector of size \f$n\f$.
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/
template <typename T>
inline T* hermite_multidimensional_cpp(const std::vector<T> &R, const std::vector<T> &y, const int &resolution) {
inline T* hermite_multidimensional_cpp(const std::vector<T> &R, const std::vector<T> &y, const int &cutoff) {
int dim = std::sqrt(static_cast<double>(R.size()));

ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);
Expand All @@ -116,10 +116,10 @@ inline T* hermite_multidimensional_cpp(const std::vector<T> &R, const std::vecto

for (ullint jj = 0; jj < Hdim-1; jj++) {

k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);

nextCoordinate = vec2index(nextPos, resolution);
fromCoordinate = vec2index(jumpFrom, resolution);
nextCoordinate = vec2index(nextPos, cutoff);
fromCoordinate = vec2index(jumpFrom, cutoff);

H[nextCoordinate] = H[fromCoordinate] * y[k];

Expand All @@ -130,7 +130,7 @@ inline T* hermite_multidimensional_cpp(const std::vector<T> &R, const std::vecto
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] - (static_cast<T>(jumpFrom[ii]-1))*(R[dim*k+ii])*H[prevCoordinate];

}
Expand All @@ -152,20 +152,20 @@ inline T* hermite_multidimensional_cpp(const std::vector<T> &R, const std::vecto
* @param R a flattened vector of size \f$n^2\f$, representing a
* \f$n\times n\f$ symmetric matrix.
* @param y a flattened vector of size \f$n\f$.
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/
template <typename T>
inline T* renorm_hermite_multidimensional_cpp(const std::vector<T> &R, const std::vector<T> &y, const int &resolution) {
inline T* renorm_hermite_multidimensional_cpp(const std::vector<T> &R, const std::vector<T> &y, const int &cutoff) {
int dim = std::sqrt(static_cast<double>(R.size()));

ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);
H[0] = 1;
std::vector<double> intsqrt(resolution+1, 0);
for (int ii = 0; ii<=resolution; ii++) {
std::vector<double> intsqrt(cutoff+1, 0);
for (int ii = 0; ii<=cutoff; ii++) {
intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
}
std::vector<int> nextPos(dim, 1);
Expand All @@ -174,10 +174,10 @@ inline T* renorm_hermite_multidimensional_cpp(const std::vector<T> &R, const st
int k;
ullint nextCoordinate, fromCoordinate;
for (ullint jj = 0; jj < Hdim-1; jj++) {
k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);

nextCoordinate = vec2index(nextPos, resolution);
fromCoordinate = vec2index(jumpFrom, resolution);
nextCoordinate = vec2index(nextPos, cutoff);
fromCoordinate = vec2index(jumpFrom, cutoff);

H[nextCoordinate] = H[fromCoordinate] * y[k];

Expand All @@ -188,7 +188,7 @@ inline T* renorm_hermite_multidimensional_cpp(const std::vector<T> &R, const st
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] - intsqrt[jumpFrom[ii]-1]*(R[k*dim+ii])*H[prevCoordinate];
}
}
Expand All @@ -202,22 +202,22 @@ inline T* renorm_hermite_multidimensional_cpp(const std::vector<T> &R, const st
*
* @param R a flattened vector of size \f$n^2\f$, representing a
* \f$n\times n\f$ symmetric matrix.
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/
template <typename T>
inline T* interferometer_cpp(const std::vector<T> &R, const int &resolution) {
inline T* interferometer_cpp(const std::vector<T> &R, const int &cutoff) {

int dim = std::sqrt(static_cast<double>(R.size()));
assert(dim % 2 == 0);
int num_modes = dim/2;
ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);
H[0] = 1;
std::vector<double> intsqrt(resolution+1, 0);
for (int ii = 0; ii<=resolution; ii++) {
std::vector<double> intsqrt(cutoff+1, 0);
for (int ii = 0; ii<=cutoff; ii++) {
intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
}
std::vector<int> nextPos(dim, 1);
Expand All @@ -227,7 +227,7 @@ inline T* interferometer_cpp(const std::vector<T> &R, const int &resolution) {
ullint nextCoordinate, fromCoordinate;

for (ullint jj = 0; jj < Hdim-1; jj++) {
k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);
int bran = 0;
for (int ii=0; ii < num_modes; ii++) {
bran += nextPos[ii];
Expand All @@ -238,8 +238,8 @@ inline T* interferometer_cpp(const std::vector<T> &R, const int &resolution) {
ketn += nextPos[ii];
}
if (bran == ketn) {
nextCoordinate = vec2index(nextPos, resolution);
fromCoordinate = vec2index(jumpFrom, resolution);
nextCoordinate = vec2index(nextPos, cutoff);
fromCoordinate = vec2index(jumpFrom, cutoff);
std::vector<int> tmpjump(dim, 0);
int low_lim;
int high_lim;
Expand All @@ -258,7 +258,7 @@ inline T* interferometer_cpp(const std::vector<T> &R, const int &resolution) {
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] - (intsqrt[jumpFrom[ii]-1])*(R[k*dim+ii])*H[prevCoordinate];
}
}
Expand All @@ -273,19 +273,19 @@ inline T* interferometer_cpp(const std::vector<T> &R, const int &resolution) {
*
* @param R a flattened vector of size 4, representing a
* \f$2\times 2\f$ symmetric matrix.
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/
template <typename T>
inline T* squeezing_cpp(const std::vector<T> &R, const int &resolution) {
inline T* squeezing_cpp(const std::vector<T> &R, const int &cutoff) {
int dim = std::sqrt(static_cast<double>(R.size()));
ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);
H[0] = std::sqrt(-R[1]);
std::vector<double> intsqrt(resolution+1, 0);
for (int ii = 0; ii<=resolution; ii++) {
std::vector<double> intsqrt(cutoff+1, 0);
for (int ii = 0; ii<=cutoff; ii++) {
intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
}
assert(dim == 2);
Expand All @@ -297,21 +297,21 @@ inline T* squeezing_cpp(const std::vector<T> &R, const int &resolution) {
int k;
ullint nextCoordinate, fromCoordinate;
for (ullint jj = 0; jj < Hdim-1; jj++) {
k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);

int bran = nextPos[0];
int ketn = nextPos[1];
if (bran % 2 == ketn % 2) {
ullint nextCoordinate = vec2index(nextPos, resolution);
ullint fromCoordinate = vec2index(jumpFrom, resolution);
ullint nextCoordinate = vec2index(nextPos, cutoff);
ullint fromCoordinate = vec2index(jumpFrom, cutoff);

std::vector<int> tmpjump(dim, 0);
for (int ii = 0; ii < dim; ii++) {
if (jumpFrom[ii] > 1) {
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] - (intsqrt[jumpFrom[ii]-1])*(R[k*dim+ii])*H[prevCoordinate];
}
}
Expand All @@ -324,23 +324,23 @@ inline T* squeezing_cpp(const std::vector<T> &R, const int &resolution) {


/**
* Returns the matrix elements of a displacement operation parametrized in terms of its double vector y
* Returns the matrix elements of a displacement operation
*
* @param y a flattened vector of size \f$2\f$, represeting the displacement via \f$\alpha, \alpha^*\f$
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/
template <typename T>
inline T* displacement_cpp(const std::vector<T> &y, const int &resolution) {
inline T* displacement_cpp(const std::vector<T> &y, const int &cutoff) {
int dim = 2;

ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);
H[0] = std::exp(-0.5*std::pow(std::abs(y[0]),2));
std::vector<double> intsqrt(resolution+1, 0);
for (int ii = 0; ii<=resolution; ii++) {
std::vector<double> intsqrt(cutoff+1, 0);
for (int ii = 0; ii<=cutoff; ii++) {
intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
}
std::vector<int> nextPos(dim, 1);
Expand All @@ -349,10 +349,10 @@ inline T* displacement_cpp(const std::vector<T> &y, const int &resolution) {
int k;
ullint nextCoordinate, fromCoordinate;
for (ullint jj = 0; jj < Hdim-1; jj++) {
k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);

nextCoordinate = vec2index(nextPos, resolution);
fromCoordinate = vec2index(jumpFrom, resolution);
nextCoordinate = vec2index(nextPos, cutoff);
fromCoordinate = vec2index(jumpFrom, cutoff);

H[nextCoordinate] = H[fromCoordinate] * y[k];

Expand All @@ -368,7 +368,7 @@ inline T* displacement_cpp(const std::vector<T> &y, const int &resolution) {
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] + (intsqrt[jumpFrom[ii]-1])*H[prevCoordinate];

}
Expand All @@ -384,24 +384,24 @@ inline T* displacement_cpp(const std::vector<T> &y, const int &resolution) {
*
* @param R a flattened vector of size \f$n^2\f$, representing a
* \f$n\times n\f$ symmetric matrix.
* @param resolution highest number of photons to be resolved.
* @param cutoff highest number of photons to be resolved.
*
*/

template <typename T>
inline T* two_mode_squeezing_cpp(const std::vector<T> &R, const int &resolution) {
inline T* two_mode_squeezing_cpp(const std::vector<T> &R, const int &cutoff) {
int dim = std::sqrt(static_cast<double>(R.size()));
assert(dim == 4);

ullint Hdim = pow(resolution, dim);
ullint Hdim = pow(cutoff, dim);
T *H;
H = (T*) malloc(sizeof(T)*Hdim);
memset(&H[0],0,sizeof(T)*Hdim);


H[0] = -R[2];
std::vector<double> intsqrt(resolution+1, 0);
for (int ii = 0; ii<=resolution; ii++) {
std::vector<double> intsqrt(cutoff+1, 0);
for (int ii = 0; ii<=cutoff; ii++) {
intsqrt[ii] = std::sqrt((static_cast<double>(ii)));
}
std::vector<int> nextPos(dim, 1);
Expand All @@ -410,12 +410,12 @@ inline T* two_mode_squeezing_cpp(const std::vector<T> &R, const int &resolution)
int k;
ullint nextCoordinate, fromCoordinate;
for (ullint jj = 0; jj < Hdim-1; jj++) {
k = update_iterator(nextPos, jumpFrom, jump, resolution, dim);
k = update_iterator(nextPos, jumpFrom, jump, cutoff, dim);
int bran = nextPos[0]-nextPos[1];
int ketn = nextPos[2]-nextPos[3];
if (bran == ketn) {
ullint nextCoordinate = vec2index(nextPos, resolution);
ullint fromCoordinate = vec2index(jumpFrom, resolution);
ullint nextCoordinate = vec2index(nextPos, cutoff);
ullint fromCoordinate = vec2index(jumpFrom, cutoff);

std::vector<int> tmpjump(dim, 0);

Expand All @@ -424,7 +424,7 @@ inline T* two_mode_squeezing_cpp(const std::vector<T> &R, const int &resolution)
std::vector<int> prevJump(dim, 0);
prevJump[ii] = 1;
std::transform(jumpFrom.begin(), jumpFrom.end(), prevJump.begin(), tmpjump.begin(), std::minus<int>());
ullint prevCoordinate = vec2index(tmpjump, resolution);
ullint prevCoordinate = vec2index(tmpjump, cutoff);
H[nextCoordinate] = H[nextCoordinate] - (intsqrt[jumpFrom[ii]-1])*(R[k*dim+ii])*H[prevCoordinate];

}
Expand Down
12 changes: 6 additions & 6 deletions thewalrus/libwalrus.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,12 +141,12 @@ cdef extern from "../include/libwalrus.hpp" namespace "libwalrus":
double complex torontonian_quad(vector[double complex] &mat)
double torontonian_fsum[T](vector[T] &mat)

T* hermite_multidimensional_cpp[T](vector[T] &mat, vector[T] &d, int &resolution)
T* renorm_hermite_multidimensional_cpp[T](vector[T] &mat, vector[T] &d, int &resolution)
T* interferometer_cpp[T](vector[T] &mat, int &resolution)
T* squeezing_cpp[T](vector[T] &mat, int &resolution)
T* displacement_cpp[T](vector[T] &y, int &resolution)
T* two_mode_squeezing_cpp[T](vector[T] &y, int &resolution)
T* hermite_multidimensional_cpp[T](vector[T] &mat, vector[T] &d, int &cutoff)
T* renorm_hermite_multidimensional_cpp[T](vector[T] &mat, vector[T] &d, int &cutoff)
T* interferometer_cpp[T](vector[T] &mat, int &cutoff)
T* squeezing_cpp[T](vector[T] &mat, int &cutoff)
T* displacement_cpp[T](vector[T] &y, int &cutoff)
T* two_mode_squeezing_cpp[T](vector[T] &y, int &cutoff)


# ==============================================================================
Expand Down

0 comments on commit b8c9059

Please sign in to comment.