diff --git a/.gitmodules b/.gitmodules index aea68447b..f5f807184 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,7 +1,3 @@ -[submodule "externals/marray"] - path = externals/marray - url = https://github.com/DerThorsten/marray - [submodule "externals/pybind11"] path = externals/pybind11 url = https://github.com/pybind/pybind11 diff --git a/CMakeLists.txt b/CMakeLists.txt index 193685675..3bbf54a2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules) project(nifty) set (NIFTY_VERSION_MAJOR 0) -set (NIFTY_VERSION_MINOR 14) +set (NIFTY_VERSION_MINOR 15) set (NIFTY_VERSION_PATCH 0) @@ -124,7 +124,7 @@ endif() #------------------------------------------------------------------------------------------------------------------- include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/qpbo") -include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/marray/include") +#include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/marray/include") include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/pybind11/include") include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/vigra/include") include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/externals/graph/include") diff --git a/deprecated/cremi/malas.py b/deprecated/cremi/malas.py deleted file mode 100644 index 3eed9c953..000000000 --- a/deprecated/cremi/malas.py +++ /dev/null @@ -1,48 +0,0 @@ -import nifty -import nifty.graph -import nifty.graph.agglo - -import numpy - - - -G = nifty.graph.UndirectedGraph -nagglo = nifty.graph.agglo - - -def testMala(): - - - - s = 1000 - g = nifty.graph.UndirectedGraph(s*s) - - def node(x,y): - return y+s*x - - for x in range(s): - for y in range(s): - if x+1 < s: - g.insertEdge(node(x,y),node(x+1,y)) - if y+1 < s: - g.insertEdge(node(x,y),node(x,y+1)) - - - - edgeIndicators = numpy.random.rand(g.numberOfEdges) - edgeSizes = numpy.ones(shape=[g.numberOfEdges]) - nodeSizes = numpy.ones(shape=[g.numberOfNodes]) - - clusterPolicy = nagglo.malaClusterPolicy( - graph=g, edgeIndicators=edgeIndicators, - edgeSizes=edgeSizes, nodeSizes=nodeSizes,verbose=True) - - - agglomerativeClustering = nagglo.agglomerativeClustering(clusterPolicy) - agglomerativeClustering.run(verbose=True) - - seg = agglomerativeClustering.result()#out=[1,2,3,4]) - - - -testMala() diff --git a/deprecated/cremi/train.py b/deprecated/cremi/train.py deleted file mode 100644 index 0196b23b1..000000000 --- a/deprecated/cremi/train.py +++ /dev/null @@ -1,99 +0,0 @@ -import nifty -import nifty.deep_learning - -import keras - -from nifty.deep_learning.targets import seg_batch_to_affinity_edges,cremi_z5_edges -from nifty.deep_learning.data_loader import MultiresCremiLoader -from nifty.deep_learning.models.homer_net import HomerNet -from nifty.deep_learning.loss.pixelwise_weighted_loss import LiftedEdgeLoss - -import numpy -import random - -# numpy.random.seed(7) -# random.seed(7) - -params = { - 'cremi_data_folder' : - # raw folder - '/media/tbeier/4D11469249B96057/work/cremi/data/raw_padded', - - # data loader - 'data_loader' : { - 'neuron_ids' : True, - 'clefts' : False, - 'patch_sizes' : [200,200,200], - 'n_z' : 5, - 'z_slice_ranges' : [[0,15 ], [50,55], [110,120]], - 'bad_slices' : [[],[],[14]] - }, - - # the net - 'net' : { - 'n_resolutions' : 3, - 'input_patch_shapes' : [[200,200,5]]*3 - } - -} - - -train_batch_size = 1 - -# the data source -loader = MultiresCremiLoader(root_folder=params['cremi_data_folder'],**params['data_loader']) - -# the targets -edges = cremi_z5_edges(r=1, atrous_rate_xy=1, add_local_edges=False) -edge_priors = numpy.ones(len(edges))*0.2 - -# the training data gen -def gen(): - while True: - raw_batch_list, neuron_ids_batch = loader(batch_size=train_batch_size) - gt_and_weights = seg_batch_to_affinity_edges(neuron_ids_batch, edges=edges, edge_priors=edge_priors) - yield raw_batch_list, gt_and_weights - - -for x,y in gen(): - break -# make the loss function(s) -input_shape = [200,200,len(edges)*2] -loss = LiftedEdgeLoss(input_shape=input_shape, n_channels=len(edges)) - - - -# make the net -net_param = params['net'] -homer_net = HomerNet(loss_function=loss, **params['net']) - - - -class Trainer(object): - def __init__(self, net, loss_functions): - - self.net = net - self.loss_functions = loss_functions - self.model = self.net.model - opt = keras.optimizers.Adam() - self.model.compile(optimizer=opt, loss=self.loss_functions) - - def train(self, train_gen, valid_gen,train_batch_size, n_val=30, epochs=100000, s_mult=50, callbacks=None): - - internal_callbacks = [] - - self.model.fit_generator(generator=train_gen(), - steps_per_epoch=s_mult*train_batch_size, - epochs=epochs, - validation_data=valid_gen(), - validation_steps=n_val, - verbose=True, - callbacks=internal_callbacks - ) - - -trainer = Trainer(net=homer_net, loss_functions=loss) -trainer.train(train_gen=gen, valid_gen=gen, train_batch_size=train_batch_size) - - - diff --git a/deprecated/learnit.py b/deprecated/learnit.py deleted file mode 100644 index 9c6ccdfe4..000000000 --- a/deprecated/learnit.py +++ /dev/null @@ -1,165 +0,0 @@ -from __future__ import print_function -from __future__ import division - -import os -import h5py -import threading - -import numpy -import vigra - -import nifty -import nifty.hdf5 -import nifty.graph.rag - - - -blockLabelsIlpFile = "/home/tbeier/block_labels.ilp" - -rawFile = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/datasets/hhess_2nm/raw_2k_2nm.h5" -pmapFile = "/home/tbeier/prediction_semantic_binary_full.h5" -oversegFile = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxels/2nm/aggloseg_0.3_50.h5" -oversegFile = "/home/tbeier/Desktop/aggloseg_0.3_50.h5" -workDir = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/opt2k" - -def getLabels(blockLabelsIlpFile): - - allLabels = dict() - blockLabelsIlpFileH5 = h5py.File(blockLabelsIlpFile,'r') - edgeLabels = blockLabelsIlpFileH5['Training and Multicut']['EdgeLabelsDict'] - - for key in edgeLabels.keys(): - el = edgeLabels[key] - labels = el['labels'][:] - spIds = el['sp_ids'][:,:] - - for l, spId in zip(labels, spIds): - if l == 1 or l == 2: - uv = long(spId[0]),long(spId[1]) - if uv in allLabels: - allLabels[uv] = max(l, allLabels[uv]) - else: - allLabels[uv] = l - blockLabelsIlpFileH5.close() - - - l = numpy.array(allLabels.values()) - uv = numpy.array(allLabels.keys()) - return l-1, uv - - - -def h5Max(dset, blockShape=None): - - shape = dset.shape - ndim = len(shape) - assert ndim == 3 - if blockShape is None: - blockShape = [100]*len(shape) - - blocking = nifty.tools.blocking(roiBegin=(0,0,0), roiEnd=shape, blockShape=list(blockShape)) - - - lock = threading.Lock() - done = [0] - gMax = [0] - - with nifty.tools.progressBar(size=blocking.numberOfBlocks) as bar: - - def f(blockIndex): - - - - - block = blocking.getBlock(blockIndex) - b,e = block.begin, block.end - - - d = dset[b[0]:e[0], b[1]:e[1], b[2]:e[2]] - - with lock: - done[0] += 1 - bar.update(done[0]) - - gMax[0] = max(d.max(),gMax[0]) - - nifty.tools.parallelForEach(range(blocking.numberOfBlocks), f=f, nWorkers=1) - - return gMax[0] - - - -def computeRag(oversegFile, dset): - - oversegFileH5 = h5py.File(oversegFile) - oversegDset = oversegFileH5[dset] - try: - maxLabel = oversegDset.attrs['maxLabel'] - #print("load max label") - except: - print("compute max label") - maxLabel = h5Max(oversegDset) - oversegDset.attrs['maxLabel'] = maxLabel - - oversegFileH5.close() - - - - print("max label",maxLabel) - - - nLabels = maxLabel + 1 - - cs = nifty.hdf5.CacheSettings() - cs.hashTabelSize = 977 - cs.nBytes = 990000000 - cs.rddc = 0.5 - - h5File = nifty.hdf5.openFile(oversegFile, cs) - labelsArray = nifty.hdf5.Hdf5ArrayUInt32(h5File, dset) - - ragFile = os.path.join(workDir,'rag.h5') - if not os.path.isfile(ragFile): - - with vigra.Timer("rag"): - gridRag = nifty.graph.rag.gridRagHdf5(labelsArray, nLabels, blockShape=[150,150,150], numberOfThreads=20) - - print("serialize") - serialization = gridRag.serialize() - - print("save serialization") - ragH5File = h5py.File(ragFile) - ragH5File['data'] = serialization - ragH5File.close() - else: - ragH5File = h5py.File(ragFile,'r') - serialization = ragH5File['data'] - gridRag = nifty.graph.rag.gridRagHdf5(labelsArray, nLabels, serialization=serialization) - ragH5File.close() - - - return gridRag - - -uv,l = getLabels(blockLabelsIlpFile) -oversegFileH5 = h5py.File(oversegFile) - - - -# rag -rag = computeRag(oversegFile,'data') - - - -h5File = nifty.hdf5.openFile(rawFile,) -rawArray = nifty.hdf5.Hdf5ArrayUInt8(h5File, 'data') - - -# accumulate features -rawFeat = nifty.graph.rag.accumulateStandartFeatures(rag=rag, data=rawArray, minVal=0.0, maxVal=255.0, blockShape=[75,75,75], numberOfThreads=-1) - - - - - - diff --git a/deprecated/norag.py b/deprecated/norag.py deleted file mode 100644 index 4056f1357..000000000 --- a/deprecated/norag.py +++ /dev/null @@ -1,540 +0,0 @@ -from __future__ import print_function -from __future__ import division - -from filters import * -import fastfilters -import numpy -import h5py -import os -import vigra -import nifty -import threading -import nifty.pipelines.neuro_seg as nseg -import nifty.graph.agglo as nagglo - -from get_ilp_labels import * - - - - -def singleBlockPipeline(rawFile, pmapFile, oversegFile, ilpFile, outFolder, outName, sigmas=(2.0,4.0,8.0), filtOnRaw=True, filtOnPmap=True, blockShape=(254,)*3, - minSegSize=None ,beta=0.5, stages=[1,1,1]): - - - - if not os.path.exists(outFolder): - os.makedirs(outFolder) - - - # save stuff - def saveRes(data,name): - f5 = h5py.File(os.path.join(outFolder,'%s_%s.h5'%(outName,name)),'w') - f5['data'] = data - f5.close() - - # load stuff - def loadRes(name): - f5 = h5py.File(os.path.join(outFolder,'%s_%s.h5'%(outName,name)),'r') - data = f5['data'][:] - f5.close() - return data - - - - - - - - - rawH5File = h5py.File(rawFile,'r') - osegH5File = h5py.File(oversegFile,'r') - pmapH5File = h5py.File(pmapFile,'r') - - - # dsets - rawDset = rawH5File['data'] - osegDset = osegH5File['data'] - pmapDset = pmapH5File['data'] - - print("raw",rawDset.shape) - print("oseg",osegDset.shape) - print("pmap",pmapDset.shape) - - shape = rawDset.shape - - #load labels - uvIdsTrain, labelsTrain = getIlpLabels(ilpFile) - - - blocking = nifty.tools.blocking(roiBegin=[0]*3, roiEnd=shape, - blockShape=blockShape) - nBlocks = blocking.numberOfBlocks - lock = threading.Lock() - lock2 = threading.Lock() - blockRes = [None]*nBlocks - - - if len(sigmas) > 0: - maxSigma = max(sigmas) - - - - - if stages[0] == 1: - print("STAGE 0") - filtList = [] - filtListPmap = [] - - bs = 100 - hs = [s//2 for s in shape] - raw = rawDset[hs[0] : hs[0] + bs, hs[1] : hs[1] + bs, hs[2] : hs[2] + bs].astype('float32') - pm = pmapDset[hs[0] : hs[0] + bs, hs[1] : hs[1] + bs, hs[2] : hs[2] + bs].astype('float32') - - for sigma in sigmas: - - if filtOnRaw: - filtList.append(GaussianGradientMagnitude(raw=raw, sigma=sigma)) - filtList.append(LaplacianOfGaussian(raw=raw, sigma=sigma)) - filtList.append(HessianOfGaussianEv(raw=raw, sigma=sigma)) - filtList.append(GaussianSmoothing(raw=raw, sigma=sigma)) - - if filtOnPmap: - filtListPmap.append(GaussianGradientMagnitude(raw=pm, sigma=sigma)) - filtListPmap.append(LaplacianOfGaussian(raw=pm, sigma=sigma)) - filtListPmap.append(HessianOfGaussianEv(raw=pm, sigma=sigma)) - filtListPmap.append(GaussianSmoothing(raw=pm, sigma=sigma)) - - - - - def f(blockIndex): - - # labels halo - blockWithBorder = blocking.getBlockWithHalo(blockIndex, [0,0,0],[1,1,1]) - outerBlock = blockWithBorder.outerBlock - bAcc, eAcc = outerBlock.begin, outerBlock.end - - - # filters - - filterRes= [] - - if len(sigmas) > 0 and (filtOnRaw or filtOnPmap): - halo = [int(round(3.5*maxSigma + 3.0))] * 3 - else: - halo = [0]*3 - - blockWithBorder = blocking.addHalo(outerBlock, halo) - outerBlock = blockWithBorder.outerBlock - innerBlockLocal = blockWithBorder.innerBlockLocal - bFilt, eFilt = outerBlock.begin, outerBlock.end - bFiltCore, eFiltCore = innerBlockLocal.begin, innerBlockLocal.end - - with lock: - rawFilt = rawDset[bFilt[0]:eFilt[0], bFilt[1]:eFilt[1], bFilt[2]:eFilt[2]].astype('float32') - pmFilt = pmapDset[bFilt[0]:eFilt[0], bFilt[1]:eFilt[1], bFilt[2]:eFilt[2]].astype('float32') - oseg = osegDset[bAcc[0]:eAcc[0], bAcc[1]:eAcc[1], bAcc[2]:eAcc[2]] - - raw = rawFilt[bFiltCore[0]:eFiltCore[0], bFiltCore[1]:eFiltCore[1], bFiltCore[2]:eFiltCore[2]]/255.0 - pmap = pmFilt[bFiltCore[0]:eFiltCore[0], bFiltCore[1]:eFiltCore[1], bFiltCore[2]:eFiltCore[2]]/255.0 - - - filterRes.append(raw[...,None]) - filterRes.append(pmap[...,None]) - - for filt in filtList: - res = filt(rawFilt) - res = res[bFiltCore[0]:eFiltCore[0], bFiltCore[1]:eFiltCore[1], bFiltCore[2]:eFiltCore[2],:] - filterRes.append(res) - - for filt in filtListPmap: - res = filt(pmFilt) - res = res[bFiltCore[0]:eFiltCore[0], bFiltCore[1]:eFiltCore[1], bFiltCore[2]:eFiltCore[2],:] - filterRes.append(res) - - - - - #for vd in filterRes: - # print(vd.shape,vd.dtype) - # - voxelData = numpy.concatenate(filterRes, axis=3) - voxelData = numpy.require(voxelData,requirements=['C_CONTIGUOUS'],dtype='float32') - - #print(voxelData.strides) - #last axis is continous - assert voxelData.strides[3] == 4 - - blockData = nseg.BlockData(blocking=blocking, blockIndex=blockIndex, - numberOfChannels=voxelData.shape[3], - numberOfBins=20) - - blockData.accumulate(oseg, voxelData) - - blockRes[blockIndex] = blockData - - nifty.tools.parallelForEach(iterable=range(nBlocks), f=f, - nWorkers=-1, showBar=True, - size=nBlocks, name="ComputeFeatures") - - - def mergePairwise(toMerge): - - while(len(toMerge) != 1): - newList = [] - - l = len(toMerge) - pairs = (l+1)//2 - print(l) - for p in range(pairs): - p0 = p*2 - p1 = p*2 + 1 - - if(p10.7)[0] - - - f5 = h5py.File(liftedProbsFile, 'w') - f5['data'] = probs - f5.close() - - -def runLiftedMc(dataDict, settings): - - - rawData = dataDict['raw'] - outDir = dataDict['outDir'] - localRfPredictDir = os.path.join(outDir,'localClfProbs') - ragsAndSuperpixels = getRagsAndSuperpixels(dataDict, settings) - liftedFeaturesDir = os.path.join(outDir,'liftedFeatures') - liftedProbsDir = os.path.join(outDir,'liftedProbs') - - ragsAndSuperpixels = getRagsAndSuperpixels(dataDict, settings) - - for sliceIndex in range(rawData.shape[0]): - - - # local probs - predictionFile = os.path.join(localRfPredictDir,'rag_pred_clf%s_%d.h5'%('0', sliceIndex)) - localProbs = h5Read(predictionFile)[:,1] - - # lifted probs - liftedProbsFile = os.path.join(liftedProbsDir,'lifted_probs_%d.h5'%(sliceIndex)) - liftedProbs = h5Read(liftedProbsFile) - - # set up the lifted objective - rag, sp = ragsAndSuperpixels[sliceIndex] - obj = nifty.graph.lifted_multicut.liftedMulticutObjective(rag) - liftedGraph = obj.liftedGraph - distance = obj.insertLiftedEdgesBfs(5, returnDistance=True).astype('float32') - liftedUvIds = obj.liftedUvIds() - - # numeric factor to not get to small weights - # might not be necessary - C = 100.0 - - # local weights - eps = 0.0001 - clipped = numpy.clip(localProbs, eps, 1.0-eps) - beta = settings['betaLocal'] - wLocal = numpy.log((1.0-clipped)/(clipped)) + numpy.log((1.0-beta)/(beta)) - # normalize by number of local edges length - wLocal *= C/len(wLocal) - wLocal *= settings['gamma'] - - # non local weights - eps = 0.0001 - clipped = numpy.clip(liftedProbs, eps, 1.0-eps) - beta = settings['betaLifted'] - wLifted = numpy.log((1.0-clipped)/(clipped)) + numpy.log((1.0-beta)/(beta)) - - - # normalize by the distance (give close one more weight) - distanceWeight = 1.0 / (distance-1.0) - #wLifted *= C/distanceWeight.sum() - wLifted *= C/len(wLifted) - - print numpy.abs(wLocal).sum(), numpy.abs(wLifted).sum() - - # write the weighs into the objective - obj.setLiftedEdgesCosts(wLifted, overwrite=True) - obj.setGraphEdgesCosts(wLocal, overwrite=True) - - - - # warm start with normal multicut - mcObj = nifty.graph.multicut.multicutObjective(rag, wLocal) - solverFactory = mcObj.multicutIlpCplexFactory() - solver = solverFactory.create(mcObj) - visitor = mcObj.multicutVerboseVisitor() - argMc = solver.optimize(visitor) - emc = obj.evalNodeLabels(argMc) - - - # finaly optimize it - solverFactory = obj.liftedMulticutAndresGreedyAdditiveFactory() - solver = solverFactory.create(obj) - visitor = obj.verboseVisitor() - arg = solver.optimize(visitor) - eg = obj.evalNodeLabels(arg) - - - solverFactory = obj.liftedMulticutAndresKernighanLinFactory() - solver = solverFactory.create(obj) - visitor = obj.verboseVisitor() - arg = solver.optimize(visitor, arg.copy()) - ekl = obj.evalNodeLabels(arg) - - - print "e",emc,eg,ekl - - #projectToPixels - pixelData = nifty.graph.rag.projectScalarNodeDataToPixels(rag, arg.astype('uint32')) - - vigra.segShow(rawData[sliceIndex,:,:], pixelData) - vigra.show() diff --git a/deprecated/segmentation_pipeline/gt.py b/deprecated/segmentation_pipeline/gt.py deleted file mode 100644 index 65ab64a84..000000000 --- a/deprecated/segmentation_pipeline/gt.py +++ /dev/null @@ -1,78 +0,0 @@ -import vigra -import nifty -import nifty.graph -import nifty.graph.rag -import nifty.graph.agglo -import nifty.ground_truth -import numpy -import h5py - -from reraise import * - - - -@reraise_with_stack -def makeGt(rawData, binaryGt, settings): - - - seeds = vigra.analysis.labelImageWithBackground(binaryGt) - - edgeIndicator = vigra.filters.hessianOfGaussianEigenvalues(rawData, 3.0)[:,:,0] - seg, nseg = vigra.analysis.watershedsNew(edgeIndicator, seeds=seeds) - - - return seg - - - - -def getTrainingData(rag, sp, pixelGt, features, settings): - - thresholds = settings['fuztGtThreshold'] - t0,t1 = thresholds - - # compute overlap - overlap = nifty.ground_truth.Overlap(rag.numberOfNodes-1, - sp, pixelGt - ) - - fuzzyEdgeGt = overlap.differentOverlaps(rag.uvIds()) - - - where0 = numpy.where(fuzzyEdgeGtt1) - - f0 = features[where0] - f1 = features[where1] - - feat = numpy.concatenate([f0,f1],axis=0) - labels = numpy.ones(feat.shape[0],dtype='uint32') - labels[0:f0.shape[0]] = 0 - - return feat,labels - - -def getLiftedTrainingData(rag, sp, pixelGt, liftedUvIds, features, settings): - - thresholds = settings['fuztGtThreshold'] - t0,t1 = thresholds - - # compute overlap - overlap = nifty.ground_truth.Overlap(rag.numberOfNodes-1, - sp, pixelGt - ) - - fuzzyEdgeGt = overlap.differentOverlaps(liftedUvIds) - - - where0 = numpy.where(fuzzyEdgeGtt1) - - f0 = features[where0] - f1 = features[where1] - - feat = numpy.concatenate([f0,f1],axis=0) - labels = numpy.ones(feat.shape[0],dtype='uint32') - labels[0:f0.shape[0]] = 0 - - return feat,labels \ No newline at end of file diff --git a/deprecated/segmentation_pipeline/lifted_features.py b/deprecated/segmentation_pipeline/lifted_features.py deleted file mode 100644 index d07d668c6..000000000 --- a/deprecated/segmentation_pipeline/lifted_features.py +++ /dev/null @@ -1,299 +0,0 @@ -import vigra -import nifty -import nifty.graph -import nifty.graph.rag -import nifty.graph.agglo -import numpy -import h5py -import sys - -nrag = nifty.graph.rag -nagglo = nifty.graph.agglo - -from reraise import * -from tools import * - - - - -@reraise_with_stack -def multicutFromLocalProbs(raw, rag, localProbs, liftedEdges): - - - u = liftedEdges[:,0] - v = liftedEdges[:,1] - - # accumulate length (todo, implement function to just accumulate length) - eFeatures, nFeatures = nrag.accumulateMeanAndLength(rag, raw, [100,100],1) - eSizes = eFeatures[:,1] - eps = 0.0001 - clipped = numpy.clip(localProbs, eps, 1.0-eps) - - features = [] - - for beta in (0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7): - - for powers in (0.0, 0.005, 0.1, 0.15, 0.2): - - # the weight of the weight - wWeight = eSizes**powers - - print "\n\nBETA ",beta - w = numpy.log((1.0-clipped)/(clipped)) + numpy.log((1.0-beta)/(beta)) * wWeight - obj = nifty.graph.multicut.multicutObjective(rag, w) - - - factory = obj.multicutIlpCplexFactory() - solver = factory.create(obj) - visitor = obj.multicutVerboseVisitor() - #ret = solver.optimize(visitor=visitor) - ret = solver.optimize() - - res = ret[u] != ret[v] - - features.append(res[:,None]) - - features = numpy.concatenate(features, axis=1) - mean = numpy.mean(features, axis=1)[:,None] - features = numpy.concatenate([features, mean], axis=1) - - return features - - -@reraise_with_stack -def ucmFromLocalProbs(raw, rag, localProbs, liftedEdges, liftedObj): - - - u = liftedEdges[:,0] - v = liftedEdges[:,1] - - # accumulate length (todo, implement function to just accumulate length) - eFeatures, nFeatures = nrag.accumulateMeanAndLength(rag, raw, [100,100],1) - eSizes = eFeatures[:,1] - nSizes = eFeatures[:,1] - - - feats = nifty.graph.lifted_multicut.liftedUcmFeatures( - objective=liftedObj, - edgeIndicators=localProbs, - edgeSizes=eSizes, - nodeSizes=nSizes, - sizeRegularizers=[0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, - 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, - 0.9, 0.95] - ) - - return feats - - - -@reraise_with_stack -def ucmFromHessian(raw, rag, liftedEdges, liftedObj): - - - u = liftedEdges[:,0] - v = liftedEdges[:,1] - - feats = [] - - - - for sigma in [1.0, 2.0, 3.0, 4.0, 5.0]: - pf = vigra.filters.hessianOfGaussianEigenvalues(raw, sigma)[:,:,0] - eFeatures, nFeatures = nrag.accumulateMeanAndLength(rag, pf, [100,100],1) - edgeIndicator = eFeatures[:,0] - eSizes = eFeatures[:,1] - nSizes = eFeatures[:,1] - - featsB = nifty.graph.lifted_multicut.liftedUcmFeatures( - objective=liftedObj, - edgeIndicators=edgeIndicator, - edgeSizes=eSizes, - nodeSizes=nSizes, - sizeRegularizers=[0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, - 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, - 0.9, 0.95] - ) - - feats.append(featsB) - - # different axis ordering as usual - return numpy.concatenate(feats,axis=0) - - - - - - - - -@reraise_with_stack -def liftedFeaturesFromLocalProbs(raw, rag, localProbs, liftedEdges, liftedObj, featureFile): - - - mcFeatureFile = featureFile + "mc.h5" - if not hasH5File(mcFeatureFile): - mcFeatures = multicutFromLocalProbs(raw=raw, rag=rag, localProbs=localProbs, - liftedEdges=liftedEdges) - - f5 = h5py.File(mcFeatureFile, 'w') - f5['data'] = mcFeatures - f5.close() - else: - mcFeatures = h5Read(mcFeatureFile) - - - ucmFeatureFile = featureFile + "ucm.h5" - if not hasH5File(ucmFeatureFile): - ucmFeatures = ucmFromLocalProbs(raw=raw, rag=rag, localProbs=localProbs, - liftedEdges=liftedEdges, - liftedObj=liftedObj) - - f5 = h5py.File(ucmFeatureFile, 'w') - f5['data'] = ucmFeatures - f5.close() - else: - ucmFeatures = h5Read(ucmFeatureFile) - - - - # combine - features = numpy.concatenate([mcFeatures, ucmFeatures.swapaxes(0,1)],axis=1) - - f5 = h5py.File(featureFile, 'w') - f5['data'] = features - f5.close() - - - - - -@reraise_with_stack -def accumlatedLiftedFeatures(raw, pmap, rag, liftedEdges, liftedObj): - - - uv = liftedEdges - u = uv[:,0] - v = uv[:,1] - - - # geometric edge features - geometricFeaturs = nifty.graph.rag.accumulateGeometricNodeFeatures(rag, - blockShape=[75, 75], - numberOfThreads=1) - - nodeSize = geometricFeaturs[:,0] - nodeCenter = geometricFeaturs[:,1:2] - nodeAxisA = geometricFeaturs[:,2:4] - nodeAxisA = geometricFeaturs[:,4:6] - - diff = (nodeCenter[u,:] - nodeCenter[v,:])**2 - diff = numpy.sum(diff,axis=1) - - - - - allEdgeFeat = [ - # sizes - (nodeSize[u] + nodeSize[v])[:,None], - (numpy.abs(nodeSize[u] - nodeSize[v]))[:,None], - (nodeSize[u] * nodeSize[v])[:,None], - (numpy.minimum(nodeSize[u] , nodeSize[v]))[:,None], - (numpy.maximum(nodeSize[u] , nodeSize[v]))[:,None], - diff[:,None] - ] - - pixelFeats = [ - raw[:,:,None], - ] - if pmap is not None: - pixelFeats.append(pmap[:,:,None]) - - for sigma in (1.0, 2.0, 4.0): - pf = [ - vigra.filters.hessianOfGaussianEigenvalues(raw, 1.0*sigma), - vigra.filters.structureTensorEigenvalues(raw, 1.0*sigma, 2.0*sigma), - vigra.filters.gaussianGradientMagnitude(raw, 1.0*sigma)[:,:,None], - vigra.filters.gaussianSmoothing(raw, 1.0*sigma)[:,:,None] - ] - pixelFeats.extend(pf) - - if pmap is not None: - pixelFeats.append(vigra.filters.gaussianSmoothing(pmap, 1.0*sigma)[:,:,None]) - - pixelFeats = numpy.concatenate(pixelFeats, axis=2) - - - - for i in range(pixelFeats.shape[2]): - - pixelFeat = pixelFeats[:,:,i] - - nodeFeat = nifty.graph.rag.accumulateNodeStandartFeatures( - rag=rag, data=pixelFeat.astype('float32'), - minVal=pixelFeat.min(), - maxVal=pixelFeat.max(), - blockShape=[75, 75], - numberOfThreads=10 - ) - - uFeat = nodeFeat[u,:] - vFeat = nodeFeat[v,:] - - - fList =[ - uFeat + vFeat, - uFeat * vFeat, - numpy.abs(uFeat-vFeat), - numpy.minimum(uFeat,vFeat), - numpy.maximum(uFeat,vFeat), - ] - - edgeFeat = numpy.concatenate(fList, axis=1) - allEdgeFeat.append(edgeFeat) - - allEdgeFeat = numpy.concatenate(allEdgeFeat, axis=1) - - return allEdgeFeat - - - -@reraise_with_stack -def liftedFeatures(raw, pmap, rag, liftedEdges, liftedObj, distances, featureFile): - - - - - ucmFeatureFile = featureFile + "ucm.h5" - if not hasH5File(ucmFeatureFile): - ucmFeatures = ucmFromHessian(raw=raw, rag=rag, - liftedEdges=liftedEdges, - liftedObj=liftedObj) - - f5 = h5py.File(ucmFeatureFile, 'w') - f5['data'] = ucmFeatures - f5.close() - else: - ucmFeatures = h5Read(ucmFeatureFile) - - - accFeatureFile = featureFile + "acc.h5" - if not hasH5File(accFeatureFile): - accFeatrues = accumlatedLiftedFeatures(raw=raw, pmap=pmap, - rag=rag, - liftedEdges=liftedEdges, - liftedObj=liftedObj) - - f5 = h5py.File(accFeatureFile, 'w') - f5['data'] = accFeatrues - f5.close() - else: - accFeatrues = h5Read(accFeatureFile) - - - # combine - features = numpy.concatenate([accFeatrues,distances[:,None], ucmFeatures.swapaxes(0,1)],axis=1) - - f5 = h5py.File(featureFile, 'w') - f5['data'] = features - f5.close() \ No newline at end of file diff --git a/deprecated/segmentation_pipeline/rag.py b/deprecated/segmentation_pipeline/rag.py deleted file mode 100644 index 343c01b75..000000000 --- a/deprecated/segmentation_pipeline/rag.py +++ /dev/null @@ -1,206 +0,0 @@ -import vigra -import nifty -import nifty.graph -import nifty.graph.rag -import nifty.graph.agglo -import numpy -import h5py -import sys - -nrag = nifty.graph.rag -nagglo = nifty.graph.agglo - -from reraise import * - -@reraise_with_stack -def makeRag(overseg, ragFile, settings): - - assert overseg.min() == 0 - rag = nifty.graph.rag.gridRag(overseg) - - # serialize the rag - serialization = rag.serialize() - - # save serialization - f5 = h5py.File(ragFile, 'w') - f5['rag'] = serialization - - # save superpixels - f5['superpixels'] = overseg - - f5.close() - - - - -def loadRag(ragFile): - # read serialization - f5 = h5py.File(ragFile, 'r') - serialization = f5['rag'][:] - - # load superpixels - overseg = f5['superpixels'][:] - - f5.close() - - rag = nifty.graph.rag.gridRag(overseg, serialization=serialization) - return rag, overseg - - - - -@reraise_with_stack -def localRagFeatures(raw, pmap, overseg, rag, featuresFile, settings): - - features = [] - - ucmFeat = ucmFeatures(raw=raw, pmap=pmap, - overseg=overseg,rag=rag, - settings=settings) - features.append(ucmFeat) - - accFeat = accumulatedFeatures(raw=raw, pmap=pmap, - overseg=overseg,rag=rag, - settings=settings) - features.append(accFeat) - - features = numpy.concatenate(features,axis=1) - - # save the features - f5 = h5py.File(featuresFile, 'w') - f5['data'] = features - f5.close() - - -@reraise_with_stack -def ucmFeatures(raw, pmap, overseg, rag, settings): - - features = [] - - def ucmTransform(edgeIndicator, sizeRegularizer): - vFeat = numpy.require(edgeIndicator, dtype='float32',requirements='C') - eFeatures, nFeatures = nifty.graph.rag.accumulateMeanAndLength(rag, vFeat, [100,100],1) - eMeans = eFeatures[:,0] - eSizes = eFeatures[:,1] - nSizes = nFeatures[:,1] - clusterPolicy = nagglo.edgeWeightedClusterPolicyWithUcm( - graph=rag, edgeIndicators=eMeans, - edgeSizes=eSizes, nodeSizes=nSizes, - numberOfNodesStop=1, - sizeRegularizer=float(sizeRegularizer) - ) - - #print numpy.abs(eMeans-eMeansCp).sum() - #assert not numpy.array_equal(eMeans, eMeansCp) - - agglomerativeClustering = nagglo.agglomerativeClustering(clusterPolicy) - mergeHeightR = agglomerativeClustering.runAndGetDendrogramHeight() - mergeHeight = agglomerativeClustering.ucmTransform(clusterPolicy.edgeIndicators) - mergeSize = agglomerativeClustering.ucmTransform(clusterPolicy.edgeSizes) - - - return mergeHeight,mergeHeightR,mergeSize - - for sigma in (1.0, 2.0, 4.0): - for reg in (0.001, 0.1, 0.2, 0.4, 0.8): - edgeIndicator = vigra.filters.hessianOfGaussianEigenvalues(raw, sigma)[:,:,0] - - A,B,C = ucmTransform(edgeIndicator,reg) - - features.append(A[:,None]) - features.append(B[:,None]) - features.append(C[:,None]) - - for sigma in (0.1, 1.0, 2.0): - for reg in (0.001, 0.1, 0.2, 0.4, 0.8): - - edgeIndicator = vigra.filters.gaussianSmoothing(pmap, sigma) - A,B,C = ucmTransform(edgeIndicator,reg) - - features.append(A[:,None]) - features.append(B[:,None]) - features.append(C[:,None]) - - return numpy.concatenate(features,axis=1) - - -@reraise_with_stack -def accumulatedFeatures(raw, pmap, overseg, rag, settings): - - #print "bincoutn", numpy.bincount(overseg.reshape([-1])).size,"nNodes",rag.numberOfNodes - - uv = rag.uvIds() - u = uv[:,0] - v = uv[:,1] - - - # geometric edge features - geometricFeaturs = nifty.graph.rag.accumulateGeometricEdgeFeatures(rag, - blockShape=[75, 75], - numberOfThreads=1) - - allEdgeFeat = [geometricFeaturs] - - pixelFeats = [ - raw[:,:,None], - ] - if pmap is not None: - pixelFeats.append(pmap[:,:,None]) - - for sigma in (1.0, 2.0, 4.0, 6.0 ,8.0): - pf = [ - vigra.filters.hessianOfGaussianEigenvalues(raw, 1.0*sigma), - vigra.filters.structureTensorEigenvalues(raw, 1.0*sigma, 2.0*sigma), - vigra.filters.gaussianGradientMagnitude(raw, 1.0*sigma)[:,:,None], - vigra.filters.gaussianSmoothing(raw, 1.0*sigma)[:,:,None] - ] - pixelFeats.extend(pf) - - if pmap is not None: - pixelFeats.append(vigra.filters.gaussianSmoothing(pmap, 1.0*sigma)[:,:,None]) - pixelFeats.append(vigra.filters.hessianOfGaussianEigenvalues(pmap, 1.0*sigma)), - pixelFeats.append(vigra.filters.structureTensorEigenvalues(pmap, 1.0*sigma, 2.0*sigma)) - - pixelFeats = numpy.concatenate(pixelFeats, axis=2) - - - - for i in range(pixelFeats.shape[2]): - - pixelFeat = pixelFeats[:,:,i] - - edgeFeat, nodeFeat = nifty.graph.rag.accumulateStandartFeatures( - rag=rag, data=pixelFeat.astype('float32'), - minVal=pixelFeat.min(), - maxVal=pixelFeat.max(), - blockShape=[75, 75], - numberOfThreads=10 - ) - - uFeat = nodeFeat[u,:] - vFeat = nodeFeat[v,:] - - du = numpy.abs(edgeFeat,uFeat) - dv = numpy.abs(edgeFeat,vFeat) - - - fList =[ - uFeat + vFeat, - uFeat * vFeat, - numpy.abs(uFeat-vFeat), - numpy.minimum(uFeat,vFeat), - numpy.maximum(uFeat,vFeat), - du + dv, - numpy.abs(du-dv), - numpy.minimum(du,dv), - numpy.maximum(du,dv), - edgeFeat - ] - - edgeFeat = numpy.concatenate(fList, axis=1) - allEdgeFeat.append(edgeFeat) - - allEdgeFeat = numpy.concatenate(allEdgeFeat, axis=1) - - return allEdgeFeat - diff --git a/deprecated/segmentation_pipeline/reraise.py b/deprecated/segmentation_pipeline/reraise.py deleted file mode 100644 index c8fcd33e1..000000000 --- a/deprecated/segmentation_pipeline/reraise.py +++ /dev/null @@ -1,16 +0,0 @@ -import functools -import traceback - - -def reraise_with_stack(func): - - @functools.wraps(func) - def wrapped(*args, **kwargs): - try: - return func(*args, **kwargs) - except Exception as e: - traceback_str = traceback.format_exc(e) - raise StandardError("Error occurred. Original traceback " - "is\n%s\n" % traceback_str) - - return wrapped \ No newline at end of file diff --git a/deprecated/segmentation_pipeline/superpixels.py b/deprecated/segmentation_pipeline/superpixels.py deleted file mode 100644 index 2d9925be5..000000000 --- a/deprecated/segmentation_pipeline/superpixels.py +++ /dev/null @@ -1,80 +0,0 @@ -import vigra -import nifty -import nifty.graph -import nifty.graph.rag -import nifty.graph.agglo -import numpy - - - -def makeSmallerSeg(oseg, edgeIndicator, reduceBy, sizeRegularizer): - import nifty - nrag = nifty.graph.rag - nagglo = nifty.graph.agglo - - # "overseg in c order starting at zero" - oseg = numpy.require(oseg, dtype='uint32',requirements='C') - oseg -= 1 - - # "make rag" - rag = nifty.graph.rag.gridRag(oseg) - - # "volfeatshape" - vFeat = numpy.require(edgeIndicator, dtype='float32',requirements='C') - - # "accumulate means and counts" - eFeatures, nFeatures = nrag.accumulateMeanAndLength(rag, vFeat, [100,100],1) - - eMeans = eFeatures[:,0] - eSizes = eFeatures[:,1] - nSizes = nFeatures[:,1] - - # "get clusterPolicy" - - numberOfNodesStop = int(float(rag.numberOfNodes)/float(reduceBy) + 0.5) - numberOfNodesStop = max(1,numberOfNodesStop) - numberOfNodesStop = min(rag.numberOfNodes, numberOfNodesStop) - - clusterPolicy = nagglo.edgeWeightedClusterPolicy( - graph=rag, edgeIndicators=eMeans, - edgeSizes=eSizes, nodeSizes=nSizes, - numberOfNodesStop=numberOfNodesStop, - sizeRegularizer=float(sizeRegularizer)) - - # "do clustering" - agglomerativeClustering = nagglo.agglomerativeClustering(clusterPolicy) - agglomerativeClustering.run() - seg = agglomerativeClustering.result()#out=[1,2,3,4]) - - # "make seg dense" - dseg = nifty.tools.makeDense(seg) - - # dseg.dtype, type(dseg) - - # "project to pixels" - pixelData = nrag.projectScalarNodeDataToPixels(rag, dseg.astype('uint32')) - # "done" - #pixelDataF = numpy.require(pixelData, dtype='uint32',requirements='F') - return pixelData - - - - - -def makeSupervoxels(rawData, settings): - - sigma = settings['spSigmaHessian'] - edgeIndicator = vigra.filters.hessianOfGaussianEigenvalues(rawData, sigma)[:,:,0] - seg, nseg = vigra.analysis.watershedsNew(edgeIndicator) - - # make smaller supervoxels via agglomerative clustering - seg = makeSmallerSeg(seg, edgeIndicator, settings['reduceBy'], settings['sizeRegularizer']) - - - - - - - - - return seg \ No newline at end of file diff --git a/deprecated/segmentation_pipeline/tools.py b/deprecated/segmentation_pipeline/tools.py deleted file mode 100644 index 185a56cc3..000000000 --- a/deprecated/segmentation_pipeline/tools.py +++ /dev/null @@ -1,171 +0,0 @@ -import random -import numpy -import os -import h5py -import warnings -import multiprocessing -import sys -import traceback -from concurrent.futures import ThreadPoolExecutor -from sklearn.cross_validation import KFold -import xgboost as xgb -from superpixels import * - -from colorama import init, Fore, Back, Style -init() - - - -def generateSplits(n, nSplits, frac=0.25): - res = [] - for x in range(nSplits): - indices = numpy.arange(n) - numpy.random.shuffle(indices) - - na = int(float(n)*frac) - a = indices[0:na] - b = indices[na:n] - - res.append((a, b)) - - return res - - - - - - - - - -class ThreadPoolExecutorStackTraced(ThreadPoolExecutor): - - def submit(self, fn, *args, **kwargs): - """Submits the wrapped function instead of `fn`""" - - return super(ThreadPoolExecutorStackTraced, self).submit( - self._function_wrapper, fn, *args, **kwargs) - - def _function_wrapper(self, fn, *args, **kwargs): - """Wraps `fn` in order to preserve the traceback of any kind of - raised exception - - """ - try: - return fn(*args, **kwargs) - except Exception: - raise sys.exc_info()[0](traceback.format_exc()) # Creates an - # exception of the - # same type with the - # traceback as - # message - - - -def printWarning(txt): - print(Fore.RED + Back.BLACK + txt) - print(Style.RESET_ALL) - -def ensureDir(f): - if not os.path.exists(f): - os.makedirs(f) - - -def hasH5File(f, key=None): - if os.path.exists(f): - return True - else: - return False - - -def hasFile(f): - if os.path.exists(f): - return True - else: - return False - -def isH5Path(h5path): - if isinstance(h5path, tuple): - if len(h5path) == 2: - return True - else: - return False - - -def threadExecutor(nThreads = multiprocessing.cpu_count()): - return ThreadPoolExecutorStackTraced(max_workers=nThreads) - - - -def h5Read(f, d='data'): - f5 = h5py.File(f,'r') - array = f5[d][:] - f5.close() - return array - - - - -class Classifier(object): - def __init__(self, nRounds=200, maxDepth=3, nThreads=multiprocessing.cpu_count(), silent=1): - - self.param = { - 'bst:max_depth':maxDepth, - 'bst:eta':1, - 'silent':silent, - 'num_class':2, - 'objective':'multi:softprob' - } - self.param['nthread'] = nThreads - - self.nRounds = nRounds - self.nThreads = nThreads - - def train(self, X, Y, getApproxError=False): - dtrain = xgb.DMatrix(X, label=Y) - self.bst = xgb.train(self.param, dtrain, self.nRounds) - - if getApproxError: - - e = 0.0 - c = 0.0 - - kf = KFold(Y.shape[0], n_folds=4) - for train_index, test_index in kf: - - XTrain = X[train_index, :] - XTest = X[test_index, :] - - YTrain = Y[train_index] - YTest = Y[test_index] - - dtrain2 = xgb.DMatrix(XTrain, label=YTrain) - bst = xgb.train(self.param, dtrain2, self.nRounds) - - - dtest = xgb.DMatrix(XTest) - probs = bst.predict(dtest) - ypred =numpy.argmax(probs, axis=1) - - - - error = float(numpy.sum(ypred != YTest)) - e += error - c += float(len(YTest)) - - e/=c - - return e - - - - def save(self, fname): - self.bst.save_model(fname) - - def load(self, fname): - self.bst = xgb.Booster({'nthread':self.nThreads}) #init model - self.bst.load_model(fname) # load data - - def predict(self, X): - dtest = xgb.DMatrix(X) - return self.bst.predict(dtest) \ No newline at end of file diff --git a/deprecated/supervoxel_pipeline.py b/deprecated/supervoxel_pipeline.py deleted file mode 100644 index ab44c6b36..000000000 --- a/deprecated/supervoxel_pipeline.py +++ /dev/null @@ -1,369 +0,0 @@ -from __future__ import print_function -from __future__ import division - -import os -import nifty -import nifty.viewer -import numpy -import vigra -import h5py - -import nifty.tools -from multiprocessing import cpu_count -import pylab -import scipy.ndimage -import math -import threading -import fastfilters - - - - -def oversegPipeline(pmapPath, outFolder, outName, pixelResNm=2.0 ,reduceBy=[5,10,20]): - - - scale = pixelResNm / 2.0 - - if not os.path.exists(outFolder): - os.makedirs(outFolder) - - - class DummyWithStatement: - def __enter__(self): - pass - def __exit__(self, type, value, traceback): - pass - - - - def makeBall(r): - size = 2*r + 1 - - mask = numpy.zeros([size]*3) - - for x0 in range(-1*r, r + 1): - for x1 in range(-1*r, r + 1): - for x2 in range(-1*r, r + 1): - - if math.sqrt(x0**2 + x1**2 + x2**2) <= r: - mask[x0+r, x1+r, x2+r] = 1 - - return mask, (r,r,r) - - - def membraneOverseg3D(pmapDset, heightMapDset, **kwargs): - - - - axisResolution = kwargs.get("axisResolution",['4nm']*3) - featureBlockShape = kwargs.get("featureBlockShape",['100']*3) - shape = pmapDset.shape[0:3] - - - roiBegin = kwargs.get("roiBegin", [0]*3) - roiEnd = kwargs.get("roiEnd", shape) - nWorkers = kwargs.get("nWorkers",cpu_count()) - invertPmap = kwargs.get("invertPmap",False) - - blocking = nifty.tools.blocking(roiBegin=roiBegin, roiEnd=roiEnd, blockShape=featureBlockShape) - margin = [45 ,45,45] - - - def pmapToHeightMap(pmap): - - r = int(min(3.0 * scale, 1.0) + 0.5) - footprint, origin = makeBall(r=r) - - - blurredSmall = fastfilters.gaussianSmoothing(pmap, 1.0*scale) - blurredLarge = fastfilters.gaussianSmoothing(pmap, 6.0*scale) - blurredSuperLarge = fastfilters.gaussianSmoothing(pmap, 10.0*scale) - - - combined = pmap + blurredSuperLarge*0.3 + 0.15*blurredLarge + 0.1*blurredSmall - - r = int(min(5.0 * scale, 1.0) + 0.5) - footprint, origin = makeBall(r=r) - - combined = scipy.ndimage.percentile_filter(input=combined, - #size=(20,20,20), - footprint=footprint, - #origin=origin, - mode='reflect', - percentile=50.0) - - - if False: - nifty.viewer.view3D(pmap, show=False, title='pm',cmap='gray') - nifty.viewer.view3D(medianImg, show=False, title='medianImg',cmap='gray') - nifty.viewer.view3D(combined, show=False, title='combined',cmap='gray') - pylab.show() - - return combined - - - numberOfBlocks = blocking.numberOfBlocks - lock = threading.Lock() - noLock = DummyWithStatement() - done = [0] - - - for blockIndex in range(numberOfBlocks): - blockWithHalo = blocking.getBlockWithHalo(blockIndex, margin) - block = blocking.getBlock(blockIndex) - outerBlock = blockWithHalo.outerBlock - innerBlock = blockWithHalo.innerBlock - innerBlockLocal = blockWithHalo.innerBlockLocal - - #print("B ",block.begin, block.end) - #print("O ",outerBlock.begin, outerBlock.end) - #print("I ",innerBlock.begin, innerBlock.end) - #print("IL",innerBlockLocal.begin, innerBlockLocal.end) - - - - with nifty.tools.progressBar(size=numberOfBlocks) as bar: - - def f(blockIndex): - blockWithHalo = blocking.getBlockWithHalo(blockIndex, margin, margin) - #print "fo" - outerBlock = blockWithHalo.outerBlock - outerSlicing = nifty.tools.getSlicing(outerBlock.begin, outerBlock.end) - b,e = outerBlock.begin, outerBlock.end - #print bi - # - - maybeLock = [lock, noLock][isinstance(pmapDset, numpy.ndarray)] - - with maybeLock: - if invertPmap: - outerPmap = 1.0 - pmapDset[b[0]:e[0], b[1]:e[1], b[2]:e[2]] - else: - outerPmap = pmapDset[b[0]:e[0], b[1]:e[1], b[2]:e[2]] - - heightMap = pmapToHeightMap(outerPmap) - - # - innerBlockLocal = blockWithHalo.innerBlockLocal - b,e = innerBlockLocal.begin, innerBlockLocal.end - innerHeightMap = heightMap[b[0]:e[0], b[1]:e[1], b[2]:e[2]] - - - b,e = blockWithHalo.innerBlock.begin, blockWithHalo.innerBlock.end - - if isinstance(heightMapDset,numpy.ndarray): - #print("NOT LOCKED") - heightMapDset[b[0]:e[0], b[1]:e[1], b[2]:e[2]] = innerHeightMap - with lock: - done[0] += 1 - bar.update(done[0]) - - else: - with lock: - #print("locked",b,e) - - heightMapDset[b[0]:e[0], b[1]:e[1], b[2]:e[2]] = innerHeightMap - done[0] += 1 - bar.update(done[0]) - - - nifty.tools.parallelForEach(range(blocking.numberOfBlocks), f=f, nWorkers=nWorkers) - - - def makeSmallerSegNifty(oseg, volume_feat, reduceBySetttings, wardnessSettings, baseFilename): - import nifty - import nifty.graph - import nifty.graph.rag - import nifty.graph.agglo - - nrag = nifty.graph.rag - nagglo = nifty.graph.agglo - - print("overseg in c order starting at zero") - oseg = numpy.require(oseg, dtype='uint32',requirements='C') - oseg -= 1 - - print("make rag") - rag = nifty.graph.rag.gridRag(oseg) - - print("volfeatshape") - vFeat = numpy.require(volume_feat, dtype='float32',requirements='C') - - print("accumulate means and counts") - eFeatures, nFeatures = nrag.accumulateMeanAndLength(rag, vFeat, [100,100,100],-1) - - eMeans = eFeatures[:,0] - eSizes = eFeatures[:,1] - nSizes = nFeatures[:,1] - - print("get clusterPolicy") - - - for wardness in wardnessSettings: - for reduceBy in reduceBySetttings: - - print("wardness",wardness,"reduceBy",reduceBy) - - numberOfNodesStop = int(float(rag.numberOfNodes)/float(reduceBy) + 0.5) - numberOfNodesStop = max(1,numberOfNodesStop) - numberOfNodesStop = min(rag.numberOfNodes, numberOfNodesStop) - - - clusterPolicy = nagglo.edgeWeightedClusterPolicy( - graph=rag, edgeIndicators=eMeans, - edgeSizes=eSizes, nodeSizes=nSizes, - numberOfNodesStop=numberOfNodesStop, - sizeRegularizer=float(wardness)) - - print("do clustering") - agglomerativeClustering = nagglo.agglomerativeClustering(clusterPolicy) - agglomerativeClustering.run() - seg = agglomerativeClustering.result()#out=[1,2,3,4]) - - print("make seg dense") - dseg = nifty.tools.makeDense(seg) - - print(dseg.dtype, type(dseg)) - - print("project to pixels") - pixelData = nrag.projectScalarNodeDataToPixels(rag, dseg.astype('uint32')) - print("done") - - - outFilename = baseFilename + str(wardness) + "_" + str(reduceBy) + ".h5" - - agglosegH5 = h5py.File(outFilename,'w') - agglosegDset = agglosegH5.create_dataset('data',shape=oseg.shape[0:3], chunks=(100,100,100),dtype='uint32',compression="gzip") - agglosegDset[:,:,:] = pixelData - agglosegH5.close() - - - - if True: - - pmapH5 = h5py.File(pmapPath,'r') - pmapDset = pmapH5['data'] - shape = list(pmapDset.shape[0:3]) - - subset = None - sshape = (subset,)*3 - heightMapFile = os.path.join(outFolder,"%s_heightMap.h5"%outName) - heightMapH5 = h5py.File(heightMapFile,'w') - heightMapDset = heightMapH5.create_dataset('data',shape=shape,dtype='float32',chunks=(64,)*3) - - print("bra",heightMapFile) - #sys.exit(1) - - print("load pmap in ram (since we have enough") - if subset is None: - pmapArray= pmapDset[:,:,:] - else: - pmapArray= pmapDset[0:subset,0:subset,0:subset] - pmapH5.close() - - - if shape[0] >= 2000: - featureBlockShape = [350,350,350] - elif shape[0] < 2000 and shape[0] >= 1000: - featureBlockShape = [250,250,250] - elif shape[0] < 1000 and shape[0] >= 500: - featureBlockShape = [100,100,100] - - with vigra.Timer("st"): - params = { - "axisResolution" : [2.0, 2.0, 2.0], - "featureBlockShape" : featureBlockShape, - "invertPmap": False, - #"roiBegin": [0,0,0], - #"roiEnd": sshape - #"nWorkers":1, - } - membraneOverseg3D(pmapArray,heightMapDset, **params) - - - heightMapH5.close() - - - if True: - - with vigra.Timer("read hmap"): - heightMapFile = os.path.join(outFolder,"%s_heightMap.h5"%outName) - heightMapH5 = h5py.File(heightMapFile,'r') - heightMapDset = heightMapH5['data'] - heightMap = heightMapDset[:,:,:] - shape = list(heightMap.shape) - heightMapH5.close() - - with vigra.Timer("do overseg"): - overseg, nseg = vigra.analysis.unionFindWatershed3D(heightMap.T, blockShape=(100,100,100)) - overseg = overseg.T - - - - - with vigra.Timer("write res"): - oversegFile = os.path.join(outFolder,"%s_ufd_overseg.h5"%outName) - oversegH5 = h5py.File(oversegFile,'w') - oversegDset = oversegH5.create_dataset('data',shape=shape, data=overseg, chunks=(64,)*3,compression='gzip') - oversegDset.attrs['nseg'] = nseg - oversegH5.close() - - #if False: - # import nifty.hdf5 - # with vigra.Timer("ws"): - # nifty.hdf5.unionFindWatershed(heightMapFile,'data', oversegFile,'data',[100,100,100]) - - - if True: - - heightMapFile = os.path.join(outFolder,"%s_heightMap.h5"%outName) - heightMapH5 = h5py.File(heightMapFile,'r') - heightMapDset = heightMapH5['data'] - - oversegFile = os.path.join(outFolder,"%s_ufd_overseg.h5"%outName) - oversegH5 = h5py.File(oversegFile,'r') - oversegDset = oversegH5['data'] - - - - - print("read hmap") - heightMap = heightMapDset[:,:,:] - - print("read oseg") - overseg = oversegDset[:,:,:] - - heightMapH5.close() - oversegH5.close() - - print("make smaller") - agglosegBaseFile = os.path.join(outFolder,"%s_agglo_seg"%outName) - makeSmallerSegNifty(overseg,heightMap, reduceBy, [0.3], agglosegBaseFile) - - - - - -# 2nm -if False: - assert False #deprecated - pmapPath = "/home/tbeier/prediction_binary_full.h5" - heightMapFile = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxels/2nm/heightMap.h5" - oversegFile = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxels/2nm/ufd_overseg6.h5" - agglosegBaseFile = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxels/2nm/aggloseg_" - -# 4nm -elif False: - pmapPath = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/pc_out/2nm/prediction_binary_4nm.h5" - outFolder = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxel4nm" - pixelResNm = 4 - outName = "_4nm" -elif True: - pmapPath = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/pc_out/2nm/prediction_binary_8nm.h5" - pmapPath = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/pnew/p8nm_0.h5" - outFolder = "/media/tbeier/4cf81285-be72-45f5-8c63-fb8e9ff4476c/supervoxel8nm" - pixelResNm = 8 - outName = "_8nm" - - -oversegPipeline(pmapPath=pmapPath, outFolder=outFolder, outName=outName, pixelResNm=pixelResNm, reduceBy=[20,70,80,90]) diff --git a/deprecated/test_cgp2d.py b/deprecated/test_cgp2d.py deleted file mode 100644 index aa8baee77..000000000 --- a/deprecated/test_cgp2d.py +++ /dev/null @@ -1,61 +0,0 @@ -from __future__ import print_function,division - -import unittest -import nifty -import nifty.cgp as ncgp - -import numpy - -seg = [ - [1,1,2], - [1,3,1], - [1,1,1], - [4,4,5], - [6,6,7], -] -seg = numpy.array(seg,dtype='uint32') - -tgrid = ncgp.TopologicalGrid2D(seg) -ftgrid = ncgp.FilledTopologicalGrid2D(tgrid) -numberOfCells = tgrid.numberOfCells - - - - -tGridArray = numpy.array(tgrid) -print(tGridArray) - -tGridArray = numpy.array(ftgrid) -print(tGridArray) - -cellBounds = ncgp.Bounds2D(tgrid) - -cell0Bounds = numpy.array(cellBounds[0]) -cell1Bounds = numpy.array(cellBounds[1]) - - -cellGeometry = ncgp.Geometry2D(tgrid) - -for cellType in (0,1,2): - - nCells = numberOfCells[cellType] - - print("cell-%d (#%d):"%(cellType, nCells)) - cellsGeo = cellGeometry[cellType] - for cellIndex in range(nCells): - cellGeo = numpy.array(cellsGeo[cellIndex]) - print(cellIndex,":\n",cellGeo) - - - -if False: - for celLType in (0,1): - nCells = tgrid.numberOfCells(celLType) - print("%d-Cell:"%celLType) - cellBounds = bounds[celLType] - for c in range(nCells): - print(" ",[i for i in cellBounds[c]]) - - - - diff --git a/externals/marray b/externals/marray deleted file mode 160000 index 8df23c17d..000000000 --- a/externals/marray +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 8df23c17d916ef9a0a4e9ed74eb551d882966b13 diff --git a/include/nifty/marray/andres/marray-fftw.hxx b/include/nifty/marray/andres/marray-fftw.hxx new file mode 100644 index 000000000..650d87863 --- /dev/null +++ b/include/nifty/marray/andres/marray-fftw.hxx @@ -0,0 +1,149 @@ +#pragma once +#ifndef ANDRES_MARRAY_FFTW_HXX +#define ANDRES_MARRAY_FFTW_HXX + +#include +#include +#include +#include + +#include "marray.hxx" + +#include + +namespace andres { + +template +class FFT { +public: + typedef S size_type; + + FFT(const andres::Marray&, andres::Marray >&); + ~FFT(); + void execute(); + +private: + fftw_plan plan_; +}; + +template +class IFFT { +public: + typedef S size_type; + + IFFT(const andres::Marray >& in, andres::Marray& out); + ~IFFT(); + void execute(); + +private: + fftw_plan plan_; +}; + +template +inline +FFT::FFT( + const andres::Marray& in, + andres::Marray >& out +) { + if(in.coordinateOrder() != andres::CoordinateOrder::FirstMajorOrder) { + throw std::runtime_error("Marray not in first-major order"); + } + if(out.coordinateOrder() != andres::CoordinateOrder::FirstMajorOrder) { + throw std::runtime_error("Marray not in first-major order"); + } + if(in.dimension() != 2) { + throw std::runtime_error("Marray not 2-dimensional."); + } + + const int nRows = in.shape(0); + const int nCols = in.shape(1); + { + int shape[] = {nRows, nCols / 2 + 1}; + out.resize(shape, shape + 2); + } + const int rank = 1; // 1D transform + int n[] = {nCols}; // 1D transforms of length nCols + int howmany = nRows; // nRows 1D transforms + int idist = nCols; // number of columns in input array + int odist = nCols/2+1; // number of columns in output array + int istride = 1; int ostride = 1; // array is contiguous in memory + int* inembed = n; int *onembed = n; +# pragma omp critical + { + plan_ = fftw_plan_many_dft_r2c( + rank, n, howmany, + &in(0), inembed, istride, idist, + reinterpret_cast(&out(0)), onembed, ostride, odist, + FFTW_ESTIMATE + ); + } +} + +template +inline +FFT::~FFT() { +# pragma omp critical + fftw_destroy_plan(plan_); +} + +template +inline void +FFT::execute() { + fftw_execute(plan_); +} + +template +IFFT::IFFT( + const andres::Marray >& in, + andres::Marray& out +) { + if(in.coordinateOrder() != andres::CoordinateOrder::FirstMajorOrder) { + throw std::runtime_error("Marray not in first-major order"); + } + if(out.coordinateOrder() != andres::CoordinateOrder::FirstMajorOrder) { + throw std::runtime_error("Marray not in first-major order"); + } + if(in.dimension() != 2) { + throw std::runtime_error("Marray not 2-dimensional."); + } + + const int nRows = in.shape(0); + const int nCols = in.shape(1); + { + int shape[] = {nRows, nCols}; + out.resize(shape, shape + 2); + } + const int rank = 1; // 1D transform + int n[] = {nCols}; // 1D transforms of length nCols + int howmany = nRows; // nRows 1D transforms + int idist = nCols; // number of columns in input array + int odist = static_cast(nCols); ///2+1; // number of columns in output array + int istride = 1; int ostride = 1; // array is contiguous in memory + int* inembed = n; int *onembed = n; +# pragma omp critical + { + plan_ = fftw_plan_many_dft_c2r( + rank, n, howmany, + reinterpret_cast(&in(0)), inembed, istride, idist, + &out(0), onembed, ostride, odist, + FFTW_ESTIMATE + ); + } +} + +template +inline +IFFT::~IFFT() { +# pragma omp critical + fftw_destroy_plan(plan_); +} + +template +inline void +IFFT::execute() { + fftw_execute(plan_); +} + +} // namespace andres + +#endif // #ifndef ANDRES_MARRAY_FFTW_HXX diff --git a/include/nifty/marray/andres/marray-hdf5.hxx b/include/nifty/marray/andres/marray-hdf5.hxx new file mode 100644 index 000000000..a2cf0deda --- /dev/null +++ b/include/nifty/marray/andres/marray-hdf5.hxx @@ -0,0 +1,1007 @@ +// Copyright (c) 2011-2013 by Bjoern Andres. +// +// This software was developed by Bjoern Andres. +// Enquiries shall be directed to bjoern@andres.sc. +// +// All advertising materials mentioning features or use of this software must +// display the following acknowledgement: ``This product includes +// andres::marray developed by Bjoern Andres. Please direct enquiries +// concerning andres::marray to bjoern@andres.sc''. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// - Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// - Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// - All advertising materials mentioning features or use of this software must +// display the following acknowledgement: ``This product includes +// andres::marray developed by Bjoern Andres. Please direct enquiries +// concerning andres::marray to bjoern@andres.sc''. +// - The name of the author must not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED +// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO +// EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +#pragma once +#ifndef MARRAY_HDF5_HXX +#define MARRAY_HDF5_HXX + +// compat fix for buggy hdf5 1.8 versions +#include +#if (H5_VERS_MAJOR == 1 && H5_VERS_MINOR >= 8 && defined(H5_USE_16_API_DEFAULT)) +#define H5Gcreate_vers 2 +#define H5Gopen_vers 2 +#define H5Dopen_vers 2 +#define H5Dcreate_vers 2 +#define H5Acreate_vers 2 +#endif + +#include +#include + +#include "marray.hxx" +#include "hdf5.h" + +namespace andres { +/// HDF5 import/export support. +namespace hdf5 { + +// assertion testing + +// \cond suppress doxygen +template class HandleCheck; +template<> class HandleCheck { +public: + HandleCheck() + { counter_ = H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL); } + void check() + { marray_detail::Assert( counter_ == H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL)); } +private: + ssize_t counter_; +}; +template<> class HandleCheck { +public: + void check() {} +}; +// \endcond + +// namespace variables + +const char reverseShapeAttributeName[14] = "reverse-shape"; + +// prototypes + +enum FileAccessMode {READ_ONLY, READ_WRITE}; +enum HDF5Version {DEFAULT_HDF5_VERSION, LATEST_HDF5_VERSION}; + +hid_t createFile(const std::string&, HDF5Version = DEFAULT_HDF5_VERSION); +hid_t openFile(const std::string&, FileAccessMode = READ_ONLY, HDF5Version = DEFAULT_HDF5_VERSION); +void closeFile(const hid_t&); + +hid_t createGroup(const hid_t&, const std::string&); +hid_t openGroup(const hid_t&, const std::string&); +void closeGroup(const hid_t&); + +template + void save(const hid_t&, const std::string&, const Marray&); +template + void save(const hid_t&, const std::string&, const View&); +template + void save(const hid_t&, const std::string&, const std::vector&); +template + void saveHyperslab(const hid_t&, const std::string&, + BaseIterator, BaseIterator, ShapeIterator, const Marray&); +template + void create(const hid_t&, const std::string&, ShapeIterator, + ShapeIterator, CoordinateOrder); + +template + void load(const hid_t&, const std::string&, Marray&); +// TODO: implement the following functions that throw an error +// if the dimension is not correct +// template +// void load(const hid_t&, const std::string&, Vector&); +// template +// void load(const hid_t&, const std::string&, Matrix&); +template + void load(const hid_t&, const std::string&, std::vector&); + +template + void load(const std::string&, const std::string&, Marray&, HDF5Version = DEFAULT_HDF5_VERSION); +// TODO: implement the following functions that throw an error +// if the dimension is not correct +// template +// void load(const std::string&, const std::string&, Vector&, HDF5Version = DEFAULT_HDF5_VERSION); +// template +// void load(const std::string&, const std::string&, Matrix&, HDF5Version = DEFAULT_HDF5_VERSION); +template + void load(const std::string&, const std::string&, std::vector&, HDF5Version = DEFAULT_HDF5_VERSION); + +template + void loadShape(const hid_t&, const std::string&, std::vector&); +template + void loadHyperslab(const hid_t&, const std::string&, + BaseIterator, BaseIterator, ShapeIterator, Marray&); + +// type conversion from C++ to HDF5 + +// \cond suppress doxygen +template +inline hid_t uintTypeHelper() { + switch(sizeof(T)) { + case 1: + return H5T_STD_U8LE; + case 2: + return H5T_STD_U16LE; + case 4: + return H5T_STD_U32LE; + case 8: + return H5T_STD_U64LE; + default: + throw std::runtime_error("No matching HDF5 type."); + } +} + +template +inline hid_t intTypeHelper() { + switch(sizeof(T)) { + case 1: + return H5T_STD_I8LE; + case 2: + return H5T_STD_I16LE; + case 4: + return H5T_STD_I32LE; + case 8: + return H5T_STD_I64LE; + default: + throw std::runtime_error("No matching HDF5 type."); + } +} + +template +inline hid_t floatingTypeHelper() { + switch(sizeof(T)) { + case 4: + return H5T_IEEE_F32LE; + case 8: + return H5T_IEEE_F64LE; + default: + throw std::runtime_error("No matching HDF5 type."); + } +} + +template +inline hid_t hdf5Type(); + +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } + +template<> inline hid_t hdf5Type() + { return intTypeHelper(); } +template<> inline hid_t hdf5Type() + { return uintTypeHelper(); } +template<> inline hid_t hdf5Type() + { return intTypeHelper(); } +template<> inline hid_t hdf5Type() + { return intTypeHelper(); } +template<> inline hid_t hdf5Type() + { return intTypeHelper(); } +template<> inline hid_t hdf5Type() + { return intTypeHelper(); } + +template<> inline hid_t hdf5Type() + { return floatingTypeHelper(); } +template<> inline hid_t hdf5Type() + { return floatingTypeHelper(); } +// \endcond + +// implementation + +/// Create and close an HDF5 dataset to store Marray data. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param begin Iterator to the beginning of a sequence that determines the shape of the dataset. +/// \param end Iterator to the end of a sequence that determines the shape of the dataset. +/// \param coordinateOrder Coordinate order of the Marray. +/// +/// \sa save(), saveHyperslab() +/// +template +void create( + const hid_t& groupHandle, + const std::string& datasetName, + ShapeIterator begin, + ShapeIterator end, + CoordinateOrder coordinateOrder +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + // build dataspace + hid_t datatype = H5Tcopy(hdf5Type()); + std::size_t dimension = std::distance(begin, end); + std::vector shape(dimension); + if(coordinateOrder == FirstMajorOrder) { + // copy shape as is + for(std::size_t j=0; j +void save( + const hid_t& groupHandle, + const std::string& datasetName, + const Marray& in +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + // build dataspace + hid_t datatype = H5Tcopy(hdf5Type()); + std::vector shape(in.dimension()); + if(in.coordinateOrder() == FirstMajorOrder) { + // copy shape as is + for(std::size_t j=0; j +inline void save( + const hid_t& groupHandle, + const std::string& datasetName, + const View& in +) { + Marray m = in; + save(groupHandle, datasetName, m); +} + +/// Save an std::vector as an HDF5 dataset. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param in std::vector. +/// +/// \sa saveHyperslab() +/// +template +void save( + const hid_t& groupHandle, + const std::string& datasetName, + const std::vector& in +) +{ + std::size_t shape[] = {in.size()}; + andres::Marray mvector(shape, shape + 1); + for(size_t j=0; j +inline void +load( + const std::string& filename, + const std::string& datasetName, + Marray& out, + HDF5Version hdf5version +) { + hid_t file = openFile(filename, READ_ONLY, hdf5version); + load(file, datasetName, out); + closeFile(file); +} + +/// Open an HDF5 file (read only), load an std::vector from an HDF5 dataset, and close the file. +/// +/// \param filename Name of the file. +/// \param datasetName Name of the HDF5 dataset. +/// \param out std::vector. +/// \param hdf5version HDF5 version tag. +/// +/// \sa load(), loadHyperslab() +/// +/// TODO: write a unit test for this function +template +inline void +load( + const std::string& filename, + const std::string& datasetName, + std::vector& out, + HDF5Version hdf5version +) { + hid_t file = openFile(filename, READ_ONLY, hdf5version); + load(file, datasetName, out); + closeFile(file); +} + +/// Load an std::vector from an HDF5 dataset. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param out std::vector. +/// +/// \sa load(), loadHyperslab() +/// +/// TODO: write a unit test for this function +template +inline void +load( + const hid_t& groupHandle, + const std::string& datasetName, + std::vector& out +) { + Marray in; + load(groupHandle, datasetName, in); + if(in.dimension() != 1) { + throw std::runtime_error("dataset is not 1-dimensional."); + } + out.assign(in.begin(), in.end()); +} + +/// Load an Marray from an HDF5 dataset. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param out Marray. +/// +/// \sa loadHyperslab() +/// +template +void load( + const hid_t& groupHandle, + const std::string& datasetName, + Marray& out +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + hid_t dataset = H5Dopen(groupHandle, datasetName.c_str(), H5P_DEFAULT); + if(dataset < 0) { + throw std::runtime_error("Marray cannot open dataset."); + } + hid_t filespace = H5Dget_space(dataset); + hid_t type = H5Dget_type(dataset); + hid_t nativeType = H5Tget_native_type(type, H5T_DIR_DESCEND); + if(!H5Tequal(nativeType, hdf5Type())) { + H5Dclose(dataset); + H5Tclose(nativeType); + H5Tclose(type); + H5Sclose(filespace); + throw std::runtime_error("Data types not equal error."); + } + int dimension = H5Sget_simple_extent_ndims(filespace); + std::vector shape(dimension); + herr_t status = H5Sget_simple_extent_dims(filespace, &shape[0], NULL); + if(status < 0) { + H5Dclose(dataset); + H5Tclose(nativeType); + H5Tclose(type); + H5Sclose(filespace); + throw std::runtime_error("H5Sget_simple_extent_dims error."); + } + hid_t memspace = H5Screate_simple(dimension, &shape[0], NULL); + + // resize marray + std::vector newShape((std::size_t)(dimension)); + for(std::size_t j=0; j 0) { + // reverse shape + out = Marray(SkipInitialization, newShape.rbegin(), + newShape.rend(), LastMajorOrder); + } + else { + // don't reverse shape + out = Marray(SkipInitialization, newShape.begin(), + newShape.end(), FirstMajorOrder); + } + + // read + status = H5Dread(dataset, nativeType, memspace, filespace, + H5P_DEFAULT, &out(0)); + H5Dclose(dataset); + H5Tclose(nativeType); + H5Tclose(type); + H5Sclose(memspace); + H5Sclose(filespace); + if(status < 0) { + throw std::runtime_error("Marray cannot read from dataset."); + } + + handleCheck.check(); +} + +/// Load the shape of an HDF5 dataset. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param out Shape. +/// +/// \sa load() +/// +template +void loadShape( + const hid_t& groupHandle, + const std::string& datasetName, + std::vector& out +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + // load shape from HDF5 file + hid_t dataset = H5Dopen(groupHandle, datasetName.c_str(), H5P_DEFAULT); + if(dataset < 0) { + throw std::runtime_error("Marray cannot open dataset."); + } + hid_t filespace = H5Dget_space(dataset); + hsize_t dimension = H5Sget_simple_extent_ndims(filespace); + hsize_t* shape = new hsize_t[(std::size_t)(dimension)]; + herr_t status = H5Sget_simple_extent_dims(filespace, shape, NULL); + if(status < 0) { + H5Dclose(dataset); + H5Sclose(filespace); + delete[] shape; + throw std::runtime_error("Marray cannot get extension of dataset."); + } + + // write shape to out + out = std::vector(dimension); + if(H5Aexists(dataset, reverseShapeAttributeName) > 0) { + for(std::size_t j=0; j +void loadHyperslab( + const hid_t& groupHandle, + const std::string& datasetName, + BaseIterator baseBegin, + BaseIterator baseEnd, + ShapeIterator shapeBegin, + Marray& out +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + // open dataset + hid_t dataset = H5Dopen(groupHandle, datasetName.c_str(), H5P_DEFAULT); + if(dataset < 0) { + throw std::runtime_error("Marray cannot open dataset."); + } + + // determine shape of hyperslab and array + std::size_t size = std::distance(baseBegin, baseEnd); + std::vector offset(size); + std::vector slabShape(size); + std::vector marrayShape(size); + CoordinateOrder coordinateOrder; + if(H5Aexists(dataset, reverseShapeAttributeName) > 0) { + // reverse base and shape + coordinateOrder = LastMajorOrder; + std::size_t j = size-1; + std::size_t k = 0; + for(;;) { + offset[j] = hsize_t(*baseBegin); + slabShape[j] = hsize_t(*shapeBegin); + marrayShape[k] = slabShape[j]; + if(j == 0) { + break; + } + else { + ++baseBegin; + ++shapeBegin; + ++k; + --j; + } + } + } + else { + // don't reverse base and shape + coordinateOrder = FirstMajorOrder; + for(std::size_t j=0; j())) { + throw std::runtime_error("data type of stored hdf5 dataset and passed array do not match in loadHyperslab"); + } + + hid_t dataspace = H5Dget_space(dataset); + herr_t status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, + &offset[0], NULL, &slabShape[0], NULL); + if(status < 0) { + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); + throw std::runtime_error("Marray cannot select hyperslab. Check offset and shape!"); + } + + // select memspace hyperslab + hid_t memspace = H5Screate_simple(int(size), &marrayShape[0], NULL); + std::vector offsetOut(size, 0); // no offset + status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, &offsetOut[0], + NULL, &marrayShape[0], NULL); + if(status < 0) { + H5Sclose(memspace); + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); + throw std::runtime_error("Marray cannot select hyperslab. Check offset and shape!"); + } + + // read from dataspace into memspace + out = Marray(SkipInitialization, &marrayShape[0], + (&marrayShape[0])+size, coordinateOrder); + status = H5Dread(dataset, datatype, memspace, dataspace, + H5P_DEFAULT, &(out(0))); + + // clean up + H5Sclose(memspace); + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); + if(status < 0) { + throw std::runtime_error("Marray cannot read from dataset."); + } + handleCheck.check(); +} + +/// Save an Marray as a hyperslab into an HDF5 dataset. +/// +/// \param groupHandle Handle of the parent HDF5 file or group. +/// \param datasetName Name of the HDF5 dataset. +/// \param baseBegin Iterator to the beginning of the sequence that determines the first coordinate of the hyperslab. +/// \param baseEnd Iterator to the end of the sequence that determines the first coordinate of the hyperslab. +/// \param shapeBegin Iterator to the beginning of the sequence that determines the shape of the hyperslab. +/// \param in Marray. +/// +/// \sa loadHyperslab(), create() +/// +template +void +saveHyperslab( + const hid_t& groupHandle, + const std::string& datasetName, + BaseIterator baseBegin, + BaseIterator baseEnd, + ShapeIterator shapeBegin, + const Marray& in +) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || groupHandle >= 0); + HandleCheck handleCheck; + + // open dataset + hid_t dataset = H5Dopen(groupHandle, datasetName.c_str(), H5P_DEFAULT); + if(dataset < 0) { + throw std::runtime_error("Marray cannot open dataset."); + } + + // determine hyperslab shape + std::vector memoryShape(in.dimension()); + for(std::size_t j=0; j offset(size); + std::vector slabShape(size); + bool reverseShapeAttribute = + (H5Aexists(dataset, reverseShapeAttributeName) > 0); + if(reverseShapeAttribute && in.coordinateOrder() == LastMajorOrder) { + // reverse base and shape + std::size_t j = size-1; + for(;;) { + offset[j] = hsize_t(*baseBegin); + slabShape[j] = hsize_t(*shapeBegin); + if(j == 0) { + break; + } + else { + ++baseBegin; + ++shapeBegin; + --j; + } + } + } + else if(!reverseShapeAttribute && in.coordinateOrder() == FirstMajorOrder) { + for(std::size_t j=0; j memoryOffset(int(in.dimension()), 0); // no offset + status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, &memoryOffset[0], NULL, + &memoryShape[0], NULL); + if(status < 0) { + H5Sclose(memspace); + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); + throw std::runtime_error("Marray cannot select hyperslab. Check offset and shape!"); + } + + // write from memspace to dataspace + status = H5Dwrite(dataset, datatype, memspace, dataspace, H5P_DEFAULT, &(in(0))); + + // clean up + H5Sclose(memspace); + H5Tclose(datatype); + H5Sclose(dataspace); + H5Dclose(dataset); + if(status < 0) { + throw std::runtime_error("Marray cannot write to dataset."); + } + handleCheck.check(); +} + +/// Create an HDF5 file. +/// +/// \param filename Name of the file. +/// \param hdf5version HDF5 version tag. +/// +/// \returns HDF5 handle +/// +/// \sa openFile(), closeFile() +/// +inline hid_t +createFile +( + const std::string& filename, + HDF5Version hdf5version +) +{ + hid_t version = H5P_DEFAULT; + if(hdf5version == LATEST_HDF5_VERSION) { + version = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_libver_bounds(version, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST); + } + + hid_t fileHandle = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, version); + if(fileHandle < 0) { + throw std::runtime_error("Could not create HDF5 file: " + filename); + } + + return fileHandle; +} + +/// Open an HDF5 file. +/// +/// \param filename Name of the file. +/// \param fileAccessMode File access mode. +/// \param hdf5version HDF5 version tag. +/// +/// \returns HDF5 handle +/// +/// \sa closeFile(), createFile() +/// +inline hid_t +openFile +( + const std::string& filename, + FileAccessMode fileAccessMode, + HDF5Version hdf5version +) +{ + hid_t access = H5F_ACC_RDONLY; + if(fileAccessMode == READ_WRITE) { + access = H5F_ACC_RDWR; + } + + hid_t version = H5P_DEFAULT; + if(hdf5version == LATEST_HDF5_VERSION) { + version = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_libver_bounds(version, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST); + } + + hid_t fileHandle = H5Fopen(filename.c_str(), access, version); + if(fileHandle < 0) { + throw std::runtime_error("Could not open HDF5 file: " + filename); + } + + return fileHandle; +} + +/// Close an HDF5 file +/// +/// \param handle Handle to the HDF5 file. +/// +/// \sa openFile(), createFile() +/// +inline void closeFile +( + const hid_t& handle +) +{ + H5Fclose(handle); +} + +/// Create an HDF5 group. +/// +/// \param parentHandle HDF5 handle on the parent group or file. +/// \param groupName Name of the group. +/// \returns HDF5 handle on the created group +/// +/// \sa openGroup(), closeGroup() +/// +inline hid_t +createGroup +( + const hid_t& parentHandle, + const std::string& groupName +) +{ + hid_t groupHandle = H5Gcreate(parentHandle, groupName.c_str(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + if(groupHandle < 0) { + throw std::runtime_error("Could not create HDF5 group."); + } + return groupHandle; +} + +/// Open an HDF5 group. +/// +/// \param parentHandle HDF5 handle on the parent group or file. +/// \param groupName Name of the group. +/// \returns HDF5 handle on the opened group. +/// +/// \sa createGroup(), closeGroup() +/// +inline hid_t +openGroup +( + const hid_t& parentHandle, + const std::string& groupName +) +{ + hid_t groupHandle = H5Gopen(parentHandle, groupName.c_str(), H5P_DEFAULT); + if(groupHandle < 0) { + throw std::runtime_error("Could not open HDF5 group."); + } + return groupHandle; +} + +/// Close an HDF5 group. +/// +/// \param handle HDF5 handle on group to close. +/// +/// \sa openGroup(), createGroup() +/// +inline void +closeGroup( + const hid_t& handle +) { + H5Gclose(handle); +} + +} // namespace hdf5 +} // namespace andres + +#endif // #ifndef MARRAY_HDF5_HXX + diff --git a/include/nifty/marray/andres/marray.hxx b/include/nifty/marray/andres/marray.hxx new file mode 100644 index 000000000..f291175c3 --- /dev/null +++ b/include/nifty/marray/andres/marray.hxx @@ -0,0 +1,6046 @@ +/// \mainpage +/// Marray: Fast Runtime-Flexible Multi-dimensional Arrays and Views in C++. +/// \newline +/// +/// Copyright (c) 2013 by Bjoern Andres, bjoern@andres.sc +/// +/// \section section_abstract Short Description +/// Marray is a single header file for fast multi-dimensional arrays and views +/// in C++. Unlike in other implementations such as boost MultiArray and +/// Blitz++, the dimension of Marray views and arrays can be set and changed at +/// runtime. Dimension is not a template parameter in Marray. Arrays and views +/// that have the same type of entries but different dimension are therefore of +/// the same C++ type. In conjunction with the comprehensive and +/// convenient Marray interface, this brings some of the flexibility known from +/// high-level languages such as Python, R and MATLAB to C++. +/// +/// \section section_features Features +/// - Multi-dimensional arrays and views whose dimension, shape, size and +/// indexing order (first or last coordinate major order) can be set and +/// changed at runtime. +/// - Access to entries via coordinates, scalar indices, STL-compliant random +/// access iterators and C++11 initializer lists. +/// - Arithmetic operators with expression templates and automatic type +/// promotion. +/// - Support for STL-compliant allocators. +/// +/// \section section_tutorial Tutorial +/// - An introductory tutorial can be found at src/tutorial/tutorial.cxx +/// +/// \section section_cpp0x C++11 Extensions +/// - C++11 extensions are enabled by defining +/// - HAVE_CPP11_VARIADIC_TEMPLATES +/// - HAVE_CPP11_INITIALIZER_LISTS +/// - HAVE_CPP11_TEMPLATE_ALIASES +/// - HAVE_CPP11_STD_ARRAY +/// . +/// +/// \section section_license License +/// Copyright (c) 2013 by Bjoern Andres. +/// +/// This software was developed by Bjoern Andres. +/// Enquiries shall be directed to bjoern@andres.sc. +/// +/// Redistribution and use in source and binary forms, with or without +/// modification, are permitted provided that the following conditions are met: +/// - Redistributions of source code must retain the above copyright notice, +/// this list of conditions and the following disclaimer. +/// - Redistributions in binary form must reproduce the above copyright notice, +/// this list of conditions and the following disclaimer in the documentation +/// and/or other materials provided with the distribution. +/// - The name of the author must not be used to endorse or promote products +/// derived from this software without specific prior written permission. +/// +/// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED +/// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +/// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO +/// EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +/// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +/// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; +/// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +/// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR +/// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +/// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +/// +#pragma once +#ifndef ANDRES_MARRAY_HXX +#define ANDRES_MARRAY_HXX + +#include +#include +#include // runtime_error +#include +#include +#include +#include // memcpy +#include // reverse_iterator, distance +#include +#include +#include // cout +#include // allocator +#include // accumulate +#include // std::multiplies +#ifdef HAVE_CPP11_INITIALIZER_LISTS + #include +#endif +#ifdef HAVE_CPP11_STD_ARRAY + #include +#endif +/// The public API. +namespace andres { + +enum StringStyle {TableStyle, MatrixStyle}; ///< Flag to be used with the member function asString() of View. +enum CoordinateOrder {FirstMajorOrder, LastMajorOrder}; ///< Flag setting the order of coordinate tuples. +struct InitializationSkipping { }; ///< Flag to indicate initialization skipping. + +static const bool Const = true; ///< Flag to be used with the template parameter isConst of View and Iterator. +static const bool Mutable = false; ///< Flag to be used with the template parameter isConst of View and Iterator. +static const CoordinateOrder defaultOrder = FirstMajorOrder; ///< Default order of coordinate tuples. +static const InitializationSkipping SkipInitialization = InitializationSkipping(); ///< Flag to indicate initialization skipping. + +template + class ViewExpression; +// \cond suppress_doxygen +template + class UnaryViewExpression; +template + class BinaryViewExpression; +template + class BinaryViewExpressionScalarFirst; +template + class BinaryViewExpressionScalarSecond; +// \endcond suppress_doxygen +template > + class View; +#ifdef HAVE_CPP11_TEMPLATE_ALIASES + template using ConstView = View; +#endif +template > + class Iterator; +template > class Marray; + +// assertion testing +// #ifdef NDEBUG + static const bool MARRAY_NO_DEBUG = true; ///< General assertion testing disabled. + static const bool MARRAY_NO_ARG_TEST = true; ///< Argument testing disabled. +// #else +// static const bool MARRAY_NO_DEBUG = false; ///< General assertion testing enabled. +// static const bool MARRAY_NO_ARG_TEST = false; ///< Argument testing enabled. +//#endif + +// \cond suppress_doxygen +namespace marray_detail { + // meta-programming + template + struct IfBool; + template + struct IfBool + { typedef TRUECASE type; }; + template + struct IfBool + { typedef FALSECASE type; }; + + template + struct IsEqual + { static const bool type = false; }; + template + struct IsEqual + { static const bool type = true; }; + + template struct TypeTraits + { static const unsigned char position = 255; }; + template<> struct TypeTraits + { static const unsigned char position = 0; }; + template<> struct TypeTraits + { static const unsigned char position = 1; }; + template<> struct TypeTraits + { static const unsigned char position = 2; }; + template<> struct TypeTraits + { static const unsigned char position = 3; }; + template<> struct TypeTraits + { static const unsigned char position = 4; }; + template<> struct TypeTraits + { static const unsigned char position = 5; }; + template<> struct TypeTraits + { static const unsigned char position = 6; }; + template<> struct TypeTraits + { static const unsigned char position = 7; }; + template<> struct TypeTraits + { static const unsigned char position = 8; }; + template<> struct TypeTraits + { static const unsigned char position = 9; }; + template<> struct TypeTraits + { static const unsigned char position = 10; }; + template struct PromoteType + { typedef typename IfBool::position >= TypeTraits::position, A, B>::type type; }; + + // assertion testing + template inline void Assert(A assertion) { + if(!assertion) throw std::runtime_error("Assertion failed."); + } + + // geometry of views + template > class Geometry; + template + inline void stridesFromShape(ShapeIterator, ShapeIterator, + StridesIterator, const CoordinateOrder& = defaultOrder); + + // operations on entries of views + template + inline void operate(View&, Functor); + template + inline void operate(View&, const T&, Functor); + template + inline void operate(View&, const View&, Functor); + template + inline void operate(View& v, const ViewExpression& expression, Functor f); + + // helper classes + template + struct OperateHelperUnary; + template + struct OperateHelperBinaryScalar; + template + struct OperateHelperBinary; + template + struct AssignmentOperatorHelper; + template + struct AccessOperatorHelper; + + // unary in-place functors + template + struct Negative { void operator()(T& x) { x = -x; } }; + template + struct PrefixIncrement { void operator()(T& x) { ++x; } }; + template + struct PostfixIncrement { void operator()(T& x) { x++; } }; + template + struct PrefixDecrement { void operator()(T& x) { --x; } }; + template + struct PostfixDecrement { void operator()(T& x) { x--; } }; + + // binary in-place functors + template + struct Assign { void operator()(T1& x, const T2& y) { x = y; } }; + template + struct PlusEqual { void operator()(T1& x, const T2& y) { x += y; } }; + template + struct MinusEqual { void operator()(T1& x, const T2& y) { x -= y; } }; + template + struct TimesEqual { void operator()(T1& x, const T2& y) { x *= y; } }; + template + struct DividedByEqual { void operator()(T1& x, const T2& y) { x /= y; } }; + + // unary functors + template + struct Negate { T operator()(const T& x) const { return -x; } }; + + // binary functors + template + struct Plus { U operator()(const T1& x, const T2& y) const { return x + y; } }; + template + struct Minus { U operator()(const T1& x, const T2& y) const { return x - y; } }; + template + struct Times { U operator()(const T1& x, const T2& y) const { return x * y; } }; + template + struct DividedBy { U operator()(const T1& x, const T2& y) const { return x / y; } }; +} +// \endcond suppress_doxygen + +/// Array-Interface to an interval of memory. +/// +/// A view makes a subset of memory look as if it was stored in an +/// Marray. With the help of a view, data in a subset of memory can +/// be accessed and manipulated conveniently. In contrast to arrays +/// which allocate and de-allocate their own memory, views only +/// reference memory that has been allocated by other means. +/// Perhaps the simplest and most important use of views is to +/// read and manipulate sub-arrays. +/// +/// Notes on arithmetic operators of View: +/// - Only the pre-fix operators ++ and -- and not the corresponding post-fix +/// operators are implemented for View because the return value of the +/// post-fix operators would have to be the View as it is prior to the +/// operator call. However, the data under the view cannot be preserved when +/// incrementing or decrementing. Some compilers might accept the post-fix +/// operators, use the pre-fix implementation implicitly and issue a warning. +/// +template +class View +: public ViewExpression, T> +{ +public: + typedef std::size_t size_type; + typedef T value_type; + typedef typename marray_detail::IfBool::type pointer; + typedef const T* const_pointer; + typedef typename marray_detail::IfBool::type reference; + typedef const T& const_reference; + typedef Iterator iterator; + typedef Iterator const_iterator; + typedef std::reverse_iterator reverse_iterator; + typedef std::reverse_iterator const_reverse_iterator; + typedef ViewExpression, T> base; + typedef typename A::template rebind::other allocator_type; + + // construction + View(const allocator_type& = allocator_type()); + View(pointer, const allocator_type& = allocator_type()); + View(const View&); + template + View(ShapeIterator, ShapeIterator, pointer, + const CoordinateOrder& = defaultOrder, + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + template + View(ShapeIterator, ShapeIterator, StrideIterator, + pointer, const CoordinateOrder&, + const allocator_type& = allocator_type()); + #ifdef HAVE_CPP11_INITIALIZER_LISTS + View(std::initializer_list, pointer, + const CoordinateOrder& = defaultOrder, + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + View(std::initializer_list, std::initializer_list, + pointer, const CoordinateOrder&, + const allocator_type& = allocator_type()); + #endif + + // assignment + View& operator=(const T&); + View& operator=(const View&); // over-write default + View& operator=(const View&); // over-write default + template + View& operator=(const View&); + template + View& + operator=(const ViewExpression&); + + void assign(const allocator_type& = allocator_type()); + template + void assign(ShapeIterator, ShapeIterator, pointer, + const CoordinateOrder& = defaultOrder, + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + template + void assign(ShapeIterator, ShapeIterator, StrideIterator, + pointer, const CoordinateOrder&, + const allocator_type& = allocator_type()); + #ifdef HAVE_CPP11_INITIALIZER_LISTS + void assign(std::initializer_list, pointer, + const CoordinateOrder& = defaultOrder, + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + void assign(std::initializer_list, + std::initializer_list, pointer, + const CoordinateOrder&, + const allocator_type& = allocator_type()); + #endif + + // query + const std::size_t dimension() const; + const std::size_t size() const; + const std::size_t shape(const std::size_t) const; + const std::size_t* shapeBegin() const; + const std::size_t* shapeEnd() const; + const std::size_t strides(const std::size_t) const; + const std::size_t* stridesBegin() const; + const std::size_t* stridesEnd() const; + const CoordinateOrder& coordinateOrder() const; + const bool isSimple() const; + template + bool overlaps(const View&) const; + + + + // element access + #ifdef HAVE_CPP11_STD_ARRAY + + template + reference operator()(const std::array & c){ + return this->operator()(c[0]); + } + template + reference operator()(const std::array & c){ + return this->operator()(c[0], c[1]); + } + template + reference operator()(const std::array & c){ + return this->operator()(c[0], c[1], c[2]); + } + template + reference operator()(const std::array & c){ + return this->operator()(c[0], c[1], c[2], c[3]); + } + template + reference operator()(const std::array & c){ + return this->operator()(c[0], c[1], c[2], c[3], c[4]); + } + template + reference operator()(const std::array & c){ + return this->operator()(c.begin()); + } + + template + reference operator()(const std::array & c)const{ + return this->operator()(c[0]); + } + template + reference operator()(const std::array & c)const{ + return this->operator()(c[0], c[1]); + } + template + reference operator()(const std::array & c)const{ + return this->operator()(c[0], c[1], c[2]); + } + template + reference operator()(const std::array & c)const{ + return this->operator()(c[0], c[1], c[2], c[3]); + } + template + reference operator()(const std::array & c)const{ + return this->operator()(c[0], c[1], c[2], c[3], c[4]); + } + template + reference operator()(const std::array & c)const{ + return this->operator()(c.begin()); + } + #endif + + + + template reference operator()(U); + template reference operator()(U) const; + #ifndef HAVE_CPP11_VARIADIC_TEMPLATES + reference operator()(const std::size_t, const std::size_t); + reference operator()(const std::size_t, const std::size_t) const; + reference operator()(const std::size_t, const std::size_t, const std::size_t); + reference operator()(const std::size_t, const std::size_t, const std::size_t) const; + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t); + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t) const; + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t); + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t) const; + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t, const std::size_t); + reference operator()(const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t, const std::size_t, const std::size_t, + const std::size_t, const std::size_t, const std::size_t) const; + #else + reference operator()(const std::size_t); + reference operator()(const std::size_t) const; + template + reference operator()(const std::size_t, const Args...); + template + reference operator()(const std::size_t, const Args...) const; + private: + std::size_t elementAccessHelper(const std::size_t, const std::size_t); + std::size_t elementAccessHelper(const std::size_t, const std::size_t) const; + template + std::size_t elementAccessHelper(const std::size_t, const std::size_t, + const Args...); + template + std::size_t elementAccessHelper(const std::size_t, const std::size_t, + const Args...) const; + public: + #endif + + // sub-views + template + void view(BaseIterator, ShapeIterator, View&) const; + template + void view(BaseIterator, ShapeIterator, const CoordinateOrder&, + View&) const; + template + View view(BaseIterator, ShapeIterator) const; + template + View view(BaseIterator, ShapeIterator, + const CoordinateOrder&) const; + template + void constView(BaseIterator, ShapeIterator, View&) const; + template + void constView(BaseIterator, ShapeIterator, const CoordinateOrder&, + View&) const; + template + View constView(BaseIterator, ShapeIterator) const; + template + View constView(BaseIterator, ShapeIterator, + const CoordinateOrder&) const; + #ifdef HAVE_CPP11_INITIALIZER_LISTS + void view(std::initializer_list, + std::initializer_list, View&) const; + void view(std::initializer_list, + std::initializer_list, const CoordinateOrder&, + View&) const; + void constView(std::initializer_list, + std::initializer_list, View&) const; + void constView(std::initializer_list, + std::initializer_list, const CoordinateOrder&, + View&) const; + #endif + + // iterator access + iterator begin(); + iterator end(); + const_iterator begin() const; + const_iterator end() const; + reverse_iterator rbegin(); + reverse_iterator rend(); + const_reverse_iterator rbegin() const; + const_reverse_iterator rend() const; + + // coordinate transformation + template + void reshape(ShapeIterator, ShapeIterator); + template + void permute(CoordinateIterator); + void transpose(const std::size_t, const std::size_t); + void transpose(); + void shift(const int); + void squeeze(); + + template + View reshapedView(ShapeIterator, ShapeIterator) const; + template + View permutedView(CoordinateIterator) const; + View transposedView(const std::size_t, const std::size_t) const; + View transposedView() const; + View shiftedView(const int) const; + View boundView(const std::size_t, const std::size_t = 0) const; + View squeezedView() const; + + #ifdef HAVE_CPP11_INITIALIZER_LISTS + void reshape(std::initializer_list); + void permute(std::initializer_list); + + View reshapedView(std::initializer_list) const; + View permutedView(std::initializer_list) const; + #endif + + // conversion between coordinates, index and offset + template + void coordinatesToIndex(CoordinateIterator, std::size_t&) const; + template + void coordinatesToOffset(CoordinateIterator, std::size_t&) const; + template + void indexToCoordinates(std::size_t, CoordinateIterator) const; + void indexToOffset(std::size_t, std::size_t&) const; + #ifdef HAVE_CPP11_INITIALIZER_LISTS + void coordinatesToIndex(std::initializer_list, + std::size_t&) const; + void coordinatesToOffset(std::initializer_list, + std::size_t&) const; + #endif + + // output as string + std::string asString(const StringStyle& = MatrixStyle) const; + +private: + typedef typename marray_detail::Geometry geometry_type; + + View(pointer, const geometry_type&); + + void updateSimplicity(); + void testInvariant() const; + + // unsafe direct memory access + const T& operator[](const std::size_t) const; + T& operator[](const std::size_t); + + // data and memory address + pointer data_; + geometry_type geometry_; + +template + friend class View; +template + friend class Marray; +// \cond suppress_doxygen +template + friend struct marray_detail::AssignmentOperatorHelper; +friend struct marray_detail::AccessOperatorHelper; +friend struct marray_detail::AccessOperatorHelper; + +template + friend class ViewExpression; +template + friend class UnaryViewExpression; +template + friend class BinaryViewExpression; +template + friend class BinaryViewExpressionScalarFirst; +template + friend class BinaryViewExpressionScalarSecond; + +template + friend void marray_detail::operate(View& v, const ViewExpression& expression, Functor f); +// \endcond end suppress_doxygen +}; + +/// STL-compliant random access iterator for View and Marray. +/// +/// In addition to the STL iterator interface, the member functions +/// hasMore(), index(), and coordinate() are defined. +/// +template +class Iterator +{ +public: + // STL random access iterator typedefs + typedef typename std::random_access_iterator_tag iterator_category; + typedef T value_type; + typedef ptrdiff_t difference_type; + typedef typename marray_detail::IfBool::type pointer; + typedef typename marray_detail::IfBool::type reference; + + // non-standard typedefs + typedef typename marray_detail::IfBool*, + View*>::type view_pointer; + typedef typename marray_detail::IfBool&, + View&>::type view_reference; + + // construction + Iterator(); + Iterator(const View&, const std::size_t = 0); + Iterator(View&, const std::size_t = 0); + Iterator(const View&, const std::size_t = 0); + Iterator(const Iterator&); + // conversion from mutable to const + + // STL random access iterator operations + reference operator*() const; + pointer operator->() const; + reference operator[](const std::size_t) const; + Iterator& operator+=(const difference_type&); + Iterator& operator-=(const difference_type&); + Iterator& operator++(); // prefix + + Iterator& operator--(); // prefix + Iterator operator++(int); // postfix + Iterator operator--(int); // postfix + Iterator operator+(const difference_type&) const; + Iterator operator-(const difference_type&) const; + template + difference_type operator-(const Iterator&) const; + template + bool operator==(const Iterator&) const; + template + bool operator!=(const Iterator&) const; + template + bool operator<(const Iterator&) const; + template + bool operator>(const Iterator&) const; + template + bool operator<=(const Iterator&) const; + template + bool operator>=(const Iterator&) const; + + // operations beyond the STL standard + bool hasMore() const; + std::size_t index() const; + template + void coordinate(CoordinateIterator) const; + +private: + void testInvariant() const; + + // attributes + view_pointer view_; + pointer pointer_; + std::size_t index_; + std::vector coordinates_; + +friend class Marray; +friend class Iterator; // for comparison operators +}; + +/// Runtime-Flexible multi-dimensional array. +template +class Marray +: public View +{ +public: + typedef View base; + typedef typename base::value_type value_type; + typedef typename base::pointer pointer; + typedef typename base::const_pointer const_pointer; + typedef typename base::reference reference; + typedef typename base::const_reference const_reference; + typedef typename base::iterator iterator; + typedef typename base::reverse_iterator reverse_iterator; + typedef typename base::const_iterator const_iterator; + typedef typename base::const_reverse_iterator const_reverse_iterator; + typedef typename A::template rebind::other allocator_type; + + // constructors and destructor + Marray(const allocator_type& = allocator_type()); + Marray(const T&, const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + template + Marray(ShapeIterator, ShapeIterator, const T& = T(), + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + template + Marray(const InitializationSkipping&, ShapeIterator, ShapeIterator, + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + #ifdef HAVE_CPP11_INITIALIZER_LISTS + Marray(std::initializer_list, const T& = T(), + const CoordinateOrder& = defaultOrder, + const allocator_type& = allocator_type()); + #endif + Marray(const Marray&); + template + Marray(const ViewExpression&, + const allocator_type& = allocator_type()); + template + Marray(const View&); + ~Marray(); + + // assignment + Marray& operator=(const T&); + Marray& operator=(const Marray&); // over-write default + template + Marray& operator=(const View&); + template + Marray& operator=(const ViewExpression&); + void assign(const allocator_type& = allocator_type()); + + // resize + template + void resize(ShapeIterator, ShapeIterator, const T& = T()); + template + void resize(const InitializationSkipping&, ShapeIterator, ShapeIterator); + #ifdef HAVE_CPP11_INITIALIZER_LISTS + void resize(std::initializer_list, const T& = T()); + void resize(const InitializationSkipping&, std::initializer_list); + #endif + +private: + typedef typename base::geometry_type geometry_type; + + void testInvariant() const; + template + void resizeHelper(ShapeIterator, ShapeIterator, const T& = T()); + + allocator_type dataAllocator_; +}; + +// implementation of View + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Compute the index that corresponds to a sequence of coordinates. +/// +/// \param coordinate Coordinate given as initializer list. +/// \param out Index (output) +/// \sa coordinatesToOffset(), indexToCoordinates(), and indexToOffset() +/// +template +inline void +View::coordinatesToIndex +( + std::initializer_list coordinate, + std::size_t& out +) const +{ + coordinatesToIndex(coordinate.begin(), out); +} +#endif + +/// Compute the index that corresponds to a sequence of coordinates. +/// +/// \param it An iterator to the beginning of the coordinate sequence. +/// \param out Index (output) +/// \sa coordinatesToOffset(), indexToCoordinates(), and indexToOffset() +/// +template +template +inline void +View::coordinatesToIndex +( + CoordinateIterator it, + std::size_t& out +) const +{ + testInvariant(); + out = 0; + for(std::size_t j=0; jdimension(); ++j, ++it) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast(*it) < shape(j)); + out += static_cast(*it) * geometry_.shapeStrides(j); + } +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Compute the offset that corresponds to a sequence of coordinates. +/// +/// \param it An iterator to the beginning of the coordinate sequence. +/// \param out Index (output) +/// \sa coordinatesToIndex(), indexToCoordinates(), and indexToOffset() +/// +template +inline void +View::coordinatesToOffset +( + std::initializer_list coordinate, + std::size_t& out +) const +{ + coordinatesToOffset(coordinate.begin(), out); +} +#endif + +/// Compute the offset that corresponds to a sequence of coordinates. +/// +/// \param it An iterator to the beginning of the coordinate sequence. +/// \param out Index (output) +/// \sa coordinatesToIndex(), indexToCoordinates(), and indexToOffset() +/// +template +template +inline void +View::coordinatesToOffset +( + CoordinateIterator it, + std::size_t& out +) const +{ + testInvariant(); + out = 0; + for(std::size_t j=0; jdimension(); ++j, ++it) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast(*it) < shape(j)); + out += static_cast(*it) * strides(j); + } +} + +/// Compute the coordinate sequence that corresponds to an index. +/// +/// \param index Index +/// \param outit An iterator into a container into which the coordinate +/// sequence is to be written (output). +/// \sa coordinatesToIndex(), coordinatesToOffset(), and indexToOffset() +/// +template +template +inline void +View::indexToCoordinates +( + std::size_t index, // copy to work on + CoordinateIterator outit +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || this->dimension() > 0); + marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size()); + if(coordinateOrder() == FirstMajorOrder) { + for(std::size_t j=0; jdimension(); ++j, ++outit) { + *outit = std::size_t(index / geometry_.shapeStrides(j)); + index = index % geometry_.shapeStrides(j); + } + } + else { // last major order + std::size_t j = this->dimension()-1; + outit += j; + for(;;) { + *outit = std::size_t(index / geometry_.shapeStrides(j)); + index = index % geometry_.shapeStrides(j); + if(j == 0) { + break; + } + else { + --outit; + --j; + } + } + } +} + +/// Compute the offset that corresponds to an index. +/// +/// \param index Index. +/// \param out Offset (output). +/// \sa coordinatesToIndex(), coordinatesToOffset(), and indexToCoordinates() +/// +template +inline void +View::indexToOffset +( + std::size_t index, + std::size_t& out +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size()); + if(isSimple()) { + out = index; + } + else { + out = 0; + if(coordinateOrder() == FirstMajorOrder) { + for(std::size_t j=0; jdimension(); ++j) { + out += geometry_.strides(j) * (index / geometry_.shapeStrides(j)); + index = index % geometry_.shapeStrides(j); + } + } + else { // last major order + if(this->dimension() == 0) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0); + return; + } + else { + std::size_t j = this->dimension()-1; + for(;;) { + out += geometry_.strides(j) * (index / geometry_.shapeStrides(j)); + index = index % geometry_.shapeStrides(j); + if(j == 0) { + break; + } + else { + --j; + } + } + } + } + } +} + +/// Empty constructor. +/// +/// The empty constructor sets the data pointer to 0. +/// It does not allocate memory for a scalar. +/// +/// \param allocator Allocator. +/// +template +inline +View::View +( + const allocator_type& allocator +) +: data_(0), + geometry_(geometry_type(allocator)) +{ + testInvariant(); +} + +// private constructor +template +inline +View::View +( + pointer data, + const geometry_type& geometry +) +: data_(data), + geometry_(geometry) +{ + testInvariant(); +} + +/// Construct View from a scalar. +/// +/// \param data Pointer to data. +/// \param allocator Allocator. +/// +template +inline +View::View +( + pointer data, + const allocator_type& allocator +) +: data_(data), + geometry_(geometry_type(0, defaultOrder, 1, true, allocator)) +{ + testInvariant(); +} + +/// Construct View from a View on mutable data. +/// +/// \param in View on mutable data. +/// +template +inline +View::View +( + const View& in +) +: data_(in.data_), + geometry_(in.geometry_) +{ + testInvariant(); +} + +/// Construct unstrided View +/// +/// \param begin Iterator to the beginning of a sequence that +/// defines the shape. +/// \param end Iterator to the end of this sequence. +/// \param data Pointer to data. +/// \param externalCoordinateOrder Flag specifying the order +/// of coordinates based on which the strides are computed. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// + +template +template +inline +View::View +( + ShapeIterator begin, + ShapeIterator end, + pointer data, + const CoordinateOrder& externalCoordinateOrder, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +: data_(data), + geometry_(begin, end, externalCoordinateOrder, + internalCoordinateOrder, allocator) +{ + testInvariant(); +} + +/// Construct strided View +/// +/// \param begin Iterator to the beginning of a sequence that +/// defines the shape. +/// \param end Iterator to the end of this sequence. +/// \param it Iterator to the beginning of a sequence that +/// defines the strides. +/// \param data Pointer to data. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +template +inline +View::View +( + ShapeIterator begin, + ShapeIterator end, + StrideIterator it, + pointer data, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +: data_(data), + geometry_(begin, end, it, internalCoordinateOrder, allocator) +{ + testInvariant(); +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Construct unstrided View +/// +/// \param shape Shape initializer list. +/// \param data Pointer to data. +/// \param externalCoordinateOrder Flag specifying the order +/// of coordinates based on which the strides are computed. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +inline +View::View +( + std::initializer_list shape, + pointer data, + const CoordinateOrder& externalCoordinateOrder, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +: data_(data), + geometry_(shape.begin(), shape.end(), externalCoordinateOrder, + internalCoordinateOrder, allocator) +{ + testInvariant(); +} + +/// Construct strided View +/// +/// \param shape Shape initializer list. +/// \param strides Strides initializer list. +/// \param data Pointer to data. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// +template +inline +View::View +( + std::initializer_list shape, + std::initializer_list strides, + pointer data, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +: data_(data), + geometry_(shape.begin(), shape.end(), strides.begin(), + internalCoordinateOrder, allocator) +{ + testInvariant(); +} +#endif + +/// Clear View. +/// +/// Leaves the View in the same state as if the empty constructor +/// had been called. +/// +/// \sa View() +/// +template +inline void +View::assign +( + const allocator_type& allocator +) +{ + geometry_ = geometry_type(allocator); + data_ = 0; + testInvariant(); +} + +/// Initialize unstrided View +/// +/// \param begin Iterator to the beginning of a sequence that +/// defines the shape. +/// \param end Iterator to the end of this sequence. +/// \param data Pointer to data. +/// \param externalCoordinateOrder Flag specifying the order +/// of coordinates based on which the strides are computed. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +template +inline void +View::assign +( + ShapeIterator begin, + ShapeIterator end, + pointer data, + const CoordinateOrder& externalCoordinateOrder, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +{ + // the invariant is not tested as a pre-condition of this + // function to allow for unsafe manipulations prior to its + // call + geometry_ = typename marray_detail::Geometry(begin, end, + externalCoordinateOrder, internalCoordinateOrder, allocator); + data_ = data; + testInvariant(); +} + +/// Initialize strided View +/// +/// \param begin Iterator to the beginning of a sequence that +/// defines the shape. +/// \param end Iterator to the end of this sequence. +/// \param it Iterator to the beginning of a sequence that +/// defines the strides. +/// \param data Pointer to data. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +template +inline void +View::assign +( + ShapeIterator begin, + ShapeIterator end, + StrideIterator it, + pointer data, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +{ + // the invariant is not tested as a pre-condition of this + // function to allow for unsafe manipulations prior to its + // call + geometry_ = typename marray_detail::Geometry(begin, end, + it, internalCoordinateOrder, allocator); + data_ = data; + testInvariant(); +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Initialize unstrided View +/// +/// \param shape Shape initializer list. +/// \param data Pointer to data. +/// \param externalCoordinateOrder Flag specifying the order +/// of coordinates based on which the strides are computed. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +inline void +View::assign +( + std::initializer_list shape, + pointer data, + const CoordinateOrder& externalCoordinateOrder, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +{ + assign(shape.begin(), shape.end(), data, externalCoordinateOrder, + internalCoordinateOrder, allocator); +} + +/// Initialize strided View +/// +/// \param shape Shape initializer list. +/// \param strides Strides initialier list. +/// defines the strides. +/// \param data Pointer to data. +/// \param internalCoordinateOrder Flag specifying the order +/// of coordinates used for scalar indexing and iterators. +/// \param allocator Allocator. +/// +template +inline void +View::assign +( + std::initializer_list shape, + std::initializer_list strides, + pointer data, + const CoordinateOrder& internalCoordinateOrder, + const allocator_type& allocator +) +{ + assign(shape.begin(), shape.end(), strides.begin(), data, + internalCoordinateOrder, allocator); +} +#endif + +/// Reference data. +/// +/// \param u If u is an integer type, scalar indexing is performed. +/// Otherwise, it is assumed that u is an iterator to the beginning +/// of a coordinate sequence. +/// \return Reference to the entry at u. +/// +template +template +inline typename View::reference +View::operator() +( + U u +) +{ + return marray_detail::AccessOperatorHelper::is_integer>::execute(*this, u); +} + +/// Reference data. +/// +/// \param u If u is an integer type, scalar indexing is performed. +/// Otherwise, it is assumed that u is an iterator to the beginning +/// of a coordinate sequence. +/// \return Reference to the entry at u. +/// +template +template +inline typename View::reference +View::operator() +( + U u +) const +{ + return marray_detail::AccessOperatorHelper::is_integer>::execute(*this, u); +} + +#ifndef HAVE_CPP11_VARIADIC_TEMPLATES + +/// Reference data in a 2-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 2-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1))); + return data_[c0*strides(0) + c1*strides(1)]; +} + +/// Reference data in a 2-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 2-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1 +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1))); + return data_[c0*strides(0) + c1*strides(1)]; +} + +/// Reference data in a 3-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 3-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)]; +} + +/// Reference data in a 3-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 3-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2 +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)]; +} + +/// Reference data in a 4-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 4-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3)]; +} + +/// Reference data in a 4-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 4-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3 +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3)]; +} + +/// Reference data in a 5-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 5-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// \param c4 5th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3, + const std::size_t c4 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3) && c4 < shape(4))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3) + c4*strides(4)]; +} + +/// Reference data in a 5-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 5-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// \param c4 5th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3, + const std::size_t c4 +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3) && c4 < shape(4))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3) + c4*strides(4)]; +} + +/// Reference data in a 10-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not + +/// 5-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// \param c4 5th coordinate. +/// \param c5 6th coordinate. +/// \param c6 7th coordinate. +/// \param c7 8th coordinate. +/// \param c8 9th coordinate. +/// \param c9 10th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3, + const std::size_t c4, + const std::size_t c5, + const std::size_t c6, + const std::size_t c7, + const std::size_t c8, + const std::size_t c9 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3) && c4 < shape(4) + && c5 < shape(5) && c6 < shape(6) && c7 < shape(7) + && c8 < shape(8) && c9 < shape(9))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6) + + c7*strides(7) + c8*strides(8) + c9*strides(9)]; +} + +/// Reference data in a 10-dimensional View by coordinates. +/// +/// This function issues a runtime error if the View is not +/// 5-dimensional. +/// +/// \param c0 1st coordinate. +/// \param c1 2nd coordinate. +/// \param c2 3rd coordinate. +/// \param c3 4th coordinate. +/// \param c4 5th coordinate. +/// \param c5 6th coordinate. +/// \param c6 7th coordinate. +/// \param c7 8th coordinate. +/// \param c8 9th coordinate. +/// \param c9 10th coordinate. +/// +template +inline typename View::reference +View::operator() +( + const std::size_t c0, + const std::size_t c1, + const std::size_t c2, + const std::size_t c3, + const std::size_t c4, + const std::size_t c5, + const std::size_t c6, + const std::size_t c7, + const std::size_t c8, + const std::size_t c9 +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10)); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1) + && c2 < shape(2) && c3 < shape(3) && c4 < shape(4) + && c5 < shape(5) && c6 < shape(6) && c7 < shape(7) + && c8 < shape(8) && c9 < shape(9))); + return data_[c0*strides(0) + c1*strides(1) + c2*strides(2) + + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6) + + c7*strides(7) + c8*strides(8) + c9*strides(9)]; +} + +#else + +template +inline std::size_t +View::elementAccessHelper +( + const std::size_t Dim, + const std::size_t value +) +{ + marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) ); + return strides(Dim-1) * value; +} + +template +inline std::size_t +View::elementAccessHelper +( + const std::size_t Dim, + const std::size_t value +) const +{ + marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) ); + return strides(Dim-1) * value; +} + +template +template +inline std::size_t +View::elementAccessHelper +( + const std::size_t Dim, + const std::size_t value, + const Args... args +) +{ + marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) ); + return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...); +} + +template +template +inline std::size_t +View::elementAccessHelper +( + const std::size_t Dim, + const std::size_t value, + const Args... args +) const +{ + marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) ); + return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...); +} + +template +inline typename View::reference +View::operator() +( + const std::size_t value +) +{ + testInvariant(); + if(dimension() == 0) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0); + return data_[0]; + } + else { + std::size_t offset; + indexToOffset(value, offset); + return data_[offset]; + } +} + +template +inline typename View::reference +View::operator() +( + const std::size_t value +) const +{ + testInvariant(); + if(dimension() == 0) { + marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0); + return data_[0]; + } + else { + std::size_t offset; + indexToOffset(value, offset); + return data_[offset]; + } +} + +template +template +inline typename View::reference +View::operator() +( + const std::size_t value, + const Args... args +) +{ + testInvariant(); + marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) ); + return data_[strides(0)*value + elementAccessHelper(sizeof...(args)+1, args...) ]; +} + +template +template +inline typename View::reference +View::operator() +( + const std::size_t value, + const Args... args +) const +{ + testInvariant(); + marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) ); + return data_[ strides(0) * static_cast(value) + + static_cast(elementAccessHelper(sizeof...(args)+1, args...)) ]; +} + +#endif // #ifndef HAVE_CPP11_VARIADIC_TEMPLATES + +/// Get the number of data items. +/// +/// \return Size. +/// +template +inline const std::size_t +View::size() const +{ + return geometry_.size(); +} + +/// Get the dimension. +/// +/// Not well-defined if the data pointer is 0. +/// +/// \return Dimension. +/// +template +inline const std::size_t +View::dimension() const +{ + marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0); + return geometry_.dimension(); +} + +/// Get the shape in one dimension. +/// +/// \param dimension Dimension +/// \return Shape in that dimension. +/// +template +inline const std::size_t +View::shape +( + const std::size_t dimension +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension()); + return geometry_.shape(dimension); +} + +/// Get a constant iterator to the beginning of the shape vector. +/// +/// \return iterator. +/// \sa shapeEnd() +/// +template +inline const std::size_t* +View::shapeBegin() const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + return geometry_.shapeBegin(); +} + +/// Get a constant iterator to the end of the shape vector. +/// +/// \return iterator. +/// \sa shapeBegin() +/// +template +inline const std::size_t* +View::shapeEnd() const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + return geometry_.shapeEnd(); +} + +/// Get the strides in one dimension. +/// +/// \param dimension Dimension +/// \return Stride in that dimension. +/// +template +inline const std::size_t +View::strides +( + const std::size_t dimension +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension()); + return geometry_.strides(dimension); +} + +/// Get a constant iterator to the beginning of the strides vector. +/// +/// \return iterator. +/// \sa stridesEnd() +/// +template +inline const std::size_t* +View::stridesBegin() const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + return geometry_.stridesBegin(); +} + +/// Get a constant iterator to the end of the strides vector. +/// +/// \return iterator. +/// \sa stridesBegin() +/// +template +inline const std::size_t* +View::stridesEnd() const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + return geometry_.stridesEnd(); +} + +/// Get the coordinate order used for scalar indexing and iterators. +/// +/// \return CoordinateOrder. enum: FirstMajorOrder, LastMajorOrder +/// +template +inline const CoordinateOrder& +View::coordinateOrder() const +{ + testInvariant(); + return geometry_.coordinateOrder(); +} + +/// Determine whether the shape strides equal the strides of the View. +/// +/// \return bool. +/// +template +inline const bool +View::isSimple() const +{ + testInvariant(); + return geometry_.isSimple(); +} + +/// Assignment. +/// +/// operator= (the assignment operator) has a non-trivial behavior. +/// In most cases, it will work as most programmers will expect. +/// Here's a complete description of the semantics of to.operator=(from) +/// or equivalently, to = from. +/// +/// Consider the following cases: +/// (1) 'to' is mutable (isConst == false) +/// (a) 'from' is mutable (isConst == false) +/// (i) 'to' is initialized (data_ != 0) +/// (ii) 'to' is un-initialized (data_ == 0) +/// (b) 'from' is constant (isConst == true) +/// (2) 'to' is constant (isConst == true) +/// +/// (i) The operator attempts to copy the data under view 'b' to +/// the memory under view 'a'. This works if both views have the +/// same size, regardless of their dimension and shape. Equality +/// of sizes is checked by an assertion. +/// +/// (ii) Unless &a == &b (self-assignment), the operator copies +/// the (data) pointer of view 'b' to view 'a', without copying +/// the data itself. In addition, all the properties of view 'b' +/// are copied to view 'a'. +/// +/// (b) The operator attempts to copy the data under view 'b' to +/// the memory under view 'a'. This works if both views have the +/// same size, regardless of their dimension and shape. Equality +/// of sizes is checked by an assertion. If 'a' is un-initialized +/// the assertion fails (because the size of a will be zero). +/// Unlike in (ii), the pointer is not copied in this case. +/// Thus, a conversion from mutable to const is prevented. +/// +/// (2) Unless &a == &b (self-assignment), the operator copies +/// the (data) pointer of view 'b' to view 'a', without copying +/// the data itself. In addition, all the properties of view 'b' +/// are copied to view 'a'. Note that changing the data under +/// a constant view would be counter-intuitive. +/// +template +inline View& +View::operator= +( + const View& in +) +{ + testInvariant(); + marray_detail::AssignmentOperatorHelper::execute(in, *this); + testInvariant(); + return *this; +} + +/// Assignment. +/// +template +inline View& +View::operator= +( + const View& in +) +{ + testInvariant(); + marray_detail::AssignmentOperatorHelper::execute(in, *this); + testInvariant(); + return *this; +} + +/// Assignment. +/// +template +template +inline View& +View::operator= +( + const View& in +) +{ + testInvariant(); + marray_detail::AssignmentOperatorHelper::execute(in, *this); + testInvariant(); + return *this; +} + +/// Assignment. +/// +/// \param value Value. +/// +/// All entries are set to value. +/// +template +inline View& +View::operator= +( + const T& value +) +{ + marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0); + if(isSimple()) { + for(std::size_t j=0; j, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 2) + marray_detail::OperateHelperBinaryScalar<2, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 3) + marray_detail::OperateHelperBinaryScalar<3, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 4) + marray_detail::OperateHelperBinaryScalar<4, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 5) + marray_detail::OperateHelperBinaryScalar<5, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 6) + marray_detail::OperateHelperBinaryScalar<6, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 7) + marray_detail::OperateHelperBinaryScalar<7, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 8) + marray_detail::OperateHelperBinaryScalar<8, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 9) + marray_detail::OperateHelperBinaryScalar<9, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else if(dimension() == 10) + marray_detail::OperateHelperBinaryScalar<10, marray_detail::Assign, T, T, A>::operate(*this, value, marray_detail::Assign(), data_); + else { + for(iterator it = begin(); it.hasMore(); ++it) { + *it = value; + } + } + return *this; +} + +template +template +inline View& +View::operator= +( + const ViewExpression& expression +) +{ + marray_detail::operate(*this, expression, marray_detail::Assign()); + return *this; +} + +/// Get a sub-view with the same coordinate order. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param out Sub-View (output). +/// +template +template +inline void +View::view +( + BaseIterator bit, + ShapeIterator sit, + View& out +) const +{ + view(bit, sit, coordinateOrder(), out); +} + +/// Get a sub-view. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// \param out Sub-View (output). +/// +template +template +inline void +View::view +( + BaseIterator bit, + ShapeIterator sit, + const CoordinateOrder& internalCoordinateOrder, + View& out +) const +{ + testInvariant(); + std::size_t offset = 0; + coordinatesToOffset(bit, offset); + out.assign(sit, sit+dimension(), geometry_.stridesBegin(), + data_+offset, internalCoordinateOrder); +} + +/// Get a sub-view with the same coordinate order. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \return Sub-View. +/// +template +template +inline View +View::view +( + BaseIterator bit, + ShapeIterator sit +) const +{ + View v; + this->view(bit, sit, v); + return v; +} + +/// Get a sub-view. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// \return Sub-View. +/// +template +template +inline View +View::view +( + BaseIterator bit, + ShapeIterator sit, + const CoordinateOrder& internalCoordinateOrder +) const +{ + View v; + this->view(bit, sit, internalCoordinateOrder, v); + return v; +} + +/// Get a sub-view to constant data with the same coordinate +/// order. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param out Sub-View (output). +/// +template +template +inline void +View::constView +( + BaseIterator bit, + ShapeIterator sit, + View& out +) const +{ + constView(bit, sit, coordinateOrder(), out); +} + +/// Get a sub-view to constant data. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// \param out Sub-View (output). +/// +template +template +inline void +View::constView +( + BaseIterator bit, + ShapeIterator sit, + const CoordinateOrder& internalCoordinateOrder, + View& out +) const +{ + testInvariant(); + std::size_t offset = 0; + coordinatesToOffset(bit, offset); + out.assign(sit, sit+dimension(), + geometry_.stridesBegin(), + static_cast(data_) + offset, + internalCoordinateOrder); +} + +/// Get a sub-view to constant data with the same coordinate +/// order. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \return Sub-View. +/// +template +template +inline View +View::constView +( + BaseIterator bit, + ShapeIterator sit +) const +{ + View v; + this->constView(bit, sit, v); + return v; +} + +/// Get a sub-view to constant data. +/// +/// \param bit Iterator to the beginning of a coordinate sequence +/// that determines the start position of the sub-view. +/// \param sit Iterator to the beginning of a sequence +/// that determines the shape of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// \return Sub-View. +/// +template +template +inline View +View::constView +( + BaseIterator bit, + ShapeIterator sit, + const CoordinateOrder& internalCoordinateOrder +) const +{ + View v; + this->constView(bit, sit, internalCoordinateOrder, v); + return v; +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Get a sub-view. +/// +/// \param b Initializer list defining the coordinate sequence +/// that determines the start position of the sub-view. +/// \param s Initializer list defining the coordinate sequence +/// that determines the stop position of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// +template +inline void +View::view +( + std::initializer_list b, + std::initializer_list s, + const CoordinateOrder& internalCoordinateOrder, + View& out +) const +{ + view(b.begin(), s.begin(), internalCoordinateOrder, out); +} + +/// Get a sub-view with the same coordinate order. +/// +/// \param b Initializer list coordinate sequence +/// that determines the start position of the sub-view. +/// \param s Initializer list coordinate sequence +/// that determines the stop position of the sub-view. +/// \param out Sub-View (output). +/// +template +inline void +View::view +( + std::initializer_list b, + std::initializer_list s, + View& out +) const +{ + view(b.begin(), s.begin(), coordinateOrder(), out); +} + +/// Get a sub-view to constant data. +/// +/// \param b Initializer list coordinate sequence +/// that determines the start position of the sub-view. +/// \param s Initializer list coordinate sequence +/// that determines the stop position of the sub-view. +/// \param internalCoordinateOrder Flag to set the coordinate order +/// for scalar indexing and iterators of the sub-view. +/// +template +inline void +View::constView +( + std::initializer_list b, + std::initializer_list s, + const CoordinateOrder& internalCoordinateOrder, + View& out +) const +{ + constView(b.begin(), s.begin(), internalCoordinateOrder, out); +} + +/// Get a sub-view to constant data with the same coordinate +/// order. +/// +/// \param b Initializer list coordinate sequence +/// that determines the start position of the sub-view. +/// \param s Initializer list coordinate sequence +/// that determines the stop position of the sub-view. +/// \param out Sub-View (output). +/// +template +inline void +View::constView +( + std::initializer_list b, + std::initializer_list s, + View& out +) const +{ + constView(b.begin(), s.begin(), coordinateOrder(), out); +} +#endif + +/// Reshape the View. +/// +/// Two conditions have to be fulfilled in order for reshape to work: +/// - The new and the old shape must have the same size. +/// - The view must be simple, cf. isSimple(). +/// . +/// +/// \param begin Iterator to the beginning of a sequence that determines +/// the new shape. +/// \param end Iterator to the end of that sequence. +/// +/// \sa reshapedView(), isSimple() +/// +template +template +inline void +View::reshape +( + ShapeIterator begin, + ShapeIterator end +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || isSimple()); + if(!MARRAY_NO_ARG_TEST) { + std::size_t size = std::accumulate(begin, end, static_cast(1), + std::multiplies()); + marray_detail::Assert(size == this->size()); + } + assign(begin, end, data_, coordinateOrder(), coordinateOrder()); + testInvariant(); +} + +/// Get a reshaped View. +/// +/// Two conditions have to be fulfilled: +/// - The new and the old shape must have the same size. +/// - The view must be simple, cf. isSimple(). +/// . +/// +/// \param begin Iterator to the beginning of a sequence that determines +/// the new shape. +/// \param end Iterator to the end of that sequence. +/// +/// \sa reshape(), isSimple() +/// +template +template +inline View +View::reshapedView +( + ShapeIterator begin, + ShapeIterator end +) const +{ + View out = *this; + out.reshape(begin, end); + return out; +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Reshape the View. +/// +/// Two conditions have to be fulfilled in order for reshape to work: +/// - The new and the old shape must have the same size. +/// - The view must be simple, cf. isSimple(). +/// . +/// +/// \param shape Initializer list defining the new shape. +/// +/// \sa reshapedView(), isSimple() +/// +template +inline void +View::reshape +( + std::initializer_list shape +) +{ + reshape(shape.begin(), shape.end()); +} + +/// Get a reshaped View. +/// +/// Two conditions have to be fulfilled: +/// - The new and the old shape must have the same size. +/// - The view must be simple, cf. isSimple(). +/// . +/// +/// \param shape Initializer list defining the new shape. +/// +/// \sa reshape(), isSimple() +/// +template +inline View +View::reshapedView +( + std::initializer_list shape +) const +{ + return reshapedView(shape.begin(), shape.end()); +} +#endif + +/// Get a View where one coordinate is bound to a value. +/// +/// Binds one coordinate to a certain value. This reduces the +/// dimension by 1. +/// +/// \param dimension Dimension of the coordinate to bind. +/// \param value Value to assign to the coordinate. +/// \return The bound view. +/// \sa squeeze(), squeezeView() +/// +template +View +View::boundView +( + const std::size_t dimension, + const std::size_t value +) const +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_ARG_TEST || (dimension < this->dimension() + && value < shape(dimension))); + if(this->dimension() == 1) { + View v(&((*this)(value))); + v.geometry_.coordinateOrder() = coordinateOrder(); + return v; + } + else { + View v; + v.geometry_.resize(this->dimension()-1); + v.geometry_.coordinateOrder() = coordinateOrder(); + v.geometry_.size() = size() / shape(dimension); + for(std::size_t j=0, k=0; jdimension(); ++j) { + if(j != dimension) { + v.geometry_.shape(k) = shape(j); + v.geometry_.strides(k) = strides(j); + ++k; + } + } + marray_detail::stridesFromShape(v.geometry_.shapeBegin(), v.geometry_.shapeEnd(), + v.geometry_.shapeStridesBegin(), v.geometry_.coordinateOrder()); + v.data_ = data_ + strides(dimension) * value; + v.updateSimplicity(); + v.testInvariant(); + return v; + } +} + +/// Remove singleton dimensions by setting their coordinates to zero. +/// +/// \sa squeezedView(), boundView() +/// +template +void +View::squeeze() +{ + testInvariant(); + if(dimension() != 0) { + std::size_t newDimension = dimension(); + for(std::size_t j=0; j +inline View +View::squeezedView() const +{ + View v = *this; + v.squeeze(); + return v; +} + +#ifdef HAVE_CPP11_INITIALIZER_LISTS +/// Permute dimensions. +/// +/// \param begin Iterator to the beginning of a sequence which +/// has to contain the integers 0, ..., dimension()-1 in any +/// order. Otherwise, a runtime error is thrown. +/// \sa permutedView(), transpose(), transposedView(), shift(), +/// shiftedView() +/// +template +void +View::permute +( + std::initializer_list permutation +) +{ + permute(permutation.begin()); +} +#endif + +/// Permute dimensions. +/// +/// \param begin Iterator to the beginning of a sequence which +/// has to contain the integers 0, ..., dimension()-1 in any +/// order. Otherwise, a runtime error is thrown. +/// \sa permutedView(), transpose(), transposedView(), shift(), +/// shiftedView() +/// +template +template +void +View::permute +( + CoordinateIterator begin +) +{ + testInvariant(); + if(!MARRAY_NO_ARG_TEST) { + marray_detail::Assert(dimension() != 0); + std::set s1, s2; + CoordinateIterator it = begin; + for(std::size_t j=0; j newShape = std::vector(dimension()); + std::vector newStrides = std::vector(dimension()); + for(std::size_t j=0; j(*begin)); + newStrides[j] = geometry_.strides(static_cast(*begin)); + ++begin; + } + for(std::size_t j=0; j +template +inline View +View::permutedView +( + CoordinateIterator begin +) const +{ + View out = *this; + out.permute(begin); + return out; +} + +/// Exchange two dimensions. +/// +/// \param c1 Dimension +/// \param c2 Dimension +/// \sa permute(), permutedView(), shift(), shiftedView() +/// +template +void +View::transpose +( + const std::size_t c1, + const std::size_t c2 +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_ARG_TEST || + (dimension() != 0 && c1 < dimension() && c2 < dimension())); + + std::size_t j1 = c1; + std::size_t j2 = c2; + std::size_t c; + std::size_t d; + + // transpose shape + c = geometry_.shape(j2); + geometry_.shape(j2) = geometry_.shape(j1); + geometry_.shape(j1) = c; + + // transpose strides + d = geometry_.strides(j2); + geometry_.strides(j2) = geometry_.strides(j1); + geometry_.strides(j1) = d; + + // update shape strides + marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(), + geometry_.shapeStridesBegin(), geometry_.coordinateOrder()); + + updateSimplicity(); + testInvariant(); +} + +/// Reverse dimensions. +/// +/// \sa transposedView(), permute(), permutedView(), shift(), +/// shiftedView() +/// +template +void +View::transpose() +{ + testInvariant(); + for(std::size_t j=0; j +inline View +View::transposedView +( + const std::size_t c1, + const std::size_t c2 +) const +{ + View out = *this; + out.transpose(c1, c2); + return out; +} + +/// Get a View with dimensions reversed. +/// +/// \return View with dimensions reversed. +/// \sa transpose(), permute(), permutedView(), shift(), +/// shiftedView() +/// +template +inline View +View::transposedView() const +{ + View out = *this; + out.transpose(); + return out; +} + +/// Cycle shift dimensions. +/// +/// \param n Number of positions to shift +/// \sa shiftedView(), permute(), permutedView(), transpose(), +/// transposedView() +/// +template +inline void +View::shift +( + const int n +) +{ + testInvariant(); + marray_detail::Assert(MARRAY_NO_DEBUG || dimension() != 0); + if(n <= -static_cast(dimension()) || n >= static_cast(dimension())) { + shift(n % static_cast(dimension())); + } + else { + if(n > 0) { + shift(n - static_cast(dimension())); + } + else { + std::vector p(dimension()); + for(std::size_t j=0; j((static_cast(j) - n)) % dimension(); + } + permute(p.begin()); + } + } + testInvariant(); +} + +/// Get a View which dimensions cycle shifted. +/// +/// \param n Number of positions to shift +/// \sa shift(), permute(), permutedView(), transpose(), transposedView() +/// +template +inline View +View::shiftedView +( + const int n +) const +{ + View out = *this; + out.shift(n); + return out; +} + +/// Get an iterator to the beginning. +/// +/// \return Iterator. +/// \sa end() +/// + +template +inline typename View::iterator +View::begin() +{ + testInvariant(); + return Iterator(*this, 0); +} + +/// Get the end-iterator. +/// +/// \return Iterator. +/// \sa begin() +/// +template +inline typename View::iterator +View::end() +{ + testInvariant(); + return Iterator(*this, geometry_.size()); +} + +/// Get an iterator to the beginning. +/// +/// \return Iterator. +/// \sa end() +/// +template +inline typename View::const_iterator +View::begin() const +{ + testInvariant(); + return Iterator(*this, 0); +} + +/// Get the end-iterator. +/// +/// \return Iterator. +/// \sa begin() +/// +template +inline typename View::const_iterator + +View::end() const +{ + testInvariant(); + return Iterator(*this, geometry_.size()); +} + +/// Get a reserve iterator to the beginning. +/// +/// \return Iterator. +/// \sa rend() +/// +template +inline typename View::reverse_iterator +View::rbegin() +{ + return reverse_iterator(end()); +} + +/// Get the reverse end-iterator. +/// +/// \return Iterator. +/// \sa rbegin() +/// +template +inline typename View::reverse_iterator +View::rend() +{ + return reverse_iterator(begin()); +} + +/// Get a reserve iterator to the beginning. +/// +/// \return Iterator. +/// \sa rend() +/// +template +inline typename View::const_reverse_iterator +View::rbegin() const +{ + return const_reverse_iterator(end()); +} + +/// Get the reverse end-iterator. +/// +/// \return Iterator. +/// \sa rbegin() +/// +template +inline typename View::const_reverse_iterator +View::rend() const +{ + return const_reverse_iterator(begin()); +} + +/// Update Simplicity. +/// +/// This function sets the redundant boolean attribute isSimple_. +/// isSimple_ is set to true if the shape strides equal the +/// strides. +/// +template +inline void +View::updateSimplicity() +{ + // no invariant test here because this function + // is called during unsafe updates of a view + geometry_.updateSimplicity(); +} + +/// Unsafe direct memory access. +/// +/// This function provides direct access to the data items under the view +/// in the order in which these items reside in memory. +/// +/// \param offset offset to be added to the data pointer. +/// \return constant reference to the data item. +/// +template +inline const T& +View::operator[] +( + const std::size_t offset +) +const +{ + return data_[offset]; +} + +/// Unsafe direct memory access. +/// +/// This function provides direct access to the data items under the view +/// in the order in which these items reside in memory. +/// +/// \param offset offset to be added to the data pointer. +/// \return reference to the data item. +/// +template +inline T& +View::operator[] +( + const std::size_t offset +) +{ + return data_[offset]; +} + +/// Test invariant. +/// +/// This function tests the invariant of View and thus the consistency +/// of redundant information. +/// +template +inline void +View::testInvariant() const +{ + /* + if(!MARRAY_NO_DEBUG) { + std::cout<<"JOOO\n"; + if(geometry_.dimension() == 0) { + marray_detail::Assert(geometry_.isSimple() == true); + if(data_ != 0) { // scalar + marray_detail::Assert(geometry_.size() == 1); + } + } + else { + marray_detail::Assert(data_ != 0); + + // test size_ to be consistent with shape_ + std::size_t testSize = 1; + for(std::size_t j=0; j +template +inline bool View::overlaps +( + const View& v +) const +{ + testInvariant(); + if(!MARRAY_NO_ARG_TEST) { + v.testInvariant(); + } + if(data_ == 0 || v.data_ == 0) { + return false; + } + else { + const void* dataPointer_ = data_; + const void* vDataPointer_ = v.data_; + const void* maxPointer = & (*this)(this->size()-1); + const void* maxPointerV = & v(v.size()-1); + if( (dataPointer_ <= vDataPointer_ && vDataPointer_ <= maxPointer) + || (vDataPointer_ <= dataPointer_ && dataPointer_ <= maxPointerV) ) + { + return true; + } + } + return false; +} + +/// Output as string. +/// +template +std::string +View::asString +( + const StringStyle& style +) const +{ + testInvariant(); + std::ostringstream out(std::ostringstream::out); + if(style == MatrixStyle) { + if(dimension() == 0) { + // scalar + out << "A = " << (*this)(0) << std::endl; + } + else if(dimension() == 1) { + // vector + out << "A = ("; + for(std::size_t j=0; jsize(); ++j) { + out << (*this)(j) << ", "; + } + out << "\b\b)" << std::endl; + } + else if(dimension() == 2) { + // matrix + if(coordinateOrder() == FirstMajorOrder) { + out << "A(r,c) =" << std::endl; + for(std::size_t y=0; yshape(0); ++y) { + for(std::size_t x=0; xshape(1); ++x) { + out << (*this)(y, x) << ' '; + } + out << std::endl; + } + } + else { + out << "A(c,r) =" << std::endl; + for(std::size_t y=0; yshape(1); ++y) { + for(std::size_t x=0; xshape(0); ++x) { + out << (*this)(x, y) << ' '; + } + out << std::endl; + } + } + } + else { + // higher dimensional + std::vector c1(dimension()); + std::vector c2(dimension()); + unsigned short q = 2; + if(coordinateOrder() == FirstMajorOrder) { + q = static_cast(dimension() - 3); + } + for(const_iterator it = this->begin(); it.hasMore(); ++it) { + it.coordinate(c2.begin()); + if(it.index() == 0 || c2[q] != c1[q]) { + if(it.index() != 0) { + out << std::endl << std::endl; + } + if(coordinateOrder() == FirstMajorOrder) { + out << "A("; + for(std::size_t j=0; j -#include -#include + + +#include "nifty/marray/marray.hxx" +#include "nifty/marray/andres/marray-hdf5.hxx" +#include "nifty/marray/andres/marray.hxx" namespace nifty{ namespace marray{ using namespace andres; - } } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f554f095c..92084fd3b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,6 +16,6 @@ if(BUILD_CPP_EXAMPLES) endif() #add_subdirectory(sandbox) -add_subdirectory(benchmark) +#add_subdirectory(benchmark) diff --git a/src/benchmark/CMakeLists.txt b/src/benchmark/CMakeLists.txt deleted file mode 100644 index c3f48bb0d..000000000 --- a/src/benchmark/CMakeLists.txt +++ /dev/null @@ -1,2 +0,0 @@ -add_executable(marray marray.cxx ) -add_executable(vs_varray vs_varray.cxx ) \ No newline at end of file diff --git a/src/benchmark/marray.cxx b/src/benchmark/marray.cxx deleted file mode 100644 index 5c6898538..000000000 --- a/src/benchmark/marray.cxx +++ /dev/null @@ -1,129 +0,0 @@ -#include -#include "nifty/tools/timer.hxx" -#include "nifty/marray/marray.hxx" -#include "vigra/multi_array.hxx" - - - -int main( int argc , char *argv[] ){ - - - - std::vector shape({100000000}); - nifty::marray::Marray a(shape.begin(), shape.end(),0, nifty::marray::LastMajorOrder); - nifty::marray::View v = a; - - nifty::tools::VerboseTimer timer(true); - - - { - const auto p = &v(0); - timer.startAndPrint("ptr"); - float s = 0; - for(auto i=0; i vv(vShape, p); - - timer.startAndPrint("vigra"); - float s = 0; - for(auto i=0; i vv(vShape, p); - - timer.startAndPrint("vigra"); - float s = 0; - for(auto i=0; i -#include - -#include -#include - - - -inline void rowMajorAcces(andres::View & view) { - for( int i = 0; i < view.shape(0); ++i ) { - for( int j = 0; j < view.shape(1); ++j ) { - view(j,i) = i + j; - } - } -} - -template -inline void marray1D(andres::View & view) { - for( int i = 0; i < view.shape(0); ++i ) { - view(i) = i; - } -} - -inline void columnMajorAcces(andres::View & view) { - for( int j = 0; j < view.shape(1); ++j ) { - for( int i = 0; i < view.shape(0); ++i ) { - view(j,i) = i + j; - } - } -} - - -inline void rowMajorAcces(vigra::MultiArray<2,int> & varray) { - for( int i = 0; i < varray.shape(0); ++i ) { - for( int j = 0; j < varray.shape(1); ++j ) { - varray(j,i) = i + j; - } - } -} - -template -inline void vigra1D(vigra::MultiArrayView<1,T> & varray) { - for( int i = 0; i < varray.shape(0); ++i ) { - varray(i) = i; - } -} - - - - -inline void columnMajorAcces(vigra::MultiArray<2,int> & varray) { - for( int j = 0; j < varray.shape(0); ++j ) { - for( int i = 0; i < varray.shape(1); ++i ) { - varray(j,i) = i + j; - } - } -} - - -void time2D(const int N) { - - size_t shape[] = {3000,3000}; - andres::Marray array(shape, shape+2); - - std::cout << "Timeing rowMajorAcces..." << std::endl; - auto t0 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - rowMajorAcces(array); - auto t1 = std::chrono::high_resolution_clock::now(); - - auto dur0 = std::chrono::duration_cast(t1 - t0); - std::cout << "... in total: " << dur0.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur0.count() / N << " ms" << std::endl; - - vigra::Shape2 vshape(3000,3000); - vigra::MultiArray<2,int> varray(vshape); - - std::cout << "Timeing rowMajorAcces Vigra..." << std::endl; - auto t00 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - rowMajorAcces(varray); - auto t10 = std::chrono::high_resolution_clock::now(); - - auto dur00 = std::chrono::duration_cast(t10 - t00); - std::cout << "... in total: " << dur00.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur00.count() / N << " ms" << std::endl; - - //std::cout << "Timeing columnMajorAcces..." << std::endl; - //auto t2 = std::chrono::high_resolution_clock::now(); - //for(int _ = 0; _ < N; ++_) - // columnMajorAcces(array); - //auto t3 = std::chrono::high_resolution_clock::now(); - - //auto dur1 = std::chrono::duration_cast(t3 - t2); - //std::cout << "... in total: " << dur1.count() << " ms" << std::endl; - //std::cout << "... per itearion " << dur1.count() / N << " ms" << std::endl; - // - //std::cout << "Timeing columnMajorAcces Vigra..." << std::endl; - //auto t20 = std::chrono::high_resolution_clock::now(); - //for(int _ = 0; _ < N; ++_) - // columnMajorAcces(varray); - //auto t30 = std::chrono::high_resolution_clock::now(); - - //auto dur10 = std::chrono::duration_cast(t30 - t20); - //std::cout << "... in total: " << dur10.count() << " ms" << std::endl; - //std::cout << "... per itearion " << dur10.count() / N << " ms" << std::endl; -} - - -void time1D(const int N) { - - size_t shape[] = {long(1e8)}; - andres::Marray array(shape, shape+1); - - std::cout << "Timeing 1d accces..." << std::endl; - auto t0 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - marray1D(array); - auto t1 = std::chrono::high_resolution_clock::now(); - - auto dur0 = std::chrono::duration_cast(t1 - t0); - std::cout << "... in total: " << dur0.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur0.count() / N << " ms" << std::endl; - - vigra::Shape1 vshape(1e8); - vigra::MultiArray<1,int> varray(vshape); - - std::cout << "Timeing 1D acces Vigra..." << std::endl; - auto t00 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - vigra1D(varray); - auto t10 = std::chrono::high_resolution_clock::now(); - - auto dur00 = std::chrono::duration_cast(t10 - t00); - std::cout << "... in total: " << dur00.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur00.count() / N << " ms" << std::endl; -} - -template -void time1DSameData(const int N) { - - // data array - size_t shape[] = {long(1e8)}; - T * data = new T[shape[0]]; - - - andres::View array(shape, shape+1, data); - - std::cout << "Timeing 1d accces..." << std::endl; - auto t0 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - marray1D(array); - auto t1 = std::chrono::high_resolution_clock::now(); - - auto dur0 = std::chrono::duration_cast(t1 - t0); - std::cout << "... in total: " << dur0.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur0.count() / N << " ms" << std::endl; - - vigra::Shape1 vshape(1e8); - vigra::MultiArrayView<1,T> varray(vshape, data); - - std::cout << "Timeing 1D acces Vigra..." << std::endl; - auto t00 = std::chrono::high_resolution_clock::now(); - for(int _ = 0; _ < N; ++_) - vigra1D(varray); - auto t10 = std::chrono::high_resolution_clock::now(); - - auto dur00 = std::chrono::duration_cast(t10 - t00); - std::cout << "... in total: " << dur00.count() << " ms" << std::endl; - std::cout << "... per itearion " << dur00.count() / N << " ms" << std::endl; - - delete[] data; -} - -int main() { - time1DSameData(25); - //time2D(25); -}