Fixes deprecated cluster_bingham.m and bingham_mlesac.m #3

madratman · 2015-07-02T07:21:34Z

Please review.
(I just browsed the deprecated folder today and realized I didn't have to write from scratch..)
Jared's paper with description of BMMs See pg 6, section E

madratman · 2015-07-02T07:40:17Z

Notes:

Look into bingham_mix_t as well. Define the data type.
More methods concerning mixtures:

~~bingham_mixture_sample~~
~~bingham_mixture_pdf~~
bingham_mixture_sample_ridge
bingham_mixture_add
bingham_mixture_mult
bingham_mixture_peak
bingham_mixture_thresh_peaks
bingham_mixture_thresh_weights
Strikethroughs are already there in matlab

madratman · 2015-07-07T21:00:41Z

This is a bit funny.
It doesn't perform as good as the one in C (based on synthetic data for the orientation of a cube in two face). Ideally, both the binghams would have concentration parameters like [-900. -900, ___] but only one of the distributions was "good" and fitted a surface.

Looking into this or worst case - wrap the C code.

EDIT: tries again. For once c and Matlab return pretty much the same output. Why is this happening?!

madratman · 2015-07-07T21:47:41Z

matlab/bingham_fit_mlesac.m

+   % update the threshold
+   if logp > pmax
+      pmax = logp;
+   end


Add condition for first iteration

first = 1 if (first || logarithm_p > pmax) pmax = logarithm_p; if first first = 0 else bing_test = bing_X_combined; end end

SebastianRiedel · 2015-07-15T15:25:57Z

Hey Ratnesh, two questions:

What is the current status on the comparison between the cluster fitting in C vs. Matlab using the same quaternion data? How big is the difference from your experience?
Could you provide me with a minimal example of Matlab Code to test the new functioniality (e.g. define two binghams, sample them & cluster the samples again). That way I can test your code before merging it in and add an example to the Matlab tutorial.

Thanks for your contribution, btw :).

madratman · 2015-07-20T16:57:55Z

Hey, I ll get back to this soon. I stopped using it as I think there's a bug in here.
I instead used the C code and wrote a script to parse the BMX files in matlab, because I needed the clusters asap and changed my priorities heh. But I completed bingham_mixture_sample.m and will add that in this PR once I am also done with cluster_bingham.m

What is weird is something I just found out an hour ago. (Maybe I should create an issue, but I want to affirm I am not making a blunder first)

I think there is a bug in bingham_sample, both in C function and Matlab code. You get repeated samples when you take more than a 100 at a time. (Jared has conditions for >100 samples)

It is okay for no_of_samples<100 because look here
which is the following snippet

void bingham_mixture_sample(double **X, bingham_mix_t *BM, int n)
{
  int i, j;

  if (n == 1) {
    i = pmfrand(BM->w, BM->n);
    bingham_sample(X, &BM->B[i], 1);
  }
  else if (n < 100) {
    for (i = 0; i < n; i++)
      bingham_mixture_sample(&X[i], BM, 1);
  }

If you just sample one at a time, you can escape this (That's what I'll do for n>100 from now on). (Same goes for bingham_mixture_sample.m )

To verify I ran ./bingham_sample in C (https://github.com/SebastianRiedel/bingham/blob/master/c/bingham_sample.c)
with Z = [-900 -900 0] and no of samples were 10. Here's what I got:

ratnesh:$ ./bingham_sample -900 -900 0 10
********* seed = 1437411427
-0.786728 -0.002284 -0.018014 -0.617033
-0.233186 0.000858 0.027903 -0.972031
-0.795369 0.012314 0.059786 -0.603045
-0.795369 0.012314 0.059786 -0.603045
-0.795369 0.012314 0.059786 -0.603045
0.119820 -0.009728 0.048072 -0.991583
-0.108677 0.013239 -0.014057 -0.993890
0.998119 -0.021786 0.030581 -0.048466
0.998119 -0.021786 0.030581 -0.048466
0.794630 -0.030609 -0.011896 0.606205

madratman · 2015-07-20T18:18:09Z

The problem is at L38-43 in bingham/matlab/bingham_sample_nd.m

    if a > rand()
        x = x2;
        p = p2;
        t = t2;
        num_accepts = num_accepts + 1;
    end

If a< rand() and i > burn_in+1,
the sample is repeated.

For the condition on burn_in+1 see #L53 (https://github.com/SebastianRiedel/bingham/blob/master/matlab/bingham_sample_nd.m#L53)

I don't know how to correct this rightaway. I ll have to look into the theory behind sampling from Binghams

Or maybe, we just need to discard if a<rand() which is not happening.

EDIT
So, it should be something on the lines of

num_accepts = 0;
X=[];
for i=1:n*sample_rate+burn_in
    x2 = X2(i,:)
    x2 = x2/norm(x2);
    t2 = T2(i); %bingham_pdf_unnormalized(x2,B);
    p2 = P2(i); %acgpdf_pcs(x2,z,V);
    a1 = t2 / t;
    a2 = p / p2;
    a = a1*a2;
    if a > rand()
        x = x2;
        p = p2;
        t = t2;
        num_accepts = num_accepts + 1;
        X = [X;x];  
   end
end

But it is not guaranteed if the condition a>rand() will be satisfied enough times so that we get back the same number of samples we passed in the function's argument. So, we increase the iterations of i to infinity and break the loop when the size of X is equal to the number of samples required.

P.S. I am just going by intuition here. Not looking at maths behind. :\

EDIT:
If you add an infinite loop..
You can't because:
Attempted to access X2(16,:); index out of bounds because size(X2)=[15,4].

So you need a lot of "proposals" - or whatever acgpdf_pcs.m returns - is called. So, I can hardcode them to say no_of_samples*1000 or something? This is so hackish

EDIT:
Hack works. :D. Hardcoded to 1000*no_of_samples.

madratman · 2015-07-20T19:34:09Z

One last bug which I forgot to correct. If you fit a BMM to some data, the last component is (generally) a uniform component, which in the C code has eigenvalues :

     0     0     0
     0     0     0
     0     0     0
     0     0     0

And the conc params will be `0 0 0' which is fine.

Though eigenvalues should be (Or could be any three pure quaternions representing orthogonal axis in 3d, right?)

     0     0     0
     1     0     0
     0     1     0
     0     0     1

This can also be seen in https://github.com/SebastianRiedel/bingham/blob/master/matlab/bingham_new_uniform.m
B.V = eye(d);
B.V = B.V(:,2:end);
I had looked at this earlier but forgot. Got to it as I got similar(I say similar, not same each time) samples for the last component. Say if 30 samples from 1000 in a BMM belong to the uniform component. They were like:

   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
    0.5237    0.5110    0.3740    0.5699
   -0.5237   -0.5110   -0.3740   -0.5699
    0.5237    0.5110    0.3740    0.5699

So these lines should be changed in bingham.c

/*
 * Create a new uniform Bingham distribution.
 */
void bingham_new_uniform(bingham_t *B, int d)
{
  bingham_alloc(B, d);
  memset(B->V[0], 0, d*(d-1)*sizeof(double));
  memset(B->Z, 0, (d-1)*sizeof(double));
  B->F = surface_area_sphere(d);
}

I don't have a complete understanding of why this happens. (similar samples)
Also, isn't such an eigenvector void. [0 0 0 0] corresponds to no unit quaternion.
Should there be an error/warning for this? (say if someone tries to come up with such a Bingham dist in their code)

Tags @SebastianRiedel @jglov3id for review.

madratman · 2015-07-20T21:23:03Z

An issue is now that the samples coming it this way are a bit noisy

SebastianRiedel · 2015-07-26T14:59:58Z

Hi Ratnesh, please use separate issues (in the issue section) for different topics otherwise we lose track of what needs to be fixed and what not.

1.) Regarding the "same samples returned" issue: Afaik there is no easy way to sample from the Bingham distribution directly, that is why in the c-function

void bingham_sample(double **X, bingham_t *B, int n)

Jared implemented a Metropolis-Hastings sampler (https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) using a multivariate Gaussian distribution as proposal distirbution for the samples. The fact that sometimes the same samples are returned is part of the algorithm and only this way the samples returned approximate the Bingham distribution.

2). "Regarding sampling from uniform binghams issue:" I can’t reproduce the uniform sampling problem, I played around with bingham_sample.c, printing the bingham distirbution before sampling using.

print_bingham(&B);

I got the following, qualitatively identical results.
With unit quaternion directional vectors:

~/software/bingham/c [master *]$ ./bingham_sample 0 0 0 10
B->F = 19.739210
B->Z = [ 0.000000 0.000000 0.000000 ]
B->V[0] = [ 0.000000 1.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 1.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 1.000000 ]
********* seed = 1437922090
0.571308 -0.236021 0.273719 0.736871
-0.342151 -0.371761 0.122001 0.854308
-0.456641 0.128070 0.874794 -0.099054
0.083509 -0.926933 0.077068 -0.357605
0.525685 -0.358047 0.764337 -0.106047
0.376756 0.212938 -0.855404 0.284598
0.416751 0.620587 0.596920 -0.291336
-0.873784 -0.004586 0.151851 0.461977
0.590044 -0.275812 0.436065 -0.620985
0.163787 -0.219224 -0.914115 0.299181

With all-zero directional vectors:

~/software/bingham/c [master *]$ ./bingham_sample 0 0 0 10
B->F = 19.739210
B->Z = [ 0.000000 0.000000 0.000000 ]
B->V[0] = [ 0.000000 0.000000 0.000000 0.000000 ]
B->V[1] = [ 0.000000 0.000000 0.000000 0.000000 ]
B->V[2] = [ 0.000000 0.000000 0.000000 0.000000 ]
********* seed = 1437921971
-0.719918 0.410494 0.376537 -0.414044
-0.414388 -0.837541 0.121301 0.334804
0.479978 0.128918 -0.237437 -0.834641
-0.208197 0.493726 0.074746 -0.841012
0.330121 0.191013 0.832717 -0.401393
0.104286 -0.082295 -0.031868 -0.990624
0.878819 0.473154 -0.060594 -0.011456
-0.877109 0.149065 0.339895 -0.304847
-0.988131 -0.057621 0.117940 0.079794
0.032479 -0.803551 -0.390414 -0.448138

Of course, all-zero directional vectors should not occur, but in this case (sampling from a uniform distribution) it doesn’t harm either.

malipoor · 2016-09-21T16:36:16Z

Hi there, trying to contact Sebastian, got confused and ended up writing here which seems not to be right place! If a newbie has some basic questions how is it handled here?! I found several definitions of Bingham distr. in literature and could not figure out how your version related to others. Also some said concentration parameters should be negative and some say positive!! Would be nice to reference a paper or tutorial with basics conforming to your notation. Thanks.

gerhardkurz · 2016-09-26T17:16:41Z

The library was originally developed by Jared Glover and is based on the conventions in the paper

Jared Glover, Leslie Pack Kaelbling, "Tracking 3-D rotations with the quaternion Bingham filter", MIT-CSAIL-TR-2013-005, 2013

It is available at
https://dspace.mit.edu/bitstream/handle/1721.1/78248/MIT-CSAIL-TR-2013-005.pdf

madratman added 2 commits July 2, 2015 03:18

fixed deprecated bingham_fit_mlesac.m

dc8304d

fixed deprecated bingham_cluster.m

2214522

madratman added 4 commits July 2, 2015 04:15

adds comments, fixes upper bound back to 100

29938dd

comments and cleanup for mlesac.m

5b0ed1d

changes directory to undeprecate

3d5eff6

accomodates previous PR, minor changes

b29e566

madratman reviewed Jul 7, 2015
View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixes deprecated cluster_bingham.m and bingham_mlesac.m #3

Fixes deprecated cluster_bingham.m and bingham_mlesac.m #3

madratman commented Jul 2, 2015

madratman commented Jul 2, 2015

madratman commented Jul 7, 2015

madratman Jul 7, 2015

SebastianRiedel commented Jul 15, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

SebastianRiedel commented Jul 26, 2015

malipoor commented Sep 21, 2016

gerhardkurz commented Sep 26, 2016

Fixes deprecated cluster_bingham.m and bingham_mlesac.m #3

Are you sure you want to change the base?

Fixes deprecated cluster_bingham.m and bingham_mlesac.m #3

Conversation

madratman commented Jul 2, 2015

madratman commented Jul 2, 2015

madratman commented Jul 7, 2015

madratman Jul 7, 2015

Choose a reason for hiding this comment

SebastianRiedel commented Jul 15, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

madratman commented Jul 20, 2015

SebastianRiedel commented Jul 26, 2015

malipoor commented Sep 21, 2016

gerhardkurz commented Sep 26, 2016