Hi, a quick search in c++ documentation did not help: is there an MPS sampling function in c++ like in julia the “sample” command?
Thanks for the information!
Ceren.
1 Like
I don’t think there’s something built in but there’s some code in the METTS tutorial to sample in the Z or X basis:
1 Like
This sampling code does not work for Boson site when I set it to hard-core boson limit such that the code would be applicable. I get this error which I think originates from Real prob_up = elt(dag(prime(psi(j),“Site”))PUppsi(j)); line in the code:
IndexSet 1 =
(dim=1|id=767|“l=1,Link”)
1: 1 QN({“Nb”,15})
(dim=2|id=979|“n=1,Site,Boson”)
1: 1 QN({“Nb”,0})
2: 1 QN({“Nb”,1})
IndexSet 2 =
(dim=1|id=767|“l=1,Link”)
1: 1 QN({“Nb”,15})
(dim=2|id=979|“n=1,Site,Boson”)
1: 1 QN({“Nb”,0})
2: 1 QN({“Nb”,1})
Mismatched QN Index from set 1 (dim=2|id=979|“n=1,Site,Boson”)
1: 1 QN({“Nb”,0})
2: 1 QN({“Nb”,1})
Mismatched QN Index from set 2 (dim=2|id=979|“n=1,Site,Boson”)
1: 1 QN({“Nb”,0})
2: 1 QN({“Nb”,1})
From line 912, file itensor.cc
Mismatched QN Index arrows
Mismatched QN Index arrows
Do you know why this might be the case? Thanks for your comments!
I think the error is coming from the line:
auto PUp = ITensor(sj,prime(sj));
which should be changed to be
auto PUp = ITensor(dag(sj),prime(sj));
to make the arrows correct for the quantum-number-conserving case.
Also note how that code is fairly specific to the case of spin 1/2 degrees of freedom. A better approach to that code would be to write it just to collapse or sample in the “computational basis” and then make the code generic to any type of site that way. Then one can make an outer wrapper that applies single-site unitary ‘gates’ to change the basis for cases where one wants to sample in a different basis.
1 Like
Thank you Miles, it worked!
I modified this function such that it works for any d local hilbert space dimension, and sampling is done with a projection to the computational basis. So it works like a snapshot function for Boson site for example. Here’s the function together with its utility functions, in case someone else needs it:
#include <itensor/all.h>
#include "itensor/util/print_macro.h"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <cstdlib>
#include <random>
#include <algorithm>
#include <vector>
#include <time.h>
#include <string>
#include <cmath>
#include <math.h>
#include <complex>
#include <iomanip>
#include <filesystem>
#include <ctime>
#include <numeric>
using namespace itensor;
using namespace std;
// Define a struct to hold the value and its original index
struct IndexedValue {
float value;
int index;
};
// Comparator function to sort IndexedValue by the value
bool compareByValueDescending(const IndexedValue& a, const IndexedValue& b) {
return a.value > b.value; // Sort in decreasing order
}
void sortWithIndices(std::vector<float>& arr, std::vector<int>& sortedIndices) {
// Create a vector of IndexedValue to keep track of values and their indices
std::vector<IndexedValue> indexedArr(arr.size());
for (size_t i = 0; i < arr.size(); ++i) {
indexedArr[i] = {arr[i], static_cast<int>(i)};
}
// Sort indexedArr by value
std::sort(indexedArr.begin(), indexedArr.end(), compareByValueDescending);
// Extract the sorted indices and values
for (size_t i = 0; i < indexedArr.size(); ++i) {
arr[i] = indexedArr[i].value; // Update original array to sorted values
sortedIndices[i] = indexedArr[i].index; // Update indices
}
}
float sumSubset(const std::vector<float>& arr, const std::vector<int>& indices) {
float sum = 0.0f;
for (int index : indices) {
if (index >= 0 && index < arr.size()) {
sum += arr[index];
} else {
std::cerr << "Index " << index << " is out of range.\n";
}
}
return sum;
}
std::vector<int> generateArray(int minValue, int maxValue, int step) {
std::vector<int> array;
if (step <= 0) {
std::cerr << "Step size must be greater than 0.\n";
return array;
}
if (minValue > maxValue) {
std::cerr << "Minimum value cannot be greater than maximum value.\n";
return array;
}
for (int value = minValue; value <= maxValue; value += step) {
array.push_back(value);
}
return array;
}
vector<int>
collapse(MPS & psi, int dim, int progressPrint,
Args const& args = Args::global())
{
// dim: local hilbert space dimension
// progressPrint = 1: details of the computation to keep track.
auto sites = siteInds(psi);
auto N = length(psi);
auto cps = vector<int>(N+1);
auto probArray = vector<float>(dim);
std::vector<int> sortedIndices(probArray.size());
/*
std::random_device rd; // Obtain a seed from the random_device (non-deterministic)
std::mt19937 generator(rd()); // Initialize the Mersenne Twister engine with the seed
generator.seed(seedNo); // Reseed the generator with the new seed*/
// orthogonal form:
psi.position(1);
for(int j = 1; j <= N; ++j)
{
Index sj = sites(j);
for(int k = 1; k < dim; ++k)
{
auto PUp = ITensor(dag(sj),prime(sj));
PUp.set(k,k,1.0);
probArray[k-1] = elt(dag(prime(psi(j),"Site"))*PUp*psi(j));
if(progressPrint == 1) std::cout << probArray[k-1] << " ";
}
std::vector<int> indices = generateArray(0, dim-2, 1);
// Calculate the sum of the specified subset
float sumProb = sumSubset(probArray,indices);
probArray[dim-1] = 1 - sumProb;
if(progressPrint == 1) std::cout << probArray[dim-1] << " ";
// sort the array in decreasing order
sortWithIndices(probArray, sortedIndices);
if(progressPrint == 1)
{
// Output the results
std::cout << "Sorted array: ";
for (float value : probArray) {
std::cout << value << " ";
}
std::cout << "\n";
std::cout << "Original indices of sorted elements: ";
for (int index : sortedIndices) {
std::cout << index << " ";
}
std::cout << "\n";
}
float randomN = Global::random();
if(progressPrint == 1)
{
std::cout << "random number generated: ";
std::cout << randomN << "\n";
}
int st = sortedIndices[0]+1;
if(randomN > probArray[0]) st = sortedIndices[1]+1;
if(progressPrint == 1)
{
std::cout << "bit determined:\n";
std::cout << st << "\n";
}
cps.at(j) = st;
auto State = ITensor(sj);
for(int k = 1; k <= dim; ++k)
{
if(k == st) State.set(st,1.0);
}
//Project state
ITensor jstate = State;
if(j < N)
{
auto newA = psi(j+1)*(dag(jstate)*psi(j));
newA /= norm(newA);
psi.set(j+1,newA);
}
//Set site j tensor
psi.set(j,jstate);
}
return cps;
}
Best,
Ceren.
1 Like
This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.
Glad to hear it! Also thanks for posting your modified code – that should be very helpful for others to have.
1 Like