Hi. I am try to define import torch class LuminDist(torch.nn.Module): for data fitting with pymc3.
here is my code
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
import pymc3 as pm
from scipy.integrate import quad
import theano.tensor as tt
import theano
import pdb
from pymc3.backends import Text
import torch
def integrand(z, H0, Om):
h = H0/100.0
Or = 4.16e-5/h**2
# For LCDM, putting following line.
fz = 1
# fz = (1 + z)**(3 * (1 + w0 + wa)) * np.exp(-3 * wa * z/(1+z))
Hz = (Or*(1+z)**4 + Om*(1+z)**3 + (1 - Om - Or)*fz)**0.5
return 1.0 / Hz
def integral(z, H0, Om,M):
como_dist= quad(integrand(z,H0,Om), 0, z, args=(H0, Om))[0]
D_L = (1 + z) * como_dist
return 5*np.log10(D_L) + M
import torch
class LuminDist(torch.nn.Module):
def __init__(self):
super(LuminDist, self).__init__()
self.H0 = torch.nn.Parameter(torch.tensor(0.0, dtype=torch.float64))
self.Om = torch.nn.Parameter(torch.tensor(0.0, dtype=torch.float64))
self.M = torch.nn.Parameter(torch.tensor(0.0, dtype=torch.float64))
self.logs = torch.nn.Parameter(torch.tensor(0.0, dtype=torch.float64))
def forward(self, x, y):
mean = integral(z,H0,Om,M)
loglike = -0.5 * torch.sum((y - mean) ** 2 * torch.exp(-2 * self.logs) + 2 * self.logs)
return loglike
import theano
import theano.tensor as tt
class TorchOp(tt.Op):
def __init__(self, module, params, args=None):
self.module = module
self.params = list(params)
if args is None:
self.args = tuple()
else:
self.args = tuple(args)
def make_node(self, *args):
if len(args) != len(self.params):
raise ValueError("dimension mismatch")
args = [tt.as_tensor_variable(a) for a in args]
return theano.gof.Apply(self, args, [tt.dscalar().type()] + [a.type() for a in args])
def infer_shape(self, node, shapes):
return tuple([()] + list(shapes))
def perform(self, node, inputs, outputs):
for p, value in zip(self.params, inputs):
p.data = torch.tensor(value)
if p.grad is not None:
p.grad.detach_()
p.grad.zero_()
result = self.module(*self.args)
result.backward()
outputs[0][0] = result.detach().numpy()
for i, p in enumerate(self.params):
outputs[i+1][0] = p.grad.numpy()
def grad(self, inputs, gradients):
for i, g in enumerate(gradients[1:]):
if not isinstance(g.type, theano.gradient.DisconnectedType):
raise ValueError(
"can't propagate gradients wrt parameter {0}".format(i + 1)
)
return [gradients[0] * d for d in self(*inputs)[1:]]
I do not know how define mean=integral(z,H0,Om,M) correctly in this class.
Any advice and suggestions will be greatly appreciated.