arcanus/koala2

1import torch
2from torch.nn import functional as F
3
4import numpy as np
5
6
7DEFAULT_MIN_BIN_WIDTH = 1e-3
8DEFAULT_MIN_BIN_HEIGHT = 1e-3
9DEFAULT_MIN_DERIVATIVE = 1e-3
10
11
12def piecewise_rational_quadratic_transform(
13    inputs,
14    unnormalized_widths,
15    unnormalized_heights,
16    unnormalized_derivatives,
17    inverse=False,
18    tails=None,
19    tail_bound=1.0,
20    min_bin_width=DEFAULT_MIN_BIN_WIDTH,
21    min_bin_height=DEFAULT_MIN_BIN_HEIGHT,
22    min_derivative=DEFAULT_MIN_DERIVATIVE,
23):
24    if tails is None:
25        spline_fn = rational_quadratic_spline
26        spline_kwargs = {}
27    else:
28        spline_fn = unconstrained_rational_quadratic_spline
29        spline_kwargs = {"tails": tails, "tail_bound": tail_bound}
30
31    outputs, logabsdet = spline_fn(
32        inputs=inputs,
33        unnormalized_widths=unnormalized_widths,
34        unnormalized_heights=unnormalized_heights,
35        unnormalized_derivatives=unnormalized_derivatives,
36        inverse=inverse,
37        min_bin_width=min_bin_width,
38        min_bin_height=min_bin_height,
39        min_derivative=min_derivative,
40        **spline_kwargs
41    )
42    return outputs, logabsdet
43
44
45def searchsorted(bin_locations, inputs, eps=1e-6):
46    bin_locations[..., -1] += eps
47    return torch.sum(inputs[..., None] >= bin_locations, dim=-1) - 1
48
49
50def unconstrained_rational_quadratic_spline(
51    inputs,
52    unnormalized_widths,
53    unnormalized_heights,
54    unnormalized_derivatives,
55    inverse=False,
56    tails="linear",
57    tail_bound=1.0,
58    min_bin_width=DEFAULT_MIN_BIN_WIDTH,
59    min_bin_height=DEFAULT_MIN_BIN_HEIGHT,
60    min_derivative=DEFAULT_MIN_DERIVATIVE,
61):
62    inside_interval_mask = (inputs >= -tail_bound) & (inputs <= tail_bound)
63    outside_interval_mask = ~inside_interval_mask
64
65    outputs = torch.zeros_like(inputs)
66    logabsdet = torch.zeros_like(inputs)
67
68    if tails == "linear":
69        unnormalized_derivatives = F.pad(unnormalized_derivatives, pad=(1, 1))
70        constant = np.log(np.exp(1 - min_derivative) - 1)
71        unnormalized_derivatives[..., 0] = constant
72        unnormalized_derivatives[..., -1] = constant
73
74        outputs[outside_interval_mask] = inputs[outside_interval_mask]
75        logabsdet[outside_interval_mask] = 0
76    else:
77        raise RuntimeError("{} tails are not implemented.".format(tails))
78
79    (
80        outputs[inside_interval_mask],
81        logabsdet[inside_interval_mask],
82    ) = rational_quadratic_spline(
83        inputs=inputs[inside_interval_mask],
84        unnormalized_widths=unnormalized_widths[inside_interval_mask, :],
85        unnormalized_heights=unnormalized_heights[inside_interval_mask, :],
86        unnormalized_derivatives=unnormalized_derivatives[inside_interval_mask, :],
87        inverse=inverse,
88        left=-tail_bound,
89        right=tail_bound,
90        bottom=-tail_bound,
91        top=tail_bound,
92        min_bin_width=min_bin_width,
93        min_bin_height=min_bin_height,
94        min_derivative=min_derivative,
95    )
96
97    return outputs, logabsdet
98
99
100def rational_quadratic_spline(
101    inputs,
102    unnormalized_widths,
103    unnormalized_heights,
104    unnormalized_derivatives,
105    inverse=False,
106    left=0.0,
107    right=1.0,
108    bottom=0.0,
109    top=1.0,
110    min_bin_width=DEFAULT_MIN_BIN_WIDTH,
111    min_bin_height=DEFAULT_MIN_BIN_HEIGHT,
112    min_derivative=DEFAULT_MIN_DERIVATIVE,
113):
114    if torch.min(inputs) < left or torch.max(inputs) > right:
115        raise ValueError("Input to a transform is not within its domain")
116
117    num_bins = unnormalized_widths.shape[-1]
118
119    if min_bin_width * num_bins > 1.0:
120        raise ValueError("Minimal bin width too large for the number of bins")
121    if min_bin_height * num_bins > 1.0:
122        raise ValueError("Minimal bin height too large for the number of bins")
123
124    widths = F.softmax(unnormalized_widths, dim=-1)
125    widths = min_bin_width + (1 - min_bin_width * num_bins) * widths
126    cumwidths = torch.cumsum(widths, dim=-1)
127    cumwidths = F.pad(cumwidths, pad=(1, 0), mode="constant", value=0.0)
128    cumwidths = (right - left) * cumwidths + left
129    cumwidths[..., 0] = left
130    cumwidths[..., -1] = right
131    widths = cumwidths[..., 1:] - cumwidths[..., :-1]
132
133    derivatives = min_derivative + F.softplus(unnormalized_derivatives)
134
135    heights = F.softmax(unnormalized_heights, dim=-1)
136    heights = min_bin_height + (1 - min_bin_height * num_bins) * heights
137    cumheights = torch.cumsum(heights, dim=-1)
138    cumheights = F.pad(cumheights, pad=(1, 0), mode="constant", value=0.0)
139    cumheights = (top - bottom) * cumheights + bottom
140    cumheights[..., 0] = bottom
141    cumheights[..., -1] = top
142    heights = cumheights[..., 1:] - cumheights[..., :-1]
143
144    if inverse:
145        bin_idx = searchsorted(cumheights, inputs)[..., None]
146    else:
147        bin_idx = searchsorted(cumwidths, inputs)[..., None]
148
149    input_cumwidths = cumwidths.gather(-1, bin_idx)[..., 0]
150    input_bin_widths = widths.gather(-1, bin_idx)[..., 0]
151
152    input_cumheights = cumheights.gather(-1, bin_idx)[..., 0]
153    delta = heights / widths
154    input_delta = delta.gather(-1, bin_idx)[..., 0]
155
156    input_derivatives = derivatives.gather(-1, bin_idx)[..., 0]
157    input_derivatives_plus_one = derivatives[..., 1:].gather(-1, bin_idx)[..., 0]
158
159    input_heights = heights.gather(-1, bin_idx)[..., 0]
160
161    if inverse:
162        a = (inputs - input_cumheights) * (
163            input_derivatives + input_derivatives_plus_one - 2 * input_delta
164        ) + input_heights * (input_delta - input_derivatives)
165        b = input_heights * input_derivatives - (inputs - input_cumheights) * (
166            input_derivatives + input_derivatives_plus_one - 2 * input_delta
167        )
168        c = -input_delta * (inputs - input_cumheights)
169
170        discriminant = b.pow(2) - 4 * a * c
171        assert (discriminant >= 0).all()
172
173        root = (2 * c) / (-b - torch.sqrt(discriminant))
174        outputs = root * input_bin_widths + input_cumwidths
175
176        theta_one_minus_theta = root * (1 - root)
177        denominator = input_delta + (
178            (input_derivatives + input_derivatives_plus_one - 2 * input_delta)
179            * theta_one_minus_theta
180        )
181        derivative_numerator = input_delta.pow(2) * (
182            input_derivatives_plus_one * root.pow(2)
183            + 2 * input_delta * theta_one_minus_theta
184            + input_derivatives * (1 - root).pow(2)
185        )
186        logabsdet = torch.log(derivative_numerator) - 2 * torch.log(denominator)
187
188        return outputs, -logabsdet
189    else:
190        theta = (inputs - input_cumwidths) / input_bin_widths
191        theta_one_minus_theta = theta * (1 - theta)
192
193        numerator = input_heights * (
194            input_delta * theta.pow(2) + input_derivatives * theta_one_minus_theta
195        )
196        denominator = input_delta + (
197            (input_derivatives + input_derivatives_plus_one - 2 * input_delta)
198            * theta_one_minus_theta
199        )
200        outputs = input_cumheights + numerator / denominator
201
202        derivative_numerator = input_delta.pow(2) * (
203            input_derivatives_plus_one * theta.pow(2)
204            + 2 * input_delta * theta_one_minus_theta
205            + input_derivatives * (1 - theta).pow(2)
206        )
207        logabsdet = torch.log(derivative_numerator) - 2 * torch.log(denominator)
208
209        return outputs, logabsdet
210