New sampling mechanisms (Sobol, Giustino)

[diegomartinez2]

Getting ready the new branch "Diegom_sobol" for CellConstructor with the Sobol sampling mechanism.

[diegomartinez2]

Geting ready also the new branch "Diegom_sobol" in python_sscha for calling the Sobol sampling from CellConstructor "Diegom_sobol" branch. Also trying to maintain retro-compatibility, so the sobol flag is False by default.

[diegomartinez2]

Testing "Diegom_sobol" branches (for python-ssha and CellConstructor). It doesn't look like the Sobol sampling improves the calculations. As I implemented it, the Sobol function always gives the same sequence. Maybe it will need improvements or some random "salt".

[mesonepigreco]

Hi Diego, Thanks for your efforts. Indeed, it is quite difficult to estimate how much this can improve. I would say we can observe several aspects: 1) In a long minimization, how many populations are required to converge with the same number of configurations? 2) When we are close to convergence, what is the minimum number of configurations required to have the gradient to go to zero without the ensemble going outside the statistical criterion?

If you can try to plot somethibg like these two values it would be important to understand if sobol improves, does nothing or worsen the convergence.

I would use the SnTe force field, as it is superfast to compute to make these tests, but if you have better ideas, please, go on!

[diegomartinez2]

Found a error in the way I introduced the Sobol sequence. I was following the random generation, but in SSCHA the random generator is later divided between the modes. This is not a problem with random numbers but is a problem with Sobol sequence. I'm changing the code to create one Sobol sequence per mode. So in the end is one Sobol sequence per mode per configuration.

[diegomartinez2]

There is another issue with Sobol. In this case is about transforming the Sobol sequence into a gaussian shape. In the process of transforming the sequence into gaussian, only the points inside the unit radius circle are taken so there are points of the Sobol sequence that are eliminated. This is not a problem due to the fact that an uniform set of points in a square are still uniform if we cut the corners, but the number of points at the end will be less and the idea of making the points powers of 2 gets more complicated. In other words we can generate 2**n points and cut the corners but at the end there will be less points and in extreme cases of low numbers there may be almost nothing at the end. I'll point out that it was not a problem for the random generator because if a point gets outside the unit circle you just cast another random point until you have the number of points you want.

[mesonepigreco]

Hi Diego. I think it is a problem eliminating points. With the SOBOL we do not want to get a random sequence, we want the "best" possible sequence of points. If you eliminate points, I think this will spoil the advantage of using Sobol.

In any case, why do you have to eliminate points? You can always transform random points on a square into 2D Gaussian distributed random variables without discarding anything

[diegomartinez2]

I'm using this code for the gaussian it uses two dimension (originally a couple of random numbers) to calculate one dimension gaussian (modified from from random_gaussian(rand) from python-sscha/SSCHAmodules/module_stochastic.f90, which comes from "Numerical Recipes"):

         v1=2.0*sample[i][0]-1.0
         v2=2.0*sample[i][1]-1.0
         Riq=v1**2+v2**2
         if(0 < Riq <= 1):
             data3=np.sqrt(-2*np.log(Riq)/Riq)
             data1.append(v1*data3)

[mesonepigreco]

No, this way is super inefficient:

Simply use the standard algorithm for gaussian numbers:

# u1 and u2 are uniform in [0,1)
r = -np.sqrt(-2 * np.log(u1))
theta = 2*np.pi * u2
random1 = r * np.cos(theta)
random2 = r*np.sin(theta)

So you get 2 random numbers for each pair of Sobol numbers. Then it is ok since we will generate powers of 2, the number of Sobol numbers are even and you do not discard anything.

[diegomartinez2]

The first Sobol number is zero, so the log is not defined. The new method requires half the Sobol numbers and looks different. Still it may be better than random generation for small configurations.

[mesonepigreco]

You can shift the origin of the uniform distribution by adding a random number between 0, 1 and then do modulus 1 (you should apply this operation to the whole Sobol sequence with the same random number so that the spacing remains fixed, is only the origin that shifts) . In this way, you will have no 0 values in the sequence.

I think this could be of some interest: https://www.sciencedirect.com/science/article/pii/S0895717710005935

Essentially it describes the best way to transform a low-discrepancy uniform sequence into a normal distribution.

The algorithm I proposed here above is the Box-Muller, the alternative is to use the direct inversion of the error function, which is available in scipy.stat (but also in that case you will have a problem if your first number is 0).

[diegomartinez2]

Yes I see the problem with zero in the direct inversion and also with the Beasley-Springer-Moro method. This method keeps the Sobol sequence properties. From this data (512 points): to this: Not so visible the gaussian but still works with 32 points: