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

dot-product, inner-product, outer-product, and many more structure functions #193

Merged
merged 13 commits into from
Dec 11, 2020

Conversation

sritchie
Copy link
Member

@sritchie sritchie commented Dec 8, 2020

Our structure operations have a few incompatibilities with the refman / scmutils implementation that I thought we might get dialed before release.

As of this PR, we have some new generic operations:

  • g/dot-product, for scalars, differentials, structures, functions and
    row/column matrices
  • g/inner-product for scalars, structures, functions and row/column matrices
  • g/outer-product for functions, structures of length 3 and matrices, between
    a row and a column only
  • g/cross-product now works for row/column matrices in addition to structures
    (and functions that accept these)

Structures got a lot of love too:

  • Structures and matrices both gain the ability to do native get-in,
    assoc-in and empty. These work as expected, like a potentially nested
    vector. (#193)

  • matrix.cljc gains up->row-matrix, up->column-matrix, row-matrix->up,
    column-matrix->up
    (#193)

  • structure.cljc gains many features in
    (#193):

    • kronecker and basis-unit for generating potentially infinite basis
      vectors
    • the ability to conj new items onto a structure: (conj (up 1 2) 3) => (up 1 2 3)
    • The structure-preserving map-chain takes a 2-arg function and presents it
      with each element of a deeply nested structure, along with a vector of its
      "chain", the path into its location. The fn's return becomes the new item at
      that location.
    • structure->prototype generates a same-shape structure as its argument,
      with symbolic entries that display their location.
    • typical-object returns a structure of the same shape and orientation as
      s, generated by substituting gensymmed symbols in for each entry.
    • compatible-zero returns a structure compatible for multiplication with s
      down to 0.
    • transpose-outer returns a new structure with the same orientation as the
      first element of s, filled with elements of the same orientation as s.
      Each element is generating by taking the first element of each entry in s,
      the the second, etc... In that sense this is similar to a traditional matrix
      transpose.
    • dot-product takes the dot product of two structures. They must be the same
      top-level orientation and dimension; beyond that, their entries are
      pairwise-multiplied and summed.
    • inner-product is the same, but the left structure is conjugated first.
    • outer-product now works multiple levels deep.
    • vector-outer-product and vector-inner-product are similar, but only
      enforce the top-level length; all internal structures are NOT flattened and
      must be compatible for g/*.
    • compatible-for-contraction? now searches recursively down into a
      structure; previously it only checked the top level.
    • The new *allow-incompatible-multiplication* dynamic variable is set to
      true by default. Set it false to force a setting where, when you multiply
      structures, they must be:
      • opposite orientation
      • every element of the right entry must be compatible for contraction with
        the left
    • structure multiplication with scalars, etc now respects ordering, just in
      case any multiplication is not commutative.

@codecov-io
Copy link

codecov-io commented Dec 9, 2020

Codecov Report

Merging #193 (b1d87e4) into master (1aa8154) will decrease coverage by 0.18%.
The diff coverage is 84.27%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #193      +/-   ##
==========================================
- Coverage   84.94%   84.76%   -0.19%     
==========================================
  Files          74       74              
  Lines        7042     7206     +164     
  Branches      386      401      +15     
==========================================
+ Hits         5982     6108     +126     
- Misses        674      697      +23     
- Partials      386      401      +15     
Impacted Files Coverage Δ
src/sicmutils/env.cljc 89.28% <ø> (ø)
src/sicmutils/generic.cljc 76.80% <0.00%> (-2.04%) ⬇️
src/sicmutils/matrix.cljc 79.65% <80.95%> (+0.12%) ⬆️
src/sicmutils/structure.cljc 84.76% <87.80%> (-2.51%) ⬇️
src/sicmutils/calculus/derivative.cljc 92.78% <100.00%> (+0.07%) ⬆️
src/sicmutils/examples/central_potential.cljc 66.12% <100.00%> (+2.33%) ⬆️
src/sicmutils/function.cljc 89.70% <100.00%> (+0.15%) ⬆️
src/sicmutils/numsymb.cljc 92.35% <100.00%> (-1.28%) ⬇️
src/sicmutils/numerical/unimin/golden.cljc 91.89% <0.00%> (-2.71%) ⬇️
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1aa8154...641ba7d. Read the comment docs.

@@ -382,11 +382,10 @@
(fn [i v_i]
(structural-derivative
(fn [w]
(g (struct/structure-assoc-in v [i] w)))
Copy link
Member Author

Choose a reason for hiding this comment

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

woohoo, I figured out the interfaces we needed to make the native thing work.

Copy link
Collaborator

Choose a reason for hiding this comment

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

nice. It is a noticeable improvement.

@@ -475,6 +476,8 @@
(define-unary-operation g/square #(diff-* % %))
(define-unary-operation g/cube #(diff-* % (diff-* % %)))

(define-binary-operation g/dot-product diff-*)
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 is the only one that scmutils implements - differential feels like it should implement conjugate and inner product, yeah?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I agree. scm punts on conjugate. inner-product is not exposed in the version I'm using. dot-product and outer-product work, so we should do all the ones we can.

src/sicmutils/matrix.cljc Show resolved Hide resolved
src/sicmutils/matrix.cljc Show resolved Hide resolved
src/sicmutils/matrix.cljc Show resolved Hide resolved
@@ -85,7 +85,7 @@
" var _000e = Math.cos(_000a);\n"
" var _0010 = Math.sin(_000a);\n"
" var _0011 = Math.pow(_0010, 2);\n"
" return [1, [thetadot, phidot], [(- l1 * m2 * _0006 * _0010 * _000e - l2 * m2 * _0005 * _0010 + g * m2 * _000e * _0008 - g * m1 * _0009 - g * m2 * _0009) / (l1 * m2 * _0011 + l1 * m1), (l2 * m2 * _0005 * _0010 * _000e + l1 * m1 * _0006 * _0010 + l1 * m2 * _0006 * _0010 + g * m1 * _0009 * _000e + g * m2 * _0009 * _000e - g * m1 * _0008 - g * m2 * _0008) / (l2 * m2 * _0011 + l2 * m1)]];\n"
Copy link
Member Author

Choose a reason for hiding this comment

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

that structure change caused some of these tests to switch variable order.

test/sicmutils/generators.cljc Show resolved Hide resolved
(is (= (if (v/one? x) 1 x)
(g/* x 1)) "First unity gets returned.")
(is (= x (g/* 1 x)) "Anything times a 1 returns itself.")))
(checking "g/*" 100 [x gen/any]
Copy link
Member Author

Choose a reason for hiding this comment

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

v/= replaces =, for the rarish case of (= 0 0.0) turning up false.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Huh! I forgot about that. (Clojure also has == that does this).

Copy link
Member Author

Choose a reason for hiding this comment

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

re: ==: I think that can ONLY do numbers, so when I generate via gen/any it solves our problem but breaks on non-numeric things.

test/sicmutils/structure_test.cljc Show resolved Hide resolved
@sritchie
Copy link
Member Author

sritchie commented Dec 9, 2020

@littleredcomputer , I ended up spending another day fleshing out the structure.cljc namespace with ~everything that GJS has in scmutils (no abstract structures yet).

This uncovered a few bugs with the way we were handling deeply nested structures. It didn't matter as much (except the double pendulum equations of motion had an issue that you'll see), but the behavior of the outer-product and inner-product that we had named in structure.cljc did not match what we want to see out of the tests.

I think this brings us up to as much parity with the refman as I could hope for.

The PR is pretty big because of a ton of tests. I am happy to split this up - most of the reading is in structure.cljc and its tests, and that is all fairly smooth to follow.

Would love feedback here, since I know you've been deep in the Structure Mines and have thoughts here!

@sritchie sritchie changed the title new structure, matrix operations [in progress] dot-product, inner-product, outer-product, and many more structure functions Dec 9, 2020
@sritchie sritchie force-pushed the sritchie/tidying_before_release branch from 641ba7d to 3d2a3c7 Compare December 9, 2020 07:17
@@ -475,6 +476,8 @@
(define-unary-operation g/square #(diff-* % %))
(define-unary-operation g/cube #(diff-* % (diff-* % %)))

(define-binary-operation g/dot-product diff-*)
Copy link
Collaborator

Choose a reason for hiding this comment

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

I agree. scm punts on conjugate. inner-product is not exposed in the version I'm using. dot-product and outer-product work, so we should do all the ones we can.

(pairs))]
(reduce - 0 (map (fn [[[m1 p1] [m2 p2]]]
(/ (* m1 m2)
(sqrt (square (- p1 p2)))))
Copy link
Collaborator

Choose a reason for hiding this comment

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

You could use abs here

Copy link
Member Author

Choose a reason for hiding this comment

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

will do in my followup!

src/sicmutils/matrix.cljc Show resolved Hide resolved
(is (= (if (v/one? x) 1 x)
(g/* x 1)) "First unity gets returned.")
(is (= x (g/* 1 x)) "Anything times a 1 returns itself.")))
(checking "g/*" 100 [x gen/any]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Huh! I forgot about that. (Clojure also has == that does this).

@@ -382,11 +382,10 @@
(fn [i v_i]
(structural-derivative
(fn [w]
(g (struct/structure-assoc-in v [i] w)))
Copy link
Collaborator

Choose a reason for hiding this comment

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

nice. It is a noticeable improvement.

If `n` is not supplied returns an infinite sequence."
([i] (map (partial kronecker i)
(range)))
([n i] (into [] (take n) (basis-unit i))))
Copy link
Collaborator

Choose a reason for hiding this comment

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

why not let them both return a seq

Copy link
Member Author

Choose a reason for hiding this comment

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

good idea, I for some reason thought this was the defensive way to go, but it looks like any place this is used makes a vector if it needs to.

:else false)
:else false))

(defn structure?
Copy link
Collaborator

Choose a reason for hiding this comment

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

Did you consider leaning on our keyword-based type system for these decisions?

Copy link
Member Author

Choose a reason for hiding this comment

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

Good idea. I'll do that in the same followup as the other suggestion.

(let [path (join ":" chain)]
(symbol
(str name path))))
(structure->access-chains s)))
Copy link
Collaborator

Choose a reason for hiding this comment

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

we could keep compatibility with scmutils of course. We could also, instead of using :, use and _. Then it would be the basis for a literal structure. WDYT?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, this is awesome. I'm going to do this in a separate PR. I can do it be returning a sequence of [::up ::down ::down] etc attached as metadata to the chain... then anyone can get access to this info if they want, and we can use it here.

That will, as you say, give us the property that this structure will match our literal structures AND our literal matrix notation.

test/sicmutils/generators.cljc Show resolved Hide resolved
row (m/up->row-matrix s)]
(is (= (g/outer-product s s)
(g/outer-product col row))))))

Copy link
Collaborator

Choose a reason for hiding this comment

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

You could have (unless I overlooked it) a test that asserts that (* (transpose v) v) is square for any structure v

Copy link
Member Author

Choose a reason for hiding this comment

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

is that true? I think any structure now will collapse to a number...

That IS true for transpose-outer with 2-tensors. Good test, will add.

@sritchie sritchie merged commit a00c78d into master Dec 11, 2020
sritchie pushed a commit that referenced this pull request Dec 14, 2020
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.

None yet

3 participants