Tophat sets PI in the RG to be the inner distance size, but the SAM spec
    states should be the insert size. This fixes the RG in the alignment
    file generated by Tophat header to match the spec
2270
def fix_insert_size(in_bam, config):
    
    fixed_file = os.path.splitext(in_bam)[0] + ".pi_fixed.bam"
    if file_exists(fixed_file):
        return fixed_file
    header_file = os.path.splitext(in_bam)[0] + ".header.sam"
    read_length = bam.estimate_read_length(in_bam)
    bam_handle= bam.open_samfile(in_bam)
    header = bam_handle.header.copy()
    rg_dict = header['RG'][0]
    if 'PI' not in rg_dict:
        return in_bam
    PI = int(rg_dict.get('PI'))
    PI = PI + 2*read_length
    rg_dict['PI'] = PI
    header['RG'][0] = rg_dict
    with pysam.Samfile(header_file, "wb", header=header) as out_handle:
        with bam.open_samfile(in_bam) as in_handle:
            for record in in_handle:
                out_handle.write(record)
    shutil.move(header_file, fixed_file)
    return fixed_file
567139
def googlenet_resize(im, targ, min_area_frac, min_aspect_ratio, max_aspect_ratio, flip_hw_p, interpolation=cv2.INTER_AREA):
    
    h,w,*_ = im.shape
    area = h*w
    for _ in range(10):
        targetArea = random.uniform(min_area_frac, 1.0) * area
        aspectR = random.uniform(min_aspect_ratio, max_aspect_ratio)
        ww = int(np.sqrt(targetArea * aspectR) + 0.5)
        hh = int(np.sqrt(targetArea / aspectR) + 0.5)
        if flip_hw_p:
            ww, hh = hh, ww
        if hh <= h and ww <= w:
            x1 = 0 if w == ww else random.randint(0, w - ww)
            y1 = 0 if h == hh else random.randint(0, h - hh)
            out = im[y1:y1 + hh, x1:x1 + ww]
            out = cv2.resize(out, (targ, targ), interpolation=interpolation)
            return out
    out = scale_min(im, targ, interpolation=interpolation)
    out = center_crop(out)
    return out
59672
private void setInUseAndPreviousInUse(int p, int s) {
    setSizeAndPreviousInUseOfInUseChunk(p, s);
    head(p + s, head(p + s) | PINUSE_BIT);
  }
1330086
private void initTimeBgRect(float topPaddingSize) {
        float mHourLeft;
        float mMinuteLeft;
        float mSecondLeft;
        boolean isInitHasBackgroundTextBaseY = false;

        if (isShowDay) {
            // initialize day background and border rectF
            if (isShowTimeBgBorder) {
                mDayBgBorderRectF = new RectF(mLeftPaddingSize, topPaddingSize, mLeftPaddingSize + mDayTimeBgWidth + (mTimeBgBorderSize * 2), topPaddingSize + mTimeBgSize + (mTimeBgBorderSize * 2));
                mDayBgRectF = new RectF(mLeftPaddingSize + mTimeBgBorderSize, topPaddingSize + mTimeBgBorderSize, mLeftPaddingSize + mDayTimeBgWidth + mTimeBgBorderSize, topPaddingSize + mTimeBgSize + mTimeBgBorderSize);
            } else {
                mDayBgRectF = new RectF(mLeftPaddingSize, topPaddingSize, mLeftPaddingSize + mDayTimeBgWidth, topPaddingSize + mTimeBgSize);
            }
            // hour left point
            mHourLeft = mLeftPaddingSize + mDayTimeBgWidth + mSuffixDayTextWidth + mSuffixDayLeftMargin + mSuffixDayRightMargin + (mTimeBgBorderSize * 2);

            if (!isShowHour && !isShowMinute && !isShowSecond) {
                isInitHasBackgroundTextBaseY = true;
                initHasBackgroundTextBaseY(mDayBgRectF);
            }
        } else {
            // hour left point
            mHourLeft = mLeftPaddingSize;
        }

        if (isShowHour) {
            // initialize hour background border rectF
            if (isShowTimeBgBorder) {
                mHourBgBorderRectF = new RectF(mHourLeft, topPaddingSize, mHourLeft + mTimeBgSize + (mTimeBgBorderSize * 2), topPaddingSize + mTimeBgSize + (mTimeBgBorderSize * 2));
                mHourBgRectF = new RectF(mHourLeft + mTimeBgBorderSize, topPaddingSize + mTimeBgBorderSize, mHourLeft + mTimeBgSize + mTimeBgBorderSize, topPaddingSize + mTimeBgSize + mTimeBgBorderSize);
            } else {
                mHourBgRectF = new RectF(mHourLeft, topPaddingSize, mHourLeft + mTimeBgSize, topPaddingSize + mTimeBgSize);
            }
            // minute left point
            mMinuteLeft = mHourLeft + mTimeBgSize + mSuffixHourTextWidth + mSuffixHourLeftMargin + mSuffixHourRightMargin + (mTimeBgBorderSize * 2);

            if (!isInitHasBackgroundTextBaseY) {
                isInitHasBackgroundTextBaseY = true;
                initHasBackgroundTextBaseY(mHourBgRectF);
            }
        } else {
            // minute left point
            mMinuteLeft = mHourLeft;
        }

        if (isShowMinute) {
            // initialize minute background border rectF
            if (isShowTimeBgBorder) {
                mMinuteBgBorderRectF = new RectF(mMinuteLeft, topPaddingSize, mMinuteLeft + mTimeBgSize + (mTimeBgBorderSize * 2), topPaddingSize + mTimeBgSize + (mTimeBgBorderSize * 2));
                mMinuteBgRectF = new RectF(mMinuteLeft + mTimeBgBorderSize, topPaddingSize + mTimeBgBorderSize, mMinuteLeft + mTimeBgSize + mTimeBgBorderSize, topPaddingSize + mTimeBgSize + mTimeBgBorderSize);
            } else {
                mMinuteBgRectF = new RectF(mMinuteLeft, topPaddingSize, mMinuteLeft + mTimeBgSize, topPaddingSize + mTimeBgSize);
            }
            // second left point
            mSecondLeft = mMinuteLeft + mTimeBgSize + mSuffixMinuteTextWidth + mSuffixMinuteLeftMargin + mSuffixMinuteRightMargin + (mTimeBgBorderSize * 2);

            if (!isInitHasBackgroundTextBaseY) {
                isInitHasBackgroundTextBaseY = true;
                initHasBackgroundTextBaseY(mMinuteBgRectF);
            }
        } else {
            // second left point
            mSecondLeft = mMinuteLeft;
        }

        if (isShowSecond) {
            // initialize second background border rectF
            if (isShowTimeBgBorder) {
                mSecondBgBorderRectF = new RectF(mSecondLeft, topPaddingSize, mSecondLeft + mTimeBgSize + (mTimeBgBorderSize * 2), topPaddingSize + mTimeBgSize + (mTimeBgBorderSize * 2));
                mSecondBgRectF = new RectF(mSecondLeft + mTimeBgBorderSize, topPaddingSize + mTimeBgBorderSize, mSecondLeft + mTimeBgSize + mTimeBgBorderSize, topPaddingSize + mTimeBgSize + mTimeBgBorderSize);
            } else {
                mSecondBgRectF = new RectF(mSecondLeft, topPaddingSize, mSecondLeft + mTimeBgSize, topPaddingSize + mTimeBgSize);
            }

            if (isShowMillisecond) {
                // millisecond left point
                float mMillisecondLeft = mSecondLeft + mTimeBgSize + mSuffixSecondTextWidth + mSuffixSecondLeftMargin + mSuffixSecondRightMargin + (mTimeBgBorderSize * 2);

                // initialize millisecond background border rectF
                if (isShowTimeBgBorder) {
                    mMillisecondBgBorderRectF = new RectF(mMillisecondLeft, topPaddingSize, mMillisecondLeft + mTimeBgSize + (mTimeBgBorderSize * 2), topPaddingSize + mTimeBgSize + (mTimeBgBorderSize * 2));
                    mMillisecondBgRectF = new RectF(mMillisecondLeft + mTimeBgBorderSize, topPaddingSize + mTimeBgBorderSize, mMillisecondLeft + mTimeBgSize + mTimeBgBorderSize, topPaddingSize + mTimeBgSize + mTimeBgBorderSize);
                } else {
                    mMillisecondBgRectF = new RectF(mMillisecondLeft, topPaddingSize, mMillisecondLeft + mTimeBgSize, topPaddingSize + mTimeBgSize);
                }
            }

            if (!isInitHasBackgroundTextBaseY) {
                initHasBackgroundTextBaseY(mSecondBgRectF);
            }
        }
    }
204181
def _fix_offset(self, state, size, arch=None):
        
        if state is not None:
            arch = state.arch

        if arch is None:
            raise ValueError('Either "state" or "arch" must be specified.')

        offset = arch.registers[self.reg_name][0]
        if size in self.alt_offsets:
            return offset + self.alt_offsets[size]
        elif size < self.size and arch.register_endness == 'Iend_BE':
            return offset + (self.size - size)
        return offset
959244
def CircularMaskGrid(gridSize, diameter=None):
    
    if diameter is None: diameter=gridSize
    return np.less_equal(RadiusGrid(gridSize),diameter/2.0)
555665
function (pixmapWidth, pixmapHeight) {
                // Find out if the mask extends beyond the visible pixels
                var paddingWanted;
                ["top", "left", "right", "bottom"].forEach(function (key) {
                    if (paddedInputBounds[key] !== visibleInputBounds[key]) {
                        paddingWanted = true;
                        return false;
                    }
                });

                // When Photoshop produces inaccurate results, the padding is adjusted to compensate
                // When no padding is requested, this may be unwanted, so return a padding of 0px
                if (!paddingWanted) {
                    return { left: 0, top: 0, right: 0, bottom: 0 };
                }

                var // How much padding is necessary in both dimensions
                    missingWidth  = finalOutputWidth  - pixmapWidth,
                    missingHeight = finalOutputHeight - pixmapHeight,
                    // How of the original padding was on which side (default 0)
                    leftRatio     = paddingInputWidth === 0 ? 0 :
                        ((visibleInputBounds.left - paddedInputBounds.left) / paddingInputWidth),
                    topRatio      = paddingInputHeight === 0 ? 0 :
                        ((visibleInputBounds.top  - paddedInputBounds.top)  / paddingInputHeight),
                    // Concrete padding size on one side so the other side can use the rest
                    leftPadding   = Math.round(leftRatio * missingWidth),
                    topPadding    = Math.round(topRatio  * missingHeight);

                // Padding: how many transparent pixels to add on which side
                return {
                    left:   Math.max(0, leftPadding),
                    top:    Math.max(0, topPadding),
                    right:  Math.max(0, missingWidth  - leftPadding),
                    bottom: Math.max(0, missingHeight - topPadding)
                };
            }
429922
def fit_size(min_length: int = 0, max_length: int = None,
             message=None) -> Filter_T:
    

    def validate(value):
        length = len(value) if value is not None else 0
        if length < min_length or \
                (max_length is not None and length > max_length):
            _raise_failure(message)
        return value

    return validate
107577
def unit_tophat_ee (x):
    
    x = np.asarray (x)
    x1 = np.atleast_1d (x)
    r = ((0 < x1) & (x1 < 1)).astype (x.dtype)
    if x.ndim == 0:
        return np.asscalar (r)
    return r
533115
def _set_topMargin(self, value):
        
        bounds = self.bounds
        if bounds is None:
            self.height = value
        else:
            xMin, yMin, xMax, yMax = bounds
            self.height = yMax + value
240813
def _fix_sitk_bug(path, metadata):
        
        ds = dicom.read_file(path)
        try:
            metadata["voxelsize_mm"][0] = ds.SpacingBetweenSlices
        except Exception as e:
            logger.warning("Read dicom 'SpacingBetweenSlices' failed: ", e)
        return metadata
58724
def _match_spec(self, nugget):
        
        # TODO: handle positional modifiers and other similar format string tricks.
        all_spec = self._all_spec

        # iterate through nugget throwing away anything which is an int
        # TODO store this in a size variable

        original_nugget = nugget
        length_str = [ ]
        length_spec = None
        length_spec_str_len = 0

        if nugget.startswith(b".*"):
            # ".*": precision is specified as an argument
            nugget = nugget[2:]
            length_spec = b".*"
            length_spec_str_len = 2

        for j, c in enumerate(nugget):
            if c in ascii_digits:
                length_str.append(c)
            else:
                nugget = nugget[j:]
                if length_spec is None:
                    length_spec = None if len(length_str) == 0 else int(bytes(length_str))
                break

        # we need the length of the format's length specifier to extract the format and nothing else
        if length_spec_str_len == 0 and length_str:
            length_spec_str_len = len(length_str)
        # is it an actual format?
        for spec in all_spec:
            if nugget.startswith(spec):
                # this is gross coz sim_type is gross..
                nugget = nugget[:len(spec)]
                original_nugget = original_nugget[:(length_spec_str_len + len(spec))]
                nugtype = all_spec[nugget]
                try:
                    typeobj = sim_type.parse_type(nugtype).with_arch(self.state.arch)
                except:
                    raise SimProcedureError("format specifier uses unknown type '%s'" % repr(nugtype))
                return FormatSpecifier(original_nugget, length_spec, typeobj.size // 8, typeobj.signed)

        return None
497382
func Pi(z *decimal.Big) *decimal.Big {
	return pi(z, decimal.Context{Precision: precision(z)})
}
1321970
def _dsinit(
       whichconst,
       cosim,  emsq,   argpo,   s1,     s2,
       s3,     s4,     s5,      sinim,  ss1,
       ss2,    ss3,    ss4,     ss5,    sz1,
       sz3,    sz11,   sz13,    sz21,   sz23,
       sz31,   sz33,   t,       tc,     gsto,
       mo,     mdot,   no,      nodeo,  nodedot,
       xpidot, z1,     z3,      z11,    z13,
       z21,    z23,    z31,     z33,    ecco,
       eccsq,  em,    argpm,  inclm, mm,
       nm,    nodem,
       irez,
       atime, d2201, d2211,  d3210, d3222,
       d4410, d4422, d5220,  d5232, d5421,
       d5433, dedt,  didt,   dmdt,
       dnodt, domdt, del1,   del2,  del3,
       xfact, xlamo, xli,    xni,
     ):

     q22    = 1.7891679e-6;
     q31    = 2.1460748e-6;
     q33    = 2.2123015e-7;
     root22 = 1.7891679e-6;
     root44 = 7.3636953e-9;
     root54 = 2.1765803e-9;
     rptim  = 4.37526908801129966e-3; # equates to 7.29211514668855e-5 rad/sec
     root32 = 3.7393792e-7;
     root52 = 1.1428639e-7;
     x2o3   = 2.0 / 3.0;
     znl    = 1.5835218e-4;
     zns    = 1.19459e-5;

     #  sgp4fix identify constants and allow alternate values
     xke = whichconst.xke

     #  -------------------- deep space initialization ------------
     irez = 0;
     if 0.0034906585 < nm < 0.0052359877:
         irez = 1;
     if 8.26e-3 <= nm <= 9.24e-3 and em >= 0.5:
         irez = 2;

     #  ------------------------ do solar terms -------------------
     ses  =  ss1 * zns * ss5;
     sis  =  ss2 * zns * (sz11 + sz13);
     sls  = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
     sghs =  ss4 * zns * (sz31 + sz33 - 6.0);
     shs  = -zns * ss2 * (sz21 + sz23);
     #  sgp4fix for 180 deg incl
     if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2:
       shs = 0.0;
     if sinim != 0.0:
       shs = shs / sinim;
     sgs  = sghs - cosim * shs;

     #  ------------------------- do lunar terms ------------------
     dedt = ses + s1 * znl * s5;
     didt = sis + s2 * znl * (z11 + z13);
     dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
     sghl = s4 * znl * (z31 + z33 - 6.0);
     shll = -znl * s2 * (z21 + z23);
     #  sgp4fix for 180 deg incl
     if inclm < 5.2359877e-2 or inclm > pi - 5.2359877e-2:
         shll = 0.0;
     domdt = sgs + sghl;
     dnodt = shs;
     if sinim != 0.0:

         domdt = domdt - cosim / sinim * shll;
         dnodt = dnodt + shll / sinim;


     #  ----------- calculate deep space resonance effects --------
     dndt   = 0.0;
     theta  = (gsto + tc * rptim) % twopi
     em     = em + dedt * t;
     inclm  = inclm + didt * t;
     argpm  = argpm + domdt * t;
     nodem  = nodem + dnodt * t;
     mm     = mm + dmdt * t;
     

     #  -------------- initialize the resonance terms -------------
     if irez != 0:

         aonv = pow(nm / xke, x2o3);

         #  ---------- geopotential resonance for 12 hour orbits ------
         if irez == 2:

             cosisq = cosim * cosim;
             emo    = em;
             em     = ecco;
             emsqo  = emsq;
             emsq   = eccsq;
             eoc    = em * emsq;
             g201   = -0.306 - (em - 0.64) * 0.440;

             if em <= 0.65:

                 g211 =    3.616  -  13.2470 * em +  16.2900 * emsq;
                 g310 =  -19.302  + 117.3900 * em - 228.4190 * emsq +  156.5910 * eoc;
                 g322 =  -18.9068 + 109.7927 * em - 214.6334 * emsq +  146.5816 * eoc;
                 g410 =  -41.122  + 242.6940 * em - 471.0940 * emsq +  313.9530 * eoc;
                 g422 = -146.407  + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
                 g520 = -532.114  + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;

             else:

                 g211 =   -72.099 +   331.819 * em -   508.738 * emsq +   266.724 * eoc;
                 g310 =  -346.844 +  1582.851 * em -  2415.925 * emsq +  1246.113 * eoc;
                 g322 =  -342.585 +  1554.908 * em -  2366.899 * emsq +  1215.972 * eoc;
                 g410 = -1052.797 +  4758.686 * em -  7193.992 * emsq +  3651.957 * eoc;
                 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
                 if em > 0.715:
                     g520 =-5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
                 else:
                     g520 = 1464.74 -  4664.75 * em +  3763.64 * emsq;

             if em < 0.7:

                 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21  * eoc;
                 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
                 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4  * eoc;

             else:

                 g533 =-37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
                 g521 =-51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
                 g532 =-40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;

             sini2=  sinim * sinim;
             f220 =  0.75 * (1.0 + 2.0 * cosim+cosisq);
             f221 =  1.5 * sini2;
             f321 =  1.875 * sinim  *  (1.0 - 2.0 * cosim - 3.0 * cosisq);
             f322 = -1.875 * sinim  *  (1.0 + 2.0 * cosim - 3.0 * cosisq);
             f441 = 35.0 * sini2 * f220;
             f442 = 39.3750 * sini2 * sini2;
             f522 =  9.84375 * sinim * (sini2 * (1.0 - 2.0 * cosim- 5.0 * cosisq) +
                     0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq) );
             f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim +
                    10.0 * cosisq) + 6.56250012 * (1.0+2.0 * cosim - 3.0 * cosisq));
             f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim+cosisq *
                    (-12.0 + 8.0 * cosim + 10.0 * cosisq));
             f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim+cosisq *
                    (12.0 + 8.0 * cosim - 10.0 * cosisq));
             xno2  =  nm * nm;
             ainv2 =  aonv * aonv;
             temp1 =  3.0 * xno2 * ainv2;
             temp  =  temp1 * root22;
             d2201 =  temp * f220 * g201;
             d2211 =  temp * f221 * g211;
             temp1 =  temp1 * aonv;
             temp  =  temp1 * root32;
             d3210 =  temp * f321 * g310;
             d3222 =  temp * f322 * g322;
             temp1 =  temp1 * aonv;
             temp  =  2.0 * temp1 * root44;
             d4410 =  temp * f441 * g410;
             d4422 =  temp * f442 * g422;
             temp1 =  temp1 * aonv;
             temp  =  temp1 * root52;
             d5220 =  temp * f522 * g520;
             d5232 =  temp * f523 * g532;
             temp  =  2.0 * temp1 * root54;
             d5421 =  temp * f542 * g521;
             d5433 =  temp * f543 * g533;
             xlamo =  (mo + nodeo + nodeo-theta - theta) % twopi
             xfact =  mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - no;
             em    = emo;
             emsq  = emsqo;

         #  ---------------- synchronous resonance terms --------------
         if irez == 1:

             g200  = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
             g310  = 1.0 + 2.0 * emsq;
             g300  = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
             f220  = 0.75 * (1.0 + cosim) * (1.0 + cosim);
             f311  = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
             f330  = 1.0 + cosim;
             f330  = 1.875 * f330 * f330 * f330;
             del1  = 3.0 * nm * nm * aonv * aonv;
             del2  = 2.0 * del1 * f220 * g200 * q22;
             del3  = 3.0 * del1 * f330 * g300 * q33 * aonv;
             del1  = del1 * f311 * g310 * q31 * aonv;
             xlamo = (mo + nodeo + argpo - theta) % twopi
             xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no;

         #  ------------ for sgp4, initialize the integrator ----------
         xli   = xlamo;
         xni   = no;
         atime = 0.0;
         nm    = no + dndt;

     return (
       em,    argpm,  inclm, mm,
       nm,    nodem,
       irez, atime,
       d2201, d2211,  d3210, d3222,
       d4410, d4422, d5220,  d5232,
       d5421, d5433, dedt,  didt,
       dmdt,  dndt, dnodt, domdt,
       del1,   del2,  del3, xfact,
       xlamo, xli,    xni,
       )
708082
def overwrite_fits(hdulist, filename):
    

    assert isinstance(hdulist, (fits.HDUList, fits.PrimaryHDU))
    temp_name = None
    flag_delete_temp = False
    if os.path.isfile(filename):
        # PyFITS does not overwrite file
        temp_name = a99.rename_to_temp(filename)
    try:
        hdulist.writeto(filename, output_verify='warn')
        flag_delete_temp = temp_name is not None
    except:
        # Writing failed, reverts renaming
        os.rename(temp_name, filename)
        raise

    if flag_delete_temp:
        os.unlink(temp_name)
1184161
def _make_scaled_srcmap(self):
        """"""

        self.logger.info('Computing scaled source map.')

        bexp0 = fits.open(self.files['bexpmap_roi'])
        bexp1 = fits.open(self.config['gtlike']['bexpmap'])
        srcmap = fits.open(self.config['gtlike']['srcmap'])

        if bexp0[0].data.shape != bexp1[0].data.shape:
            raise Exception('Wrong shape for input exposure map file.')

        bexp_ratio = bexp0[0].data / bexp1[0].data

        self.logger.info(
            'Min/Med/Max exposure correction: %f %f %f' % (np.min(bexp_ratio),
                                                           np.median(
                                                               bexp_ratio),
                                                           np.max(bexp_ratio)))

        for hdu in srcmap[1:]:

            if hdu.name == 'GTI':
                continue
            if hdu.name == 'EBOUNDS':
                continue
            hdu.data *= bexp_ratio

        srcmap.writeto(self.files['srcmap'], overwrite=True)
18283
def PC_PI_calc(P, TOP, POP):
    
    try:
        result = 0
        for i in P.keys():
            result += ((P[i] + TOP[i]) / (2 * POP[i]))**2
        return result
    except Exception:
        return "None"