Adds @compound macro

[LalitChauhan56]

[WIP] Extension to this PR which creates the @compound macro with the updated input.. Also adds the iscompound() function to check if a species is a compound or not.

julia>@variables t
@parameters k
@species C(t) H(t) O(t)

julia>@compound C6H12O2(t) 6C 12H 2O 
3-element Vector{Num}:
  6C(t)
 12H(t)
  2O(t)

julia> iscompound(C6H12O2)
true

[isaacsas]

Can you add a corresponding tests file with some tests?

[TorkelE]

Could you do a test that checks that the components of a compound is the expected one? Like that H2O indeed consist of 2 H and 1 O?

Also, for the test file to be run, you need to add it to the runtests.jl file. You could add it to this section:

    # Runs various miscellaneous tests.
    @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end
    @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end
    @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end
    @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end

Maybe like:

    # Runs various miscellaneous tests.
    @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end
    @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end
    @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end
    @time @safetestset "Compound species" begin include("miscellaneous_tests/compound_macro.jl") end
    @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end

I'd probably also move it to misc tests (?) that it does not necessarily require the DSL.

Otherwise looks good :)

[LalitChauhan56]

From what I understand the CompoundSpecies struct is not being imported when the tests run? I included the compound.jl file using include("compound.jl") and exported the @compound macro along with the iscompound and components functions. Is there something else I am missing?

[isaacsas]

Can you store the component information in a more easy to query manner. You could either have two metadata, CompoundComponent and CompoundCoefficients or have CompoundComponent return a vector of pairs of the form [A => 3, B => 1, C => 2] for A3BC2?

[isaacsas]

@TorkelE and I spoke, and I think we prefer the approach of storing the symbolic components and their coefficients separately.

[LalitChauhan56]

Okay, will implement this in the next commit.

[isaacsas]

You should test equality of the numeric coefficients and the symbolic variables separately. For symbolic variables you can use isequal(A,B) to test if they are the same symbolic. You shouldn't ever be converting to strings when working with symbolic variables except for display reasons.

[isaacsas]

You probably need to use Catalyst.CompoundSpecies. Take a look at the @species macro definition.

[TorkelE]

Ok, so some specifications to this.

1. Storage. For each compound, let's store 1 vector components and vectorcoefficients`. Here, e.g. for

@variables t
@species C(t) H(t) O(t)
@compound C6H12O2(t) 6C 12H 2O 

we would store components as [C, H, O], where C is the species C, and coefficients as [6, 12, 2]. Create a test and ensure that this works:

let
  @variables t
  @species C(t) H(t) O(t)
  @compound C6H12O2(t) 6C 12H 2O 

  @test isequal([C, H, O], components(C6H12O2))
  @test isequal([6, 12, 2], coefficients(C6H12O2))
end

Here, components and coefficients are getter functions that retrieve the corresponding fields.

2. It would be good if the @compound command threw an error if one of the components where not described. E.g. make a test to ensure that

let
  @variables t
  @species C(t) H(t)
  @compound C6H12O2(t) 6C 12H 2O 
end

throws an error.

3. Nested components should work. E.g.

  @variables t
  @species C(t) H(t) O
  @compound OH O H
  @compound C3H5OH3 3C 5H 3OH

  @test isequal(components(components(C3H5OH3)[3]), [O, H])
  @test isequal(coefficients(components(C3H5OH3)[3]), [1, 1])

4. Interpolation. One should be able to do e.g.

@variables t
@species C(t) H(t) O(t)
s = C
@compound C6H12O2(t) 6$s 12H 2O 

and get the same result as

@variables t
@species C(t) H(t) O(t)
@compound C6H12O2(t) 6C 12H 2O 

I will try to write some proper test cases for you to add to your file, which you can use as specifications, when I get back to the office tomorrow.

[LalitChauhan56]

Understood. Will implement these in the following commits. Thank you for the detailed explanation :)

[TorkelE]

Ok, here is a slightly more comprehensive test file you can use to ensure that the functionality works. If there are any oddities, just report it and I will have a look.

using Catalyst, Test

# Test base funcationality in two cases.
let 
    @variables t
    @species C(t) H(t) O(t)
    @compound C6H12O2(t) 6C 12H 2O 
    
    @test iscompound(C6H12O2)
    @test isspecies(C6H12O2)
    @test !iscompound(C)
    @test !iscompound(H)
    @test !iscompound(O)

    @test isequal([C, H, O], components(C6H12O2))
    @test isequal([6, 12, 2], coefficients(C6H12O2))
    @test isequal([C => 6, H => 12, O => 2], component_coefficients(C6H12O2)) # Maybe chaneg "component_coefficients" name soem time.
    @test all(.! iscompound.components(C6H12O2))    
end
let 
    @variables t
    @species O(t)
    @compound O2(t) 2O
    
    @test iscompound(O2)
    @test isspecies(O2)
    @test !iscompound(O)

    @test isequal([O], components(O2))
    @test isequal([2], coefficients(O2))
    @test isequal([O => 2], component_coefficients(O2)) # Maybe chaneg "component_coefficients" name soem time.
    @test all(.! iscompound.components(O2))    
end

# Checks that compounds cannot be created from non-existing species.
let
    @variables t
    @species C(t) H(t)
    @test_throws Exception @compound C6H12O2(t) 6C 12H 2O    # @test_throws might behave weirdly here, if you think this fails as it should but the test comes out starneg, just report and we will figure it out.
end
let
    @variables t
    @test_throws Exception @compound O2(t) 2O    # @test_throws might behave weirdly here, if you think this fails as it should but the test comes out starneg, just report and we will figure it out.
end

# Checks that nested components works as expected.
let 
    @variables t
    @species C(t) H(t) O
    @compound OH O H
    @compound C3H5OH3 3C 5H 3OH

    @test !iscompound(O)
    @test !iscompound(H)
    @test iscompound(OH)
    @test iscompound(C3H5OH3)
    
    @test !iscompound(components(C3H5OH3)[1])
    @test !iscompound(components(C3H5OH3)[2])
    @test iscompound(components(C3H5OH3)[3])

    @test isequal([C, H, OH], components(C3H5OH3))
    @test isequal([O, H], components(components(C3H5OH3)[3]))
    @test isequal([3, 5, 3], coefficients(C3H5OH3))
end

# Checks that interpolation works.
let 
    @variables t
    @species C(t) H(t) O(t)
    s = C
    @compound C6H12O2_1(t) 6$s 12H 2O 
    @compound C6H12O2_2(t) 6C 12H 2O 

    @test isequal(components(C6H12O2_1), components(C6H12O2_2))
    @test isequal(coefficients(C6H12O2_1), coefficients(C6H12O2_2))
end
let 
    @variables t
    @species C(t) H(t)
    @compound Cyclopentadiene 5C 6H
    C5H6 = Cyclopentadiene
    @compound C10H12 2$(C5H6)

    @test iscompound(C10H12)
    @test iscompound(components(C10H12)[1])
    
    @test isequal(components(C10H12)[1], C5H6)
    @test isequal(components(C10H12)[1], Cyclopentadiene)
    @test isequal(coefficients(C10H12)[1], 2)
end

[LalitChauhan56]

Thank you! Will test this and let you know.

[LalitChauhan56]

After numerous attempts of trying to store the components and coefficients separately within the macro itself, I ran into hygiene and scoping issues with the macro. As a workaround, I switched back to the previous approach and instead made changes to the retrieval functions components and coefficients themselves. They now retrieve the metadata array from the compound and returns the needed values. All the tests pass now for the components, coefficients and iscompound and nested functions.

Also, the @test all(.! iscompound.components(C6H12O2)) type lines in the test file gave error iscompound has no field components (also added the error as comments in the test file) so instead I switched it with @test all(!iscompound(i) for i in components(C6H12O2)) which does the same thing (correct me if I'm wrong).

Also, the interpolation can be implemented with this approach without the use of $.

julia>  @variables t
    @species C(t) H(t)
    @compound Cyclopentadiene(t) 5C 6H
    C5H6 = Cyclopentadiene
    @compound C10H12(t) 2C5H6
    
  julia> iscompound(C10H12)
true

julia> iscompound(components(C10H12)[1])
true

julia> components(C10H12)
1-element Vector{Num}:
 Cyclopentadiene(t)
 
 julia> coefficients(C10H12)
1-element Vector{Int64}:
 2

[LalitChauhan56]

Removed $ from the tests as mentioned above. All tests pass now.
@isaacsas @TorkelE

[LalitChauhan56]

Also, how can I fix the format check tests? :D

[ErikQQY]

To fix the format check, you will need to run using JuliaFormatter, Catalyst; format(joinpath(dirname(pathof(Catalyst)), "..")) in the package developing environment before commit.

[isaacsas]

Don't reformat in this PR. The style recently changed, so it will polute the PR by modifying like every Catalyst file...

There should be a new style update in the near future that will also require a complete repo reformat, and we can rerun formatting at that time.

[LalitChauhan56]

From what I understand GCD is actually creating issues in some cases. For example -
When balancing the reaction Reaction(k,[CO2,H2O],[C6H12O6,O2])
the final stoich values come out to be

4-element Vector{Float64}:
 1.0
 1.0
 0.16666666666666666
 1.0

I am using this as a source to get the answers. Here too, x3 comes out to be 1/6 = 0.16666.... which is rounded off to 0.2, so I think some amount of rounding will be required.

[TorkelE]

Are you adding balance of equations and such stuff here as well?

Generally, it is better to restrict 1 PR to 1 function. E.g. do one PR to add the @compound functionality, and then another for balancing equations and what not.

[isaacsas]

@LalitChauhan56 will try to take a closer look in the next day or two, but shouldn't there be a way to do this that avoids floating point? Also, is your back substitution step really doing backsub, or is it doing a floating point least squares solve? (Is the matrix non-square, if so I believe it is likely doing least squares?)

[isaacsas]

Also, @TorkelE is right that it probably makes sense to do the balancing as a separate PR so we can merge the compound stuff now.

[LalitChauhan56]

Understood, I'll create a separate PR for the balancing part.

Also, yes there should be a way to avoid floating point numbers. I'll try to figure it out.

[isaacsas]

BTW, it would be fine to drop into rational numbers if really needed, though probably slow. It would just be nice to avoid floating point if possible. (That said, if we need to have the first version use floating point and can then refine it over time to be better that is ok too.)

I'm not up on the latest methods for balancing, so it may be there is no "standard" integer-based approach. I'd have to look at the literature.

[LalitChauhan56]

@isaacsas This commit uses rational numbers and skips floating point numbers entirely.

Also, once you have reviewed this I can delete the balancing code from this PR so that the @compound part can be merged and then I can make a new PR for the balancing part as that depends on the @compound part to work.

[isaacsas]

Why don’t you move the balancing stuff out of this PR now so we can get this merged, and then we can review it in the new PR you open. That’ll be easier for us all I think. Thanks!

[isaacsas]

It is ok if the new PR has the code for this PR too; once this is merged you can just merge master into your new PR and the redundant code will get dropped.

[LalitChauhan56]

@isaacsas I removed the balancing code from this PR. Also can you please take a look at the formatting test failing? I formatted the file using Julia Formatter 1.0.32 but it still fails.

[isaacsas]

You are exporting non-existent functions here? (i.e. balance and get_stoich?) Also, let's not export component_coefficients.

[isaacsas]

What are these commented lines about?

[isaacsas]

There don't seem to be any tests allowing using a parameter or its value as a coefficient?

i.e.

alpha = 2
@compound 2H $alpha*H

or

@parameter alpha = 2
@compound 2H $alpha*H

or something similar (dropping the $ if you want).

Similarly are there tests that allow using a variable like

@species A(t)
B = A
@compound A2 2*B   # or @compound 2*$B

[isaacsas]

Basically, I want to ensure someone can use this programmatically where coefficients and/or species may be stored in a variable.

[LalitChauhan56]

I added more tests for the functionalities you suggested.

Also, regarding this issue, I am unable to make this work. In cases where there is a coefficient in the value of the variable itself such as S = 2C , in the macro I am left with two arguments. First, which is the original coefficient defined when using the @compound macro and the second is the one which is present in the value of the variable (2 in this case). To simplify this I tried using the simplify() function but it does not help and eval() runs into scoping issues because the macro cannot access the previously defined species/variables ('C' is not defined in this case). esc() also does not work as I already need to escape the expressions to set the metadata for components and coefficients and I cannot escape the same variable twice. Is there a way to make it work? If yes, then should I make the changes here itself or open a new PR later once this is merged as this enables basic functionality for @compound?

Letting it work the way it currently does leaves me incorrect components and coefficients. For example -

  @variables t
  @species C(t) H(t) O(t)
  s = 2C
  @compound C6H12O2_1(t) 6s 12H 2O
    
  julia> components(C6H12O2_1)
  3-element Vector{Num}:
   2C(t)  #Should be C
    H(t)
    O(t)

  julia> coefficients(C6H12O2_1)
  3-element Vector{Int64}:
    6 #Should be 2*6 = 12
   12
    2

[LalitChauhan56]

However cases such as S = C and the ones you have mentioned above work perfectly fine.

[isaacsas]

That sounds ok. I think it is fine if the right most term in the product has to be just a variable and not a general expression. We just want to make sure the coefficients can be expressions (and eventually document all this behavior).

[isaacsas]

Looks good, we can clean up later if needed. Thanks, nice work on this PR!!!

[isaacsas]

I'm not sure about the formatter issue. We can figure it out later.

I'll be traveling the next two weeks, but try to look at your balancing PR when I get some time. @TorkelE said he will also try to work with you on any issues and reviewing that PR.

The next step is for you to now update that PR to master (either rebase master onto it, or merge master into it -- I usually do the latter).

[LalitChauhan56]

Understood. Thank you for guiding me throughout the PR :)

Also, I have updated the balancing PR.