jeanluct/braidlab

Matlab package for analyzing data using braids

braidlab's Issues

one-strand braid of data has serious problems!

Erick Fredj found a bug when he ran this code on his data:

#!matlab

LonLat(:,1,:)=east';
LonLat(:,2,:)=north';
XY=LonLat(:,:,2000);
b = braid (XY);
b = compact(b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Braid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(b)

This creates the horrible figure:

Here's a simpler code that breaks everything:

#!matlab

>> XY=randomwalk(1,10,.1);  % a braid of one strand!
>> b=braid(XY)
Error using braidlab.braid/set.n (line 196)
Too few strings for generators.

Error in braidlab.braid (line 183)
          br.n = max(abs(b))+1;
>> b

b = 

Error using horzcat
Dimensions of matrices being concatenated are not consistent.

Error in braidlab.braid/char (line 368)
      str = ['< ' str ' >'];

Error in braidlab.braid/disp (line 377)
        c = char(b);
>> plot(b) 

The output of plot(b) has serious issues!

The fix should be, I think, to make sure a braid of one strand is always the identity? Shouldn't be too hard to implement. Make sure plot(b) handles things ok for the identity braid, even for one strand.

Plotting Loops in Subplot Axis

In an attempt to plot multiple loops in subfigures, the function loop.plot fails to plot in the subplot axis. For example

#!matlab

l = braidlab.loop([1 1 0 0; 1 2 1 0])
subplot(2,1,1); l(1).plot; 
subplot(2,1,2); l(2).plot;

Will remove the subplot axis and plot as if there were no subplots.

Move TODO contents to issues.

The file TODO contains some ideas for improvements or possible problems: these should be converted to issues on BitBucket.

Empty array of braids

Initializing an empty array of braids results in incorrect dimensions.

E.g.

a(1,10) = braid;

creates a 10x10 matrix of braids, instead of 1x10 matrix.
The same happens if a(10) and a(10,1) are issued.

color_braiding: option to minimize crossings over projection line?

Is this a good idea?

How to best distribute?

What's the best way to distribute braidlab? The LCS toolbox is not quite ready, so might be best to create a tarball of the rest. Make script to tar the relevant bits.

What about compilation of MEX files? Could provide binaries for OSX, Windows, and Linux (32, 64). Would that be portable enough? Is static compilation an option?

Asymptotic spectral gap

It's probably easy to derive the asymptotics of the spectral gap for the psi braids. This could be used instead of calling psiroots in braid.entropy.

characteristic polynomial of a braid

This may be tough: return the characteristic polynomial of pA braids. This can be done iteratively if the braid is truly pA, I think. The problem might be to decide if the braid is pA or not. Perhaps limit the number of iterations until the min/max settle down?

Could also extract from Toby's code for shorter braids.

Bug in braid creation when some strands are equal.

Bug discovered while writing taffy_xrods.m (talk), found that if z is zero get 'somehow there are still coincident coordinates' error.

The problem is that the code only looks for coincident particles in one coordinate (along the projection line). If particles are coincident in both coordinates, need to issue an error since it's impossible to define a braid. Adding noise, as is done currently, will give a braid, but it should't.

Here's a first code that should return a more proper error:

#!matlab

>> braid(zeros(10,2,4))

This returns the error

#!matlab

Warning: Coincident coordinates... adding a bit of noise. 
> In @braid/private/color_braiding at 50
  In braid.braid>braid.braid at 126 
Error using color_braiding (line 95)
Somehow there are still coincident trajectories...

Error in braidlab.braid (line 126)
        br = color_braiding(b,1:size(b,1),secnd);

But the error "Somehow there are still coincident trajectories..." should never actually happen. It is a safety check.

The second, related bug is when two strands are equal:

#!matlab

npts = 10; n = 4; XY = zeros(npts,2,n);
% Create constant braid (fixed strands).
for i = 1:n, XY(:,:,i) = i*ones(npts,2); end
XY(:,:,1) = XY(:,:,2); % Make two strands equals.
% Bug: nontrivial braid returned, since routine adds noise.  Should
% return error instead.
braid(XY)

This returns

#!matlab

Warning: Coincident coordinates... adding a bit of noise. 
> In @braid/private/color_braiding at 50
  In braid.braid>braid.braid at 126
  In bug_color_braiding_coincident at 9 

ans = 

< -1  1 -1  1  1 -1 -1 >

The braid will vary, since it is random. The warning is appropriate, since the coordinates are indeed coincident, but it should be an error about coincident particles.

This should be easy to fix: before even computing anything else, check if any coordinates are coincident. Maybe this is a bit difficult to do fast?

Improve loop.plot

The options are a mess. I had to hack in a few changes for plot_loop_iter (part of my braid-talk), but it's really kludgy.

Make it such that we pass a list of options of the form 'Name',value.

jL

compact_helper: group same generators together

i.e., prefer < 1 1 1 3 > to < 1 1 3 1 >.

Is this worth doing?

Name for "open braid"

What's a good mathematical name for an "open braid", that is, a braid that doesn't necessarily close. These are the main braid objects used by braidlab. The form a groupoid, a group-like structure where not all the elements can be multiplied together. So maybe just call them "elements of the braid groupoid"? Kinda long. (Tangle is already taken.)

How about the process of turning an open braid into a true braid? "Closure" already means something different, but maybe that's the right thing. "Completion"? "Closed extension"? "Braid extension"? Start writing notes formalizing this. Can still define topological equivalence rel the punctures. Form a groupoid, a group-like structure where not all the elements can be multiplied together.

Implement Lawrence-Krammer representation

This is not so necessary, but could be fun. I've already done this in Mathematica somewhere (Braid.m, I think), so dig up the code first. But first solve issue #9, so that the LK rep can also have polynomial entries.

mpz_set_str in cyclemat

#!matlab
>> b = braid([4  3 -1  4  1  4  2 -3 -4  3 -2 -1  4  2  4 -3]);
>> b.cyclemat
Unexpected Standard exception from MEX file.
What() is:mpz_set_str
..

Error in loopsigma (line 56)
      [varargout{1:nargout}] = loopsigma_helper(ii,ustr);

Error in braidlab.braid/mtimes (line 330)
        [varargout{1:nargout}] = loopsigma(b1.word,vertcat(b2.coords));

Error in braidlab.braid/cycle (line 109)
  [l,pn] = b*l;

Error in braidlab.braid/cyclemat (line 56)
[pn,it] = cycle(b,varargin{:});

Bug when changing projection line in color_braiding

Bug discovered while writing taffy_xrods.m (talk). (This was originally filed with issue #2, but has been split into a separae issue.)

This is a very 'symmetric' case where noise is needed for the trajectories. By changing the projection angle it should have fixed it? But we still get an error.

#!matlab

npts = 200; r = .75; rodr = .05; n = 5;

th = linspace(0,2*pi,npts); th = th(end:-1:1);

z = zeros(npts,n);
z(:,1) = 0 + r*exp(1i*(th-pi));
z(:,2) = 0;
z(:,3) = 1 + r*exp(1i*th);
z(:,4) = 1;
z(:,5) = .5*ones(size(z(:,1)));

XY = zeros(npts,2,n); XY(:,1,:) = real(z); XY(:,2,:) = imag(z);

% Bug? If we specify a projection angle we sometimes get an error.
b = braid(XY,.1)

This returns the error

#!matlab

Error using color_braiding (line 181)
crossdat inconsistency at crossing 6, time 124.110878, index 2, with permutation [2
5  3  1  4].

Error in braidlab.braid (line 126)
        br = color_braiding(b,1:size(b,1),secnd);

The longer version attached also plots the trajectories and the braid.

cyclic equality of braids

When dealing with, say, periodic orbits, braids could be checked for "cyclic equality", meaning they are the same no matter where you start in the cycle. It would be nice to implement this in a "cycliceq" method.

example: [1 -2 3 -4] is cyclically equal to [1 -2 -4 3]:

[1 -2 3 -4] = [-2 3 -4 1] = [-2 3 1 -4] = [-2 1 3 -4] = [1 3 -2 -4]

[1 -2 -4 3] = [-2 -4 3 1] = [-2 -4 1 3] = [3 -2 -4 1] = [1 3 -2 -4]

This is the same as saying their closure gives the same knot?

Doesn't seem to be a matter of just permuting and trying possibilities:

#!matlab

>> a=braid([1 -2 3 -4]);
>> b=braid([1 -2 -4 3]); a==b
ans = 0
>> b=braid([-2 -4 3 1]); a==b
ans = 0
>> b=braid([-4 3 1 -2]); a==b
ans = 0
>> b=braid([3 1 -2 -4]); a==b
ans = 0

Why? Because we may have to "commute periodically" several times. Is this a hard problem? Or maybe they really aren't the same and I'm cheating somehow by using the periodicity several times.

Maybe extend (repeat) the braids until all the punctures return to their initial position. Would this help?

Polynomial in Burau matrices

The Burau representation of a braid right now only allows numerical specification of the parameter t. It would be better to have the option to return the Laurent polynomial. Could then implement Alexander polynomial.

braidlab.pdf - font error

Evince and xpdf on my Debian distribution do not display shaded lines containing code in the documentation braidlab.pdf

E.g. evince outputs:
"Error: could not create type1 face
some font thing failed"

Can't plot output of loopcoords?

For any loop obtained by loopcoords:

#!matlab
>> plot(loopcoords(braid([1])))
Error using linspace (line 21)
Inputs must be floats, namely single or double.

Error in braidlab.loop/plot (line 251)
    loop_curve_x = sign(nl)*linspace(0,rad,50);

subbraid method for databraid class

Currently databraid just uses its inherited superclass method braid.subbraid, which returns a braid object. It should really return a databraid, but this would mean preserving only the crossing times of the subbraid.

Check for existence of MEX files, use backup code

Braidlab should gracefully handle the absence of a compile MEX file. In most cases it should be possible to substitute Matlab code (i.e., compact) or at least return an error explaining things.

Should we also try to compile the file if it exists?

Check types when invoking MEX files.

Calling the MEX helper files with the wrong argument type can cause a hard crash (say, double vs int loops). Make sure this is checked better.

Extract reducing curve from Toby Hall's trains code?

Is it easy to parse the output in trains, and convert to Dynnikov coordinates? That would be very useful.

See the doc for Toby's code

Use GMP when finding entropy?

The current method braid.entropy uses the iterative method, or optionally the train-track method. Recently (d766e7a) I wrote a unit test for this method, and I find it hard to get accurate results for low-entropy braids (braids such as braid('psi',20)). Even increasing the number of iterations doesn't help much.

Could this be resolved by using arbitrary precision arithmetic (like GNU's GMP library) within a MEX file? It would be good to be able to theoretically get any precision we want, up to machine accuracy.

braid.compact may overflow precision for long lists of generators

Before compactifying the braid,
braid.compact calls istrivial, which in turns calls loopcoords, which has 'int64' as its default data type. For long strings of generators, this means that just checking whether the braid is trivial might overflow the int64 and cause an error.

This should be handled more elegantly. istrivial was originally introduced to make sure that empty braids do not crash the compactification process.

• So perhaps just short-circuiting this test would be sufficient.
• Alternatively, perhaps using 'double' precision explicitly here would work as a bandaid solution.

Check that MEX is on the path in Makefiles

MRA had problem compiling and it took us a while to realize that mex wasn't on the path. Makefile should check that.

insertboiler somehow modified .mat file

In rev 2881524bf14a, the file testdata.mat was corrupted by running devel/inserboiler. I've reverted it in 4785de700025 but I'd like to understand how that happened.

Asymptotic entropy for one large braid?

Given a braid of n strands, obtained from a smooth dynamical system, can we predict the asymptotic entropy by sampling substrands? That is, there are nchoosek braids of size k, which is a huge number for most k. If we randomly pick and average the entropy, can we get a plot like Marko's and predict the asymptote, or does the fact that we're working from a single braid kill us?

Bibliography

Eventually start adding a list of projects that use braidlab, or perhaps Dynnikov coordinates in general. Maybe could put these in the braidlab guide.

Warn of overflow of loops

Matlab is not great for manipulating types. Margaux was confused by something she rightfully expected to work. Load the attached MAT file, then do

#!matlab

b = compact(braid(XY));
entropy(b)
   ans = 5.1058
entropy(b,'trains')
   ans = 5.1058    % answers agree
l = loopcoords(b);
log(minlength(b^100*l)/minlength(l)) / 100   % find entropy by iteration
   ans = 0.3556

These are radically different, and they shouldn't be. The problem is that since loopcoords tries to use int, it can overflow.

Compare this to using way fewer iterations:

#!matlab

log(minlength(b^6*l)/minlength(l)) / 6
  ans = 5.0872

or using a loop object directly:

#!matlab

l = loop(b.n);
log(minlength(b^100*l)/minlength(l)) / 100
  ans = 5.1082

which are both close to the real answer.

When I wrote this I looked for a way to check for overflow: it's not trivial. It's also possible to force use of exact types, but there are no built-in ones. I tried the package VariablePrecisionIntegers.

But in any case, there is a flag 'checkoverflow' that is set to true, so something should have complained. The problem, I think, is that the check is done in loopcoords but not in the multiplication.

Implement braid tensor product

This is simply the lateral joining of two diagrams, with no extra crossings. Call it tensor(a,b)? Not a method but a friend?

Empty braid acting on loop

Problem:

#!matlab
>> b=braid([],3);
>> l = loop(3);
>> b^0*l
Index exceeds matrix dimensions.

Error in braidlab.braid/mtimes (line 331)
        varargout{1} = braidlab.loop(varargout{1});

compact fails on very large braids

Works for random braid of length 1000, but not 10000. I don't see why not.

#!matlab

>> rng('default'); b=braid('random',10,1000); b2 = b*inv(b); b3 = compact(b2); isempty(b3.word)

ans =

     1

>> rng('default'); b=braid('random',10,10000); b2 = b*inv(b); b3 = compact(b2); isempty(b3.word)

ans =

     0

Entropy overflow computation

I have been looking into the entropy overflow computation, so in the commit 3de3219 I have added a function that computes longer and longer random braids in an effort to zero-in on the point at which entropy computation overflow.

Entropy is computed as

entropy(B, [], 1)

to minimize any overflow due to iterative computation.

This is the graph and it appears to show that the more strands take longer to overflow. It goes against my experience from dynamical system simulations, where braids with more strands ended up running into overflow sooner. However, more strands usually means that there are more generators per unit of physical time so it's not clear if the correlation is expected or not.

Even though generators are random, every run re-seeds the RNG so the list of integers specifying the generators is fixed for each given number of strands.

I am not sure whether this is expected behavior or not, but if we're working towards solutions that increase number of generators allowed before the overflow happens, then this could be the benchmark that we work off of.

Add testsuite: testLoop.m

The test in changeset 42e4e2c should be moved to make a new loop-testing file in the testsuite.

Clean up README

Right now the README is being displayed in the welcome page, but it's not very helpful. Move the contents of that file somewhere else (or get rid of most of it), then make a proper README.

loopcoords: automatically set type?

In loopcoords.m, currently the user can specify the data type for loops. The default is int64. Maybe better to do like braid.eq and catch an overflow, and use VPI in those cases (print warning?)? This makes the function simpler. Could leave the option to specify a type anyways, for the power user.

Wrap long braids and loops

Currently long braids are not wrapped:

#!matlab

>> b=braid('random',10,100)

b = 

< 2  8 -3  5  2 -6 -3  6  7  7 -5  1  3 -9  2 -8 -5  9  1  4  1 -9  1  7  8  8  1 -4  3 -8  4 -9 -2 -3 -2  2  8 -6  5 -2 -8 -6  4  5 -4 -1 -3  2  2  3 -4  1  9 -9 -5 -5  4 -9  4 -2  8 -4  3  4  1  2  9 -9 -6 -1  3 -4  8 -1  1 -2 -6 -7  6 -5 -5  3  7  2  7 -2  4  6 -8  1 -9 -7  5  4  5 -3  5 -5 -8 -8 >
>> l=loop(100)

l = 

(( 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ))

Probably they should be? The output should follow the setting of format loose or format compact. We need to know the display width to do that.

We can get the current window size with

#!matlab

>> get(0, 'CommandWindowSize')
ans =
    88    37

See also the command textwrap.

Fix equality for long braids

Now that we have VPI in place (see discussion in issue #11), we can address a problem in the equality of long braids. The current code for eq in braid.m reads:

#!matlab

      % Check if the loop coordinates are the same.
      % This can fail if the braids are too long, since the coordinates
      % overflow.  Check for that.
      ee = all(loopcoords(b1) == loopcoords(b2));
    end

The problem is that for long braids the int64 used by default in loopcoords can easily overflow.

#!matlab

>> b=braid([1 -2])

b = 

< 1 -2 >
>> b^50 == b^50
Error using sumg (line 28)
Summation of 7540113804746346429 and 7540113804746346428 has overflowed.

Error in sumg (line 19)
  out = sumg( varargin{1:floor(end/2)}, sumg(varargin{floor(end/2)+1:end}) );

Error in sumg (line 19)
  out = sumg( varargin{1:floor(end/2)}, sumg(varargin{floor(end/2)+1:end}) );

Error in loopsigma (line 56)
      d = sumg(a(:,i-1), -a(:,i), pos(b(:,i)), -neg(b(:,i-1)));

Error in braidlab.braid/loopcoords (line 89)
l = braidlab.loop(loopsigma(w,htyp(l.coords)));

Error in braidlab.braid/eq (line 182)
      ee = all(loopcoords(b1) == loopcoords(b2));

Use try-catch block and retry with vpi?

The closure issue

The function braid.closure closes a braid, by default using the most direct method (no new crossings in projection), but also by minimizing L^2. Are these the best choices? Why not, say, the lowest entropy?

Maybe find some cases where we have a target entropy, that is, we can generate a very long braid so we can compute the entropy, but we'd like to estimate the entropy from a subset of that. (Use Duffing in the chaotic sea?).

How do we avoid the (n!) scaling of such methods?

Plotting braids issue

This is similar to issue #43:

#!matlab

imagesc([1 1]); plot(braid([1 2]))

doesn't clear the axes properly. Why not? I tried copying the cla from loop.plot, but it's not helping.

Note that

#!matlab

imagesc([1 1]); plot(braid([1 2])); plot(braid([1 2]))

works, that is, the second time around it clears the plot. However,

#!matlab

imagesc([1 1]); plot(braid); plot(braid)

doesn't work, but plot(braid) by itself does.

Number of connected components in a loop

Write a method to find this.

Alexander polynomial

If/when issue #9 is addressed, can use to computer the Alexander polynomial of a braid.

Array display

When an empty array of braids is created, Matlab fails in 'disp' method of the braid as it immediately tries to access .word property, even though an array has been passed to it.

E.g.

a(1,10) = braid

fails.

compact_helper.cpp unsigned expression warning

I've just compiled braidlab on Mac and it goes through fine, but its compiler throws a warning midway. Is this a potential bug?

#!shell

compact_helper.cpp:153:23: warning: comparison of unsigned expression >= 0 is always true [-Wtautological-compare]
          if (i+2*dir >= 0 && i+2*dir <= b.size()-1)
              ~~~~~~~ ^  ~
compact_helper.cpp:221:10: note: in instantiation of function template specialization
      'commute_and_cancel<std::__1::vector<int, std::__1::allocator<int> > >' requested here
  while (commute_and_cancel(bw,1) || commute_and_cancel(bw,-1)) {}
         ^
compact_helper.cpp:182:24: warning: comparison of unsigned expression >= 0 is always true [-Wtautological-compare]
        } while (i+dir >= 0 && i+dir <= b.size()-1);
                 ~~~~~ ^  ~

Compacting a trivial braid does nothing

a = braidlab.braid([1,-2])

compact(a*inv(a)) should return , but it returns <1 -2 2 1>

As seen in BIRS Tutorial.

Braid arrays vs. cell arrays

Issue #30 raised the problem of whether we should allow for Matlab arrays of braids.

Cell array syntax for initialization

a{10} = braid(1);

creates a cell array with 9 empty cells, with the 10th containing the braid.
Conversely,

a(10) = braid(1),

creates an array with 9 identity (default) braids, with the 10th element initialized to the generator <1>

A user might expect that array syntax implicitly allows vector operations, e.g., element-wise multiplication of two braid vectors, or operation by a set of braids on a single loop. We should decide whether to allow this and how to implement it/make it apparent to users.

Note that just forbidding creation of braid arrays might be difficult from Matlab syntax.

color_braiding: perturb projection line rather than trajectories

color_braiding: perturb projection line rather than trajectories when coincidence found? Or just report error message. Maybe better to just report error and let user deal with it by changing the projection line.

Warning label typo?

In some .m files, it appears that error/warning labels have a typo. Running

find . -name *.m | xargs grep BRAIDLAD

in braidlab/ directory reveals them.

This is trivial to fix, namely run

find . -name *.m | xargs sed --in-place 's/BRAIDLAD/BRAIDLAB/g'

in braidlab/ directory, but I'm not sure if this is really a typo, or is there a reason for these labels to be different.

Geometric intersection number of two loops

Is this hard to do?