run AFP chaining allowing up to maxTra flexible regions.
Input is original coordinates.
775745
private static Group[] chainAfp(FatCatParameters params,AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException{

		// we don;t want to rotate input atoms, do we?
		Atom[] ca2clone = StructureTools.cloneAtomArray(ca2);

		List<AFP> afpSet = afpChain.getAfpSet();

		if (debug)
			System.out.println("entering chainAfp");
		int afpNum = afpSet.size();

		if ( afpNum < 1)
			return new Group[0];

		long bgtime = System.currentTimeMillis();
		if(debug)    {
			System.out.println(String.format("total AFP %d\n", afpNum));
		}

		//run AFP chaining

		AFPChainer.doChainAfp(params,afpChain ,ca1,ca2);

		int afpChainLen = afpChain.getAfpChainLen();

		if(afpChainLen < 1)     {

			afpChain.setShortAlign(true);
			return new Group[0];
		} //very short alignment

		long chaintime = System.currentTimeMillis();
		if(debug)    {

			System.out.println("Afp chaining: time " + (chaintime-bgtime));
		}

		// do post processing

		AFPPostProcessor.postProcess(params, afpChain,ca1,ca2);

		// Optimize the final alignment

		AFPOptimizer.optimizeAln(params, afpChain,ca1,ca2);

		AFPOptimizer.blockInfo( afpChain);

		AFPOptimizer.updateScore(params,afpChain);

		AFPAlignmentDisplay.getAlign(afpChain,ca1,ca2);

		Group[] twistedPDB = AFPTwister.twistPDB(afpChain, ca1, ca2clone);

		SigEva sig =  new SigEva();
		double probability = sig.calSigAll(params, afpChain);
		afpChain.setProbability(probability);
		double normAlignScore = sig.calNS(params,afpChain);
		afpChain.setNormAlignScore(normAlignScore);
		double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2,false);
		afpChain.setTMScore(tmScore);

		/*

		SIGEVA  sig;
		probability = sig.calSigAll(maxTra, sparse, pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);
		normAlignScore = sig.calNS(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);

		 */

		//if(maxTra == 0)       probability = sig.calSigRigid(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength);
		//else  probability = sig.calSigFlexi(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);

		if(debug)    {

			long nowtime = System.currentTimeMillis();
			long diff = nowtime - chaintime;
			System.out.println("Alignment optimization: time "+ diff);

			System.out.println("score:      " + afpChain.getAlignScore());
			System.out.println("opt length: " + afpChain.getOptLength());
			System.out.println("opt rmsd:   "+ afpChain.getTotalRmsdOpt());

		}
		return twistedPDB;

	}
726859
def _parse_run_control(block):
    """"""
    _, max_iterations, strain_ratio, _, _ = parse_fixed_width(
        2 * [(5, int)] + [(10, float)] + 2 * [(5, int)], block)

    return propagation.EquivalentLinearCalculation(
        strain_ratio, max_iterations, tolerance=10.)
1087237
def chain(args):
    
    valid_merge_op = ('sum', 'min', 'max', 'mean', 'collapse')

    p = OptionParser(chain.__doc__)
    p.add_option("--key", dest="attrib_key", default=None,
                help="Attribute to use as `key` for chaining operation")
    p.add_option("--chain_ftype", default="cDNA_match",
                help="GFF feature type to use for chaining operation")
    p.add_option("--parent_ftype", default=None,
                help="GFF feature type to use for the chained coordinates")
    p.add_option("--break", dest="break_chain", action="store_true",
                help="Break long chains which are non-contiguous")
    p.add_option("--transfer_attrib", dest="attrib_list",
                help="Attributes to transfer to parent feature; accepts comma" + \
                " separated list of attribute names [default: %default]")
    p.add_option("--transfer_score", dest="score_merge_op", choices=valid_merge_op,
                help="Transfer value stored in score field to parent feature." + \
                " Score is reported based on chosen operation")
    p.set_outfile()

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    gffile, = args
    attrib_key = opts.attrib_key
    attrib_list = opts.attrib_list
    score_merge_op = opts.score_merge_op
    break_chain = opts.break_chain

    chain_ftype = opts.chain_ftype
    parent_ftype = opts.parent_ftype if opts.parent_ftype else chain_ftype

    gffdict = {}
    fw = must_open(opts.outfile, "w")
    gff = Gff(gffile)
    if break_chain:
        ctr, prev_gid = dict(), None
    for g in gff:
        if g.type != chain_ftype:
            print(g, file=fw)
            continue

        id = g.accn
        gid = id
        if attrib_key:
            assert attrib_key in g.attributes.keys(), \
                "Attribute `{0}` not present in GFF3".format(attrib_key)
            gid = g.get_attr(attrib_key)
        curr_gid = gid
        if break_chain:
            if prev_gid != curr_gid:
                if curr_gid not in ctr:
                    ctr[curr_gid] = 0
                else:
                    ctr[curr_gid] += 1
            gid = "{0}:{1}".format(gid, ctr[curr_gid])
        gkey = (g.seqid, gid)
        if gkey not in gffdict:
            gffdict[gkey] = { 'seqid': g.seqid,
                            'source': g.source,
                            'strand': g.strand,
                            'type': parent_ftype,
                            'coords': [],
                            'children': [],
                            'score': [],
                            'attrs': DefaultOrderedDict(set)
                          }
            gffdict[gkey]['attrs']['ID'].add(gid)

        if attrib_list:
            for a in attrib_list.split(","):
                if a in g.attributes:
                    [gffdict[gkey]['attrs'][a].add(x) for x in g.attributes[a]]
                    del g.attributes[a]

        if break_chain:
            _attrib = "Alias" if attrib_list and ("Name" not in attrib_list) else "Name"
            gffdict[gkey]['attrs'][_attrib].add(curr_gid)

        gffdict[gkey]['coords'].append((g.start, g.end))
        if score_merge_op:
            if is_number(g.score):
                gffdict[gkey]['score'].append(float(g.score))
                g.score = "."

        g.attributes["Parent"] = [gid]
        g.attributes["ID"] = ["{0}-{1}".\
                format(gid, len(gffdict[gkey]['children']) + 1)]
        g.type = valid_gff_parent_child[g.type]
        g.update_attributes()
        gffdict[gkey]['children'].append(g)
        if break_chain:
            prev_gid = curr_gid

    for gkey, v in sorted(gffdict.items()):
        gseqid, key = gkey
        seqid = v['seqid']
        source = v['source']
        type = v['type']
        strand = v['strand']
        start, stop = range_minmax(gffdict[gkey]['coords'])

        score = "."
        if score_merge_op:
            v['score'].sort()
            if score_merge_op == "sum":
                score = sum(v['score'])
            elif score_merge_op == "min":
                score = min(v['score'])
            elif score_merge_op == "max":
                score = max(v['score'])
            elif score_merge_op == "mean":
                score = sum(v['score'], 0.0)/len(v['score'])
            elif score_merge_op == "collapse":
                score = ",".join((str(x) for x in v['score']))

        g = GffLine("\t".join(str(x) for x in [seqid, source, type, start, stop, \
            score, strand, ".", None]))
        g.attributes = v['attrs']
        g.update_attributes()

        print(g, file=fw)

        for child in gffdict[gkey]['children']:
            print(child, file=fw)

    fw.close()
697645
function constrainedMDSGradient(x, fxprime, distances, constraints) {
        var loss = 0, i;
        for (i = 0; i < fxprime.length; ++i) {
            fxprime[i] = 0;
        }

        for (i = 0; i < distances.length; ++i) {
            var xi = x[2 * i], yi = x[2 * i + 1];
            for (var j = i + 1; j < distances.length; ++j) {
                var xj = x[2 * j], yj = x[2 * j + 1],
                    dij = distances[i][j],
                    constraint = constraints[i][j];

                var squaredDistance = (xj - xi) * (xj - xi) + (yj - yi) * (yj - yi),
                    distance$$1 = Math.sqrt(squaredDistance),
                    delta = squaredDistance - dij * dij;

                if (((constraint > 0) && (distance$$1 <= dij)) ||
                    ((constraint < 0) && (distance$$1 >= dij))) {
                    continue;
                }

                loss += 2 * delta * delta;

                fxprime[2*i]     += 4 * delta * (xi - xj);
                fxprime[2*i + 1] += 4 * delta * (yi - yj);

                fxprime[2*j]     += 4 * delta * (xj - xi);
                fxprime[2*j + 1] += 4 * delta * (yj - yi);
            }
        }
        return loss;
    }
822430
def _cont_norm_running_quantile_regions_mp(wl, fluxes, ivars, q, delta_lambda,
                                           ranges, n_proc=2, verbose=False):
    
    print("contnorm.py: continuum norm using running quantile")
    print("Taking spectra in %s chunks" % len(ranges))
    # nstars = fluxes.shape[0]
    nchunks = len(ranges)
    norm_fluxes = np.zeros(fluxes.shape)
    norm_ivars = np.zeros(ivars.shape)
    for i in xrange(nchunks):
        chunk = ranges[i, :]
        start = chunk[0]
        stop = chunk[1]
        if verbose:
            print('@Bo Zhang: Going to normalize Chunk [%d/%d], pixel:[%d, %d] ...'
                  % (i+1, nchunks, start, stop))
        output = _cont_norm_running_quantile_mp(
            wl[start:stop], fluxes[:, start:stop],
            ivars[:, start:stop], q, delta_lambda,
            n_proc=n_proc, verbose=verbose)
        norm_fluxes[:, start:stop] = output[0]
        norm_ivars[:, start:stop] = output[1]
    return norm_fluxes, norm_ivars
1231564
def _array_io(self, action, array, frames):
        """"""
        if (array.ndim not in (1, 2) or
                array.ndim == 1 and self.channels != 1 or
                array.ndim == 2 and array.shape[1] != self.channels):
            raise ValueError("Invalid shape: {0!r}".format(array.shape))
        if not array.flags.c_contiguous:
            raise ValueError("Data must be C-contiguous")
        ctype = self._check_dtype(array.dtype.name)
        assert array.dtype.itemsize == _ffi.sizeof(ctype)
        cdata = _ffi.cast(ctype + '*', array.__array_interface__['data'][0])
        return self._cdata_io(action, cdata, ctype, frames)
861033
def chain(args):
    
    p = OptionParser(chain.__doc__)
    p.add_option("--dist", dest="dist",
            default=100, type="int",
            help="extent of flanking regions to search [default: %default]")

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    blastfile, = args
    dist = opts.dist
    assert dist > 0

    blast = BlastSlow(blastfile)
    logging.debug("A total of {} records imported".format(len(blast)))
    chained_hsps = chain_HSPs(blast, xdist=dist, ydist=dist)
    logging.debug("A total of {} records after chaining".format(len(chained_hsps)))

    for b in chained_hsps:
        print(b)
459855
def _offset_from_spaces(dom, ran):
    """"""
    affected = np.not_equal(dom.shape, ran.shape)
    diff_l = np.abs(ran.grid.min() - dom.grid.min())
    offset_float = diff_l / dom.cell_sides
    offset = np.around(offset_float).astype(int)
    for i in range(dom.ndim):
        if affected[i] and not np.isclose(offset[i], offset_float[i]):
            raise ValueError('in axis {}: range is shifted relative to domain '
                             'by a non-multiple {} of cell_sides'
                             ''.format(i, offset_float[i] - offset[i]))
    offset[~affected] = 0
    return tuple(offset)
104604
def offsets(self):
      """"""
      for path, tailedfile in self._tailedfiles.iteritems():
         yield path, tailedfile._offset
844523
private static Group[] rChainAfp(FatCatParameters params, AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException{
		params.setMaxTra(0);
		afpChain.setMaxTra(0);
		return chainAfp(params,afpChain,ca1,ca2);
	}
1331024
private static void permuteOptAln(AFPChain afpChain, int cp)
	{
		int ca2len = afpChain.getCa2Length();

		if( ca2len <= 0) {
			throw new IllegalArgumentException("No Ca2Length specified in "+afpChain);
		}

		// Allow negative cp points for convenience.
		if(cp == 0) {
			return;
		}
		if(cp <= -ca2len || cp >= ca2len) {
			// could just take cp%ca2len, but probably its a bug if abs(cp)>=ca2len
			throw new ArrayIndexOutOfBoundsException( String.format(
					"Permutation point %d must be between %d and %d for %s",
					cp, 1-ca2len,ca2len-1, afpChain.getName2() ) );
		}
		if(cp < 0) {
			cp = cp + ca2len;
		}

		// the unprocessed alignment
		int[][][] optAln = afpChain.getOptAln();
		int[] optLen = afpChain.getOptLen();

		// the processed alignment
		List<List<List<Integer>>> blocks = new ArrayList<List<List<Integer>>>(afpChain.getBlockNum()*2);

		//Update residue indices
		// newi = (oldi-cp) % N
		for(int block = 0; block < afpChain.getBlockNum(); block++) {
			if(optLen[block]<1)
				continue;

			// set up storage for the current block
			List<List<Integer>> currBlock = new ArrayList<List<Integer>>(2);
			currBlock.add( new ArrayList<Integer>());
			currBlock.add( new ArrayList<Integer>());
			blocks.add(currBlock);

			// pos = 0 case
			currBlock.get(0).add( optAln[block][0][0] );
			currBlock.get(1).add( (optAln[block][1][0]+cp ) % ca2len);

			for(int pos = 1; pos < optLen[block]; pos++) {
				//check if we need to start a new block
				//this happens when the new alignment crosses the protein terminus
				if( optAln[block][1][pos-1]+cp<ca2len &&
						optAln[block][1][pos]+cp >= ca2len) {
					currBlock = new ArrayList<List<Integer>>(2);
					currBlock.add( new ArrayList<Integer>());
					currBlock.add( new ArrayList<Integer>());
					blocks.add(currBlock);
				}
				currBlock.get(0).add( optAln[block][0][pos] );
				currBlock.get(1).add( (optAln[block][1][pos]+cp ) % ca2len);
			}
		}

		// save permuted blocks to afpChain
		assignOptAln(afpChain,blocks);
	}
697833
private void trackFeatures(Image input, RectangleRotate_F64 region) {
		pairs.reset();

		currentImage.process(input);
		for( int i = 0; i < currentImage.getNumLayers(); i++ ) {
			Image layer = currentImage.getLayer(i);
			gradient.process(layer,currentDerivX[i],currentDerivY[i]);
		}

		// convert to float to avoid excessive conversions from double to float
		float cx = (float)region.cx;
		float cy = (float)region.cy;

		float height = (float)(region.height);
		float width = (float)(region.width);

		float c = (float)Math.cos(region.theta);
		float s = (float)Math.sin(region.theta);

		float p = 1.0f/(config.numberOfSamples-1);
		for( int i = 0; i < config.numberOfSamples; i++ ) {
			float y = (p*i-0.5f)*height;
			for( int j = 0; j < config.numberOfSamples; j++ ) {
				float x = (p*j - 0.5f)*width;
				float xx = cx + x*c - y*s;
				float yy = cy + x*s + y*c;

				// track in the forward direction
				track.x = xx;
				track.y = yy;

				klt.setImage(previousImage,previousDerivX,previousDerivY);
				if( !klt.setDescription(track) ) {
					continue;
				}

				klt.setImage(currentImage,currentDerivX,currentDerivY);
				KltTrackFault fault = klt.track(track);
				if( fault != KltTrackFault.SUCCESS ) {
					continue;
				}

				float xc = track.x;
				float yc = track.y;

				// validate by tracking backwards
				if( !klt.setDescription(track) ) {
					continue;
				}
				klt.setImage(previousImage,previousDerivX,previousDerivY);
				fault = klt.track(track);
				if( fault != KltTrackFault.SUCCESS ) {
					continue;
				}

				float error = UtilPoint2D_F32.distanceSq(track.x, track.y, xx, yy);

				if( error > maximumErrorFB ) {
					continue;
				}

				// create a list of the observations
				AssociatedPair a = pairs.grow();

				a.p1.x = xx;
				a.p1.y = yy;
				a.p2.x = xc;
				a.p2.y = yc;
			}
		}
	}
541475
def affine(sx, rx, ry, sy, tx, ty)
      primitive 'affine ' + format('%g,%g,%g,%g,%g,%g', sx, rx, ry, sy, tx, ty)
    end
543482
function(){

            var tr
                ,intr
                ,aabb
                ,list = this.tracked
                ,i = list.length
                ;

            // for all tracked bodies
            while ( (--i) >= 0 ){

                tr = list[ i ];
                intr = tr.interval;
                aabb = tr.body.aabb();

                // copy the position (plus or minus) the aabb half-dimensions
                // into the min/max intervals
                intr.min.val.clone( aabb ).sub( aabb.hw, aabb.hh );
                intr.max.val.clone( aabb ).add( aabb.hw, aabb.hh );
            }
        }
668687
def _nanstd(array, axis=None, ddof=0):
    """"""

    if isinstance(axis, tuple):
        array = _move_tuple_axes_first(array, axis=axis)
        axis = 0
    return bottleneck.nanstd(array, axis=axis, ddof=ddof)
1200814
private void applyHorizontalLimit()
    {
        if (mover.getX() < limitLeft && limitLeft != Integer.MIN_VALUE)
        {
            mover.teleportX(limitLeft);
        }
        else if (mover.getX() > limitRight && limitRight != Integer.MAX_VALUE)
        {
            mover.teleportX(limitRight);
        }
    }
531030
def _calc_strain(self, loc_input, loc_layer, motion, *args):
        """"""
        strain_max = self._calc_strain_max(loc_input, loc_layer, motion, *args)
        return self.strain_ratio * strain_max