SymPy and the GSoC

Google has announced the results for Google Summer of Code. The following projects have been accepted for SymPy:

(Project, Student, Mentor, Link to proposal on the wiki)
- Category Theory Module, Sergiu Ivanov, Tom Bachmann
- Density Operators for Quantum Module in sympy.physics.quantum, Guru
Devanla, Brian Granger (co-mentor Sean Vig)
- Enhancements to sympy.physics.mechanics, Angadh Nanjangud, Gilbert Gede
- Group Theory, Aleksandar Makelov, David Joyner (Aaron Meurer co-mentor)
- Implicit Plotting Module, Bharath M R, Aaron Meurer
- Vector Analysis, Stefan Krastanov, Matthew Rocklin

I will help mentor Aleksandar Makelov’s work on group theory. He is a freshman at Harvard.

Bifid cipher and Sage

The Bifid cipher was invented around 1901 by Felix Delastelle. It is a “fractional substitution” cipher, where letters are replaced by pairs of symbols from a smaller alphabet. The cipher uses a 5×5 square filled with some ordering of the alphabet, except that “i”‘s and “j”‘s are identified (this is a so-called Polybius square; there is a 6×6 analog if you add back in “j” and also append onto the usual 26 letter alphabet, the digits 0, 1, …, 9). According to Helen Gaines’ book “Cryptanalysis”, this type of cipher was used in the field by the German army during World War I.

The Bifid cipher was discusses in Alasdair McAndrew’s book on Cryptography and Sage. We shall follow his discussion. As an example of a Polybius square for the Bifid cipher, pick the key to be “encrypt” (as Alasdair does). In that case, the Polybius square is BTW, the analog is: .

Here is Sage code to produce the case (the case is in Alasdair’s book):

def bifid(pt, key):
    """
    INPUT:
        pt - plaintext string      (digits okay)
        key - short string for key (no repetitions, digits okay)
    
    OUTPUT:
        ciphertext from Bifid cipher (all caps, no spaces)

    This is the version of the Bifid cipher that uses the 6x6
    Polybius square.

    AUTHOR:
        Alasdair McAndrew

    EXAMPLES:
        sage: key = "encrypt"
        sage: pt = "meet me on monday at 8am"
        sage: bifid(pt, key)
        [[2, 5], [0, 0], [0, 0], [1, 0], [2, 5], [0, 0], [3, 0], 
         [0, 1], [2, 5], [3, 0], [0, 1], [1, 3], [1, 1], [0, 4], 
         [1, 1], [1, 0], [5, 4], [1, 1], [2, 5]]
        'HNHOKNTA5MEPEGNQZYG'

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [x.capitalize() for x in key if not(x.isspace())]
    pt0 = [x.capitalize() for x in pt if not(x.isspace())]
    # create long key
    long_key = key0+[x for x in A if not(x in key0)]
    n = len(pt0)
    # the fractionalization
    pairs = [[long_key.index(x)//6, long_key.index(x)%6] for x in pt0]
    print pairs
    tmp_cipher = flatten([x[0] for x in pairs]+[x[1] for x in pairs])
    ct = join([long_key[6*tmp_cipher[2*i]+tmp_cipher[2*i+1]] for i in range(n)], sep="")
    return ct

def bifid_square(key):
    """
    Produced the Polybius square for the 6x6 Bifid cipher.

    EXAMPLE:
        sage: key = "encrypt"
        sage: bifid_square(key)
        [E N C R Y P]
        [T A B C D E]
        [F G H I J K]
        [L M N O P Q]
        [R S T U V W]
        [X Y Z 0 1 2]

    """
    AS = AlphabeticStrings()
    A = list(AS.alphabet())+[str(x) for x in range(10)]
    # first make sure the letters are capitalized
    # and text has no spaces
    key0 = [SR(x.capitalize()) for x in key if not(x.isspace())]
    # create long key
    long_key = key0+[SR(x) for x in A if not(x in key0)]
    # the fractionalization
    pairs = [[long_key.index(SR(x))//6, long_key.index(SR(x))%6] for x in A]
    f = lambda i,j: long_key[6*i+j] 
    M = matrix(SR, 6, 6, f)
    return M

Have fun!

uniform matroids and MDS codes

It is known that uniform (resp. paving) matroids correspond to MDS (resp. “almost MDS” codes). This post explains this connection.

An MDS code is an linear error correcting block code which meets the Singleton bound, . A uniform matroid is a matroid for which all circuits are of size , where is the rank of M. Recall, a circuit in a matroid M=(E,J) is a minimal dependent subset of E — that is, a dependent set whose proper subsets are all independent (i.e., all in J).

Consider a linear code whose check matrix is an matrix . The vector matroid M=M[H] is a matroid for which the smallest sized dependency relation among the columns of H is determined by the check relations , where is a codeword (in C which has minimum dimension d). Such a minimum dependency relation of H corresponds to a circuit of M=M[H].

"""

Implements a Matroid class and some methods.

A finite matroid M is a pair (E, I), where E is a finite set

and I is a collection of subsets of E (called the independent sets)

with the following properties:
 * The empty set is independent.

 * Every subset of an independent set is independent.

 * If A and B are two independent sets of I and A has

   more elements than B, then there exists an element in A

   which is not in B that when added to B still gives an independent set.
A subset of E that is not independent is called dependent.

A maximal independent set is called a basis for the matroid.

Any two bases of a matroid M have the same number of elements.

This number is called the rank of M. A circuit in a matroid M is

a minimal dependent subset of E.
REFERENCES:
http://en.wikipedia.org/wiki/Matroid
"""
class Matroid_generic():

    """

    Implements a generic matroid.
    """

    def __init__(self, X, IndepSets = None, Bases = None):

        """

        A matroid can either be specified by its independent

        sets, its bases or its circuits (probably not both).

        This class provides methods for converting one type

        into the others.
        TODO:

           Add "Circuits = None" option.
        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0,1],[0,1,0,1,0,1,1],[0,0,1,0,1,1,1]])

            sage: Fano = vector_matroid(A)

            sage: Fano

            Matroid([0, 1, 2, 3, 4, 5, 6], Bases = [[0, 2, 3],

            [1, 3, 6], [1, 2, 4], [0, 3, 5], [1, 5, 6], [0, 1, 5],

            [2, 3, 4], [0, 1, 4], [0, 1, 2], [3, 5, 6], [3, 4, 6],

            [0, 1, 6], [1, 4, 5], [1, 3, 4], [1, 3, 5], [0, 3, 6],

            [0, 2, 6], [0, 4, 6], [0, 2, 5], [2, 5, 6], [2, 4, 6],

            [1, 2, 3], [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 6],

            [2, 4, 5], [4, 5, 6]])

            sage: Fano.circuits()

            [[0, 1, 3], [0, 2, 4], [0, 5, 6], [1, 2, 5], [1, 4, 6],

             [2, 3, 6], [3, 4, 5], [0, 1, 2, 6], [0, 1, 4, 5],

             [0, 2, 3, 5], [0, 3, 4, 6], [1, 2, 3, 4], [1, 3, 5, 6],

             [2, 4, 5, 6]]

        """

        self.ground_set = X

        self.bases = Bases

        if Bases != None:

            S = {}

            for B in Bases:

                SB = Set(B).subsets()

                S.union(SB)

            self.indep_sets = [list(x) for x in Set(S)]

            J = Set([Set(x) for x in self.indep_sets])

            D = [list(x) for x in SX if not(x in J)]

            self.dep_sets = D

        if IndepSets != None:

            self.indep_sets = IndepSets

            J = Set([Set(x) for x in self.indep_sets])

            M = max([x.cardinality() for x in J])

            B = [list(x) for x in J if x.cardinality() == M]

            self.bases = B

            SX = Set(X).subsets()

            D = [list(x) for x in SX if not(x in J)]

            self.dep_sets = D

        # now create the circuits - BUG

        C = []

        for x in D:

            addx = [Set([b for b in x if b != a]) in J for a in x]

            if not(False in addx):

                C.append(x)

        self.circuts = [list(x) for x in C]
    def __repr__(self):

        return "Matroid(%s, Bases = %s)"%(self.base_set(), self.bases())
    def __str__(self):

        return "Matroid on base set %s of rank %s"%(self.base_set(), self.rank())
    def basis(self):

        """

        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

            sage: M = vector_matroid(A)

            sage: M.basis()

            [[0, 2, 3], [1, 2, 4], [0, 3, 5], [0, 1, 5], [2, 3, 4],

             [0, 1, 4], [0, 1, 2], [0, 2, 5], [1, 4, 5], [1, 3, 5],

             [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 3], [1, 3, 4],

             [2, 4, 5]]

        """

        return self.bases
    def circuits(self):

        """

        Returns the minimal dependent sets.
        EXAMPLES:

            sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

            sage: M = vector_matroid(A)

            sage: M.circuits()

            [[0, 1, 3], [0, 2, 4], [1, 2, 5], [3, 4, 5], [0, 1, 4, 5],

             [0, 2, 3, 5], [1, 2, 3, 4]]

        """

        return self.circuts
    def dependent_sets(self):

        return self.dep_sets
    def base_set(self):

        return self.ground_set
    def independent_sets(self):

        S = self.indep_sets

        S.sort()

        return S
    def rank(self):

        r = len(self.bases()[0])

        return r
def vector_matroid(mat):

    """

    Returns the vector matroid defined by the matrix M
    EXAMPLES:

        sage: A = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: M = vector_matroid(A)

        sage: M.rank()

        3

        sage: print M

        Matroid on base set [0, 1, 2, 3, 4, 5] of rank 3

        sage: M.base_set()

        [0, 1, 2, 3, 4, 5]

        sage: M.basis()

        [[0, 2, 3], [1, 2, 4], [0, 3, 5], [0, 1, 5], [2, 3, 4],

         [0, 1, 4], [0, 1, 2], [0, 2, 5], [1, 4, 5], [1, 3, 5],

         [0, 3, 4], [0, 4, 5], [2, 3, 5], [1, 2, 3], [1, 3, 4], [2, 4, 5]]

        sage: M.circuits()

        [[1], [2], [2], [4], [1, 3], [1, 4], [1, 5], [1, 4], [1, 5],

         [1, 5], [2, 4], [2, 5], [2, 5], [3, 5], [2, 4], [2, 5], [2, 5],

         [3, 5], [3, 5], [1, 3], [1, 3], [1, 4], [1, 4], [2, 4], [2, 4], [1, 3, 5]]

    """

    F = mat.base_ring()

    n = len(mat.columns())

    k = len(mat.rows())

    J = Combinations(n,k)

    indE = []

    for x in J:

        M = matrix([mat.column(i) for i in x])

        if k == M.rank():     # all indep sets of max size

            indE.append(x)

            for y in powerset(x): # all smaller indep sets

                if not(y in indE):

                    indE.append(y)

    return Matroid_generic(range(n), IndepSets=indE)
def bases(mat, invariant = False):

    """

    Find all bases in a vector/representable matroid.
    EXAMPLES:

        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: L = [x[1] for x in bases(M, invariant=True)]; L.sort(); L

        [5, 16, 17, 33, 37, 44, 45, 48, 81, 84, 92, 93, 96, 112, 113, 116]

        sage: C = ReedSolomonCode(6,4,GF(7))

        sage: G = C.gen_mat(); G

        [1 1 1 1 1 1]

        [0 1 2 3 4 5]

        [0 1 4 2 2 4]

        [0 1 1 6 1 6]

        sage: B = bases(G); B

        [[0, 1, 2, 3], [0, 1, 2, 4], [0, 1, 2, 5], [0, 1, 3, 4],

         [0, 1, 3, 5], [0, 1, 4, 5], [0, 2, 3, 4], [0, 2, 3, 5],

         [0, 2, 4, 5], [0, 3, 4, 5], [1, 2, 3, 4], [1, 2, 3, 5],

         [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]

        sage: len(B)

        15

        sage: binomial(6,4)

        15

    """

    r = mat.rank()

    ind = independent_sets(mat)

    B = []

    for x in ind:

        if len(x) == r:

            B.append(x)

    B.sort()

    if invariant == True:

        B = [(b,sum([k*2^(k-1) for k in b])) for b in B]

    return B
def circuits(mat):

    """

    Finds all circuits in a vector matroid.
    EXAMPLES:

        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: M

        [1 0 0 1 1 0]

        [0 1 0 1 0 1]

        [0 0 1 0 1 1]

        sage: circuits(M)

        [[0, 1, 3], [0, 2, 4], [1, 2, 3, 4], [1, 2, 5], [0, 2, 3, 5],

        [0, 1, 4, 5], [3, 4, 5]]

        sage: M = matrix(GF(3), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])

        sage: circuits(M)

        [[0, 1, 3], [0, 2, 4], [1, 2, 3, 4], [1, 2, 5], [0, 2, 3, 5],

        [0, 1, 4, 5], [0, 3, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]
    """

    J = independent_sets(mat)

    n = len(mat.columns())

    E = set(range(n))

    subsetsE = powerset(E)

    C = []

    for S in subsetsE:

        is_circuit = True

        if not(S in J):

            for x in powerset(S):

                if xS and not(x in J):

                    is_circuit = False

        if is_circuit == True:

            if not(S in J):

                C.append(S)

    return C

                

             

Floyd-Warshall-Roy, 3

As in the previous post, let be a weighted digraph having vertices and let denote its adjacency matrix. We identify the vertices with the set .

The previous post discussed the following result, due to Sturmfels et al.

Theorem: The entry of the matrix A ^n-1 in row i and column j equals the length of a shortest path from vertex i to vertex j in G. (Here A ^n-1 denotes the n-1-st tropical power of A.)

This post discusses an implementation in Python/Sage.

Consider the following class definition.

class TropicalNumbers:
    """
    Implements the tropical semiring.

    EXAMPLES:
        sage: T = TropicalNumbers()
        sage: print T
        Tropical Semiring

    """
    def __init__(self):
        self.identity = Infinity

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:
            sage: TropicalNumbers()
            TropicalNumbers()

        """
        return "TropicalNumbers()"

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T
            Tropical Semiring

        """
        return "Tropical Semiring"

    def __call__(self, a):
        """
        Coerces a into the tropical semiring.

        EXAMPLES:
            sage: T(10)
            TropicalNumber(10)
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return TropicalNumber(a)

    def __contains__(self, a):
        """
        Implements "in".

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True

        """
        if a in RR or a == Infinity:
            return a==Infinity or (RR(a) in RR)
        else:
            return a==Infinity or (RR(a.element) in RR)

class TropicalNumber:
    def __init__(self, a):
        self.element = a
        self.base_ring = TropicalNumbers()

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalNumber(%s)"%self.element

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: print T(10)
            Tropical element 10 in Tropical Semiring
        """
        return "%s"%(self.number())

    def number(self):
        return self.element

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalNumber class.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a+b
            10

        """
        T = TropicalNumbers()
        return T(min(self.element,other.element))

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: T = TropicalNumbers()
            sage: a = T(10)
            sage: a in T
            True
            sage: b = T(15)
            sage: a*b
            25
        """
        T = TropicalNumbers()
        return T(self.element+other.element)

class TropicalMatrix:
    def __init__(self, A):
        T = TropicalNumbers()
        self.base_ring = T
        self.row_dimen = len(A)
        self.column_dimen = len(A[0])
        # now we coerce the entries into T
        A0 = A
        m = self.row_dimen
        n = self.column_dimen
        for i in range(m):
            for j in range(n):
                A0[i][j] = T(A[i][j])
        self.array = A0

    def matrix(self):
        """
        Returns the entries (as ordinary numbers).

        EXAMPLES:
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: AT.matrix()
            [[0, 1, 3, 7], [2, 0, 1, 3], [4, 5, 0, 1], [6, 3, 1, 0]]
        """
        m = self.row_dim()
        n = self.column_dim()
        A0 = [[0 for i in range(n)] for j in range(m)]
        for i in range(m):
            for j in range(n):
                A0[i][j] = (self.array[i][j]).number()
        return A0

    def row_dim(self):
        return self.row_dimen

    def column_dim(self):
        return self.column_dimen

    def __repr__(self):
        """
        Called to compute the "official" string representation of an object.
        If at all possible, this should look like a valid Python expression
        that could be used to recreate an object with the same value.

        EXAMPLES:

        """
        return "TropicalMatrix(%s)"%self.array

    def __str__(self):
        """
        Called to compute the "informal" string description of an object. 

        EXAMPLES:

        """
        return "Tropical matrix %s"%(self.matrix())

    def __add__(self, other):
        """
        Implements +. Assumes both self and other are instances of
        TropicalMatrix class.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: B = [[2,Infinity,1],[3,-1,1]]
            sage: AT = TropicalMatrix(A)
            sage: BT = TropicalMatrix(B)
            sage: AT
            TropicalMatrix([[TropicalNumber(1), TropicalNumber(2), TropicalNumber(+Infinity)],
              [TropicalNumber(3), TropicalNumber(+Infinity), TropicalNumber(0)]])
            sage: AT+BT
            [[TropicalNumber(1), TropicalNumber(2), TropicalNumber(1)],
             [TropicalNumber(3), TropicalNumber(-1), TropicalNumber(0)]]
        """
        A = self.array
        B = other.array
        C = []
        m = self.row_dim()
        n = self.column_dim()
        if m != other.row_dim:
            raise ValueError, "Row dimensions must be equal."
        if n != other.column_dim:
            raise ValueError, "Column dimensions must be equal."
        for i in range(m):
            row = [A[i][j]+B[i][j] for j in range(n)] # + as tropical numbers
            C.append(row)
        return C

    def __mul__(self, other):
        """
        Implements multiplication *.

        EXAMPLES:
            sage: A = [[1,2,Infinity],[3,Infinity,0]]
            sage: AT = TropicalMatrix(A)
            sage: B = [[2,Infinity],[-1,1],[Infinity,0]]
            sage: BT = TropicalMatrix(B)
            sage: AT*BT
             [[TropicalNumber(1), TropicalNumber(3)],
              [TropicalNumber(5), TropicalNumber(0)]]
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
            sage: AT = TropicalMatrix(A)
            sage: print AT*AT*AT
            Tropical matrix [[0, 1, 2, 3], [2, 0, 1, 2], [4, 4, 0, 1], [5, 3, 1, 0]]
        """
        T = TropicalNumbers()
        A = self.matrix()
        B = other.matrix()
        C = []
        mA = self.row_dim()
        nA = self.column_dim()
        mB = other.row_dim()
        nB = other.column_dim()
        if nA != mB:
            raise ValueError, "Column dimension of A and row dimension of B must be equal."
        for i in range(mA):
            row = []
            for j in range(nB):
                c = T(Infinity)
                for k in range(nA):
                    c = c+T(A[i][k])*T(B[k][j])
                row.append(c.number())
            C.append(row)
        return TropicalMatrix(C)

This shows that the shortest distances of digraph with adjacency matrix is equal to A ³, which is equal to This verifies an example given in chapter 1 of the book by Maclagan and Sturmfels, Introduction to Tropical Geometry .

Sage and the future of mathematics

I am not a biologist nor a chemist, and maybe neither are you, but suppose we were. Now, if I described a procedure, in “standard” detail, to produce a result XYZ, you would (based on your reasonably expected expertise in the field), follow the steps you believe were described and either reproduce XYZ or, if I was mistaken, not be able to reproduce XYZ. This is called scientific reproducibility. It is cructial to what I believe is one of the fundamental principles of science, namely Popper’s Falsifiability Criterion.

More and more people are arguing, correctly in my opinion, that in the computational realm, in particular in mathematical research which uses computational experiments, that much higher standards are needed. The Ars Technica article linked above suggests that “it’s time for a major revision of the scientific method.” Also, Victoria Stodden argues one must “facilitate reproducibility. Without this data may be open, but will remain de facto in the ivory tower.” The argument basically is that to reproduce computational mathematical experiments is unreasonably difficult, requiring more that a reasonable expertise. These days, it may in fact (unfortunately) require purchasing very expensive software, or possessing of very sophisticated programming skills, accessibility to special hardware, or (worse) guessing parameters and programming procedures only hinted at by the researcher.

Hopefully, Sage can play the role of a standard bearer for such computational reproducibility. Sage is free, open source and there is a publically available server it runs on (sagenb.org).

What government agencies should require such reproducibility? In my opinion, all scientific funding agencies (NSF, etc) should follow these higher standards of computational accountability.

The enlightenment of Professor Bigglesnot

An enthusiastic “Yippee!” echoed down the corridor. so loud it woke several faculty members in nearby offices. Some even got up out of their chairs and looked up and down the hallway before returning to grading or research or freecell before falling asleep again. But Bigglesnot was excited. After all, computing automorphism semigroups of quantum hyperalgebras was his life’s passion, ever since he was a graduate student. In front of him, was the latest issue of the Quantum Hyperalgebra Journal, newly released from its plastic shrink-wrap. It was opened to the article which was the focus of Bigglesnot’s attention – the esteemed Ziggotwat’s discussion of a new algorithm to compute automorphism semigroups of quantum hyperalgebras. Bigglesnot could see immediately from the tables of new data presented that Ziggotwat’s implementation was faster, better and more general than his own program. Whispering “Awesome! Awesome! So, awesome! …” under his breath, he shot off an email to Ziggotwat asking for more information, and, if at all possible, further details on the implementation. Could he please post or email the code, for others to look at? That would be awesome!

Days and weeks went by, but with no reply. One morning B found an email from Z: “…the code needs to be cleaned up first …”, “… so sorry for the late reply …”, but “.. thanks for your interest!”, was the gist. As luck would have it, B spotted Z a few weeks later at the annual meeting of the Society of Quantum Hyperalgebraists. Undaunted, one night after the talks of the meeting were finished, B bombarded Z with free beer and flattery peppered with questions about his program. “In fact”, Z finally confessed, “all the work was done by my former student Pipperpop, who has graduated and does not reply to my emails. I can send you what I have – but no promises!” Drunk, but now estatic, Bigglesnot managed to say “Awesome!” before he fell off his barstool.

In the months that followed, B pored through the incomplete, undocumented code. It was provided as a sequence of files, each one seeming dependent on another. They wouldn’t compile, no matter how B tried. Each day for a month, after attending to his classes, B would try to modify one of the files, hoping that a small change would allow him to compile the files into a functioning program. Each day, he would draft an email to Z (or to P, or to the editor of QHJ) asking what kind of LSD was he tripping on to make him high enough to think this mess of code would ever compile?
And each day, he wisely deleted the draft. Bigglesnot, usually filled to overflowing with self-confidence, was defeated.

Finally, B realized the solution! Not the solution of to how to compile Pipperpop’s poop, but the solution to the general problem. Software computations submitted to scientific journals must be “open” – if scientific data obtained as a result of a software computation is part of research paper submitted for publication then the source code for that software must also be made public and verified before publication. Enlightened, Bigglesnot was optimistic once again.

---

This postwas inspired by the excellent paper: Pederson, Empiricism is Not a Matter of Faith, Computational Linguistics, Volume 34, Number 3, pp. 465-470, September 2008.
http://www.d.umn.edu/~tpederse/full-pubs.html