RDataReporting::residuals with adjoints

[dweindl]

So far, when using adjoint sensitivities, it's only possible to get ReturnData::y with RDataReporting::full. RDataReporting::residuals is not supported:

AMICI/src/solver.cpp

Lines 738 to 741 in 2d58f54

	if (rdata_mode_ == RDataReporting::residuals
	&& sensi_meth == SensitivityMethod::adjoint)
	throw AmiException("Adjoint Sensitivity Analysis is not compatible with"
	" only reporting residuals!");

We should either allow that, and just skip the respective sensitivities (or basically just use ReturnData::initializeResidualReporting(false) in case of ASA), or introduce a new mode.

@FFroehlich : Any preference?

CC @Doresic

[FFroehlich]

I agree that RDataReporting::full should give you everything that is available. ReturnData::y and RDataReporting::residuals would be available but corresponding sensitivities are of course not.

I don't quite get the relationship to the linked piece of code, it checks for rdata_mode_ == RDataReporting::residuals, so this does not apply for RDataReporting::full.

[dweindl]

I agree that RDataReporting::full should give you everything that is available.

That's what's happening. All fine there

I don't quite get the relationship to the linked piece of code, it checks for rdata_mode_ == RDataReporting::residuals, so this does not apply for RDataReporting::full.

That is also correct.

But I want to avoid using RDataReporting::full, because I only need y and sllh. The closest thing to that would be RDataReporting::residuals. But the referenced code forbids that when using adjoint sensitivities.

[Doresic]

A nicer thing for my particular application (due to memory issues) would be to have another amici.RDataReporting.observables (naming up to change) that would give only llh, grad, y, and sigmay.

However, this differs from amici.RDataReporting.residuals only by res, if I'm not missing something. So that shouldn't be much more memory and should be a sufficient reduction from amici.RDataReporting.full.

[FFroehlich]

Yeah I agree that it sounds like introducing a new mode would be the way to go.

[dweindl]

New mode is fine for me. However, I am not sure if this makes too much sense:

AMICI/src/solver.cpp

Lines 738 to 741 in 2d58f54

	if (rdata_mode_ == RDataReporting::residuals
	&& sensi_meth == SensitivityMethod::adjoint)
	throw AmiException("Adjoint Sensitivity Analysis is not compatible with"
	" only reporting residuals!");

Because the respective sensitivity members are only initialized for FSA anyways and would just be skipped later on:

AMICI/src/rdata.cpp

Lines 81 to 96 in 2d58f54

	void ReturnData::initializeResidualReporting(bool enable_res) {
	y.resize(nt * ny, 0.0);
	sigmay.resize(nt * ny, 0.0);
	if (enable_res)
	res.resize((sigma_res ? 2 : 1) * nt * nytrue, 0.0);
	if ((sensi_meth == SensitivityMethod::forward
	&& sensi >= SensitivityOrder::first)
	|| sensi >= SensitivityOrder::second) {
	sy.resize(nt * ny * nplist, 0.0);
	ssigmay.resize(nt * ny * nplist, 0.0);
	if (enable_res)
	sres.resize((sigma_res ? 2 : 1) * nt * nytrue * nplist, 0.0);
	}
	}

Removing that check above would be sufficient in my opinion, and we wouldn't need a new mode. But maybe I am missing something.

[FFroehlich]

New mode is fine for me. However, I am not sure if this makes too much sense:

AMICI/src/solver.cpp

Lines 738 to 741 in 2d58f54

	if (rdata_mode_ == RDataReporting::residuals
	&& sensi_meth == SensitivityMethod::adjoint)
	throw AmiException("Adjoint Sensitivity Analysis is not compatible with"
	" only reporting residuals!");

Because the respective sensitivity members are only initialized for FSA anyways and would just be skipped later on:

AMICI/src/rdata.cpp

Lines 81 to 96 in 2d58f54

	void ReturnData::initializeResidualReporting(bool enable_res) {
	y.resize(nt * ny, 0.0);
	sigmay.resize(nt * ny, 0.0);
	if (enable_res)
	res.resize((sigma_res ? 2 : 1) * nt * nytrue, 0.0);
	if ((sensi_meth == SensitivityMethod::forward
	&& sensi >= SensitivityOrder::first)
	|| sensi >= SensitivityOrder::second) {
	sy.resize(nt * ny * nplist, 0.0);
	ssigmay.resize(nt * ny * nplist, 0.0);
	if (enable_res)
	sres.resize((sigma_res ? 2 : 1) * nt * nytrue * nplist, 0.0);
	}
	}

Removing that check above would be sufficient in my opinion, and we wouldn't need a new mode. But maybe I am missing something.

Ah I see, yeah sounds plausible to simply remove that check.

[dweindl]

Unfortunately, it won't work that way. sllh is not computed then. I somehow thought, RDataReporting::residuals was a superset of RDataReporting::likelihood, but it's not.

Options

1. always including sllh with RDataReporting::residuals (maybe not so desirable, although probably negligible impact)
2. including sllh with RDataReporting::residuals if ASA is used (would make sense, otherwise there is no reason to use this combination)
3. new option RDataReporting::observables_sllh
4. supporting more or less arbitrary combinations via some flags (reporting=Y|SIGMA|SLLH, reporting=SX, ...)

I think, I'd go for option 2. No strong opinion, though.

[FFroehlich]

Unfortunately, it won't work that way. sllh is not computed then. I somehow thought, RDataReporting::residuals was a superset of RDataReporting::likelihood, but it's not.

Options

1. always including sllh with RDataReporting::residuals (maybe not so desirable, although probably negligible impact)
2. including sllh with RDataReporting::residuals if ASA is used (would make sense, otherwise there is no reason to use this combination)
3. new option RDataReporting::observables_sllh
4. supporting more or less arbitrary combinations via some flags (reporting=Y|SIGMA|SLLH, reporting=SX, ...)

I think, I'd go for option 2. No strong opinion, though.

the quick and dirty option would be to simply set sllh as schi2. Otherwise I would prefer option 3 as we would really blur the meaning of RDataReporting flags.