find best magnetometer offset fit to a log file
487320
def magfit(logfile):
    ''''''

    print("Processing log %s" % filename)
    mlog = mavutil.mavlink_connection(filename, notimestamps=args.notimestamps)

    data = []

    last_t = 0
    offsets = Vector3(0,0,0)
    motor_ofs = Vector3(0,0,0)
    motor = 0.0

    # now gather all the data
    while True:
        m = mlog.recv_match(condition=args.condition)
        if m is None:
            break
        if m.get_type() == "PARAM_VALUE" and m.param_id == "RC3_MIN":
            rc3_min = float(m.param_value)
        if m.get_type() == "SENSOR_OFFSETS":
            # update current offsets
            offsets = Vector3(m.mag_ofs_x, m.mag_ofs_y, m.mag_ofs_z)
        if m.get_type() == "SERVO_OUTPUT_RAW":
            motor_pwm = m.servo1_raw + m.servo2_raw + m.servo3_raw + m.servo4_raw
            motor_pwm *= 0.25
            rc3_min = mlog.param('RC3_MIN', 1100)
            rc3_max = mlog.param('RC3_MAX', 1900)
            motor = (motor_pwm - rc3_min) / (rc3_max - rc3_min)
            if motor > 1.0:
                motor = 1.0
            if motor < 0.0:
                motor = 0.0
        if m.get_type() == "RAW_IMU":
            mag = Vector3(m.xmag, m.ymag, m.zmag)
            # add data point after subtracting the current offsets
            data.append((mag - offsets + noise(), motor))

    print("Extracted %u data points" % len(data))
    print("Current offsets: %s" % offsets)

    data = select_data(data)

    # do an initial fit with all data
    (offsets, motor_ofs, field_strength) = fit_data(data)

    for count in range(3):
        # sort the data by the radius
        data.sort(lambda a,b : radius_cmp(a,b,offsets,motor_ofs))

        print("Fit %u    : %s  %s field_strength=%6.1f to %6.1f" % (
            count, offsets, motor_ofs,
            radius(data[0], offsets, motor_ofs), radius(data[-1], offsets, motor_ofs)))

        # discard outliers, keep the middle 3/4
        data = data[len(data)/8:-len(data)/8]

        # fit again
        (offsets, motor_ofs, field_strength) = fit_data(data)

    print("Final    : %s  %s field_strength=%6.1f to %6.1f" % (
        offsets, motor_ofs,
        radius(data[0], offsets, motor_ofs), radius(data[-1], offsets, motor_ofs)))
    print("mavgraph.py '%s' 'mag_field(RAW_IMU)' 'mag_field_motors(RAW_IMU,SENSOR_OFFSETS,(%f,%f,%f),SERVO_OUTPUT_RAW,(%f,%f,%f))'" % (
        filename,
        offsets.x,offsets.y,offsets.z,
        motor_ofs.x, motor_ofs.y, motor_ofs.z))
154003
function positionToMercator(position) {
    var lng = position[0];
    var lat = Math.max(Math.min(position[1], 89.99999), -89.99999);
    return [degToRad(lng) * EarthRadius, EarthRadius/2.0 * Math.log( (1.0 + Math.sin(degToRad(lat))) / (1.0 - Math.sin(degToRad(lat))) )];
  }
297607
def lognormcdf(x, mu, tau):
    """"""
    x = np.atleast_1d(x)
    return np.array(
        [0.5 * (1 - flib.derf(-(np.sqrt(tau / 2)) * (np.log(y) - mu))) for y in x])
1063341
def mjd2gmst(mjd):
    """"""
    
    D2PI= 6.2831853071795864769252867665590057683943387987502
    DS2R= 7.2722052166430399038487115353692196393452995355905e-5
    
    tu = (mjd-51544.5)/36525.0;
    
    st = fmod(mjd,1.0) * D2PI + ( 24110.54841 + ( 8640184.812866 + ( 0.093104 - 6.2e-6 * tu ) * tu ) * tu ) * DS2R 
    
    w = fmod(st,D2PI)
    if w >= 0.0 :
        return w
    else:
        return w+D2PI
1361160
def _get_cf_grid_mapping_var(self):
        """"""
        gmaps = ['fixedgrid_projection', 'goes_imager_projection',
                 'lambert_projection', 'polar_projection',
                 'mercator_projection']
        if 'grid_mapping' in self.filename_info:
            gmaps = [self.filename_info.get('grid_mapping')] + gmaps
        for grid_mapping in gmaps:
            if grid_mapping in self.nc:
                return self.nc[grid_mapping]
        raise KeyError("Can't find grid mapping variable in SCMI file")
269346
def _offset_setup(self,sigangle,leading,deltaAngleTrack):
        """"""
        #From the progenitor orbit, determine the sigmas in J and angle
        self._sigjr= (self._progenitor.rap()-self._progenitor.rperi())/numpy.pi*self._sigv
        self._siglz= self._progenitor.rperi()*self._sigv
        self._sigjz= 2.*self._progenitor.zmax()/numpy.pi*self._sigv
        #Estimate the frequency covariance matrix from a diagonal J matrix x dOdJ
        self._sigjmatrix= numpy.diag([self._sigjr**2.,
                                      self._siglz**2.,
                                      self._sigjz**2.])
        self._sigomatrix= numpy.dot(self._dOdJp,
                                    numpy.dot(self._sigjmatrix,self._dOdJp.T))
        #Estimate angle spread as the ratio of the largest to the middle eigenvalue
        self._sigomatrixEig= numpy.linalg.eig(self._sigomatrix)
        self._sigomatrixEigsortIndx= numpy.argsort(self._sigomatrixEig[0])
        self._sortedSigOEig= sorted(self._sigomatrixEig[0])
        if sigangle is None:
            self._sigangle= self._sigv*1.8
        else:
            self._sigangle= sigangle
        self._sigangle2= self._sigangle**2.
        self._lnsigangle= numpy.log(self._sigangle)
        #Estimate the frequency mean as lying along the direction of the largest eigenvalue
        self._dsigomeanProgDirection= self._sigomatrixEig[1][:,numpy.argmax(self._sigomatrixEig[0])]
        self._progenitor_Omega_along_dOmega= \
            numpy.dot(self._progenitor_Omega,self._dsigomeanProgDirection)
        #Make sure we are modeling the correct part of the stream
        self._leading= leading
        self._sigMeanSign= 1.
        if self._leading and self._progenitor_Omega_along_dOmega < 0.:
            self._sigMeanSign= -1.
        elif not self._leading and self._progenitor_Omega_along_dOmega > 0.:
            self._sigMeanSign= -1.
        self._progenitor_Omega_along_dOmega*= self._sigMeanSign 
        self._sigomean= self._progenitor_Omega\
            +self._sigMeanOffset*self._sigMeanSign\
            *numpy.sqrt(numpy.amax(self._sigomatrixEig[0]))\
            *self._dsigomeanProgDirection
#numpy.dot(self._dOdJp,
#                          numpy.array([self._sigjr,self._siglz,self._sigjz]))
        self._dsigomeanProg= self._sigomean-self._progenitor_Omega
        self._meandO= self._sigMeanOffset\
            *numpy.sqrt(numpy.amax(self._sigomatrixEig[0]))
        #Store cholesky of sigomatrix for fast evaluation
        self._sigomatrixNorm=\
            numpy.sqrt(numpy.sum(self._sigomatrix**2.))
        self._sigomatrixinv, self._sigomatrixLogdet= \
            fast_cholesky_invert(self._sigomatrix/self._sigomatrixNorm,
                                 tiny=10.**-15.,logdet=True)
        self._sigomatrixinv/= self._sigomatrixNorm
        deltaAngleTrackLim = (self._sigMeanOffset+4.) * numpy.sqrt(
            self._sortedSigOEig[2]) * self._tdisrupt
        if (deltaAngleTrack is None):  
            deltaAngleTrack = deltaAngleTrackLim
        else:
            if (deltaAngleTrack > deltaAngleTrackLim):
                warnings.warn("WARNING: angle range large compared to plausible value.", galpyWarning)
        self._deltaAngleTrack= deltaAngleTrack
        return None
47937
def _lsm_fix_strip_offsets(self):
        
        if self.filehandle.size < 2**32:
            return

        pages = self.pages
        npages = len(pages)
        series = self.series[0]
        axes = series.axes

        # find positions
        positions = 1
        for i in 0, 1:
            if series.axes[i] in 'PM':
                positions *= series.shape[i]

        # make time axis first
        if positions > 1:
            ntimes = 0
            for i in 1, 2:
                if axes[i] == 'T':
                    ntimes = series.shape[i]
                    break
            if ntimes:
                div, mod = divmod(npages, 2*positions*ntimes)
                assert mod == 0
                shape = (positions, ntimes, div, 2)
                indices = numpy.arange(product(shape)).reshape(shape)
                indices = numpy.moveaxis(indices, 1, 0)
        else:
            indices = numpy.arange(npages).reshape(-1, 2)

        # images of reduced page might be stored first
        if pages[0]._offsetscounts[0][0] > pages[1]._offsetscounts[0][0]:
            indices = indices[..., ::-1]

        # unwrap offsets
        wrap = 0
        previousoffset = 0
        for i in indices.flat:
            page = pages[int(i)]
            dataoffsets = []
            for currentoffset in page._offsetscounts[0]:
                if currentoffset < previousoffset:
                    wrap += 2**32
                dataoffsets.append(currentoffset + wrap)
                previousoffset = currentoffset
            page._offsetscounts = dataoffsets, page._offsetscounts[1]
427524
def check_offset(self):
        """"""

        for d in range(self.dmps):
            if (self.y0[d] == self.goal[d]):
                self.goal[d] += 1e-4
459986
def _fit_position_tsmap(self, name, **kwargs):
        """"""

        prefix = kwargs.get('prefix', '')
        dtheta_max = kwargs.get('dtheta_max', 0.5)
        zmin = kwargs.get('zmin', -3.0)

        kw = {
            'map_size': 2.0 * dtheta_max,
            'write_fits':  kwargs.get('write_fits', False),
            'write_npy':  kwargs.get('write_npy', False),
            'use_pylike': kwargs.get('use_pylike', True),
            'max_kernel_radius': self.config['tsmap']['max_kernel_radius'],
            'loglevel': logging.DEBUG
        }

        src = self.roi.copy_source(name)

        if src['SpatialModel'] in ['RadialDisk', 'RadialGaussian']:
            kw['max_kernel_radius'] = max(kw['max_kernel_radius'],
                                          2.0 * src['SpatialWidth'])

        skydir = kwargs.get('skydir', src.skydir)
        tsmap = self.tsmap(utils.join_strings([prefix, name.lower().
                                               replace(' ', '_')]),
                           model=src.data,
                           map_skydir=skydir,
                           exclude=[name],
                           make_plots=False, **kw)

        # Find peaks with TS > 4
        peaks = find_peaks(tsmap['ts'], 4.0, 0.2)
        peak_best = None
        o = {}
        for p in sorted(peaks, key=lambda t: t['amp'], reverse=True):
            xy = p['ix'], p['iy']
            ts_value = tsmap['ts'].data[xy[1], xy[0]]
            posfit = fit_error_ellipse(tsmap['ts'], xy=xy, dpix=2,
                                       zmin=max(zmin, -ts_value * 0.5))
            offset = posfit['skydir'].separation(self.roi[name].skydir).deg
            if posfit['fit_success'] and posfit['fit_inbounds']:
                peak_best = p
                break

        if peak_best is None:
            ts_value = np.max(tsmap['ts'].data)
            posfit = fit_error_ellipse(tsmap['ts'], dpix=2,
                                       zmin=max(zmin, -ts_value * 0.5))

        o.update(posfit)
        pix = posfit['skydir'].to_pixel(self.geom.wcs)
        o['xpix'] = float(pix[0])
        o['ypix'] = float(pix[1])
        o['skydir'] = posfit['skydir'].transform_to('icrs')
        o['pos_offset'] = posfit['skydir'].separation(
            self.roi[name].skydir).deg
        o['loglike'] = 0.5 * posfit['zoffset']
        o['tsmap'] = tsmap['ts']

        return o
448991
def log_(t, alpha=1, beta=1):
    """"""
    return beta * np.log(alpha * np.maximum(0, t - 10 * beta) + 1)
1140279
def LogSoftmax(x, params, axis=-1, **kwargs):
  """"""
  del params, kwargs
  return x - backend.logsumexp(x, axis, keepdims=True)
1340548
def _margdist_loglr(self, mf_snr, opt_snr):
        
        mf_snr_marg = mf_snr/self._dist_array
        opt_snr_marg = opt_snr/self._dist_array**2
        return special.logsumexp(mf_snr_marg - 0.5*opt_snr_marg,
                                 b=self._deltad*self.dist_prior)
1322502
def parse_skewer_log(self, f):
        
        fh = f['f']
        regexes = {
            'fq1': "Input file:\s+(.+)",
            'fq2': "Paired file:\s+(.+)",
            'r_processed': "(\d+) read|reads pairs? processed",
            'r_short_filtered': "(\d+) \(\s*\d+.\d+%\) short read",
            'r_empty_filtered': "(\d+) \(\s*\d+.\d+%\) empty read",
            'r_avail': "(\d+) \(\s*\d+.\d+%\) read",
            'r_trimmed': "(\d+) \(\s*\d+.\d+%\) trimmed read",
            'r_untrimmed': "(\d+) \(\s*\d+.\d+%\) untrimmed read"
        }
        regex_hist = "\s?(\d+)\s+(\d+)\s+(\d+.\d+)%"

        data = dict()
        for k, v in regexes.items():
            data[k] = 0
        data['fq1'] = None
        data['fq2'] = None
        readlen_dist = OrderedDict()

        for l in fh:
            for k, r in regexes.items():
                match = re.search(r, l)
                if match:
                    data[k] = match.group(1).replace(',', '')

            match = re.search(regex_hist, l)
            if match:
                read_length = int(match.group(1))
                pct_at_rl = float(match.group(3))
                readlen_dist[read_length] = pct_at_rl

        if data['fq1'] is not None:
            s_name = self.clean_s_name(data['fq1'], f['root'])
            if s_name in self.skewer_readlen_dist:
                log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
            self.add_data_source(f, s_name)
            self.add_skewer_data(s_name, data, f)
            self.skewer_readlen_dist[s_name] = readlen_dist

        if data['fq2'] is not None:
            s_name = self.clean_s_name(data['fq1'], f['root'])
            if s_name in self.skewer_readlen_dist:
                log.debug("Duplicate sample name found in {}! Overwriting: {}".format(f['fn'], s_name))
            self.add_data_source(f, s_name)
            self.add_skewer_data(s_name, data, f)
            self.skewer_readlen_dist[s_name] = readlen_dist
1371801
def log_m_probs(self):
        """"""
        m = self.kernel.feature_log_prob_[self._match_class_pos()]
        return self._prob_inverse_transform(m)
157229
def _determine_rotated_logfile(self):
        
        rotated_filename = self._check_rotated_filename_candidates()
        if rotated_filename and exists(rotated_filename):
            if stat(rotated_filename).st_ino == self._offset_file_inode:
                return rotated_filename

            # if the inode hasn't changed, then the file shrank; this is expected with copytruncate,
            # otherwise print a warning
            if stat(self.filename).st_ino == self._offset_file_inode:
                if self.copytruncate:
                    return rotated_filename
                else:
                    sys.stderr.write(
                        "[pygtail] [WARN] file size of %s shrank, and copytruncate support is "
                        "disabled (expected at least %d bytes, was %d bytes).\n" %
                        (self.filename, self._offset, stat(self.filename).st_size))

        return None
263757
def _parse_input_loc(block, profile, **kwargs):
    """"""
    layer, wave_field = parse_fixed_width(2 * [(5, int)], block)

    return profile.location(
        motion.WaveField[wave_field],
        index=(layer - 1), )
175935
def block_offset_bounds(self, namespace):
        """"""
        cursor = self.cursor
        cursor.execute('SELECT MIN("offset"), MAX("offset") '
                       'FROM gauged_statistics WHERE namespace = %s',
                       (namespace,))
        return cursor.fetchone()