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

Support multiple input for forward #354

Open
wants to merge 12 commits into
base: main
Choose a base branch
from
Prev Previous commit
Next Next commit
update demo
  • Loading branch information
amal-ghamdi committed Jan 22, 2024
commit b35c97e65d7bc9498a1cde418e1637646e874599
70 changes: 38 additions & 32 deletions demos/howtos/multiple_input_fwd.py
Original file line number Diff line number Diff line change
@@ -1,48 +1,61 @@
# %%Imports
import numpy as np
import cuqi

# %%Create a fwd function
def fwd(a,b):
assert len(b) == 10
return a + 2*b
from cuqi.geometry import Continuous1D, Discrete, ConcatenatedGeometries,\
MappedGeometry

# %% GEOMETRIES
range_geometry = cuqi.geometry.Continuous1D(np.linspace(0, 1, 10))
domain_geometry = cuqi.geometry.GeometryTuple(
cuqi.geometry.Continuous1D(np.linspace(0, 1, 10)),
cuqi.geometry.Discrete('a'))
range_geometry = Continuous1D(np.linspace(0, 1, 10))
geom1 = MappedGeometry(Continuous1D(np.linspace(0, 1, 10)),
lambda x: x+1,
lambda x: x-1)
geom2 = Discrete(['a'])
domain_geometry = ConcatenatedGeometries( geom1, geom2)
# Note:
#1. direct product implements functions such as par2fun, fun2par, plot,
# gradient, etc, by calling each function on each component of the direct
# product
#2. another name choice could be ConcatenatedGeometries


domain_geometry.par_dim # returns 11
domain_geometry.fun_dim # returns [(10,), 1]
print(domain_geometry.par_shape) # returns [(10,), (1,)]
print(domain_geometry.fun_shape) # returns [(10,), (1,)]
print(domain_geometry.par_dim) # returns 11
print(domain_geometry.fun_dim) # returns [(10,), 1]
# how to distinguish between 2D geometry and direct product?

vec1 = np.random.randn(10)
funval = domain_geometry.par2fun(vec1, 5) # or passing tuple or list [vec1, 5]?
# funval is a tuple or list of [par2fun(vec1), par2fun(5)]
vec2 = np.random.randn(11)

funval = domain_geometry.par2fun(np.random.randn(11))
# this format is also understood by the geometry
funval = domain_geometry.par2fun(vec2)

par = domain_geometry.fun2par(funval)
# par_stacked is a numpy array of shape (11,)
# this format is used for samplers
par = domain_geometry.fun2par(funval[0], funval[1])
# par is a numpy array of shape (11,)
# this format is helpful for samplers


#%% Assert values are equal
assert np.allclose(funval[0], geom1.par2fun(vec2[:geom1.par_dim]))
assert np.allclose(funval[1], geom2.par2fun(vec2[geom1.par_dim:]))
assert np.allclose(par[:geom1.par_dim], geom1.fun2par(funval[0]))
assert np.allclose(par[geom1.par_dim:], geom2.fun2par(funval[1]))

# %% CUQIarray
a = cuqi.array.CUQIarray(vec1, 5, geometry=domain_geometry)
a_arr = cuqi.array.CUQIarray(vec2, geometry=domain_geometry)
print("a_arr")
print(repr(a_arr))
print("a_arr.funvals")
print(repr(a_arr.funvals))
# 1. Another name could be CUQIarrayList or CUQIarrayTuple or ConcatenatedCUQIArray
# 2. The object can be indexed like a[0] or a[1]
# 3. similarly here, a.funvals returns
# cuqi.array.CUQIarray(a[0].funvals, a[1].funvals, geometry=domain_geometry,
# is_par=False)

a_arr2 = cuqi.array.CUQIarray(funval, geometry=domain_geometry, is_par=False)
print("a_arr2")
print(repr(a_arr2))
print("a_arr2.parameters")
print(repr(a_arr2.parameters))

# %%PRIORS
a = cuqi.distribution.Uniform(0, 1)
Expand All @@ -52,7 +65,6 @@ def fwd(a,b):
# or possibly another name
prior = cuqi.distribution.JointPrior(a, b)


logpdf_val = prior.logpdf(vec1, 5) # returns a.logpdf(vec1) + b.logpdf(5)

samples = prior.sample(10)
Expand All @@ -74,14 +86,18 @@ def fwd(a,b):
# object

#%% Forward model
# %%Create a fwd function
def fwd(a,b):
assert len(b) == 10
return a + 2*b

F = cuqi.model.Model(fwd, range_geometry, domain_geometry)
# OR
F = cuqi.model.Model(fwd, 10, [1, 10]) #(2,2)

y_exact = F(1, np.random.randn(10))
y_exact = F(np.random.randn(11))


# In the second case, the model uses information from the geometry to
# to determine the first and the second argument

Expand Down Expand Up @@ -116,14 +132,4 @@ def grad_b(a, b, direction):
# sampler uses the par_dim to create a proposal of size (11,)
# model can accept input of size (11,)



#%% Update gradient chain rule in model to consider chain from two distributions