Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds working bitstring Sample function to MPS. #490

Merged
merged 6 commits into from
Dec 22, 2021
Merged

Adds working bitstring Sample function to MPS. #490

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f3d63cb
Adds working bitstring Sample function to MPS.
MichaelBroughton Dec 17, 2021
dd0f7f1
remove iostream include.
MichaelBroughton Dec 17, 2021
fa36fd7
remove use of gmock.
MichaelBroughton Dec 17, 2021
70f7a2b
Moved to O(n) sampling algorithm.
MichaelBroughton Dec 21, 2021
e4ef583
remove unused includes.
MichaelBroughton Dec 21, 2021
fc1fcc8
Merge branch 'master' into mps_sample
MichaelBroughton Dec 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 63 additions & 0 deletions lib/mps_statespace.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include <cstdlib>
#include <cstring>
#include <memory>
#include <random>

#include "../eigen/Eigen/Dense"
#include "../eigen/unsupported/Eigen/CXX11/Tensor"
Expand Down Expand Up @@ -372,6 +373,68 @@ class MPSStateSpace {
out = t_4d.contract(t_2d, product_dims);
}

// Draw a single bitstring sample from state using scratch and scratch2
// as working space.
static void SampleOnce(MPS& state, MPS& scratch, MPS& scratch2,
std::mt19937* random_gen, std::vector<bool>* sample) {
const auto bond_dim = state.bond_dim();
const auto num_qubits = state.num_qubits();
std::default_random_engine generator;
fp_type* scratch_raw = scratch.get();
fp_type rdm[8];

sample->reserve(num_qubits);
Copy(state, scratch);
Copy(state, scratch2);

// Sample left block.
ReduceDensityMatrix(scratch, scratch2, 0, rdm);
auto p0 = rdm[0] / (rdm[0] + rdm[6]);
std::bernoulli_distribution distribution(1 - p0);
auto bit_val = distribution(*random_gen);

sample->push_back(bit_val);
MatrixMap tensor_block((Complex*)scratch_raw, 2, bond_dim);
tensor_block.row(!bit_val).setZero();
tensor_block.imag() *= -1;

// Sample internal blocks.
for (unsigned i = 1; i < num_qubits - 1; i++) {
ReduceDensityMatrix(scratch, scratch2, i, rdm);
p0 = rdm[0] / (rdm[0] + rdm[6]);
distribution = std::bernoulli_distribution(1 - p0);
bit_val = distribution(*random_gen);

sample->push_back(bit_val);
const auto mem_start = GetBlockOffset(scratch, i);
new (&tensor_block) MatrixMap((Complex*)(scratch_raw + mem_start),
bond_dim * 2, bond_dim);
for (unsigned j = !bit_val; j < 2 * bond_dim; j += 2) {
tensor_block.row(j).setZero();
}
tensor_block.imag() *= -1;
}

// Sample right block.
ReduceDensityMatrix(scratch, scratch2, num_qubits - 1, rdm);
p0 = rdm[0] / (rdm[0] + rdm[6]);
distribution = std::bernoulli_distribution(1 - p0);
bit_val = distribution(*random_gen);
sample->push_back(bit_val);
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @MichaelBroughton for adding this! If I understand the implementation correctly, you can save time from n^2 to n where n is the number of sites by storing intermediate contractions. That is, first compute \rho_{1,...,n-1} (rdm obtained by tracing out the last qubit) and store it in memory, then from this compute \rho_{1,...,n-2} and store it in memory, then from this ... compute \rho_{1}. Then sample bits from 1 to n, again reusing these already-computed rdms.

I hope that makes sense... I can't find a reference implementation right now. Basically the idea is not to contract edges more than once - e.g., currently to get \rho_{1} and \rho_{2}, the last n - 2 edges of the MPS are individually contracted twice, but you can avoid this.

Copy link
Collaborator Author

@MichaelBroughton MichaelBroughton Dec 18, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @rmlarose I'm not sure I follow. Can you write out the approach here quickly using something like python and einsum ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, looking back I think you may already be doing what I suggested since scratch and scratch2 get modified in each call to ReduceDensityMatrix.

For an n=3 MPS, the sequence of rdms should be:

|  |  |
@--@--@  (rho_123)
|  |  |

|  |
@--@     (rho_12 obtained from rho_123)
|  |

|  
@        (rho_1 obtained from rho_12, *not* from rho_123)
|

That is, for each i you want to obtain rho_1...i from rho_1...(i + 1), not from rho_1...n. If you are already doing this (it's not 100% clear to me how scratch and scratch2 get modified in ReduceDensityMatrix), you can ignore my comment.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not entirely sure I see how this contraction process is O(n) and not O(n^2). The contraction process in this PR is this:

  1. Computes RDM on the first qubit.
  2. Toss a coin with probabilities set by the RDM.
  3. Projects the qubit to zero or one state depening on coin outcome.
  4. Repeat this process for qubits [1,....n].

Step 1 in this operation always considers tensors [0,...n]. You compute n RDM matrices at the cost of n contractions each so it's O(n^2) to draw a single bitstring. Looking at your proposal IIUC you do this:

  1. Compute RDM on the first qubit.
  2. Toss a coin with probability set by the RDM.
  3. Project the qubit to zero or one state depending on the coin outcome.
  4. As I move to contracting the next qubit after the one I just did I re-use some of the contractions I've done up until this point (on the left side only) to compute the RDM to prevent having to contract all the way from the first qubit every time.

Suppose I have this MPS:

  │     │     │     │     │
┌─┴─┐ ┌─┴─┐ ┌─┴─┐ ┌─┴─┐ ┌─┴─┐
│ 0 ├─┤ 1 ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└───┘ └───┘ └───┘ └───┘ └───┘

and part way through my sampling process I need to compute the RDM on qubit 3 (having already sampled and projected down qubits 0,1 and 2):

                    │
┌───┐ ┌───┐ ┌───┐ ┌─┴─┐ ┌───┐
│ 0 ├─┤ 1 ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└─┬─┘ └─┬─┘ └─┬─┘ └───┘ └─┬─┘
  │     │     │           │
┌─┴─┐ ┌─┴─┐ ┌─┴─┐ ┌───┐ ┌─┴─┐
│ 0 ├─┤ 1 ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└───┘ └───┘ └───┘ └─┬─┘ └───┘
                    │

With my approach I would just naively recontract everything together. IIUC your approach would keep another temporary variable containing the contractions of the projections from tensor 0 up to i - 1 if I want to compute the RDM for qubit i. This way I can use the fact that I have previously computed:

              │
┌───┐ ┌───┐ ┌─┴─┐ ┌───┐ ┌───┐
│ 0 ├─┤ 1 ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└─┬─┘ └─┬─┘ └───┘ └─┬─┘ └─┬─┘
  │     │           │     │
┌─┴─┐ ┌─┴─┐ ┌───┐ ┌─┴─┐ ┌─┴─┐
│ 0 ├─┤ 1 ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└───┘ └───┘ └─┬─┘ └───┘ └───┘
              │

Where I could store this:

        │
┌───┐ ┌─┴─┐ ┌───┐ ┌───┐
│   ├─┤ 2 ├─┤ 3 ├─┤ 4 │
│   │ └───┘ └─┬─┘ └─┬─┘
│0,1│         │     │
│   │ ┌───┐ ┌─┴─┐ ┌─┴─┐
│   ├─┤ 2 ├─┤ 3 ├─┤ 4 │
└───┘ └─┬─┘ └───┘ └───┘
        │

And then (from the diagram above) I can contract 2 into 0,1 -> 0,1,2 and have the entire left side of my 3 RDM calculation done by simply remembering some of the computations I did from the previous step. This means on iteration number:

  1. I want to compute RDM-0 I contract n-1 tensors on the right of tensor 0 and none on the left.
  2. I want to compute RDM-1 I contract n-2 tensors on the right of tensor 1 and just one on the left (creating my history variable)
  3. I want to compute RDM-2 I contract n-3 tensors on the right of tensor 2 and just one on the left (applying one contraction onto my history variable)
  4. I want to compute RDM-3 I contract n-4 tensors on the right of tensor 3 and just one on the left (applying one contraction to my history variable)
    .
    .
    .

There will be n steps where step i requires n-i+1 contractions -> \sum_i^n (n-i) -> n * (n - 1) / 2 -> O(n^2) with a prefactor of 1/2 instead of a prefactor of 1. How do I get this down to O(n) ? You can't "cache" tensor contractions on the right side of your RDM index as you go using this same approach can you ?

If we can get it down to O(n) I'd say it's worth updating in this PR, but if it's just the prefactor of 2 I'd say that's fine to come in a subsequent PR.

Copy link
Contributor

@rmlarose rmlarose Dec 20, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @MichaelBroughton for the very detailed and clear response. It seems my original comment does apply.

With my approach I would just naively recontract everything together. IIUC your approach would keep another temporary variable containing the contractions of the projections from tensor 0 up to i - 1 if I want to compute the RDM for qubit i.

Indeed this is precisely the idea. There is one important change to your sequence of diagrams to get to the O(n) algorithm. Below I sketch out the full O(n) sampling algorithm, which proceeds in two "phases".

Phase 1: Computing RDMs

Starting from the first step with an n=4 MPS, you have this tensor network:

|  |  |  |
@--@--@--@
|  |  |  |

First, you store this. (Note this is not a single-site RDM, which I think is the key difference from how you are thinking about this.) Then you connect the rightmost free legs and contract every edge on the last tensor to get the RDM on the first three sites:

|  |  | 
@--@--@
|  |  |

Store this. Then connect the rightmost free legs and contract every edge on the last tensor to get the RDM on the first two sites:

|  | 
@--@
|  |

Store this. Then connect...to get the RDM on the first site:

|
@
|

Store this. Note: There are n steps here, and each step requires contracting O(1) edges, yielding O(n) total contracted edges. The output of this step is the list of n "sub tensor networks" shown above.

Phase 2: Sampling

To sample, you iterate backwards through this list. Starting with the last "sub tensor network", namely the RDM on the first site

|
@
|

sample a bit a.

Next, move to the previous "sub tensor network", namely the RDM on the first two sites, and attach a (as a vector, |0> or |1>, corresponding to the bit) to the first site:

a
|  | 
@--@
|  | 
a

Contract this tensor network to get a 2x2 rdm, then sample a bit b.

Next, move to the previous "sub tensor network", namely the RDM on the first three sites, and attach a and b (as vectors) to the first two sites:

a  b
|  |  | 
@--@--@
|  |  |  
a  b

Contract this tensor network to get a 2x2 rdm, then sample a bit c. Note: In the previous step, the tensor network

a
|
@--
| 
a

was already contracted. To make "Phase 2" also have O(n) contractions, you need to store and reuse this result.

Finally, repeat the same process to get the last bit. Namely, form

a  b  c
|  |  |  |
@--@--@--@
|  |  |  | 
a  b  c

contract to get a 2x2 rdm, then sample a bit d. Note again that the tensor network

a  b
|  |  
@--@--
|  |  
a  b

appearing in this diagram has already been contracted in the previous step.

Output abcd. Note: This phase has n steps, each of which require O(1) contractions, yielding O(n) contractions. The total for both phase 1 and phase 2 then is O(n).

Pseudocode

If helpful, I sketched this procedure out in the following rough-around-the-edges pseudocode.

def sample(mps) -> List[int]:
  """Input mps is |psi><psi|:
  |  |  |  |  (top_free_legs)
  @--@--@--@  (tensors)
  |  |  |  |  (bottom_free_legs)
  """
  # Phase 1: Compute reduced density matrices rho_1, ..., rho_1...n.
  n = number of qubits in mps
  rdms = [mps]  # List of tensor networks corresponding to each rdm
  for i in range(n - 1, 0, -1):
    mps = connect top_free_legs[i] and bottom_free_legs[i] then contract all legs on tensors[i]
    rdms.append(mps)
    
  # Phase 2: Sample bits.
  bits: List[int] = []
  for i in range(n):
    rdm = rdms[-i - 1]
    bit = sample a bit from rdm, a 2x2 density matrix up to normalization
    bits.append(bit)

    left_tensor_network = add bit to free legs of rdm and contract
    replace i + 1 leftmost tensors in rdms[-i - 2] with left_tensor_network
  return bits

Copy link
Collaborator Author

@MichaelBroughton MichaelBroughton Dec 20, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. This makes sense. One final question: I'm still not sold on the Phase 1 portion. Looking at the pseudocode here:

rdms = [mps]  # List of tensor networks corresponding to each rdm
for i in range(n - 1, 0, -1):
  mps = connect top_free_legs[i] and bottom_free_legs[i] then contract all legs on tensors[i]
  rdms.append(mps)

rdms will contain initially an n site tensor network. After the first loop iteration it will contain an n site network and an n-1 site network etc. all the way down to a 1 site network. In terms of space requirements it feels like \sum_i^n i -> n*(n-1) / 2 -> O(n^2) has shown up again, only now in space and not time. Am I still missing something ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're absolutely right that phase 1, as I've described it, requires O(n^2) memory. However, after talking with a "famous" mps guy - who has formally kicked me out of the tensor network community for describing it this way 😄 - there's a better way to do phase 1 which only requires O(n) memory:

You still store the copy of the first tensor network (using the same n=4 example):

|  |  |  |
u--v--t--w
|  |  |  |

But now you just compute and store the following tensors by successively connecting top/bottom free indices and contracting:

--w

--(tw)

--(vtw)

Explicitly, as an example, --(tw) is

 i|  |j
--t--w
 i|  |j

(i and j indicate these legs are connected, respectively).

The time is still O(n), but now the space is also O(n).

The only difference to phase 2 is how the RDMs are obtained. For example, in the first step you form

|
u--(vtw)
|

and contract to get the 2x2 density matrix, then sample a bit a. Note that

|
u
|

is selected from the copy of the tensor network in the first step of phase 1.

Next you form

a
|  |
u--v--(tw)
|  |
a

and contract, sample, etc.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright then. I've gone ahead and implemented this. It's a lot faster on the big systems which is great. I didn't look too carefully into profiling it for things like cache friendliness, but generally speaking it looks like it's spending a lot of time doing BLAS operations and moving data:
image

This might be a good follow up for anyone interested in trying to squeeze some more performance out.


// Draw num_samples bitstring samples from state and store the result
// bit vectors in results. Uses scratch and scratch2 as workspace.
static void Sample(MPS& state, MPS& scratch, MPS& scratch2,
unsigned num_samples, unsigned seed,
std::vector<std::vector<bool>>* results) {
std::mt19937 rand_source(seed);
results->reserve(num_samples);
for (unsigned i = 0; i < num_samples; i++) {
MichaelBroughton marked this conversation as resolved.
Show resolved Hide resolved
SampleOnce(state, scratch, scratch2, &rand_source, &(*results)[i]);
}
}

// Testing only. Convert the MPS to a wavefunction under "normal" ordering.
// Requires: wf be allocated beforehand with bond_dim * 2 ^ num_qubits -1
// memory.
Expand Down
204 changes: 204 additions & 0 deletions tests/mps_statespace_test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@

#include "../lib/formux.h"
#include "gtest/gtest.h"
#include "gmock/gmock.h"

namespace qsim {

Expand Down Expand Up @@ -900,6 +901,209 @@ TEST(MPSStateSpaceTest, ReduceDensityMatrixLarge){

}

TEST(MPSStateSpaceTest, SampleOnceSimple){
95-martin-orion marked this conversation as resolved.
Show resolved Hide resolved
auto ss = MPSStateSpace<For, float>(1);
auto mps = ss.Create(3, 4);
auto scratch = ss.Create(3, 4);
auto scratch2 = ss.Create(3, 4);
std::mt19937 rand_source(1234);
std::vector<bool> results;

// Set to |100>.
results.clear();
ss.SetStateZero(mps);
mps.get()[0] = 0;
mps.get()[8] = 1;
ss.SampleOnce(mps, scratch, scratch2, &rand_source, &results);
ASSERT_THAT(results, testing::ElementsAre(1, 0, 0));
//EXPECT_EQ(1,2);

// Set to |010>.
results.clear();
ss.SetStateZero(mps);
mps.get()[16] = 0;
mps.get()[24] = 1;
ss.SampleOnce(mps, scratch, scratch2, &rand_source, &results);
ASSERT_THAT(results, testing::ElementsAre(0, 1, 0));

// Set to |001>.
results.clear();
ss.SetStateZero(mps);
mps.get()[80] = 0;
mps.get()[82] = 1;
ss.SampleOnce(mps, scratch, scratch2, &rand_source, &results);
ASSERT_THAT(results, testing::ElementsAre(0, 0, 1));

// Set to |101>.
results.clear();
ss.SetStateZero(mps);
mps.get()[0] = 0;
mps.get()[8] = 1;
mps.get()[80] = 0;
mps.get()[82] = 1;
ss.SampleOnce(mps, scratch, scratch2, &rand_source, &results);
ASSERT_THAT(results, testing::ElementsAre(1, 0, 1));
}

TEST(MPSStateSpaceTest, SampleGHZ){
const int num_samples = 10000;
auto ss = MPSStateSpace<For, float>(1);
auto mps = ss.Create(3, 4);
auto scratch = ss.Create(3, 4);
auto scratch2 = ss.Create(3, 4);
std::vector<std::vector<bool>> results(
num_samples, std::vector<bool>({}));

memset(mps.get(), 0, ss.RawSize(mps));
mps.get()[0] = 1;
mps.get()[10] = 1;
mps.get()[16] = 1;
mps.get()[42] = -1;
mps.get()[80] = 0.70710677;
mps.get()[86] = -0.70710677;

float count = 0;
ss.Sample(mps, scratch, scratch2, num_samples, 1234, &results);
for(int i = 0 ; i < num_samples; i++){
ASSERT_THAT(results[i], testing::AnyOf(testing::ElementsAre(1, 1, 1),
testing::ElementsAre(0, 0, 0)));
count += results[i][0];
}
EXPECT_NEAR(count / float(num_samples), 0.5, 1e-2);
}

TEST(MPSStateSpaceTest, SampleComplex){
const int num_samples = 10000;
auto ss = MPSStateSpace<For, float>(1);
auto mps = ss.Create(4, 4);
auto scratch = ss.Create(4, 4);
auto scratch2 = ss.Create(4, 4);
std::vector<std::vector<bool>> results(
num_samples, std::vector<bool>({}));

memset(mps.get(), 0, ss.RawSize(mps));
mps.get()[ 0 ] = 0.033688569334715854 ;
mps.get()[ 1 ] = -0.10444182602180123 ;
mps.get()[ 2 ] = 0.9076354671683359 ;
mps.get()[ 3 ] = 0.405160344657187 ;
mps.get()[ 8 ] = -0.9595253512026178 ;
mps.get()[ 9 ] = -0.25936097827312377 ;
mps.get()[ 10 ] = -0.03987001675676861 ;
mps.get()[ 11 ] = 0.10224185693597321 ;
mps.get()[ 16 ] = -0.4350591822776815 ;
mps.get()[ 17 ] = 0.22228546667942578 ;
mps.get()[ 18 ] = -0.6285732819602607 ;
mps.get()[ 19 ] = 0.5943422063507785 ;
mps.get()[ 20 ] = -0.02428345908884816 ;
mps.get()[ 21 ] = 0.026256572727652475 ;
mps.get()[ 22 ] = 0.0728063572325396 ;
mps.get()[ 23 ] = -0.07991114142962712 ;
mps.get()[ 24 ] = -0.1642035571020447 ;
mps.get()[ 25 ] = -0.8209212529030018 ;
mps.get()[ 26 ] = 0.21124207331921135 ;
mps.get()[ 27 ] = 0.4033152234452636 ;
mps.get()[ 28 ] = -0.11315780332634073 ;
mps.get()[ 29 ] = -0.18477947087021204 ;
mps.get()[ 30 ] = -0.11199707215175961 ;
mps.get()[ 31 ] = -0.17985444082650426 ;
mps.get()[ 32 ] = -0.18059162771674087 ;
mps.get()[ 33 ] = -0.12173196101857839 ;
mps.get()[ 34 ] = 0.19817171168239098 ;
mps.get()[ 35 ] = 0.063054719070231 ;
mps.get()[ 36 ] = 0.45024015745008505 ;
mps.get()[ 37 ] = 0.11068157212593255 ;
mps.get()[ 38 ] = 0.8106501683288581 ;
mps.get()[ 39 ] = 0.19287240226762353 ;
mps.get()[ 40 ] = -0.08702392413741208 ;
mps.get()[ 41 ] = 0.07370887245848737 ;
mps.get()[ 42 ] = -0.018412987786278347 ;
mps.get()[ 43 ] = 0.027921764369775018 ;
mps.get()[ 44 ] = 0.42351439662329743 ;
mps.get()[ 45 ] = -0.7466201917698305 ;
mps.get()[ 46 ] = -0.24298735917837008 ;
mps.get()[ 47 ] = 0.4359199055764641 ;
mps.get()[ 80 ] = 0.422436255430577 ;
mps.get()[ 81 ] = 0.0 ;
mps.get()[ 82 ] = 0.1211402132186689 ;
mps.get()[ 83 ] = -0.819174648113452 ;
mps.get()[ 84 ] = 0.0 ;
mps.get()[ 85 ] = -7.333691512826885e-20 ;
mps.get()[ 88 ] = -0.8676720638499252 ;
mps.get()[ 89 ] = 0.1360568551008419 ;
mps.get()[ 90 ] = 0.011793867549118184 ;
mps.get()[ 91 ] = -0.44097834083157716 ;
mps.get()[ 96 ] = -0.08978879818674973 ;
mps.get()[ 97 ] = 0.0 ;
mps.get()[ 98 ] = -0.021807957717091198 ;
mps.get()[ 99 ] = -0.05873893136151775 ;
mps.get()[ 100 ] = 0.0 ;
mps.get()[ 101 ] = -9.89518074266081e-19 ;
mps.get()[ 104 ] = 0.14307979940517454 ;
mps.get()[ 105 ] = 0.06032194765563529 ;
mps.get()[ 106 ] = 0.22589440044405648 ;
mps.get()[ 107 ] = -0.2397609424987549 ;
mps.get()[ 112 ] = 0.12734430206944722 ;
mps.get()[ 113 ] = 0.0 ;
mps.get()[ 114 ] = -0.003114595079760157 ;
mps.get()[ 115 ] = 0.06816683893204967 ;
mps.get()[ 116 ] = 2.722010100011512e-17 ;
mps.get()[ 117 ] = 1.9629880528929172e-18 ;
mps.get()[ 120 ] = 0.014022255715263434 ;
mps.get()[ 121 ] = 0.017127855001478075 ;
mps.get()[ 122 ] = 0.025812082320798548 ;
mps.get()[ 123 ] = -0.027110021000464 ;
mps.get()[ 128 ] = -0.018262196707018574 ;
mps.get()[ 129 ] = 0.0 ;
mps.get()[ 130 ] = 0.0032725358458428836 ;
mps.get()[ 131 ] = 0.0310845568816579 ;
mps.get()[ 132 ] = 1.935877805637811e-17 ;
mps.get()[ 133 ] = -1.0370773958773989e-18 ;
mps.get()[ 136 ] = -0.028632305212090994 ;
mps.get()[ 137 ] = 0.012199896816087576 ;
mps.get()[ 138 ] = 0.0009323445588941451 ;
mps.get()[ 139 ] = -0.014212789540748644 ;
mps.get()[ 144 ] = -0.07762944756130831 ;
mps.get()[ 145 ] = -0.25063255485414937 ;
mps.get()[ 146 ] = 0.515385895406013 ;
mps.get()[ 147 ] = 0.7314486807404007 ;
mps.get()[ 148 ] = 0.20689214104052 ;
mps.get()[ 149 ] = 0.2781707321332216 ;
mps.get()[ 150 ] = 0.08286244916183945 ;
mps.get()[ 151 ] = 0.05888783848647657 ;

ss.Sample(mps, scratch, scratch2, num_samples, 12345, &results);

std::vector<float> expected({
0.00467637,
0.0020386,
0.00112952,
0.0269848,
0.00704221,
0.00147802,
0.00243688,
0.0350753,
0.0324814,
0.0141599,
0.0412371,
0.0780275,
0.0644363,
0.140995,
0.0355866,
0.512215,
});
std::vector<float> hist(16, 0);
for(int i =0;i<num_samples;i++){
int index = 0;
for(int j =0;j<4;j++){
index |= (results[i][j]) << (4 - j - 1);
}
hist[index] += 1;
}
for(int i =0;i<16;i++){
EXPECT_NEAR(hist[i] / float(num_samples), expected[i], 1e-2);
}
}

} // namespace
} // namespace mps
} // namespace qsim
Expand Down