Tracks KLT features in forward/reverse direction and the tracking error metrics
497152
protected boolean trackFeature() {

		pairs.reset();
		// total number of tracks which contribute to FB error
		int numTracksFB = 0;
		// tracks which are not dropped
		int numTracksRemaining = 0;

		for( int i = 0; i < tracks.length; i++ ) {
			Track t = tracks[i];
			if( !t.active )
				continue;

			float prevX = t.klt.x;
			float prevY = t.klt.y;

			// track in forwards direction
			tracker.setImage(currentImage,currentDerivX,currentDerivY);
			KltTrackFault result = tracker.track(t.klt);
			if( result != KltTrackFault.SUCCESS ) {
				t.active = false;
				continue;
			}

			float currX = t.klt.x;
			float currY = t.klt.y;

			// track in reverse direction
			tracker.setDescription(t.klt);
			tracker.setImage(previousImage, previousDerivX, previousDerivY);
			result = tracker.track(t.klt);
			if( result != KltTrackFault.SUCCESS ) {
				t.active = false;
				continue;
			}

			// compute forward-backwards error
			double errorForwardBackwards = UtilPoint2D_F32.distanceSq(prevX,prevY,t.klt.x,t.klt.y);

			// put into lists for computing the median error
			errorsFB[numTracksFB++] = errorForwardBackwards;

			// discard if error is too large
			if( errorForwardBackwards > maxErrorFB ) {
				t.active = false;
				continue;
			}

			// create data structure used for group motion estimation
			AssociatedPair p = pairs.grow();
			p.p1.set( prevX, prevY );
			p.p2.set( currX, currY );

			numTracksRemaining++;
		}

		// if the forward-backwards error is too large, give up
		double medianFB = QuickSelect.select(errorsFB,numTracksFB/2,numTracksFB);

//		System.out.println("Median tracking error FB: "+medianFB);

		if( medianFB > maxErrorFB || numTracksRemaining < 4 )
			return false;

		return true;
	}
1266537
def doExperiment(numObjects,
                 featuresPerObject,
                 objectWidth,
                 numFeatures,
                 featureDistribution,
                 useTrace,
                 noiseFactor,
                 moduleNoiseFactor,
                 numModules,
                 thresholds,
                 inverseReadoutResolution,
                 enlargeModuleFactor,
                 bumpOverlapMethod,
                 seed1,
                 seed2):
  
  if not os.path.exists("traces"):
    os.makedirs("traces")

  if seed1 != -1:
    np.random.seed(seed1)

  if seed2 != -1:
    random.seed(seed2)

  features = [str(i) for i in xrange(numFeatures)]
  objects = generateObjects(numObjects, featuresPerObject, objectWidth,
                            numFeatures, featureDistribution)

  locationConfigs = []
  scale = 40.0

  if thresholds == -1:
    thresholds = int(math.ceil(numModules*0.8))
  elif thresholds == 0:
    thresholds = numModules
  perModRange = float(60.0 / float(numModules))
  for i in xrange(numModules):
    orientation = (float(i) * perModRange) + (perModRange / 2.0)

    locationConfigs.append({
        "scale": scale,
        "orientation": np.radians(orientation),
        "activationThreshold": 8,
        "initialPermanence": 1.0,
        "connectedPermanence": 0.5,
        "learningThreshold": 8,
        "sampleSize": 10,
        "permanenceIncrement": 0.1,
        "permanenceDecrement": 0.0,
        "inverseReadoutResolution": inverseReadoutResolution,
        "enlargeModuleFactor": enlargeModuleFactor,
        "bumpOverlapMethod": bumpOverlapMethod,
    })

  l4Overrides = {
    "initialPermanence": 1.0,
    "activationThreshold": thresholds,
    "reducedBasalThreshold": thresholds,
    "minThreshold": numModules,
    "sampleSize": numModules,
    "cellsPerColumn": 16,
  }

  column = PIUNCorticalColumn(locationConfigs, L4Overrides=l4Overrides,
                              bumpType="gaussian2")
  exp = PIUNExperiment(column, featureNames=features,
                       numActiveMinicolumns=10,
                       noiseFactor=noiseFactor,
                       moduleNoiseFactor=moduleNoiseFactor)

  for objectDescription in objects:
    exp.learnObject(objectDescription)

  convergence = collections.defaultdict(int)
  try:
    if useTrace:
      filename = os.path.join(
        SCRIPT_DIR,
        "traces/{}-resolution-{}-modules-{}-objects-{}-feats.html".format(
          inverseReadoutResolution, numModules, numObjects, numFeatures)
      )
      traceFileOut = io.open(filename, "w", encoding="utf8")
      traceHandle = trace(traceFileOut, exp, includeSynapses=True)
      print "Logging to", filename

    for objectDescription in objects:
      steps = exp.inferObjectWithRandomMovements(objectDescription)
      convergence[steps] += 1
      if steps is None:
        print 'Failed to infer object "{}"'.format(objectDescription["name"])
  finally:
    if useTrace:
      traceHandle.__exit__()
      traceFileOut.close()

  for step, num in sorted(convergence.iteritems()):
    print "{}: {}".format(step, num)

  result = {
    "convergence": convergence,
  }

  return result
1144437
def get_monitor_forward_and_backward(summary_writer, bins):
    
    def monitor_forward_and_backward(module, input, output):
        """
        Iterate over the module parameters and monitor their forward values.
        Then iterate over all of the monitored_vars, monitor their forward values
        and set their grad_hooks
        """
        # Parameters
        if module.is_monitoring:
            param_names = [ name for name, _ in module.named_parameters()]
            for name, param in zip(param_names, module.parameters()):
                if module.track_grad and param.requires_grad:
                    hook = grad_hook(module, name, summary_writer, bins)
                    module.param_hooks[name] = param.register_hook(hook)
                if module.track_data:
                    summary_writer.add_histogram('{}/data'.format(name.replace('.','/')),
                                                 param.detach().cpu().numpy(),
                                                 module.global_step,
                                                 bins=bins)
                if name in module.last_state_dict:
                    if module.track_update:
                        update = param - module.last_state_dict[name]
                        summary_writer.add_histogram('{}/update-val'.format(name.replace('.','/')),
                                                     update.detach().cpu().numpy(),
                                                     module.global_step-1,
                                                     bins=bins)
                    if module.track_update and module.track_update_ratio:
                        update_ratio = update / (module.last_state_dict[name]+1e-15)

                        summary_writer.add_histogram('{}/update-ratio'.format(name.replace('.','/')),
                                                     update_ratio.detach().cpu().numpy(),
                                                     module.global_step-1,
                                                     bins=bins)
                if module.track_update:
                    module.last_state_dict[name] = param.clone()

            # Intermediate Vars
            for prefix, mod in module.named_modules():
                for tensor_name, entry in mod.monitored_vars.items():
                    name = '{}/{}'.format(prefix, tensor_name) if prefix else tensor_name
                    tensor = entry['tensor']
                    if entry['track_grad'] and tensor.requires_grad:
                        hook = grad_hook(module, name, summary_writer, bins)
                        module.var_hooks[name] = tensor.register_hook(hook)
                    if entry['track_data']:
                        summary_writer.add_histogram('{}/data'.format(name.replace('.','/')),
                                                     tensor.detach().cpu().numpy(),
                                                     module.global_step,
                                                     bins=bins)
            # Update step
        module.global_step += 1
    return monitor_forward_and_backward
5440
private boolean estimateMotion() {
		// organize the data
		List<Stereo2D3D> data = new ArrayList<>();

		for( PointTrack l : candidates ) {
			LeftTrackInfo info = l.getCookie();
			PointTrack r = info.right;

			Stereo2D3D stereo = info.location;
			// compute normalized image coordinate for track in left and right image
			leftImageToNorm.compute(l.x,l.y,info.location.leftObs);
			rightImageToNorm.compute(r.x,r.y,info.location.rightObs);

			data.add(stereo);
		}

		// Robustly estimate left camera motion
		if( !matcher.process(data) )
			return false;

		Se3_F64 keyToCurr = matcher.getModelParameters();
		keyToCurr.invert(currToKey);

		// mark tracks that are in the inlier set
		int N = matcher.getMatchSet().size();
		for( int i = 0; i < N; i++ ) {
			int index = matcher.getInputIndex(i);
			LeftTrackInfo info = candidates.get(index).getCookie();
			info.lastInlier = tick;
		}

//		System.out.println("Inlier set size: "+N);

		return true;
	}
916952
private void saveInlierMatches(ModelMatcher<?, ?> ransac,
								   FastQueue<AssociatedIndex> matches, PairwiseImageGraph2.Motion edge) {

		int N = ransac.getMatchSet().size();
		edge.inliers.reset();
		for (int i = 0; i < N; i++) {
			int idx = ransac.getInputIndex(i);
			edge.inliers.grow().set(matches.get(idx));
		}
	}
900959
public boolean update( Image input , RectangleRotate_F64 output ) {

		if( trackLost )
			return false;

		trackFeatures(input, region);

		// See if there are enough points remaining.  use of config.numberOfSamples is some what arbitrary
		if( pairs.size() < config.numberOfSamples ) {
			System.out.println("Lack of sample pairs");
			trackLost = true;
			return false;
		}

		// find the motion using tracked features
		if( !estimateMotion.process(pairs.toList()) ) {
			System.out.println("estimate motion failed");
			trackLost = true;
			return false;
		}

		if( estimateMotion.getFitQuality() > config.robustMaxError ) {
			System.out.println("exceeded Max estimation error");
			trackLost = true;
			return false;
		}

		// update the target's location using the found motion
		ScaleTranslateRotate2D model = estimateMotion.getModelParameters();

		region.width *= model.scale;
		region.height *= model.scale;

		double c = Math.cos(model.theta);
		double s = Math.sin(model.theta);

		double x = region.cx;
		double y = region.cy;

		region.cx = (x*c - y*s)*model.scale + model.transX;
		region.cy = (x*s + y*c)*model.scale + model.transY;

		region.theta += model.theta;

		output.set(region);

		// make the current image into the previous image
		swapImages();

		return true;
	}
904406
private void extractImageFeatures(final ProgressMonitor progressMonitor, final int progress,
									  T image,
									  List<TupleDesc> descs, List<Point2D_F64> locs) {
		SwingUtilities.invokeLater(new Runnable() {
			public void run() {
				progressMonitor.setNote("Detecting");
			}
		});
		detector.detect(image);
		SwingUtilities.invokeLater(new Runnable() {
			public void run() {
				progressMonitor.setProgress(progress + 1);
				progressMonitor.setNote("Describing");
			}
		});
		describe.setImage(image);
		orientation.setImage(image);

		// See if the detector can detect the feature's scale
		if (detector.hasScale()) {
			for (int i = 0; i < detector.getNumberOfFeatures(); i++) {
				double yaw = 0;

				Point2D_F64 pt = detector.getLocation(i);
				double radius = detector.getRadius(i);
				if (describe.requiresOrientation()) {
					orientation.setObjectRadius(radius);
					yaw = orientation.compute(pt.x, pt.y);
				}

				TupleDesc d = describe.createDescription();
				if ( describe.process(pt.x, pt.y, yaw, radius, d) ) {
					descs.add(d);
					locs.add(pt.copy());
				}
			}
		} else {
			// just set the radius to one in this case
			orientation.setObjectRadius(1);
			for (int i = 0; i < detector.getNumberOfFeatures(); i++) {
				double yaw = 0;

				Point2D_F64 pt = detector.getLocation(i);
				if (describe.requiresOrientation()) {
					yaw = orientation.compute(pt.x, pt.y);
				}

				TupleDesc d = describe.createDescription();
				if (describe.process(pt.x, pt.y, yaw, 1, d)) {
					descs.add(d);
					locs.add(pt.copy());
				}
			}
		}
		SwingUtilities.invokeLater(new Runnable() {
			public void run() {
				progressMonitor.setProgress(progress + 2);
			}
		});
	}
785699
def tenengrad(img, ksize=3):
    '''''''
    Gx = cv2.Sobel(img, ddepth=cv2.CV_64F, dx=1, dy=0, ksize=ksize)
    Gy = cv2.Sobel(img, ddepth=cv2.CV_64F, dx=0, dy=1, ksize=ksize)
    FM = Gx*Gx + Gy*Gy
    mn = cv2.mean(FM)[0]
    if np.isnan(mn):
        return np.nanmean(FM)
    return mn
1092159
private void selectCandidateTracks() {
		// mark tracks in right frame that are active
		List<PointTrack> activeRight = trackerRight.getActiveTracks(null);
		for( PointTrack t : activeRight ) {
			RightTrackInfo info = t.getCookie();
			info.lastActiveList = tick;
		}

		int mutualActive = 0;
		List<PointTrack> activeLeft = trackerLeft.getActiveTracks(null);
		candidates.clear();
		for( PointTrack left : activeLeft ) {
			LeftTrackInfo info = left.getCookie();

//			if( info == null || info.right == null ) {
//				System.out.println("Oh Crap");
//			}

			// for each active left track, see if its right track has been marked as active
			RightTrackInfo infoRight = info.right.getCookie();
			if( infoRight.lastActiveList != tick ) {
				continue;
			}

			// check epipolar constraint and see if it is still valid
			if( stereoCheck.checkPixel(left, info.right) ) {
				info.lastConsistent = tick;
				candidates.add(left);
			}
			mutualActive++;
		}

//		System.out.println("Active Tracks: Left "+trackerLeft.getActiveTracks(null).size()+" right "+
//				trackerRight.getActiveTracks(null).size());
//		System.out.println("All Tracks:    Left "+trackerLeft.getAllTracks(null).size()+" right "+
//				trackerRight.getAllTracks(null).size());
//		System.out.println("Candidates = "+candidates.size()+" mutual active = "+mutualActive);
	}
1288312
def kl_error(self, input_data, targets, average=True,
                 cache=None, prediction=True):
        

        if cache is not None:
            activations = cache
        else:
            activations = \
              self.feed_forward(input_data, prediction=prediction)

        targets_non_nan = gpuarray.empty_like(targets)
        nan_to_zeros(targets, targets_non_nan)
        kl_error = gpuarray.sum(targets_non_nan *
                                (cumath.log(targets_non_nan + eps) -
                                 cumath.log(activations + eps)))
        if average:
            kl_error /= targets.shape[0]
        return kl_error.get()
468766
def _backtrace(elt, dom):
    '''
    '''
    s = ''
    while elt != dom:
        name, parent = elt.nodeName, elt.parentNode
        if parent is None: break
        matches = [ c for c in _child_elements(parent) 
                        if c.nodeName == name ]
        if len(matches) == 1:
            s = '/' + name + s
        else:
            i = matches.index(elt) + 1
            s = ('/%s[%d]' % (name, i)) + s
        elt = parent
    return s
529461
def _tracker(self, _observer_, _self_, *args, **kwds):
        
        self.track_object(_self_,
                          name=_observer_.name,
                          resolution_level=_observer_.detail,
                          keep=_observer_.keep,
                          trace=_observer_.trace)
        _observer_.init(_self_, *args, **kwds)
1129914
def inferObjectsWithRandomMovements(self, objectPlacements, maxTouches=20,
                                      settlingTime=2):
    

    for monitor in self.monitors.values():
      monitor.afterPlaceObjects(objectPlacements)

    for objectName, objectDict in self.objects.iteritems():
      self.reset()

      visitCounts = defaultdict(int)

      learnedObjectPlacement = self.learnedObjectPlacements[objectName]

      allocentricLocation = random.choice(objectDict.keys())
      nextTransitionSDR = np.empty(0, dtype="uint32")

      # Traverse the object until it is inferred.
      success = False

      for _ in xrange(maxTouches):
        featureName = objectDict[allocentricLocation]
        egocentricLocation = (allocentricLocation[0] +
                              objectPlacements[objectName][0],
                              allocentricLocation[1] +
                              objectPlacements[objectName][1])
        featureSDR = self.features[featureName]

        steps = ([nextTransitionSDR] +
                 [np.empty(0)]*settlingTime)
        for transitionSDR in steps:
          self.doTimestep(np.empty(0), transitionSDR, featureSDR,
                          egocentricLocation, learn=False)

        visitCounts[allocentricLocation] += 1

        # We should eventually infer the egocentric location where we originally
        # learned this location on the object.
        learnedEgocentricLocation = (
          allocentricLocation[0] + learnedObjectPlacement[0],
          allocentricLocation[1] + learnedObjectPlacement[1])

        if (set(self.objectLayer.getActiveCells()) ==
            set(self.objectRepresentations[objectName]) and

            set(self.inputLayer.getActiveCells()) ==
            set(self.inputRepresentations[(featureName,
                                           learnedEgocentricLocation)]) and

            set(self.locationLayer.getActiveCells()) ==
            set(self.locations[learnedEgocentricLocation])):
          success = True
          break
        else:
          allocentricLocation, transitionName = self._selectTransition(
            allocentricLocation, objectDict, visitCounts)
          nextTransitionSDR = self.transitions[transitionName]
364802
private void associateL2R( T left , T right ) {
		// make the previous new observations into the new old ones
		ImageInfo<TD> tmp = featsLeft1;
		featsLeft1 = featsLeft0; featsLeft0 = tmp;
		tmp = featsRight1;
		featsRight1 = featsRight0; featsRight0 = tmp;

		// detect and associate features in the two images
		featsLeft1.reset();
		featsRight1.reset();

//		long time0 = System.currentTimeMillis();
		describeImage(left,featsLeft1);
		describeImage(right,featsRight1);
//		long time1 = System.currentTimeMillis();

		// detect and associate features in the current stereo pair
		for( int i = 0; i < detector.getNumberOfSets(); i++ ) {
			SetMatches matches = setMatches[i];
			matches.swap();
			matches.match2to3.reset();

			FastQueue<Point2D_F64> leftLoc = featsLeft1.location[i];
			FastQueue<Point2D_F64> rightLoc = featsRight1.location[i];

			assocL2R.setSource(leftLoc,featsLeft1.description[i]);
			assocL2R.setDestination(rightLoc, featsRight1.description[i]);
			assocL2R.associate();

			FastQueue<AssociatedIndex> found = assocL2R.getMatches();

//			removeUnassociated(leftLoc,featsLeft1.description[i],rightLoc,featsRight1.description[i],found);
			setMatches(matches.match2to3, found, leftLoc.size);
		}
//		long time2 = System.currentTimeMillis();
//		System.out.println("  desc "+(time1-time0)+" assoc "+(time2-time1));
	}
1096295
def run_tracking(self, label_image_1, label_image_2):
        
        self.scale = self.parameters_tracking["avgCellDiameter"] / 35.0

        detections_1 = self.derive_detections(label_image_1)

        detections_2 = self.derive_detections(label_image_2)

        # Calculate tracking based on cell features and position.
        traces = self.find_initials_traces(detections_1, detections_2)

        # Use neighbourhoods to improve tracking.
        for _ in range(int(self.parameters_tracking["iterations"])):
            traces = self.improve_traces(detections_1, detections_2, traces)

        # Filter traces.
        return [(trace.previous_cell.number, trace.current_cell.number) for trace in traces]
1122817
protected void updateTrackLocation(T image) {
		get_subwindow(image, templateNew);

		// calculate response of the classifier at all locations
		// matlab: k = dense_gauss_kernel(sigma, x, z);
		dense_gauss_kernel(sigma, templateNew, template,k);

		fft.forward(k,kf);

		// response = real(ifft2(alphaf .* fft2(k)));   %(Eq. 9)
		DiscreteFourierTransformOps.multiplyComplex(alphaf, kf, tmpFourier0);
		fft.inverse(tmpFourier0, response);

		// find the pixel with the largest response
		int N = response.width*response.height;
		int indexBest = -1;
		double valueBest = -1;
		for( int i = 0; i < N; i++ ) {
			double v = response.data[i];
			if( v > valueBest ) {
				valueBest = v;
				indexBest = i;
			}
		}

		int peakX = indexBest % response.width;
		int peakY = indexBest / response.width;

		// sub-pixel peak estimation
		subpixelPeak(peakX, peakY);

		// peak in region's coordinate system
		float deltaX = (peakX+offX) - templateNew.width/2;
		float deltaY = (peakY+offY) - templateNew.height/2;

		// convert peak location into image coordinate system
		regionTrack.x0 = regionTrack.x0 + deltaX*stepX;
		regionTrack.y0 = regionTrack.y0 + deltaY*stepY;

		updateRegionOut();
	}
370851
def spline_backwards_hankel(ht, htarg, opt):
    """
    # Ensure ht is all lowercase
    ht = ht.lower()

    # Only relevant for 'fht' and 'hqwe', not for 'quad'
    if ht in ['fht', 'qwe', 'hqwe']:

        # Get corresponding htarg
        if ht == 'fht':
            htarg = _check_targ(htarg, ['fhtfilt', 'pts_per_dec'])
        elif ht in ['qwe', 'hqwe']:
            htarg = _check_targ(htarg, ['rtol', 'atol', 'nquad', 'maxint',
                                'pts_per_dec', 'diff_quad', 'a', 'b', 'limit'])

        # If spline (qwe, fht) or lagged (fht)
        if opt == 'spline':

            # Issue warning
            mesg = ("\n    The use of `opt='spline'` is deprecated and will " +
                    "be removed\n    in v2.0.0; use the corresponding " +
                    "setting in `htarg`.")
            warnings.warn(mesg, DeprecationWarning)

            # Reset opt
            opt = None

            # Check pts_per_dec; set to old default values if not given
            if 'pts_per_dec' not in htarg:
                if ht == 'fht':
                    htarg['pts_per_dec'] = -1  # Lagged Convolution DLF

                elif ht in ['qwe', 'hqwe']:
                    htarg['pts_per_dec'] = 80  # Splined QWE; old default value

    return htarg, opt