Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interop with SymbolicUtils #64

Merged
merged 35 commits into from
Oct 11, 2024
Merged

Interop with SymbolicUtils #64

merged 35 commits into from
Oct 11, 2024

Conversation

olynch
Copy link
Member

@olynch olynch commented Aug 14, 2024

This is a more-or-less feature-complete conversion between DecaExpr and an equivalent data structure built on SymbolicUtils, but it requires some testing/naming bikeshedding.

@lukem12345
Copy link
Member

@GeorgeR227 Looks like this should get rebased onto #54 ; there is apparently some redundancy in dicts for typing/ untyping symbols.

src/ThDEC.jl Outdated
Comment on lines 140 to 151
@nospecialize
function ★(s::Sort)
@match s begin
Scalar() => throw(SortError("Cannot take Hodge star of a scalar"))
Form(i, isdual) => Form(2 - i, !isdual)
end
end

function Base.nameof(::typeof(★), s)
inv = isdual(s) ? "⁻¹" : ""
Symbol("★$(as_sub(isdual(s) ? 2 - dim(s) : dim(s)))$(inv)")
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we need a deeper consideration of the dimension that our physics in occurring in here. I've noticed that the current code is mostly hard-coded for 2D and in the case of the Hodge star, the resulting operator would be wrong if we were working in 1D.

I can say that the way our current inference/overloading rules work is by grouping them into contexts based on the dimension, so we have separate rules to apply in the 1D case and in the 2D case. The way the dimension is figured out is by the user explicitly passing us that information, since some physics can work without modification in multiple dimensions, like the Heat equation.

Copy link
Member

@quffaro quffaro Aug 16, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think one way to work around this is to store the Space as a parameter in the Form sort. Space may be a wrapper for the dimension of Space, but can be embiggened with other attributes down the line, with the provision we have to update our signature checking.

I think this makes philosophical sense. n-forms are the same in any (n+k)-dimensional space, but since they are also determined by duality, their dual varies according to the dimensionality. So this means every n-form is determined by a different dual per dimension. Tagging this dimension Form(dim, isdual, space) allows us to track this.

@lukem12345 what do you think about this philosophically and practically?

Copy link
Member

@lukem12345 lukem12345 Aug 16, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Philosophically, here are some probing questions: How do you want to distinguish between the Heat equation in 1D, 2D, and 3D? What do morphisms between diagrams (of possibly differing dimension) look like?

Practically: the entire PR seems to be hard-coded for the 2D Discrete Exterior Calculus over a double de Rham complex. I believe more consideration needs to be taken as to how we reconcile having 3 similar algebraic theories in the same package. Is there any advantage to using parameters here besides code reduction? What is the scale of that code reduction? And is the generalization worth the mental overhead? Are we going to be auto-generating rules by doing index arithmetic on these parameters or not? Do we want it to be clear when we are programming in a 1D environment vs a 3D environment? How does this impact representations with respect to the multi-dimensional multiphysics project? Do we want to represent the various DECs over a single de Rham complex?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This all makes sense. This was supposed to be a first stab at what the API would look like. I think the next step is to add the space to the types, and part of the information about the space will be its dimension.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When we convert a DecaExpr to a DecaSymbolic, you should have to also provide a mapping from space symbol to dimension.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What would you all think of something like:

struct Space
    name::Symbol
    dim::Int
end

@data Sort begin
    Scalar()
    Form(dim::Int, isdual::Bool, space::Space)
    VField(isdual::Bool, space::Space)
end

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jpfairbanks and I discussed this exact suggestion earlier today. I'd like to implement it, but I think @lukem12345 's practical questions still remain open. @lukem12345 @GeorgeR227 for clarification, should these questions + considerations should be resolved in this PR?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nothing is immediately keeping this PR open as long as tests pass.

@GeorgeR227
Copy link
Contributor

@GeorgeR227 Looks like this should get rebased onto #54 ; there is apparently some redundancy in dicts for typing/ untyping symbols.

I'm not sure because in this case OPERATOR_LOOKUP is taking typed names down to generic names. In that PR we take those typed names down to a single representative which is still typed. It's true that they're similar but their purposes feel different.

@lukem12345
Copy link
Member

Yes the rebase is for purposes of grouping together code that accomplishes similar tasks.

@olynch olynch mentioned this pull request Aug 19, 2024
quffaro and others added 3 commits August 19, 2024 14:20
1. add show methods for DecaEquation and DecaSymbolic
2. remove spurious exports in code gen
3. fix binary form constructors from pre-space code
4. remove aqua export check because of code gen
5. add Klausmeier tests
s = ThDEC.$binop(Sort(T1), Sort(T2))
SymbolicUtils.Term{Number(s)}(ThDEC.$binop, [v, w])
end
export $binop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why do you need to export these 3 times for the different nospecializes? I replaced them with 1 export and tests pass

src/ThDEC.jl Outdated
struct SpaceLookup
default::Space
named::Dict{Symbol,Space}
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should add a default constructor that just takes the default name.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Realized this needs to be default name and dimension

src/ThDEC.jl Outdated
:DualForm0 => Form(0, true, space)
:DualForm1 => Form(1, true, space)
:DualForm2 => Form(2, true, space)
:Constant => Scalar()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

name could also be a Parameter or a Literal. Parameters are "constants" that have values that depend on time so in the solver the t variable has to get threaded to them.

src/ThDEC.jl Outdated

# TODO
function Base.nameof(::typeof(♯), s)
Symbol("♯s")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should there be a $ before s here?

src/ThDEC.jl Outdated

# TODO
function Base.nameof(::typeof(♭), s)
Symbol("♭s")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should there be a $ before s here?


# a struct carry the symbolic variables and their equations
struct DecaSymbolic
vars::Vector{Symbolic}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we call this context?

# App1(f, App1(...)
AppCirc1(fs, arg) => foldr(
# panics with constants like :k
# see test/language.jl
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is actually a huge problem. In applications, we always have to add some arbitrary nonlinear functions at model specification time.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also fits into best practice of "being liberal in what we accept" https://en.wikipedia.org/wiki/Robustness_principle

src/decasymbolic.jl Outdated Show resolved Hide resolved
SymbolicUtils.Term{Number(s)}(ThDEC.$unop, [v])
end
end
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this loop should be a function that users can eval to register new symbolic functions. We almost always need at least 1 arbitrary functions that users should be able to register with a syntax like

@register f(x::PrimalForm{1})::PrimalForm{1}

It would also be cool if they could give the implementation there too.

@register f(x::PrimalForm{1})::PrimalForm{1} = map(x) do val
  val  ^ 2
end

Copy link
Member

@quffaro quffaro Aug 23, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this comment was moved to the end of the thread for record keeping purposes.

src/decasymbolic.jl Outdated Show resolved Hide resolved
@jpfairbanks
Copy link
Member

jpfairbanks commented Aug 21, 2024

Instead of creating julia functions for the operators of the DEC, could we just create SymbolicUtils function types?

https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/df8e84d31e60f4dc70f0df2ad88ede91f6f3fd27/src/types.jl#L913

That is what you get when you do:

julia> @syms f(x)::Complex
(f(::Number)::Complex,)

You can construct function types with SymbolicUtils.Sym{SymbolicUtils.FnType{Tuple{Real,Real}, Complex}}(:h) so these are what should come out of registering a new function type and be used to represent the basic DEC operators. That way any code that expects function types will work for our DEC operators.

Copy link

codecov bot commented Aug 23, 2024

Codecov Report

Attention: Patch coverage is 86.77966% with 39 lines in your changes missing coverage. Please review.

Project coverage is 86.58%. Comparing base (55d2972) to head (df41b18).
Report is 8 commits behind head on main.

Files with missing lines Patch % Lines
src/deca/ThDEC.jl 72.72% 24 Missing ⚠️
src/SymbolicUtilsInterop.jl 86.11% 10 Missing ⚠️
src/symbolictheoryutils.jl 93.33% 3 Missing ⚠️
src/graph_traversal.jl 94.28% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main      #64      +/-   ##
==========================================
- Coverage   86.80%   86.58%   -0.22%     
==========================================
  Files          13       18       +5     
  Lines         894     1170     +276     
==========================================
+ Hits          776     1013     +237     
- Misses        118      157      +39     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@lukem12345
Copy link
Member

lukem12345 commented Aug 23, 2024

An implication of Decapode type-inference that seems baked-into this PR is that it is assumed that the choice of primal/ dual is completely inferable based on the types of the parents. However certain operations such as the musical isomorphisms or wedge products can either return a primal or dual form depending on the discrete operator chosen. Note also that for any Decapode, there are multiple means of specifying whether forms are primal or dual. Performing type inference via breadth-first or depth-first traversal at parse-time biases the system towards only representing certain discretizations. In certain cases, choosing a particular discretization for a certain variable or operation may constrain successor operations to certain discretizations, for which CombinatorialSpaces may not yet offer a method. And so a satisfiable model may not be returned, although one may exist.

The appearance of dual complexes leads to a proliferation of the operators in the discrete theory. For example there are
primal-primal, primal-dual etc. versions of many operators. This is of course unique to the discrete side.

  • Hirani, Discrete Exterior Calculus

Resolving this issue by allowing the Decapode macro to represent any formal combination of terms should simplify parsing logic. As it stands now, some algebraic rules are implicit in the code that parses the body of the macro, and other algebraic laws are (will be) encoded as rewrite rules.

@jpfairbanks
Copy link
Member

We can always do dispatch on return type to let you write sharp(Dual, x)

@lukem12345
Copy link
Member

lukem12345 commented Aug 23, 2024

Yes but that is still local. The issue I described is not with choosing which functions to emit once types have been decided. It is an issue with deciding types.

@jpfairbanks
Copy link
Member

I'm moving a comment from @quffaro from inline reply to top level comment so that it stays here after we change the code and mark that thread resolved:

Tests currently pass with exception to klausmeier. To get it working, @jpfairbanks suggested a @register macro to add function terms. However I'm not sure if I'm understanding completely, so I'm seeking clarification.

Currently (in this PR), for a given operator we define two methods:

  1. a Julia function over ThDEC for checking its type. They are aggregated in a OPERATOR_LOOKUP const. This step is analogous to SymbolicUtil's promote_symtype
  2. a Julia function over DecaSymbolic, our type for the SymbolicUtils implementation of ThDEC, which calls (1) for validation

we use the OPERATOR_LOOKUP to also check if a model has admissible terms (e.g., Δ), so we want the ability to append the OPERATOR_LOOKUP by @registering new functions.

@register foo(x::Scalar)::Scalar = x + 1
# codegens
foo(s::Scalar) = Scalar
foo(x::BasicSymbolic{T<:DECType}) = begin
    s = foo(Sort(T)) # TYPE CHECKING
    Sym{FnType{Tuple{Number(s)}, Real}(foo)
end

Another example comes with a hazard: PrimalForm{T} is a parameterized type, but PrimalForm::Sort is a function. If we want to continue typechecking with ThDEC types, do our current way of checking signatures, we need to make sure codegen for typechecking supports the engineering task to do this.

@register foo(x::PrimalFormT{1})::PrimalFormT{1}
# codegens
foo(x::Form) = @match x begin
    PrimalForm(1) => PrimalForm(1)
    _ => error
end

The steps:

Suppose we wanted to add foo to the OPERATOR_LOOKUP dictionary. Notice PrimalFormT{1} is the DecaSymbolic type, which accepts parameters.

@register foo(x::PrimalFormT{1})::PrimalFormT{1} = map(x) do val
  val  ^ 2
end

We need to define its (1) type-checking (src/ThDEC), (2) its DecaSymbolic implementaton (src/decasymbolic), and possibly (3) its implementation.

We define the implementation of a function

func(x::PrimalFormT{1})::PrimalFormT{1} = map(x) do val; val ^ 2 end

❓ Where does this live? Does this act on BasicSymbolic terms?

Next, we define it symbolically:

symfunc = Sym{FnType{Args, Real}}(ThDEC.func)

Then, we add it to OPERATOR_LOOKUP (possible haskey validations omitted)

ThDEC.OPERATOR_LOOKUP[op] = func

❓ Currently OPERATOR_LOOKUP is typed as Dict{Symbol, Function}. Is the intent to change Dict{Symbol, SymbolicFunction} and do typechecking by replacing the # TYPE CHECKING (above) line with promote_symtype?

@jpfairbanks
Copy link
Member

We define the implementation of a function

func(x::PrimalFormT{1})::PrimalFormT{1} = map(x) do val; val ^ 2 end

❓ Where does this live? Does this act on BasicSymbolic terms?

I am rethinking my idea to define the implementation here. The implementations will often need to refer to parts of the mesh, including vertex coordinates or volumes of a simplex. So we should reserve this syntax for defining a new function symbol and a definition in terms of more primitive functions that could get turned into a rewrite. So that you could register a function like

f, r = @register f(v::PrimalFormT{1},x::PrimalFormT{0})::PrimalFormT{1} = d(Δ(x)) + (v   x)

and have that expand out to

f, = @syms f(v::PrimalFormT{1},x::PrimalFormT{0})::PrimalFormT{1}
r = @rule f(~v::PrimalFormT{1},~x::PrimalFormT{0})::PrimalFormT{1} => d(Δ(~x)) + (~v   ~x)

Next, we define it symbolically:

symfunc = Sym{FnType{Args, Real}}(ThDEC.func)

Then, we add it to OPERATOR_LOOKUP (possible haskey validations omitted)

ThDEC.OPERATOR_LOOKUP[op] = func

❓ Currently OPERATOR_LOOKUP is typed as Dict{Symbol, Function}. Is the intent to change Dict{Symbol, SymbolicFunction} and do typechecking by replacing the # TYPE CHECKING (above) line with promote_symtype?

I think that we can create a function for each DEC operator, then create a symbolic function with that function as the head.

function Δ end
Δ, = @syms Δ(x::PrimalFormT{0})::PrimalFormT{0}

If the function symbol we are creating needs to do dispatch, then it should overload promote_symtype to hook that logic into SymbolicUtils the way that it is expected.

@lukem12345
Copy link
Member

#35

I think using a defintion operator := like that mentioned in this old issue would be nice.

jpfairbanks added a commit that referenced this pull request Aug 26, 2024
This was referenced Aug 30, 2024
@nospecialize
function $(esc(f))($(newargs...)) where $(constraints...)
s = $(esc(f))($(innerargs...))
SymbolicUtils.Term{Number(s)}($(esc(f)), [$(getindex.(binding, 1)...)])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I macro expanded a call to @register and saw a bunch of gensim stuff. I don't think we need to be as aggressive in esc in here.

@@ -6,6 +6,7 @@ import SymbolicUtils: promote_symtype
function rules end
export rules

# TODO: Probable piracy here
function promote_symtype(f::ComposedFunction, args)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should upstream this method JuliaSymbolics/SymbolicUtils.jl#650

GeorgeR227 and others added 4 commits October 3, 2024 11:15
* Set version to 0.1.7

* Added more exports (#44)

Added `apex` and `@relation`, `to_graphviz` from Catlab

Co-authored-by: James <[email protected]>

* Add type rules for vectorfields

* Add musical overload resolution

* Take advantage of :infer in type rules

* Initial attempt at rewriting

Converts ACSet to a series of Symbolic terms that can be rewritten with a provided rewriter

* Added proof of concept

Added a short script showcasing how rewriting could be done with the `Sort` types and a reference ACSet.

* Added ability to do through-op rewrites

This now supports the ability for ACSet intermediate expressions to be merged into one single expression upon which rewriting rules (like dd=0) may be performed.

* Added Space import

* Completed full pipeline

Can take ACSets to Symbolics back to  ACSets

* Remove metadata usage

This needs to switch to use the new type system

* Added DECQuantity types

Also switched to using SymbolicsUtils' `substitute`. Still needs tests and code needs to be cleaned up.

* Completed pipeline again

Addition now works as well but rewriting seems to be janky, unrelated to this pipeline specifically I believe.

* fixed bug where type-checking subtraction uses +(S1,S2), which is obsolete

* George and I debugged rewriting. Incorrect type passed to resulting term meant typed rewriting would fail

* Cleaning up pipeline

This black boxes the intermediate symbolic expressions to the user. The user will simply submit a rewriter that will then be applied

* Fixed order of inclusions

* adding support for Parameters and Constants

* Added tests for acset2symbolic

* etc

* Literals testing

* parameters test passing after some debugging.

* supporting Infer, better Base.nameof, better tests

* Clean out-of-order vector constructions

* Convert to symbolics inside merge_equations

* Reduce cases of topological sort

* Reify via recursive function, not lambda case

* Further improvement of acset2symbolics

Remove special DerivOp handling, fixed bug where multiple equations with the same variable result were being dropped, more tests to cover these cases and further clean up.

* Remove extraneous tangents

* Remove redundant helper functions

* Pass indexed names and types directly

* Removed extraneous d arg

* fixing work on tumor invasion

* macros which create export stmts will fail inside @testset due to JuliaLang issue #51325

* removed ghost emoji and added convenience function for rules. aqua's failing persistent tasks.

* Added more tests for acset2symbolics

* Fixed persistence issue

Also set default form dim to 2 and allowed it to vary.

* Final touches

* Remove unused fuctionality

---------

Co-authored-by: AlgebraicJulia Bot <[email protected]>
Co-authored-by: James <[email protected]>
Co-authored-by: Luke Morris <[email protected]>
Co-authored-by: Matt <[email protected]>
These tests attempt to cover all of the basic functionality of ThDEC. It also points out some strange or desired behavior.
Comment on lines 79 to 114
@testset "Nameof" begin
@test nameof(symtype(c)) == :Constant
@test nameof(symtype(t)) == :Parameter
@test nameof(symtype(a)) == :Scalar
@test nameof(symtype(ℓ)) == :Literal

@test nameof(symtype(u)) == :Form0
@test nameof(symtype(ω)) == :Form1
@test nameof(symtype(η)) == :DualForm2

# TODO: Do we want this style of typed subtraction?
@test nameof(-, symtype(u), symtype(u)) == Symbol("₀-₀")
# TODO: This breaks since this expects the types to have a `dim` function
@test_broken nameof(-, symtype(a), symtype(b)) == Symbol("-")

@test nameof(∧, symtype(u), symtype(u)) == Symbol("∧₀₀")
@test nameof(∧, symtype(u), symtype(ω)) == Symbol("∧₀₁")
@test nameof(∧, symtype(ω), symtype(u)) == Symbol("∧₁₀")

# TODO: Do we need a special designation for wedges with duals in them?
@test nameof(∧, symtype(ω), symtype(η)) == Symbol("∧₁₂")

# TODO: Why is this being named as such?
@test nameof(∂ₜ, symtype(u)) == Symbol("∂ₜ(Form0)")
@test nameof(∂ₜ, symtype(d(u))) == Symbol("∂ₜ(Form1)")

@test nameof(d, symtype(u)) == Symbol("d₀")
@test_broken nameof(d, symtype(η)) == Symbol("dual_d₂")

@test_broken nameof(Δ, symtype(u)) == Symbol("Δ₀")
@test_broken nameof(Δ, symtype(ω)) == Symbol("Δ₁")

@test nameof(★, symtype(u)) == Symbol("★₀")
@test nameof(★, symtype(ω)) == Symbol("★₁")
@test nameof(★, symtype(η)) == Symbol("★₀⁻¹")
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lukem12345 I noticed while going through the Base.nameof functionality that support for names of operators with dual inputs is not well supported. Could you advise on how we should be naming these operators to differentiate them from our common primal-primal operators?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think most notable here is the wedge and the Laplacian.

Copy link
Member

@jpfairbanks jpfairbanks Oct 4, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that either dual_ or \~ would be appropriate.

jpfairbanks and others added 5 commits October 4, 2024 16:12
Discovered that the `Base.nameof` functions require logic to support InferredTypes. If any of the args is an InferredType, then we should simply return the generic operators, e.g. d or Δ.
Added a fix to the acset2symbolics code that helps fix the issue caused there by #77. A full fix to this issue can restore the code to just use infer_terminal_names
We now better deal with being passed infer types and added support for dual exterior derivatives.
@jpfairbanks jpfairbanks merged commit 2f82754 into main Oct 11, 2024
7 of 8 checks passed
@jpfairbanks jpfairbanks deleted the symbolicutilsinterop branch October 11, 2024 17:15
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants