JuMP-ing with AlgebraicJulia I: The IJKLM model

Conjunctive Queries on Acsets

The critical thing that the JuMP devs did to speed thing up was to replace the for loops with 2 inner joins, to get the “paths” through the relations. How to do this with acsets? Well one thing to do is execute a conjunctive query on the acset to get the same thing. This is described in a C-sets for data analysis: relational data and conjunctive queries.

connected_paths_query = @relation (i=i,j=j,k=k,l=l,m=m) begin
    IJK(IJK_I=i, IJK_J=j, IJK_K=k)
    JKL(JKL_J=j, JKL_K=k, JKL_L=l)
    KLM(KLM_K=k, KLM_L=l, KLM_M=m)
end

Catlab.to_graphviz(connected_paths_query, box_labels=:name, junction_labels=:variable, graph_attrs=Dict(:dpi=>"72",:size=>"3.5",:ratio=>"expand"))

While the blog post should be consulted for a complete explanation, the conjunctive query is expressed using an undirected wiring diagram (UWD) which is visualized above. Nodes (labeled ovals) in the UWD correspond to tables (primary keys) in the acset. Junctions (labeled dots) correspond to variables. Ports, which are unlabed in this graphical depiction, are where wires connect junctions to nodes. These correspond to columns of the table they are connected to. Outer ports, which are wires that run “off the page”, are the columns of the table that will be returned as the result of the query. Conceptually, rows that are returned from the query come from filtering the Cartesian product of the tables (nodes) such that variables in columns match according to ports that share a junction.

The JuMP blog post notes that while the data frames version doesn’t resemble the nested summation it is arguably just as readable, especially if the columns were related to the process that was being modeled. We suggest that the acsets version is also just as readable, if not more, as the data schema and query diagram directly represent the data that parameterizes the optimization model. Furthermore because the schema of the acset is known at compile time, incorrect queries (or other operations) on acsets will be caught as compile time errors.

The query can then be evaluated on the specific acset instance. We can confirm that both the acsets and data frame methods return the same number of rows.

ijklm_query = query(ijklm_acs, connected_paths_query)
size(ijklm_query) == size(ijklm_df)

true

Now we can go ahead and see how fast the acsets version is. The fact that the acsets based query is right on the tails of the DataFrames version is a performance win for the acsets library, as it is usually a generic conjunctive query engine across very general data structures (i.e., acsets are in general much more complex than a single dataframe, due to presence of multiple tables connected via foreign keys).

@benchmark let
    ijklm = query(ijklm_acs, connected_paths_query)
    model = JuMP.Model()
    set_silent(model)
    ijklm[!, :x] = @variable(model, x[1:size(ijklm, 1)] >= 0)
    for df in DataFrames.groupby(ijklm, :i)
        @constraint(model, sum(df.x) >= 0)
    end
end

BenchmarkTools.Trial: 2759 samples with 1 evaluation.
 Range (min … max):  1.549 ms …   9.245 ms  ┊ GC (min … max): 0.00% … 20.42%
 Time  (median):     1.685 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.808 ms ± 446.566 μs  ┊ GC (mean ± σ):  5.04% ± 10.74%

    ▂██▆▃                                             ▁▁▁     ▁
  ▅▅██████▅▅▄▆▅▅▄▁▃▅▄▄▁▃▁▁▁▃▁▃▁▃▄▃▁▆▅▅▄▆▄▁▁▁▁▁▁▁▃▁▃▆▇█████▇▆▆ █
  1.55 ms      Histogram: log(frequency) by time      3.38 ms <

 Memory estimate: 2.74 MiB, allocs estimate: 23500.

Data Migrations between Acsets

While the query execution in the previous section is already quite useful for practical application, it does have one downside, and that is the type of the return object is a DataFrame. While this is appropriate for many cases, one of the benefits of acsets is that, from one point of view, they are in-memory relational databases, and are therefore capable of representing data more complex than can be expressed in a single table. Therefore, it would be nice if one could execute a data migration from one type of acset, where “type” means the schema, to another, that is able to carry along further data we need. For more details on data migration, please see Evan Patterson’s Topos Institute colloquium talk “Categories of diagrams in data migration and computational physics”.

In this context, if the “set” objects (I, J, etc) were further connected to other tables, maybe, say, a list of suppliers, or a list of materials, or even a process graph of downstream work, it would be inconvenient at least, if we lost that relational information during a query. In that case, we’d really want to return another acset on a different schema that is precisely the right shape for what we want to do.

In this simple case, the schema we want is shown below. We’ll think of the object IJKLM as being mapped to the subset of the n-ary product that has those “paths” through the relations we are seeking.

@present IJKLMRelSch(FreeSchema) begin
    (IJKLM,I,J,K,L,M)::Ob
    i::Hom(IJKLM,I)
    j::Hom(IJKLM,J)
    k::Hom(IJKLM,K)
    l::Hom(IJKLM,L)
    m::Hom(IJKLM,M)
end

@acset_type IJKLMRelType(IJKLMRelSch)

Catlab.to_graphviz(IJKLMRelSch, graph_attrs=Dict(:dpi=>"72",:ratio=>"expand",:size=>"3.5"))

Now we formulate the data migration, using AlgebraicJulia/DataMigrations.jl. While we will not be able to rigorously explain data migration here, if one has C-Set (instance of data on schema C) and wants to migrate it to a D-Set, a data migration functor F needs to be specified.

Here, C is our schema IJKLMSch and D is IJKLMRelSch. The functor F is a mapping from D to the category of diagrams on C; formally we denote it F:D\rightarrow \text{Diag}^{\text{op}}(C). Each object in D gets assigned a diagram into C, and morphisms in D get assigned to contravariant morphisms of diagrams.

M = @migration IJKLMRelSch IJKLMSch begin
    IJKLM => @join begin
        ijk::IJK
        jkl::JKL
        klm::KLM
        i::I
        j::J
        k::K
        l::L
        m::M
        IJK_I(ijk) == i
        IJK_J(ijk) == j
        JKL_J(jkl) == j
        IJK_K(ijk) == k
        JKL_K(jkl) == k
        KLM_K(klm) == k
        JKL_L(jkl) == l
        KLM_L(klm) == l
        KLM_M(klm) == m
    end
    I => I
    J => J
    K => K
    L => L
    M => M
    i => i
    j => j
    k => k
    l => l
    m => m
end;

A diagram is itself a functor D:J\rightarrow C, where J is (usually) a small category, and D will point at some instance of the diagram in C. We can plot what the largest diagram looks like, that which object IJKLM in D is mapped to. Note the similarity to the conjunctive query visualized as a UWD previously. In particular, note that “relation” elements must agree upon the relevant “set” elements via their morphisms. The object in C that each object in J corresponds to is given by the text after the colon in the relevant node.

F = functor(M)
to_graphviz(F.ob_map[:IJKLM],node_labels=true)

Because of the simplicity of the schema IJKLMRelSch, the contravariant morphisms of diagrams simply pick out the object in D associated with the source of the morphism. Likewise, the natural transformation part of morphisms of diagrams simply selects for each object its identity morphism.

We run the data migration to move data from the schema IJKLMSch to IJKLMRelSch using the function migrate, and check that the result has the same number of records as other methods.

ijklm_migrate_acset = migrate(IJKLMRelType, ijklm_acs, M)
nparts(ijklm_migrate_acset, :IJKLM) == size(ijklm_query,1)

true

Once again, let’s benchmark. The data migration is slightly slower than the conjunctive query method, but data migrations can express a much richer language of data manipulation than conjunctive queries are capable of.

@benchmark let
    ijklm = migrate(IJKLMRelType, ijklm_acs, M)
    model = JuMP.Model()
    set_silent(model)
    @variable(model, x[parts(ijklm, :IJKLM)] >= 0)
    for i in parts(ijklm, :I)
        @constraint(model, sum(x[incident(ijklm, i, :i)]) >= 0)
    end
end

BenchmarkTools.Trial: 1380 samples with 1 evaluation.
 Range (min … max):  3.212 ms …  10.213 ms  ┊ GC (min … max): 0.00% … 18.52%
 Time  (median):     3.347 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.621 ms ± 654.179 μs  ┊ GC (mean ± σ):  6.28% ± 10.98%

  ▃▇██▇▅▄▂                                  ▁  ▁ ▂▂▁▁         ▁
  █████████▇▅▅▄▅▅▁▁▅▁▄▅▁▅▁▄▁▁▄▁▁▄▁▁▁▄▁▁▁▄▄▅▆██▇████████▆▆█▇▇▆ █
  3.21 ms      Histogram: log(frequency) by time      5.36 ms <

 Memory estimate: 7.33 MiB, allocs estimate: 42773.

To summarize, acsets and tools for working with them provided by AlgebraicJulia can greatly ease development of complex mathematical programming models. Relations, unions, and other concepts from sets are elegantly generalized by category theory, and can help develop more correct, formal, and expressive models. While the optimization “model” examined in this article is fairly abstract, we will investigate model of more practical interest in future installments on the blog.