modsim / hopsy

Python interface to C++ sampling library HOPS
MIT License
11 stars 3 forks source link

Functions on the latest versions of hopsy #4

Open vlsena opened 1 year ago

vlsena commented 1 year ago

Hi, I'm Physics master student in Brazil and I have been using hopsy library for a few months in my research project. Basically I am trying to sort out evenly distributed points within a polytope. I used your 'test.py' file as a base to build my code and it worked fine for a while, but when I recently updated hopsy library on my computer, some of the functions like "MultivariateGaussianModel" and "Run" stopped working. Have these functions had their names changed? How do I use them exactly the same as before?

I'm attaching my code (that's not working anymore) bellow:

import hopsy
import numpy as np

def sorteio(A,b,mu,cov,n_points): 
    model = hopsy.MultivariateGaussianModel(mu, cov)

    problem = hopsy.Problem(A, b, model)

    # the run object contains and constructs the markov chains. in the default case, the
    # Run object will have a single chain using the Hit-and-Run proposal algorithm and is
    # set to produce 10,000 samples.
    run = hopsy.Run(problem)

    # we finally sample
    run.sample(n_points)

    # from the run, we can now extract the produced data
    data = run.data

    # the states is a list of lists of numpy.ndarrays, which can be casted to a numpy.ndarray
    # which then has the shape (m,n,d), where m is the number of chains, n the number of samples
    # and d the dimenion
    states  = data.states

    return np.array(states)
ripaul commented 1 year ago

Hi there! Yes, we have changed the API. The hopsy.MultivariateGaussianModel has become simply a hopsy.Gaussian and the hopsy.Run has been completely discarded. Instead you now have to setup the Markov chains yourself as

chains = [hopsy.MarkovChain(problem, starting_point=some_valid_starting_point) for i in range(n_chains)]
rngs = [hopsy.RandomNumberGenerator(seed, i) for i in range(n_chains)]

accrate, states = hopsy.sample(chains, rngs, n_samples=10_000)

Here you still have to set the starting point some_valid_starting_point, the number of parallel Markov chains n_chains you want to use and your random seed seed. If you don't have a starting point at hand, you can try to compute the Chebyshev center of the polytope and use this one as a starting point using hopsy.compute_chebyshev_center(problem). It is always recommended to use multiple parallel chains for sampling in order to estimate the convergence of your samplers, I think arviz recommends like 4? In any case, if you use parallel chains, you can also parallelize them by calling hopsy.sample(..., n_threads=64) or whatever your computer supports.

Regarding what you're trying to achieve: If you want to sample a uniform distribution on a convex polytope, then you not use the hopsy.Gaussian class at all. Instead, just set up a problem with hopsy.Problem(A, b). In hopsy, no density means uniform density and, thus, distribution.

I hope that helps! Feel free to ask if you have any more questions, I'm glad to help. Also checkout the docs. They're still pretty much WIP, but maybe they can still help you.

qacwnfq commented 1 year ago

Hey :) The changes in the api allowed us to incoroporate some cool stuff in the backend. While we strongly recommend upgrading, in case you are very limited in time you could try pinning hopsy to the old version.

Also out of curiousity: What problem in physics are you working on? Is it related to astronomy? I've only ever seen convex polytopes in the context of free form gravitational lense reconstruction.

vlsena commented 1 year ago

Hey guys, thanks a lot for the support! My code is already functional again thanks to you.

Richard, I am indeed trying to sample a uniform distribution on a convex polytope, so I must use only hopsy.Problem(A, b) , as you suggested, right? I am very grateful for the clarifications.

Johann, the problem I'm dealing with is a problem on the foundations of quantum theory and it is not related to astronomy. If you are familiar with the subject of Quantum Mechanics, you may be familiar with "Bell scenarios"; what happens is that, for a given Bell scenario, if we classify all probability distributions as vector (or points) coordinates, the union of these points form certain geometries in a hyperspace. If these probabilities represent a "local behaviour" of the system, the points are restricted to a polytope defined by the famous Bell inequalities; if the probabilities represent a "quantum behavior", the points form a larger convex set that contains the local polytope. This is the basics about the topic I'm researching, but we are trying to do more interesting things with these geometries by sorting some points inside it hehe.

Best, Vitor

On Mon, Aug 1, 2022 at 7:57 PM Johann Fredrik Jadebeck < @.***> wrote:

Hey :) The changes in the api allowed us to incoroporate some cool stuff in the backend. While we strongly recommend upgrading, in case you are very limited in time you could try pinning hopsy to the old version.

Also out of curiousity: What problem in physics are you working on? Is it related to astronomy? I've only ever seen convex polytopes in the context of free form gravitational lense reconstruction.

— Reply to this email directly, view it on GitHub https://github.com/modsim/hopsy/issues/4#issuecomment-1201815842, or unsubscribe https://github.com/notifications/unsubscribe-auth/AXO3PAVO6RKJCC4EPOMOU3LVXBI6LANCNFSM54XHFXOA . You are receiving this because you authored the thread.Message ID: @.***>