-
Notifications
You must be signed in to change notification settings - Fork 11
/
Interpolation.py
157 lines (144 loc) · 6.27 KB
/
Interpolation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
"""interpolatons"""
import numpy as np
from numpy import *
import torch
import torch.nn as nn
from ..functional.interpolation import lagrangeinterp,_ele2coe,_fix_inputs,_base
from ...utils import meshgen
__all__ = ['LagrangeInterp', 'LagrangeInterpFixInputs']
class LagrangeInterp(nn.Module):
"""
piecewise Lagrange Interpolation in R^m
Arguments:
interp_dim (int): spatial dimension, m=interp_dim
interp_coe (Tensor): DoubleTensor (cuda) or FloatTensor (cuda).
torch.size(np.array(mesh_size)*interp_degree+1)
interp_degree (int): degree of Lagrange Interpolation Polynomial
mesh_bound (tuple): ((l_1,l_2,...,l_n),(u_1,u_2,...,u_n)). mesh_bound
defines the interpolation domain. l_i,u_i is lower and upper bound
of dimension i.
mesh_size (tuple): mesh_size defines the grid number of
piecewise interpolation. mesh_size[i] is mesh num of dimension i.
Usage:
# 3rd-order polynomial Lagrange interpolation in 2d domain
# [0,1]x[1,2], with mesh_size 50x50
Approximator = LagrangeInterp(2,3,[(0,1),(1,2)],(50,50))
print(Approximator.interp_coe) # a 151x151 tensor (mesh_size*interp_degree+1)
InputCoord = torch.rand(10,10,10,10,2)
InputCoord[...,1] += 1
Approximator(InputCoord)
"""
def __init__(self, interp_dim, interp_degree, mesh_bound, mesh_size):
super(LagrangeInterp, self).__init__()
self.__m = interp_dim
self.__d = interp_degree
mesh_bound = array(mesh_bound).reshape(2,self.__m)
mesh_size = array(mesh_size).reshape(self.__m)
self.__mesh_bound = mesh_bound.copy()
self.__mesh_size = mesh_size.copy()
__ele2coe = _ele2coe(self.m, self.d)
__ele2coe = torch.from_numpy(__ele2coe).long()
# nn.Module.to(dtype) will only cast the floating point parameters
# and buffers to dtype
self.register_buffer('ele2coe', __ele2coe)
mesh_size = list(map(lambda x:int(x), list(mesh_size*self.d+1)))
interp_coe = torch.Tensor(*mesh_size).normal_()
self.interp_coe = nn.Parameter(interp_coe)
def init(self, func, is_numpy_func=False):
inputs = meshgen(self.mesh_bound, self.mesh_size*self.d, endpoint=True)
if not is_numpy_func:
inputs = torch.from_numpy(inputs).to(self.interp_coe)
self.interp_coe.data = func(inputs)
else:
self.interp_coe.data.copy_(torch.from_numpy(func(inputs)))
return None
@property
def m(self):
return self.__m
@property
def d(self):
return self.__d
@property
def mesh_bound(self):
return self.__mesh_bound
@property
def mesh_size(self):
return self.__mesh_size
def forward(self, inputs):
"""
piecewise Lagrange Interpolation in R^m
Arguments:
inputs (Tensor): DoubleTensor (cuda) or FloatTensor (cuda).
torch.size=[...,m], where m is the spatial dimension.
"""
size = inputs.size()
if self.m == 1 and size[-1] != 1:
inputs = inputs[...,newaxis]
size = inputs.size()
inputs = inputs.contiguous()
inputs = inputs.view([-1,self.m])
outputs = lagrangeinterp(inputs, self.interp_coe, self.m, self.d,
self.mesh_bound, self.mesh_size, ele2coe=self.ele2coe)
return outputs.view(size[:-1])
class LagrangeInterpFixInputs(LagrangeInterp):
"""
piecewise Lagrange Interpolation in R^m for fixed inputs.
Arguments:
inputs (Tensor): DoubleTensor (cuda) or FloatTensor (cuda).
torch.size=[...,m], where m is the spatial dimension.
interp_dim (int): spatial dimension, m=interp_dim
interp_coe (Tensor): DoubleTensor (cuda) or FloatTensor (cuda).
torch.size(np.array(mesh_size)*interp_degree+1)
interp_degree (int): degree of Lagrange Interpolation Polynomial
mesh_bound (tuple): ((l_1,l_2,...,l_n),(u_1,u_2,...,u_n)). mesh_bound
defines the interpolation domain. l_i,u_i is lower and upper bound
of dimension i.
mesh_size (tuple): mesh_size defines the grid number of
piecewise interpolation. mesh_size[i] is mesh num of dimension i.
Usage:
LagrangeInterpFixInputs is similar with LagrangeInterp except
it was optimized for interpolating fixed inputs.
"""
def __init__(self, inputs, interp_dim, interp_degree, mesh_bound, mesh_size):
super(LagrangeInterpFixInputs, self).__init__(interp_dim,
interp_degree, mesh_bound, mesh_size)
assert isinstance(inputs, torch.Tensor)
self.to(dtype=inputs.dtype, device=inputs.device)
self.register_buffer('_inputs',inputs.new(1))
self.register_buffer('flat_indices', inputs.new(1).long())
self.register_buffer('points_shift', inputs.new(1))
self.register_buffer('base', inputs.new(1))
self.update_inputs(inputs)
def update_inputs(self, inputs):
inputs = inputs.clone()
inputs.data = inputs.data.to(self._inputs.device)
inputs.data = inputs.data.contiguous()
size = inputs.size()
if self.m == 1 and size[-1] != 1:
inputs = inputs[...,newaxis]
size = inputs.size()
self.__inputs_size = size
inputs = inputs.view([-1,self.m])
self._inputs = inputs.data.view(size)
flat_indices, points_shift = _fix_inputs(inputs, self.m, self.d, \
self.mesh_bound, self.mesh_size, self.ele2coe)
self.flat_indices = flat_indices.data
self.points_shift = points_shift.data
base = _base(points_shift, self.m, self.d)
self.base = base.data
@property
def inputs(self):
return self._inputs
@inputs.setter
def inputs(self, v):
self.update_inputs(v)
@property
def inputs_size(self):
return self.__inputs_size
def forward(self):
return lagrangeinterp(self._inputs.view([-1,self.m]), self.interp_coe,
self.m, self.d, self.mesh_bound, self.mesh_size,
ele2coe=self.ele2coe, fix_inputs=True,
flat_indices=self.flat_indices,
points_shift=self.points_shift,
base=self.base).view(self.__inputs_size[:-1])