parameter substitution

added 1 commit

3938aaf2 - add tests for parameter substitution

added 7 commits

59b4195f - add parameter substituter and tests
61ab6d96 - make '_common.get_parameters' return a named tuple
319cdf35 - add generic utilities
5103c4ad - factor out common shallow copying logic from Builder.reversed
0f21dfe5 - add Builder.subs method
67e5f7b2 - refactor Builder.edges to use interleave
049a04b4 - refactor Builder.reversed to use Builder._shallow_copy

OK the functionality's there, just need to add tests now.

Note to reviewers: while implementing this change I needed 2 common pieces of functionality: interleaving elements of a list and shallow copying a Builder.

In both cases I factored these pieces out from existing functionality in the builder module (edges and Builder.reversed), but I only changed the existing implementations of edges and Builder.reversed to use this factored-out functionality in the top 2 commits of this MR.

I think this makes sense to have this change in this MR, but to keep the commits separate to make it clear for people looking through the history what's going on.

added 7 commits

2c42c3db - add parameter substituter and tests
cee758b4 - make '_common.get_parameters' return a named tuple
3d224525 - add generic utilities
7130dacd - factor out common shallow copying logic from Builder.reversed
b11fe540 - add Builder.subs method
698ab877 - refactor Builder.edges to use interleave
a4ca4fd8 - refactor Builder.reversed to use Builder._shallow_copy

Compare with previous version

added 3 commits

e01655e4 - add tests for Builder.subs
13e0772b - refactor Builder.edges to use interleave
cb1451aa - refactor Builder.reversed to use Builder._shallow_copy

Compare with previous version

added 2 commits

b93f11ba - replace 'flatten' with 'chain.from_iterable'
48e3cf0a - rename Builder._shallow_copy to Builder.__copy__

Compare with previous version

resolved all discussions

unmarked as a Work In Progress

added 5 commits

57411b54 - add tests for Builder.subs
b0074a19 - refactor Builder.edges to use interleave
0c6ce8d0 - refactor Builder.reversed to use Builder._shallow_copy
edf21dbe - replace 'flatten' with 'chain.from_iterable'
001dfa43 - rename Builder._shallow_copy to Builder.__copy__

Compare with previous version

Still need to add WhatsNew and add an example to the tutorial.

It probably belongs wherever we introduce "models" that we use to fill builders.

A relevant point that we should at the very least discuss (even if we don't address it in the MR).

This MR nicely allows to replace Hamiltonian parameters. Thanks to Python dynamic nature, passing functions as params also allows to provide different functional forms, so that the same model can correspond to different position dependence of the potential. However builder.subs won't allow replacing a constant value with a callable.

To explain with an example consider this code:

graphene = graphene_model()
scattering_region = Builder()
scattering_region.fill(graphene.subs(mu=mu_bulk), shape)
lead = Builder(TranslationalSymmetry(...))
lead.fill(graphene.subs(mu=mu_lead), shape)
scattering_region.attach_lead(lead)
...

if I want to make chemical potential site-dependent in the scattering region (relevant if it is disordered), I must define the model like this:

def graphene_model_onsite(site, mu):
    return mu(site)

Then I'd use the finalized system as

kwant.smatrix(syst, params={'mu_bulk': lambda site: tanh(site.pos[0])}

But then I am forced to

Declare all parameters that may reasonably vary as callables
Pass even the constant parameters as lambda functions.

mentioned in issue #191 (closed)

This patch series looks fine as always (up to the minor comments listed below), but I'm troubled by a development that has started with the introduction of params and continues here: hamiltonian() is getting complex internally which is bad in itself but also leaves to slowness. I think that this will continue to be a problem even once we have vecotrization.

With this MR, when hamiltonian() is called for a substituted builder, as far as I can see three (!) dictionaries will be constructed (the first one is params in _FinalizedBuilderMixin.hamiltonian(), the second and third one is in _invert_pmap) There will be also a flurry of loops and function calls. Compare this with an old version of hamiltonian().

I think that most of this could be replaced by something simple along the lines of:

def hamiltonian(self, i, j, params):
    if i == j:
        if i >= self.cell_size:
            i -= self.cell_size
        value = self.onsite_hamiltonians[i]
        if callable(value):
            args = (params[name] for name in self._param_names[value])
            if isinstance(value, Substitution):
                # The substitution has been taken account of in _param_names.
                value = value.orig
            value = value(self.symmetry.to_fd(self.sites[i]), *args)
    else:
        # Treat hopping value functions accordingly.

The dictionary param_names would be precomputed. Note how simple and efficient it is and how there's minimal extra cost of the parameter substitution! Even if I forgot to treat some extra cases, I think the idea is clear.

There's still a callable Substitution class, but it's only there so that users can get values from a builder and call them manually.

I'm for radically simplifying how parameters are treated internally even if this means for example that value functions cannot have keyword-only args. IMHO that's totally OK. We don't have to support every single feature of Python for value functions.

Here are some minor points (that partially concern code that should be removed/rewritten if the above is taken into account).

I wonder whether the name "substitute" is the best for what we're doing here. When I was explaining this MR to someone today it took a while until he understood that what we're doing here is simply renaming the parameters (of a copy of a builder). Quite often (in many CAS for example) substitution is a more powerful operation that allows to replace an expression by another one. However, I cannot think of a better alternative. rename seems confusing in other ways and rename_params is too long. So perhaps it should be subs after all.
Shouldn't ParameterSubsitution.__eq__ also check whether isinstance(other, ParameterSubsitution)?
I suggest only putting stuff in _common that is truly used in multiple modules. Otherwise _common becomes a dumping ground for assorted functionality. IMHO all the parameter substitution stuff can be moved from _common to builder.
I assume that this is an incremental patch series and that it will be rebased and partially squashed when finished. For example the three commits: (1) add shallow copy operations to builders, (2) use this inside reversed, (3) rename it to copy should be all squashed into a single one for clarity.
There should be a tutorial for this. I realize that there isn't even a tutorial for params. We should probably modernize the tutorial to use fill and subs, at least in some places.

I assume that this is an incremental patch series and that it will be rebased and partially squashed when finished. For example the three commits: (1) add shallow copy operations to builders, (2) use this inside reversed, (3) rename it to copy should be all squashed into a single one for clarity.

You know that I always do this. The only reason that I sometimes prefer to wait until after review to squash is so that reviewers can easily see if I addressed their points by just looking at commit messages.

There should be a tutorial for this

Already noted in my comment above yours

I think that most of this could be replaced by something simple along the lines of...

So when you said " This patch series looks fine as always", what you really meant was "please rewrite everything" ;)?

So when you said " This patch series looks fine as always", what you really meant was "please rewrite everything" ;)?

This remark was about the quality of the code, which is very nice as usual.

My doubts are about how hamiltonian() works internally. Most of what you wrote would remain if we adopt what I propose. Specifically, the class ParameterSubstitution and the related bits must remain, because substituted values should be callable. (I didn't look whether ParameterSubstitution.__call__ can be optimized - if not it's not a problem.)

BTW, even with my proposed design we can support keyword-only args in value functions if this is a concern. That would require replacing the dictionary _param_names of a builder instance with two: _arg_names and _kwarg_names. The latter would be only used for keyword-only args, so the only penality would be testing whether the value function at hand occurs in _kwarg_names. The kwargs dictionary for it would be constructed only if needed (when there are keyword-only args).

@cwg I don't understand how computing self._param_names[value] would work in the following pseudocode:

model = model_with_parameters_a_and_b
builder1.fill(shape1, model.subs(a='c'))
builder2.fill(shape2, model.subs(b='c'))
builder1.update(builder2)
builder1.finalized().hamiltonian_submatrix(a=1, b=2, c=3)

It appears we do need a more detailed data structure.

I'm for radically simplifying how parameters are treated internally even if this means for example that value functions cannot have keyword-only args. IMHO that's totally OK. We don't have to support every single feature of Python for value functions.

Do we care that we'd break backwards compatibility if we did this now? Probably there are not many Kwant users who are using keyword-only parameters in their value functions so maybe it's ok to do this.

Joseph wrote:

Do we care that we'd break backwards compatibility if we did this now? Probably there are not many Kwant users who are using keyword-only parameters in their value functions so maybe it's ok to do this.

I think that it's OK to require that value functions have no keyword-only args and that default parameter values are ignored (because the user would be forced to supply all the parameters). Technically that's breaking backwards compatibility, but I don't think that this will hurt many people.

The question is rather whether it's useful to have default args or keyword-only args in value functions. With default values there's the risk that some value function already uses a keyword with a default value. Now the user adds another value function with the same keyword and when the param is set this also affects the forgotten keyword. That's why I think that args with default values are dangerous for value functions and should not be accepted (i.e. value functions can have args with default values, but those would be ignored by systems.)

Keyword-only args don't seem to cause problems with value functions, but then their doesn't seem to be any point in using them, because all the parameters are set using a dictionary anyway. So I think that for simplicity we may choose not to support them.

However, if you do not agree with the above reasoning, my proposed simplification/optimization of the inner workings of hamiltonian can be done in a way that supports any kind of Python function:

We just need to store in the _param_names dictionary for each value function triples of (params, required_keywords, optional_keywords). The first item, params, would store the names of the positional parameter for that function, like in my proposal above. The other two items would be typically empty sequences, but when non-empty required_keywords would store the names of keyword-only args without a default value, and all the keyword args with a default value (both keyword-only and not) would go into the last one, optional_keywords.

In this way, it seems to me, all possible function signatures would be supported. But in the common case where there are no args with default values nor keyword-only args, no kwargs dictionary would have to be built. And if kwargs would have to be built it would be just a single dict with simple logic, and not three.

Anton wrote:

@cwg I don't understand how computing self._param_names[value] would work in the following pseudocode: (...)

Let's assume without loss of generality that model_with_parameters_a_and_b has a single value function that requires parameters a and b.

Then the finished builder1 will contain two different value functions, namely two different substitutions of the same sole value function of model. These two values are separate instances of Substitution but their orig attribute (to use the names of my snippet above) contains the same object, the value function of model. So there will be two entries in _param_names, because builder contains two different callable values, even though there is only one "real" value function.

OK, I see. Thanks for the clarification.

I have a workflow proposal based on the following observations:

The MR as it stands is not a performance regression: existing code does not have a slowdown; only the Substitution instances are called with a suboptimal performance.
The MR currently does not modify FinalizedBuilder.hamiltonian, something that @cwg's proposal would require.
@cwg's suggestion has nontrivial consequences for functionality, in particular it changes the scope of what callable values we support.
@cwg's suggestion does not impact the new API introduced by this MR.

Based on these, it seems to make most sense to limit the scope of our review, and separate the rewriting of the FinalizedBuilder.hamiltonian into a different issue/MR.

added 6 commits

3c5e292d - add parameter substituter and tests
ef0af748 - make '_common.get_parameters' return a named tuple
4031d4e2 - add generic 'interleave' utility and use it in Builder.edges
40f346d6 - factor out common shallow copying logic from Builder.reversed
ad0b4e76 - add Builder.subs method and tests
8062b696 - add whatsnew entry

Compare with previous version

I just squashed everything down.

AFAICT this MR is ready to go:

the additions to the tutorial should go in another MR (see the unresolved discussion above)
the rewrite of hamiltonian should also go in anther MR (see previous post by anton)

added 2 commits

e2c6edd0 - add Builder.subs method and tests
0ebe6076 - add whatsnew entry

Compare with previous version

As said, the code of this MR is nice and clean! I do not have a strong opinion on how modifications to Kwant are split into MRs.

That said, I have a very strong opinion that this code is not ready for release. We must fix the performance issue (and we should fix the dangerous semantics of "params") before releasing this as part of Kwant 1.4, or revert this change. Here is a real-world benchmark based on the first Kwant tutorial:

import time
import kwant


lat = kwant.lattice.square()


def make_syst_const(W, L):
    syst = kwant.Builder()
    syst[(lat(x, y) for x in range(L) for y in range(W))] = 4
    syst[lat.neighbors()] = -1
    lead = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((-1, 0))))
    lead[(lat(0, j) for j in range(W))] = 4
    lead[lat.neighbors()] = -1
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())
    return syst


def onsite(site, t, u):
    return 4 * t + u


def hopping(a, b, t, u):
    return -t


def make_syst_vfs(W, L):
    syst = kwant.Builder()
    syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
    syst[lat.neighbors()] = hopping
    lead = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((-1, 0))))
    lead[(lat(0, j) for j in range(W))] = onsite
    lead[lat.neighbors()] = hopping
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())
    return syst


def make_syst_subs(W, L):
    syst = kwant.Builder()
    syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
    syst[lat.neighbors()] = hopping
    syst = syst.subs(t='tt', u='uu')
    lead = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((-1, 0))))
    lead[(lat(0, j) for j in range(W))] = onsite
    lead[lat.neighbors()] = hopping
    lead = lead.subs(t='tt', u='uu')
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())
    return syst


def main(W=50, L=50, energies=tuple(i/100 for i in range(50))):
    for name, f, kwargs in [("consts", make_syst_const, {}),
                            ("args", make_syst_vfs,
                             dict(args=(1, 0))),
                            ("params", make_syst_vfs,
                             dict(params=dict(t=1, u=0))),
                            ("subs", make_syst_subs,
                             dict(params=dict(tt=1, uu=0)))]:
        syst = f(W, L).finalized()
        t = time.time()
        transmissions = [kwant.smatrix(syst, e, **kwargs).transmission(1, 0)
                         for e in energies]
        print(name, time.time() - t)


if __name__ == '__main__':
    main()

On my machine the output is:

consts 5.318358659744263
args 5.975970268249512
params 6.638056993484497
subs 14.160792112350464

Already "params" should be able to get way closer to "args". "params" is clearly out of proportion. There is no reason why "params" should not be almost as fast as "args".

The above comment is outside the scope of this merge request.

After a discussion with @anton-akhmerov, we think that for now we should merely add a "warning" note to the docstring of Builder.subs.

mentioned in issue #195 (closed)

Given the above I propose to merge this in. We have already added #195 (closed) where we can address @cwg's only complaint

added 1 commit

ef1eead3 - add a warning to Builder.subs about potential slowdowns

Compare with previous version

resolved all discussions

Sorry that it took me so long to have another look.

Thanks, please go ahead and merge this. Some final comments:

Couldn't the first line of _invert_map be moved to __init__? I.e. pmap would be precomputed and stored instead of substitutions. If this is possible, then I think there's not much point in not inlining the remainder of _invert_map into __call__.
Other optimizations are possible, but I think it's better to do them when we review/rework how params work.

Couldn't the first line of _invert_map be moved to __init__? I.e. pmap would be precomputed and stored instead of substitutions. If this is possible, then I think there's not much point in not inlining the remainder of _invert_map into __call__

This was for 2 reasons: to store only immutable datastructures, and to have the meaning of the stored substitutions be old → new. If we change this now we'll have to invert the logic everywhere

added 12 commits

ef1eead3...40b89650 - 5 commits from branch kwant:master
5b8675f9 - add parameter substituter and tests
cd15f40f - make '_common.get_parameters' return a named tuple
2ecd6f8f - add generic 'interleave' utility and use it in Builder.edges
7bb1bab2 - factor out common shallow copying logic from Builder.reversed
790a6f92 - add Builder.subs method and tests
7457563a - add whatsnew entry
6e4f89f7 - add a warning to Builder.subs about potential slowdowns

Compare with previous version

enabled an automatic merge when the pipeline for 6e4f89f7 succeeds

mentioned in commit bdc21039

merged

mentioned in commit jbweston/kwant@bbb864fe

parameter substitution

Merge request reports

Activity