Skip to content

Instantly share code, notes, and snippets.

@davidbau
Last active Dec 28, 2021
Embed
What would you like to do?
Running stats objects for pytorch: mean, variance, covariance, second-moment, quantiles, topk, and combinations.
'''
To use a runningstats object,
1. Create the the desired stat object, e.g., `m = Mean()`
2. Feed it batches via the add method, e.g., `m.add(batch)`
3. Repeat step 2 any number of times.
4. Read out the statistic of interest, e.g., `m.mean()`
Built-in runningstats objects include:
Mean - produces mean()
Variance - mean() and variance() and stdev()
Covariance - mean(), covariance(), correlation(), variance(), stdev()
SecondMoment - moment() is the non-mean-centered covariance, E[x x^T].
Quantile - quantile(), min(), max(), median(), mean(), variance(), stdev()
TopK - topk() returns (values, indexes)
Bincount - bincount() histograms integral data
IoU - intersection(), union(), iou() tally binary co-occurrences.
History - history() returns concatenation of data.
CrossCovariance - covariance between two signals, without self-covariance.
CrossIoU - iou between two signals, without self-IoU.
CombinedStat - aggregates any set of stats.
Add more running stats by subclassing the Stat class.
These statistics are vectorized along dim>=1, so stat.add()
should supply a two-dimensional input where the zeroth
dimension is the batch/sampling dimension and the first
dimension is the feature dimension.
The data type and device used matches the data passed to add();
for example, for higher-precision covariances, convert to double
before calling add().
It is common to want to compute and remember a statistic sampled
over a Dataset, computed in batches, possibly caching the computed
statistic in a file. The tally(stat, dataset, cache) handles
this pattern. It takes a statistic, a dataset, and a cache filename
and sets up a data loader that can be run (or not, if cached) to
compute the statistic, adopting the convention that cached stats are
saved to and loaded from numpy npz files.
'''
import torch
import numpy
import math
import random
import os
import struct
from torch.utils.data.sampler import Sampler
def tally(stat, dataset, cache=None, quiet=False, **kwargs):
'''
To use tally, write code like the following.
stat = Mean()
ds = MyDataset()
for batch in tally(stat, ds, cache='mymean.npz', batch_size=50):
stat.add(batch)
mean = stat.mean()
The first argument should be the Stat being computed. After the
loader is exhausted, tally will bring this stat to the cpu and
cache it (if a cache is specified).
The dataset can be a torch Dataset or a plain Tensor, or it can
be a callable that returns one of those.
Details on caching via the cache= argument:
If the given filename cannot be loaded, tally will leave the
statistic object empty and set up a DataLoader object so that
the loop can be run. After the last iteration of the loop, the
completed statistic will be moved to the cpu device and also
saved in the cache file.
If the cached statistic can be loaded from the given file, tally
will not set up the data loader and instead will return a fully
loaded statistic object (on the cpu device) and an empty list as
the loader.
The `with cache_load_enabled(False):` context manager can
be used to disable loading from the cache.
If needed, a DataLoader will be created to wrap the dataset:
Keyword arguments of tally are passed to the DataLoader,
so batch_size, num_workers, pin_memory, etc. can be specified.
Subsampling is supported via sample_size= and random_sample=:
If sample_size=N is specified, rather than loading the whole
dataset, only the first N items are sampled. If additionally
random_sample=S is specified, the pseudorandom seed S will be
used to select a fixed psedorandom sample of size N to sample.
'''
assert isinstance(stat, Stat)
args = {}
for k in ['sample_size']:
if k in kwargs:
args[k] = kwargs[k]
cached_state = load_cached_state(cache, args, quiet=quiet)
if cached_state is not None:
stat.load_state_dict(cached_state)
def empty_loader():
return
yield
return empty_loader()
loader = make_loader(dataset, **kwargs)
def wrapped_loader():
yield from loader
stat.to_(device='cpu')
if cache is not None:
save_cached_state(cache, stat, args)
return wrapped_loader()
class cache_load_enabled():
'''
When used as a context manager, cache_load_enabled(False) will prevent
tally from loading cached statsitics, forcing them to be recomputed.
'''
def __init__(self, enabled=True):
self.prev = False
self.enabled = enabled
def __enter__(self):
global global_load_cache_enabled
self.prev = global_load_cache_enabled
global_load_cache_enabled = self.enabled
def __exit__(self, exc_type, exc_value, traceback):
global global_load_cache_enabled
global_load_cache_enabled = self.prev
class Stat:
'''
Abstract base class for a running pytorch statistic.
'''
def __init__(self, state):
'''
By convention, all Stat subclasses can be initialized by passing
state=; and then they will initialize by calling load_state_dict.
'''
self.load_state_dict(resolve_state_dict(state))
def add(self, x, *args, **kwargs):
'''
Observes a batch of samples to be incorporated into the statistic.
Dimension 0 should be the batch dimension, and dimension 1 should
be the feature dimension of the pytorch tensor x.
'''
pass
def load_state_dict(self, d):
'''
Loads this Stat from a dictionary of numpy arrays as saved
by state_dict.
'''
pass
def state_dict(self):
'''
Saves this Stat as a dictionary of numpy arrays that can be
stored in an npz or reloaded later using load_state_dict.
'''
return {}
def save(self, filename):
'''
Saves this stat as an npz file containing the state_dict.
'''
save_cached_state(filename, self, {})
def load(self, filename):
'''
Loads this stat from an npz file containing a saved state_dict.
'''
self.load_state_dict(
load_cached_state(filename, {}, quiet=True, throw=True))
def to_(self, device):
'''
Moves this Stat to the given device.
'''
pass
def cpu_(self):
'''
Moves this Stat to the cpu device.
'''
self.to_('cpu')
def cuda_(self):
'''
Moves this Stat to the default cuda device.
'''
self.to_('cuda')
def _normalize_add_shape(self, x, attr='data_shape'):
'''
Flattens input data to 2d.
'''
if not torch.is_tensor(x):
x = torch.tensor(x)
if len(x.shape) < 1:
x = x.view(-1)
data_shape = getattr(self, attr, None)
if data_shape is None:
data_shape = x.shape[1:]
setattr(self, attr, data_shape)
else:
assert x.shape[1:] == data_shape
return x.view(x.shape[0], int(numpy.prod(data_shape)))
def _restore_result_shape(self, x, attr='data_shape'):
'''
Restores output data to input data shape.
'''
data_shape = getattr(self, attr, None)
if data_shape is None:
return x
return x.view(data_shape * len(x.shape))
class Mean(Stat):
'''
Running mean.
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self.batchcount = 0
self._mean = None
self.data_shape = None
def add(self, a):
a = self._normalize_add_shape(a)
if len(a) == 0:
return
batch_count = a.shape[0]
batch_mean = a.sum(0) / batch_count
centered = a - batch_mean
self.batchcount += 1
# Initial batch.
if self._mean is None:
self.count = batch_count
self._mean = batch_mean
return
# Update a batch using Chan-style update for numerical stability.
self.count += batch_count
new_frac = float(batch_count) / self.count
# Update the mean according to the batch deviation from the old mean.
delta = batch_mean.sub_(self._mean).mul_(new_frac)
self._mean.add_(delta)
def size(self):
return self.count
def mean(self):
return self._restore_result_shape(self._mean)
def to_(self, device):
if self._mean is not None:
self._mean = self._mean.to(device)
def load_state_dict(self, state):
self.count = state['count']
self.batchcount = state['batchcount']
self._mean = torch.from_numpy(state['mean'])
self.data_shape = None if state['data_shape'] is None else tuple(state['data_shape'])
def state_dict(self):
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
count=self.count,
data_shape=self.data_shape and tuple(self.data_shape),
batchcount=self.batchcount,
mean=self._mean.cpu().numpy())
class Variance(Stat):
'''
Running computation of mean and variance. Use this when you just need
basic stats without covariance.
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self.batchcount = 0
self._mean = None
self.v_cmom2 = None
self.data_shape = None
def add(self, a):
a = self._normalize_add_shape(a)
if len(a) == 0:
return
batch_count = a.shape[0]
batch_mean = a.sum(0) / batch_count
centered = a - batch_mean
self.batchcount += 1
# Initial batch.
if self._mean is None:
self.count = batch_count
self._mean = batch_mean
self.v_cmom2 = centered.pow(2).sum(0)
return
# Update a batch using Chan-style update for numerical stability.
oldcount = self.count
self.count += batch_count
new_frac = float(batch_count) / self.count
# Update the mean according to the batch deviation from the old mean.
delta = batch_mean.sub_(self._mean).mul_(new_frac)
self._mean.add_(delta)
# Update the variance using the batch deviation
self.v_cmom2.add_(centered.pow(2).sum(0))
self.v_cmom2.add_(delta.pow_(2).mul_(new_frac * oldcount))
def size(self):
return self.count
def mean(self):
return self._restore_result_shape(self._mean)
def variance(self, unbiased=True):
return self._restore_result_shape(self.v_cmom2
/ (self.count - (1 if unbiased else 0)))
def stdev(self, unbiased=True):
return self.variance(unbiased=unbiased).sqrt()
def to_(self, device):
if self._mean is not None:
self._mean = self._mean.to(device)
if self.v_cmom2 is not None:
self.v_cmom2 = self.v_cmom2.to(device)
def load_state_dict(self, state):
self.count = state['count']
self.batchcount = state['batchcount']
self._mean = torch.from_numpy(state['mean'])
self.v_cmom2 = torch.from_numpy(state['cmom2'])
self.data_shape = None if state['data_shape'] is None else tuple(state['data_shape'])
def state_dict(self):
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
count=self.count,
data_shape=self.data_shape and tuple(self.data_shape),
batchcount=self.batchcount,
mean=self._mean.cpu().numpy(),
cmom2=self.v_cmom2.cpu().numpy())
class Covariance(Stat):
'''
Running computation. Use this when the entire covariance matrix is needed,
and when the whole covariance matrix fits in the GPU.
Chan-style numerically stable update of mean and full covariance matrix.
Chan, Golub. LeVeque. 1983. http://www.jstor.org/stable/2683386
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self._mean = None
self.cmom2 = None
self.data_shape = None
def add(self, a):
a = self._normalize_add_shape(a)
if len(a) == 0:
return
batch_count = a.shape[0]
# Initial batch.
if self._mean is None:
self.count = batch_count
self._mean = a.sum(0) / batch_count
centered = a - self._mean
self.cmom2 = centered.t().mm(centered)
return
# Update a batch using Chan-style update for numerical stability.
self.count += batch_count
# Update the mean according to the batch deviation from the old mean.
delta = a - self._mean
self._mean.add_(delta.sum(0) / self.count)
delta2 = a - self._mean
# Update the variance using the batch deviation
self.cmom2.addmm_(mat1=delta.t(), mat2=delta2)
def to_(self, device):
if self._mean is not None:
self._mean = self._mean.to(device)
if self.cmom2 is not None:
self.cmom2 = self.cmom2.to(device)
def mean(self):
return self._restore_result_shape(self._mean)
def covariance(self, unbiased=True):
return self._restore_result_shape(
self.cmom2 / (self.count - (1 if unbiased else 0)))
def correlation(self, unbiased=True):
cov = self.cmom2 / (self.count - (1 if unbiased else 0))
rstdev = cov.diag().sqrt().reciprocal()
return self._restore_result_shape(
rstdev[:, None] * covariance * rstdev[None, :])
def variance(self, unbiased=True):
return self._restore_result_shape(
self.cmom2.diag() / (self.count - (1 if unbiased else 0)))
def stdev(self, unbiased=True):
return self.variance(unbiased=unbiased).sqrt()
def state_dict(self):
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
count=self.count,
data_shape=self.data_shape and tuple(self.data_shape),
mean=self._mean.cpu().numpy(),
cmom2=self.cmom2.cpu().numpy())
def load_state_dict(self, state):
self.count = state['count']
self._mean = torch.from_numpy(state['mean'])
self.cmom2 = torch.from_numpy(state['cmom2'])
self.data_shape = None if state['data_shape'] is None else tuple(state['data_shape'])
class SecondMoment(Stat):
'''
Running computation. Use this when the entire non-centered 2nd-moment
'covariance-like' matrix is needed, and when the whole matrix fits
in the GPU.
'''
def __init__(self, split_batch=True, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self.mom2 = None
self.split_batch = split_batch
def add(self, a):
a = self._normalize_add_shape(a)
if len(a) == 0:
return
# Initial batch reveals the shape of the data.
if self.count == 0:
self.mom2 = a.new(a.shape[1], a.shape[1]).zero_()
batch_count = a.shape[0]
# Update the covariance using the batch deviation
self.count += batch_count
self.mom2 += a.t().mm(a)
def to_(self, device):
if self.mom2 is not None:
self.mom2 = self.mom2.to(device)
def moment(self):
return self.mom2 / self.count
def state_dict(self):
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
count=self.count,
mom2=self.mom2.cpu().numpy())
def load_state_dict(self, state):
self.count = int(state['count'])
self.mom2 = torch.from_numpy(state['mom2'])
class Bincount(Stat):
'''
Running bincount. The counted array should be an integer type with
non-negative integers.
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self._bincount = None
def add(self, a, size=None):
a = a.view(-1)
bincount = a.bincount()
if self._bincount is None:
self._bincount = bincount
elif len(self._bincount) < len(bincount):
bincount[:len(self._bincount)] += self._bincount
self._bincount = bincount
else:
self._bincount[:len(bincount)] += bincount
if size is None:
self.count += len(a)
else:
self.count += size
def to_(self, device):
self._bincount = self._bincount.to(device)
def size(self):
return self.count
def bincount(self):
return self._bincount
def state_dict(self):
return dict(
constructor=self.__module__ + '.' +
self.__class__.__name__ + '()',
count=self.count,
bincount=self._bincount.cpu().numpy())
def load_state_dict(self, dic):
self.count = int(dic['count'])
self._bincount = torch.from_numpy(dic['bincount'])
class CrossCovariance(Stat):
'''
Covariance. Use this when an off-diagonal block of the covariance
matrix is needed (e.g., when the whole covariance matrix does
not fit in the GPU, this could use a quarter of the memory).
Chan-style numerically stable update of mean and full covariance matrix.
Chan, Golub. LeVeque. 1983. http://www.jstor.org/stable/2683386
'''
def __init__(self, split_batch=True, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self._mean = None
self.cmom2 = None
self.v_cmom2 = None
self.split_batch = split_batch
def add(self, a, b):
if len(a.shape) == 1:
a = a[None, :]
b = b[None, :]
assert(a.shape[0] == b.shape[0])
if len(a.shape) > 2:
a, b = [d.view(d.shape[0], d.shape[1], -1).permute(0, 2, 1)
.reshape(-1, d.shape[1]) for d in [a, b]]
batch_count = a.shape[0]
# Initial batch.
if self._mean is None:
self.count = batch_count
self._mean = [d.sum(0) / batch_count for d in [a, b]]
centered = [d - bm for d, bm in zip([a, b], self._mean)]
self.v_cmom2 = [c.pow(2).sum(0) for c in centered]
self.cmom2 = centered[0].t().mm(centered[1])
return
# Update a batch using Chan-style update for numerical stability.
self.count += batch_count
# Update the mean according to the batch deviation from the old mean.
delta = [(d - bm) for d, bm in zip([a, b], self._mean)]
for m, d in zip(self._mean, delta):
m.add_(d.sum(0) / self.count)
delta2 = [(d - bm) for d, bm in zip([a, b], self._mean)]
# Update the cross-covariance using the batch deviation
self.cmom2.addmm_(mat1=delta[0].t(), mat2=delta2[1])
# Update the variance using the batch deviation
for vc2, d, d2 in zip(self.v_cmom2, delta, delta2):
vc2.add_((d * d2).sum(0))
def mean(self):
return self._mean
def variance(self, unbiased=True):
return [vc2 / (self.count - (1 if unbiased else 0))
for vc2 in self.v_cmom2]
def stdev(self, unbiased=True):
return [v.sqrt() for v in self.variance(unbiased=unbiased)]
def covariance(self, unbiased=True):
return self.cmom2 / (self.count - (1 if unbiased else 0))
def correlation(self):
covariance = self.covariance(unbiased=False)
rstdev = [s.reciprocal() for s in self.stdev(unbiased=False)]
cor = rstdev[0][:, None] * covariance * rstdev[1][None, :]
# Remove NaNs
cor[torch.isnan(cor)] = 0
return cor
def to_(self, device):
self._mean = [m.to(device) for m in self._mean]
self.v_cmom2 = [vcs.to(device) for vcs in self.v_cmom2]
self.cmom2 = self.cmom2.to(device)
def state_dict(self):
return dict(
constructor=self.__module__ + '.' +
self.__class__.__name__ + '()',
count=self.count,
mean_a=self._mean[0].cpu().numpy(),
mean_b=self._mean[1].cpu().numpy(),
cmom2_a=self.v_cmom2[0].cpu().numpy(),
cmom2_b=self.v_cmom2[1].cpu().numpy(),
cmom2=self.cmom2.cpu().numpy())
def load_state_dict(self, state):
self.count = int(state['count'])
self._mean = [torch.from_numpy(state[f'mean_{k}']) for k in 'ab']
self.v_cmom2 = [torch.from_numpy(state[f'cmom2_{k}']) for k in 'ab']
self.cmom2 = torch.from_numpy(state['cmom2'])
def _float_from_bool(a):
'''
Since pytorch only supports matrix multiplication on float,
IoU computations are done using floating point types.
This function binarizes the input (positive to True and
nonpositive to False), and converts from bool to float.
If the data is already a floating-point type, it leaves
it keeps the same type; otherwise it uses float.
'''
if a.dtype == torch.bool:
return a.float()
if a.dtype.is_floating_point:
return a.sign().clamp_(0)
return (a > 0).float()
class IoU(Stat):
'''
Running computation of intersections and unions of all features.
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self._intersection = None
def add(self, a):
assert len(a.shape) == 2
a = _float_from_bool(a)
if self._intersection is None:
self._intersection = torch.mm(a.t(), a)
else:
self._intersection.addmm_(a.t(), a)
self.count += len(a)
def size(self):
return self.count
def intersection(self):
return self._intersection
def union(self):
total = self._intersection.diagonal(0)
return total[:, None] + total[None, :] - self._intersection
def iou(self):
return self.intersection() / (self.union() + 1e-20)
def to_(self, _device):
self._intersection = self._intersection.to(_device)
def state_dict(self):
return dict(constructor=self.__module__ + '.' +
self.__class__.__name__ + '()',
count=self.count,
intersection=self._intersection.cpu().numpy())
def load_state_dict(self, state):
self.count = int(state['count'])
self._intersection = torch.tensor(state['intersection'])
class CrossIoU(Stat):
'''
Running computation of intersections and unions of two binary vectors.
'''
def __init__(self, state=None):
if state is not None:
return super().__init__(state)
self.count = 0
self._intersection = None
self.total_a = None
self.total_b = None
def add(self, a, b):
assert len(a.shape) == 2 and len(b.shape) == 2
assert len(a) == len(b), f'{len(a)} vs {len(b)}'
a = _float_from_bool(a) # CUDA only supports mm on float...
b = _float_from_bool(b) # otherwise we would use integers.
intersection = torch.mm(a.t(), b)
asum = a.sum(0)
bsum = b.sum(0)
if self._intersection is None:
self._intersection = intersection
self.total_a = asum
self.total_b = bsum
else:
self._intersection += intersection
self.total_a += asum
self.total_b += bsum
self.count += len(a)
def size(self):
return self.count
def intersection(self):
return self._intersection
def union(self):
return self.total_a[:, None] + self.total_b[None, :] - self._intersection
def iou(self):
return self.intersection() / (self.union() + 1e-20)
def to_(self, _device):
self.total_a = self.total_a.to(_device)
self.total_b = self.total_b.to(_device)
self._intersection = self._intersection.to(_device)
def state_dict(self):
return dict(constructor=self.__module__ + '.' +
self.__class__.__name__ + '()',
count=self.count,
total_a=self.total_a.cpu().numpy(),
total_b=self.total_b.cpu().numpy(),
intersection=self._intersection.cpu().numpy())
def load_state_dict(self, state):
self.count = int(state['count'])
self.total_a = torch.tensor(state['total_a'])
self.total_b = torch.tensor(state['total_b'])
self._intersection = torch.tensor(state['intersection'])
class Quantile(Stat):
'''
Streaming randomized quantile computation for torch.
Add any amount of data repeatedly via add(data). At any time,
quantile estimates be read out using quantile(q).
Implemented as a sorted sample that retains at least r samples
(by default r = 3072); the number of retained samples will grow to
a finite ceiling as the data is accumulated. Accuracy scales according
to r: the default is to set resolution to be accurate to better than about
0.1%, while limiting storage to about 50,000 samples.
Good for computing quantiles of huge data without using much memory.
Works well on arbitrary data with probability near 1.
Based on the optimal KLL quantile algorithm by Karnin, Lang, and Liberty
from FOCS 2016. http://ieee-focs.org/FOCS-2016-Papers/3933a071.pdf
'''
def __init__(self, r=3 * 1024, buffersize=None, seed=None, state=None):
if state is not None:
return super().__init__(state)
self.depth = None
self.dtype = None
self.device = None
resolution = r * 2 # sample array is at least half full before discard
self.resolution = resolution
# Default buffersize: 128 samples (and smaller than resolution).
if buffersize is None:
buffersize = min(128, (resolution + 7) // 8)
self.buffersize = buffersize
self.samplerate = 1.0
self.data = None
self.firstfree = [0]
self.randbits = torch.ByteTensor(resolution)
self.currentbit = len(self.randbits) - 1
self.extremes = None
self.count = 0
self.batchcount = 0
def size(self):
return self.count
def _lazy_init(self, incoming):
self.depth = incoming.shape[1]
self.dtype = incoming.dtype
self.device = incoming.device
self.data = [
torch.zeros(
self.depth, self.resolution, dtype=self.dtype, device=self.device
)
]
self.extremes = torch.zeros(self.depth, 2, dtype=self.dtype, device=self.device)
self.extremes[:, 0] = float('inf')
self.extremes[:, -1] = -float('inf')
def to_(self, device):
'''Switches internal storage to specified device.'''
if device != self.device:
old_data = self.data
old_extremes = self.extremes
self.data = [d.to(device) for d in self.data]
self.extremes = self.extremes.to(device)
self.device = self.extremes.device
del old_data
del old_extremes
def add(self, incoming):
if self.depth is None:
self._lazy_init(incoming)
assert len(incoming.shape) == 2
assert incoming.shape[1] == self.depth, (incoming.shape[1], self.depth)
self.count += incoming.shape[0]
self.batchcount += 1
# Convert to a flat torch array.
if self.samplerate >= 1.0:
self._add_every(incoming)
return
# If we are sampling, then subsample a large chunk at a time.
self._scan_extremes(incoming)
chunksize = int(math.ceil(self.buffersize / self.samplerate))
for index in range(0, len(incoming), chunksize):
batch = incoming[index : index + chunksize]
sample = sample_portion(batch, self.samplerate)
if len(sample):
self._add_every(sample)
def _add_every(self, incoming):
supplied = len(incoming)
index = 0
while index < supplied:
ff = self.firstfree[0]
available = self.data[0].shape[1] - ff
if available == 0:
if not self._shift():
# If we shifted by subsampling, then subsample.
incoming = incoming[index:]
if self.samplerate >= 0.5:
# First time sampling - the data source is very large.
self._scan_extremes(incoming)
incoming = sample_portion(incoming, self.samplerate)
index = 0
supplied = len(incoming)
ff = self.firstfree[0]
available = self.data[0].shape[1] - ff
copycount = min(available, supplied - index)
self.data[0][:, ff : ff + copycount] = torch.t(
incoming[index : index + copycount, :]
)
self.firstfree[0] += copycount
index += copycount
def _shift(self):
index = 0
# If remaining space at the current layer is less than half prev
# buffer size (rounding up), then we need to shift it up to ensure
# enough space for future shifting.
while self.data[index].shape[1] - self.firstfree[index] < (
-(-self.data[index - 1].shape[1] // 2) if index else 1
):
if index + 1 >= len(self.data):
return self._expand()
data = self.data[index][:, 0 : self.firstfree[index]]
data = data.sort()[0]
if index == 0 and self.samplerate >= 1.0:
self._update_extremes(data[:, 0], data[:, -1])
offset = self._randbit()
position = self.firstfree[index + 1]
subset = data[:, offset::2]
self.data[index + 1][:, position : position + subset.shape[1]] = subset
self.firstfree[index] = 0
self.firstfree[index + 1] += subset.shape[1]
index += 1
return True
def _scan_extremes(self, incoming):
# When sampling, we need to scan every item still to get extremes
self._update_extremes(
torch.min(incoming, dim=0)[0], torch.max(incoming, dim=0)[0]
)
def _update_extremes(self, minr, maxr):
self.extremes[:, 0] = torch.min(
torch.stack([self.extremes[:, 0], minr]), dim=0
)[0]
self.extremes[:, -1] = torch.max(
torch.stack([self.extremes[:, -1], maxr]), dim=0
)[0]
def _randbit(self):
self.currentbit += 1
if self.currentbit >= len(self.randbits):
self.randbits.random_(to=2)
self.currentbit = 0
return self.randbits[self.currentbit]
def state_dict(self):
state = dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
resolution=self.resolution,
depth=self.depth,
buffersize=self.buffersize,
samplerate=self.samplerate,
sizes=numpy.array([d.shape[1] for d in self.data]),
extremes=self.extremes.cpu().detach().numpy(),
size=self.count,
batchcount=self.batchcount,
)
for i, (d, f) in enumerate(zip(self.data, self.firstfree)):
state[f'data.{i}'] = d.cpu().detach().numpy()[:, :f].T
return state
def load_state_dict(self, state):
self.resolution = int(state['resolution'])
self.randbits = torch.ByteTensor(self.resolution)
self.currentbit = len(self.randbits) - 1
self.depth = int(state['depth'])
self.buffersize = int(state['buffersize'])
self.samplerate = float(state['samplerate'])
firstfree = []
buffers = []
for i, s in enumerate(state['sizes']):
d = state[f'data.{i}']
firstfree.append(d.shape[0])
buf = numpy.zeros((d.shape[1], s), dtype=d.dtype)
buf[:, : d.shape[0]] = d.T
buffers.append(torch.from_numpy(buf))
self.firstfree = firstfree
self.data = buffers
self.extremes = torch.from_numpy((state['extremes']))
self.count = int(state['size'])
self.batchcount = int(state.get('batchcount', 0))
self.dtype = self.extremes.dtype
self.device = self.extremes.device
def min(self):
return self.minmax()[0]
def max(self):
return self.minmax()[-1]
def minmax(self):
if self.firstfree[0]:
self._scan_extremes(self.data[0][:, : self.firstfree[0]].t())
return self.extremes.clone()
def median(self):
return self.quantiles(0.5)
def mean(self):
return self.integrate(lambda x: x) / self.count
def variance(self, unbiased=True):
mean = self.mean()[:, None]
return self.integrate(lambda x:
(x - mean).pow(2)) / (self.count - (1 if unbiased else 0))
def stdev(self, unbiased=True):
return self.variance(unbiased=unbiased).sqrt()
def _expand(self):
cap = self._next_capacity()
if cap > 0:
# First, make a new layer of the proper capacity.
self.data.insert(
0, torch.zeros(self.depth, cap, dtype=self.dtype, device=self.device)
)
self.firstfree.insert(0, 0)
else:
# Unless we're so big we are just subsampling.
assert self.firstfree[0] == 0
self.samplerate *= 0.5
for index in range(1, len(self.data)):
# Scan for existing data that needs to be moved down a level.
amount = self.firstfree[index]
if amount == 0:
continue
position = self.firstfree[index - 1]
# Move data down if it would leave enough empty space there
# This is the key invariant: enough empty space to fit half
# of the previous level's buffer size (rounding up)
if self.data[index - 1].shape[1] - (amount + position) >= (
-(-self.data[index - 2].shape[1] // 2) if (index - 1) else 1
):
self.data[index - 1][:, position : position + amount] = self.data[
index
][:, :amount]
self.firstfree[index - 1] += amount
self.firstfree[index] = 0
else:
# Scrunch the data if it would not.
data = self.data[index][:, :amount]
data = data.sort()[0]
if index == 1:
self._update_extremes(data[:, 0], data[:, -1])
offset = self._randbit()
scrunched = data[:, offset::2]
self.data[index][:, : scrunched.shape[1]] = scrunched
self.firstfree[index] = scrunched.shape[1]
return cap > 0
def _next_capacity(self):
cap = int(math.ceil(self.resolution * (0.67 ** len(self.data))))
if cap < 2:
return 0
# Round up to the nearest multiple of 8 for better GPU alignment.
cap = -8 * (-cap // 8)
return max(self.buffersize, cap)
def _weighted_summary(self, sort=True):
if self.firstfree[0]:
self._scan_extremes(self.data[0][:, : self.firstfree[0]].t())
size = sum(self.firstfree)
weights = torch.FloatTensor(size) # Floating point
summary = torch.zeros(self.depth, size, dtype=self.dtype, device=self.device)
index = 0
for level, ff in enumerate(self.firstfree):
if ff == 0:
continue
summary[:, index : index + ff] = self.data[level][:, :ff]
weights[index : index + ff] = 2.0 ** level
index += ff
assert index == summary.shape[1]
if sort:
summary, order = torch.sort(summary, dim=-1)
weights = weights[order.view(-1).cpu()].view(order.shape)
summary = torch.cat(
[self.extremes[:, :1], summary, self.extremes[:, 1:]], dim=-1
)
weights = torch.cat(
[
torch.zeros(weights.shape[0], 1),
weights,
torch.zeros(weights.shape[0], 1),
],
dim=-1,
)
return (summary, weights)
def quantiles(self, quantiles):
if not hasattr(quantiles, 'cpu'):
quantiles = torch.tensor(quantiles)
qshape = quantiles.shape
if self.count == 0:
return torch.full((self.depth,) + qshape, torch.nan)
summary, weights = self._weighted_summary()
cumweights = torch.cumsum(weights, dim=-1) - weights / 2
cumweights /= torch.sum(weights, dim=-1, keepdim=True)
result = torch.zeros(
self.depth, quantiles.numel(), dtype=self.dtype, device=self.device
)
# numpy is needed for interpolation
nq = quantiles.view(-1).cpu().detach().numpy()
ncw = cumweights.cpu().detach().numpy()
nsm = summary.cpu().detach().numpy()
for d in range(self.depth):
result[d] = torch.tensor(
numpy.interp(nq, ncw[d], nsm[d]), dtype=self.dtype, device=self.device
)
return result.view((self.depth,) + qshape)
def integrate(self, fun):
result = []
for level, ff in enumerate(self.firstfree):
if ff == 0:
continue
result.append(torch.sum(fun(self.data[level][:, :ff]) * (2.0 ** level), dim=-1))
if len(result) == 0:
return None
return torch.stack(result).sum(dim=0) / self.samplerate
def readout(self, count=1001):
return self.quantiles(torch.linspace(0.0, 1.0, count))
def normalize(self, data):
'''
Given input data as taken from the training distirbution,
normalizes every channel to reflect quantile values,
uniformly distributed, within [0, 1].
'''
assert self.count > 0
assert data.shape[0] == self.depth
summary, weights = self._weighted_summary()
cumweights = torch.cumsum(weights, dim=-1) - weights / 2
cumweights /= torch.sum(weights, dim=-1, keepdim=True)
result = torch.zeros_like(data).float()
# numpy is needed for interpolation
ndata = data.cpu().numpy().reshape((data.shape[0], -1))
ncw = cumweights.cpu().numpy()
nsm = summary.cpu().numpy()
for d in range(self.depth):
normed = torch.tensor(
numpy.interp(ndata[d], nsm[d], ncw[d]),
dtype=torch.float,
device=data.device,
).clamp_(0.0, 1.0)
if len(data.shape) > 1:
normed = normed.view(*(data.shape[1:]))
result[d] = normed
return result
class TopK:
'''
A class to keep a running tally of the the top k values (and indexes)
of any number of torch feature components. Will work on the GPU if
the data is on the GPU. Tracks largest by default, but tracks smallest
if largest=False is passed.
This version flattens all arrays to avoid crashes.
'''
def __init__(self, k=100, largest=True, state=None):
if state is not None:
return super().__init__(state)
self.k = k
self.count = 0
# This version flattens all data internally to 2-d tensors,
# to avoid crashes with the current pytorch topk implementation.
# The data is puffed back out to arbitrary tensor shapes on ouput.
self.data_shape = None
self.top_data = None
self.top_index = None
self.next = 0
self.linear_index = 0
self.perm = None
self.largest = largest
def add(self, data, index=None):
'''
Adds a batch of data to be considered for the running top k.
The zeroth dimension enumerates the observations. All other
dimensions enumerate different features.
'''
if self.top_data is None:
# Allocation: allocate a buffer of size 5*k, at least 10, for each.
self.data_shape = data.shape[1:]
feature_size = int(numpy.prod(self.data_shape))
self.top_data = torch.zeros(
feature_size, max(10, self.k * 5), out=data.new())
self.top_index = self.top_data.clone().long()
self.linear_index = 0 if len(data.shape) == 1 else torch.arange(
feature_size, out=self.top_index.new()).mul_(
self.top_data.shape[-1])[:, None]
size = data.shape[0]
sk = min(size, self.k)
if self.top_data.shape[-1] < self.next + sk:
# Compression: if full, keep topk only.
self.top_data[:, :self.k], self.top_index[:, :self.k] = (
self.topk(sorted=False, flat=True))
self.next = self.k
free = self.top_data.shape[-1] - self.next
# Pick: copy the top sk of the next batch into the buffer.
# Currently strided topk is slow. So we clone after transpose.
# TODO: remove the clone() if it becomes faster.
cdata = data.reshape(size, numpy.prod(data.shape[1:])).t().clone()
td, ti = cdata.topk(sk, sorted=False, largest=self.largest)
self.top_data[:, self.next:self.next + sk] = td
if index is not None:
ti = index[ti]
else:
ti = ti + self.count
self.top_index[:, self.next:self.next + sk] = ti
self.next += sk
self.count += size
def size(self):
return self.count
def topk(self, sorted=True, flat=False):
'''
Returns top k data items and indexes in each dimension,
with channels in the first dimension and k in the last dimension.
'''
k = min(self.k, self.next)
# bti are top indexes relative to buffer array.
td, bti = self.top_data[:, :self.next].topk(k,
sorted=sorted, largest=self.largest)
# we want to report top indexes globally, which is ti.
ti = self.top_index.view(-1)[
(bti + self.linear_index).view(-1)
].view(*bti.shape)
if flat:
return td, ti
else:
return (td.view(*(self.data_shape + (-1,))),
ti.view(*(self.data_shape + (-1,))))
def to_(self, device):
if self.top_data is not None:
self.top_data = self.top_data.to(device)
if self.top_index is not None:
self.top_index = self.top_index.to(device)
if isinstance(self.linear_index, torch.Tensor):
self.linear_index = self.linear_index.to(device)
def state_dict(self):
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
k=self.k,
count=self.count,
largest=self.largest,
data_shape=self.data_shape and tuple(self.data_shape),
top_data=self.top_data.cpu().detach().numpy(),
top_index=self.top_index.cpu().detach().numpy(),
next=self.next,
linear_index=(self.linear_index.cpu().numpy()
if isinstance(self.linear_index, torch.Tensor)
else self.linear_index),
perm=self.perm)
def load_state_dict(self, state):
self.k = int(state['k'])
self.count = int(state['count'])
self.largest = bool(state.get('largest', True))
self.data_shape = None if state['data_shape'] is None else tuple(state['data_shape'])
self.top_data = torch.from_numpy(state['top_data'])
self.top_index = torch.from_numpy(state['top_index'])
self.next = int(state['next'])
self.linear_index = (torch.from_numpy(state['linear_index'])
if len(state['linear_index'].shape) > 0
else int(state['linear_index']))
class History(Stat):
'''
Accumulates the concatenation of all the added data.
'''
def __init__(self, data=None, state=None):
if state is not None:
return super().__init__(state)
self._data = data
self._added = []
def _cat_added(self):
if len(self._added):
self._data = torch.cat(
([self._data] if self._data is not None else [])
+ self._added)
self._added = []
def add(self, d):
self._added.append(d)
if len(self._added) > 100:
self._cat_added()
def history(self):
self._cat_added()
return self._data
def load_state_dict(self, state):
data = state['data']
self._data = None if data is None else torch.from_numpy(data)
self._added = []
def state_dict(self):
self._cat_added()
return dict(
constructor=self.__module__ + '.' + self.__class__.__name__ + '()',
data=None if self._data is None else self._data.cpu().numpy())
def to_(self, device):
'''Switches internal storage to specified device.'''
self._cat_added()
if self._data is not None:
self._data = self._data.to(device)
class CombinedStat(Stat):
'''
An Stat that bundles together multiple Stat objects.
Convenient for loading and saving a state_dict made up of a
hierarchy of stats, and for use with the tally() function.
Example:
cs = CombinedStat(m=Mean(), q=Quantile())
for [b] in tally(cs, MyDataSet(), cache=fn, batch_size=100):
cs.add(b)
print(cs.m.mean())
print(cs.q.median())
'''
def __init__(self, state=None, **kwargs):
self._objs = kwargs
if state is not None:
return super().__init__(state)
def __getattr__(self, k):
if k in self._objs:
return self._objs[k]
raise AttributeError()
def add(self, d, *args, **kwargs):
for obj in self._objs.values():
obj.add(d, *args, **kwargs)
def load_state_dict(self, state):
for prefix, obj in self._objs.items():
obj.load_state_dict(pull_key_prefix(prefix, state))
def state_dict(self):
result = {}
for prefix, obj in self._objs.items():
result.update(push_key_prefix(prefix, obj.state_dict()))
return result
def to_(self, device):
'''Switches internal storage to specified device.'''
for v in self._objs.values():
v.to_(device)
def push_key_prefix(prefix, d):
'''
Returns a dict with the same values as d, but where each key
adds the prefix, followed by a dot.
'''
return {prefix + '.' + k: v for k, v in d.items()}
def pull_key_prefix(prefix, d):
'''
Returns a filtered dict of all the items of d that start with
the given key prefix, plus a dot, with that prefix removed.
'''
pd = prefix + '.'
lpd = len(pd)
return {k[lpd:]: v for k, v in d.items() if k.startswith(pd)}
# We wish to be able to save None (null) values in numpy npz files,
# yet do so without setting the unsecure 'allow_pickle' flag. To do
# that, we will encode null as a special kind of IEEE 754 NaN value.
# Inspired by https://github.com/zuiderkwast/nanbox/blob/master/nanbox.h
# we follow the same Nanboxing scheme used in JavaScriptCore
# (search for JSCJSValue.h#L435), which encodes null values in NaN
# as the NaN value with hex pattern 0xfff8000000000002.
null_numpy_value = numpy.array(
struct.unpack('>d', struct.pack('>Q', 0xfff8000000000002))[0],
dtype=numpy.float64)
def is_null_numpy_value(v):
return (isinstance(v, numpy.ndarray) and numpy.ndim(v) == 0 and
v.dtype == numpy.float64 and numpy.isnan(v) and
0xfff8000000000002 == struct.unpack('>Q', struct.pack('>d', v))[0])
def box_numpy_null(d):
try:
return {k: box_numpy_null(v) for k, v in d.items()}
except:
return null_numpy_value if d is None else d
def unbox_numpy_null(d):
try:
return {k: unbox_numpy_null(v) for k, v in d.items()}
except:
return None if is_null_numpy_value(d) else d
def resolve_state_dict(s):
'''
Resolves a state, which can be a filename or a dict-like object.
'''
if isinstance(s, str):
return unbox_numpy_null(numpy.load(s))
return s
global_load_cache_enabled = True
def load_cached_state(cachefile, args, quiet=False, throw=False):
'''
Resolves a state, which can be a filename or a dict-like object.
'''
if not global_load_cache_enabled or cachefile is None:
return None
try:
if isinstance(cachefile, dict):
dat = cachefile
cachefile = 'state' # for printed messages
else:
dat = unbox_numpy_null(numpy.load(cachefile))
for a, v in args.items():
if a not in dat or dat[a] != v:
pbar.print('%s %s changed from %s to %s' % (
cachefile, a, dat[a], v))
return None
except (FileNotFoundError, ValueError) as e:
if throw:
raise e
return None
else:
if not quiet:
print('Loading cached %s' % cachefile)
return dat
def save_cached_state(cachefile, obj, args):
if cachefile is None:
return
dat = obj.state_dict()
for a, v in args.items():
if a in dat:
assert (dat[a] == v)
dat[a] = v
if isinstance(cachefile, dict):
cachefile.clear()
cachefile.update(dat)
else:
os.makedirs(os.path.dirname(cachefile), exist_ok=True)
numpy.savez(cachefile, **box_numpy_null(dat))
class FixedSubsetSampler(Sampler):
'''Represents a fixed sequence of data set indices.
Subsets can be created by specifying a subset of output indexes.
'''
def __init__(self, samples):
self.samples = samples
def __iter__(self):
return iter(self.samples)
def __len__(self):
return len(self.samples)
def __getitem__(self, key):
return self.samples[key]
def subset(self, new_subset):
return FixedSubsetSampler(self.dereference(new_subset))
def dereference(self, indices):
'''
Translate output sample indices (small numbers indexing the sample)
to input sample indices (larger number indexing the original full set)
'''
return [self.samples[i] for i in indices]
class FixedRandomSubsetSampler(FixedSubsetSampler):
'''Samples a fixed number of samples from the dataset, deterministically.
Arguments:
data_source,
sample_size,
seed (optional)
'''
def __init__(self, data_source, start=None, end=None, seed=1):
rng = random.Random(seed)
shuffled = list(range(len(data_source)))
rng.shuffle(shuffled)
self.data_source = data_source
super(FixedRandomSubsetSampler, self).__init__(shuffled[start:end])
def class_subset(self, class_filter):
'''
Returns only the subset matching the given rule.
'''
if isinstance(class_filter, int):
def rule(d): return d[1] == class_filter
else:
rule = class_filter
return self.subset([i for i, j in enumerate(self.samples)
if rule(self.data_source[j])])
def make_loader(dataset, sample_size=None, batch_size=1, sampler=None,
random_sample=None, **kwargs):
'''Utility for creating a dataloader on fixed sample subset.'''
import typing
if isinstance(dataset, typing.Callable):
# To support deferred dataset loading, support passing a factory
# that creates the dataset when called.
dataset = dataset()
if isinstance(dataset, torch.Tensor):
# The dataset can be a simple tensor.
dataset = torch.utils.data.TensorDataset(dataset)
if sample_size is not None:
assert sampler is None, 'sampler cannot be specified with sample_size'
if sample_size > len(dataset):
pbar.print('Warning: sample size %d > dataset size %d' %
(sample_size, len(dataset)))
sample_size = len(dataset)
if random_sample is None:
sampler = FixedSubsetSampler(list(range(sample_size)))
else:
sampler = FixedRandomSubsetSampler(
dataset, seed=random_sample, end=sample_size)
return torch.utils.data.DataLoader(
dataset,
sampler=sampler,
batch_size=batch_size,
**kwargs)
# Unit Tests
def _unit_test():
import warnings
warnings.filterwarnings('error')
import time
import argparse
import tempfile
import shutil
import random
parser = argparse.ArgumentParser(
description='Test things out')
parser.add_argument('--mode', default='cpu', help='cpu or cuda')
parser.add_argument('--test_size', type=int, default=1000000)
args = parser.parse_args()
testdir = tempfile.mkdtemp()
batch_size = random.randint(500, 1500)
# Test NaNboxing.
assert numpy.isnan(null_numpy_value)
assert is_null_numpy_value(null_numpy_value)
assert not is_null_numpy_value(numpy.nan)
# Test Covariance
goal = torch.tensor(numpy.random.RandomState(1).standard_normal(
10 * 10)).view(10, 10)
data = torch.tensor(numpy.random.RandomState(2).standard_normal(
args.test_size * 10)).view(args.test_size, 10).mm(goal)
data += torch.randn(1, 10) * 999
dcov = data.t().cov()
rcov = Covariance()
rcov.add(data) # All one batch
assert((rcov.covariance() - dcov).abs().max() < 1e-16)
cs = CombinedStat(cov=Covariance(), xcov=CrossCovariance())
ds = torch.utils.data.TensorDataset(data)
for [a] in tally(cs, ds, batch_size=9876):
cs.cov.add(a)
cs.xcov.add(a[:,:3], a[:,3:])
assert((data.mean(0) - cs.cov.mean()).abs().max() < 1e-12)
assert((dcov - cs.cov.covariance()).abs().max() < 2e-12)
assert((dcov[:3,3:] - cs.xcov.covariance()).abs().max() < 1e-12)
assert((dcov.diagonal() - torch.cat(
cs.xcov.variance())).abs().max() < 1e-12)
# Test CrossCovariance and CrossIoU
fn = f'{testdir}/cross_cache.npz'
ds = torch.utils.data.TensorDataset(
(torch.arange(args.test_size)[:,None] %
torch.arange(1, 6)[None, :] == 0).double(),
(torch.arange(args.test_size)[:,None] %
torch.arange(5, 8)[None, :] == 0).double(),
)
c = CombinedStat(c=CrossCovariance(), iou=CrossIoU())
riou = IoU()
count = 0
for [a, b] in tally(c, ds, cache=fn, batch_size=100):
count += 1
c.add(a, b)
riou.add(torch.cat([a, b], dim=1))
assert count == -(-args.test_size // 100)
cor = c.c.correlation()
iou = c.iou.iou()
assert cor.shape == iou.shape == (5, 3)
assert iou[4, 0] == 1.0
assert abs(iou[0, 2] + (-args.test_size // 7 / float(args.test_size))) < 1e-6
assert abs(cor[4, 0] - 1.0) < 1e-2
assert abs(cor[0, 2] - 0.0) < 1e-6
assert all((riou.iou()[:5,-3:] == iou).view(-1))
assert all(riou.iou().diagonal(0) == 1)
c = CombinedStat(c=CrossCovariance(), iou=CrossIoU())
count = 0
for [a, b] in tally(c, ds, cache=fn, batch_size=10):
count += 1
c.add(a, b)
assert count == 0
assert all((c.c.correlation() == cor).view(-1))
assert all((c.iou.iou() == iou).view(-1))
# Test Concatantaion, Mean, Bincount and tally.
fn = f'{testdir}/series_cache.npz'
count = 0
ds = torch.utils.data.TensorDataset(torch.arange(args.test_size))
c = CombinedStat(s=History(), m=Mean(), b=Bincount())
for [b] in tally(c, ds, cache=fn, batch_size=batch_size):
count += 1
c.add(b)
assert count == -(-args.test_size // batch_size)
assert len(c.s.history()) == args.test_size
assert c.s.history()[-1] == args.test_size - 1
assert all(c.s.history() == ds.tensors[0])
assert all(c.b.bincount() == torch.ones(args.test_size))
assert c.m.mean() == float(args.test_size - 1) / 2.0
c2 = CombinedStat(s=History(), m=Mean(), b=Bincount())
batches = tally(c2, ds, cache=fn)
assert len(c2.s.history()) == args.test_size
assert all(c2.s.history() == c.s.history())
assert all(c2.b.bincount() == torch.ones(args.test_size))
assert c2.m.mean() == c.m.mean()
count = 0
for b in batches:
count += 1
assert count == 0 # Shouldn't do anything when it's cached
# An adverarial case: we keep finding more numbers in the middle
# as the stream goes on.
amount = args.test_size
quantiles = 1000
data = numpy.arange(float(amount))
data[1::2] = data[-1::-2] + (len(data) - 1)
data /= 2
depth = 50
test_cuda = torch.cuda.is_available()
alldata = data[:, None] + (numpy.arange(depth) * amount)[None, :]
actual_sum = torch.FloatTensor(numpy.sum(alldata * alldata, axis=0))
amt = amount // depth
for r in range(depth):
numpy.random.shuffle(alldata[r * amt:r * amt + amt, r])
if args.mode == 'cuda':
alldata = torch.cuda.FloatTensor(alldata)
dtype = torch.float
device = torch.device('cuda')
else:
alldata = torch.FloatTensor(alldata)
dtype = torch.float
device = None
starttime = time.time()
cs = CombinedStat(
qc=Quantile(),
m=Mean(),
v=Variance(),
c=Covariance(),
s=SecondMoment(),
t=TopK(),
i=IoU())
# Feed data in little batches
i = 0
while i < len(alldata):
batch_size = numpy.random.randint(1000)
cs.add(alldata[i:i+batch_size])
i += batch_size
# Test state dict
saved = cs.state_dict()
# numpy.savez(f'{testdir}/saved.npz', **box_numpy_null(saved))
# saved = unbox_numpy_null(numpy.load(f'{testdir}/saved.npz'))
cs.save(f'{testdir}/saved.npz')
loaded = unbox_numpy_null(numpy.load(f'{testdir}/saved.npz'))
assert set(loaded.keys()) == set(saved.keys())
# Restore using state=saved in constructor.
cs2 = CombinedStat(
qc=Quantile(),
m=Mean(),
v=Variance(),
c=Covariance(),
s=SecondMoment(),
t=TopK(),
i=IoU(),
state=saved)
# saved = unbox_numpy_null(numpy.load(f'{testdir}/saved.npz'))
assert not cs2.qc.device.type == 'cuda'
cs2.to_(device)
# alldata = alldata.cpu()
cs2.add(alldata)
actual_sum *= 2
# print(abs(alldata.mean(0) - cs2.m.mean()) / alldata.mean())
assert all(abs(alldata.mean(0) - cs2.m.mean()) / alldata.mean() < 1e-5)
assert all(abs(alldata.mean(0) - cs2.v.mean()) / alldata.mean() < 1e-5)
assert all(abs(alldata.mean(0) - cs2.c.mean()) / alldata.mean() < 1e-5)
# print(abs(alldata.var(0) - cs2.v.variance()) / alldata.var(0))
assert all(abs(alldata.var(0) - cs2.v.variance()) / alldata.var(0) < 1e-3)
assert all(abs(alldata.var(0) - cs2.c.variance()) / alldata.var(0) < 1e-2)
# print(abs(alldata.std(0) - cs2.v.stdev()) / alldata.std(0))
assert all(abs(alldata.std(0) - cs2.v.stdev()) / alldata.std(0) < 1e-4)
# print(abs(alldata.std(0) - cs2.c.stdev()) / alldata.std(0))
assert all(abs(alldata.std(0) - cs2.c.stdev()) / alldata.std(0) < 2e-3)
moment = (alldata.t() @ alldata) / len(alldata)
# print(abs(moment - cs2.s.moment()) / moment.abs())
assert all((abs(moment - cs2.s.moment()) / moment.abs()).view(-1) < 1e-2)
assert all(alldata.max(dim=0)[0] == cs2.t.topk()[0][:,0])
assert cs2.i.iou()[0, 0] == 1
assert all((cs2.i.iou()[1:, 1:] == 1).view(-1))
assert all(cs2.i.iou()[1:, 0] < 1)
assert all(cs2.i.iou()[1:, 0] == cs2.i.iou()[0, 1:])
# Restore using cs.load() method.
cs = CombinedStat(
qc=Quantile(),
m=Mean(),
v=Variance(),
c=Covariance(),
s=SecondMoment(),
t=TopK(),
i=IoU())
cs.load(f'{testdir}/saved.npz')
assert not cs.qc.device.type == 'cuda'
cs.to_(device)
cs.add(alldata)
# actual_sum *= 2
# print(abs(alldata.mean(0) - cs.m.mean()) / alldata.mean())
assert all(abs(alldata.mean(0) - cs.m.mean()) / alldata.mean() < 1e-5)
assert all(abs(alldata.mean(0) - cs.v.mean()) / alldata.mean() < 1e-5)
assert all(abs(alldata.mean(0) - cs.c.mean()) / alldata.mean() < 1e-5)
# print(abs(alldata.var(0) - cs.v.variance()) / alldata.var(0))
assert all(abs(alldata.var(0) - cs.v.variance()) / alldata.var(0) < 1e-3)
assert all(abs(alldata.var(0) - cs.c.variance()) / alldata.var(0) < 1e-2)
# print(abs(alldata.std(0) - cs.v.stdev()) / alldata.std(0))
assert all(abs(alldata.std(0) - cs.v.stdev()) / alldata.std(0) < 1e-4)
# print(abs(alldata.std(0) - cs.c.stdev()) / alldata.std(0))
assert all(abs(alldata.std(0) - cs.c.stdev()) / alldata.std(0) < 2e-3)
moment = (alldata.t() @ alldata) / len(alldata)
# print(abs(moment - cs.s.moment()) / moment.abs())
assert all((abs(moment - cs.s.moment()) / moment.abs()).view(-1) < 1e-2)
assert all(alldata.max(dim=0)[0] == cs.t.topk()[0][:,0])
assert cs.i.iou()[0, 0] == 1
assert all((cs.i.iou()[1:, 1:] == 1).view(-1))
assert all(cs.i.iou()[1:, 0] < 1)
assert all(cs.i.iou()[1:, 0] == cs.i.iou()[0, 1:])
# Randomized quantile test
qc = cs.qc
ro = qc.readout(1001).cpu()
endtime = time.time()
gt = torch.linspace(0, amount, quantiles + 1)[None, :] + (
torch.arange(qc.depth, dtype=torch.float) * amount)[:, None]
maxreldev = torch.max(torch.abs(ro - gt) / amount) * quantiles
print('Randomized quantile test results:')
print('Maximum relative deviation among %d perentiles: %f' % (
quantiles, maxreldev))
minerr = torch.max(torch.abs(qc.minmax().cpu()[:, 0] -
torch.arange(qc.depth, dtype=torch.float) * amount))
maxerr = torch.max(torch.abs((qc.minmax().cpu()[:, -1] + 1) -
(torch.arange(qc.depth, dtype=torch.float) + 1) * amount))
print('Minmax error %f, %f' % (minerr, maxerr))
interr = torch.max(torch.abs(qc.integrate(lambda x: x * x).cpu()
- actual_sum) / actual_sum)
print('Integral error: %f' % interr)
medianerr = torch.max(torch.abs(qc.median() -
alldata.median(0)[0]) / alldata.median(0)[0]).cpu()
print('Median error: %f' % interr)
meanerr = torch.max(
torch.abs(qc.mean() - alldata.mean(0)) / alldata.mean(0)).cpu()
print('Mean error: %f' % meanerr)
varerr = torch.max(
torch.abs(qc.variance() - alldata.var(0)) / alldata.var(0)).cpu()
print('Variance error: %f' % varerr)
counterr = ((qc.integrate(lambda x: torch.ones(x.shape[-1]).cpu())
- qc.size()) / (0.0 + qc.size())).item()
print('Count error: %f' % counterr)
print('Time %f' % (endtime - starttime))
# Algorithm is randomized, so some of these will fail with low probability.
assert maxreldev < 1.0
assert minerr == 0.0
assert maxerr == 0.0
assert interr < 0.01
assert abs(counterr) < 0.001
shutil.rmtree(testdir, ignore_errors=True)
print('OK')
if __name__ == '__main__':
_unit_test()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment