The first time delay weighted statistical moment of the
        instantaneous unit hydrograph.
180007
def moment1(self):
        """"""
        delays, response = self.delay_response_series
        return statstools.calc_mean_time(delays, response)
842144
Stream<WavefrontMetricLineData> writeMeter(Meter meter) {
        long wallTime = clock.wallTime();
        Stream.Builder<WavefrontMetricLineData> metrics = Stream.builder();

        stream(meter.measure().spliterator(), false)
            .forEach(measurement -> {
                Meter.Id id = meter.getId().withTag(measurement.getStatistic());
                addMetric(metrics, id, null, wallTime, measurement.getValue());
            });

        return metrics.build();
    }
710025
def clusterMetrics(timeValues, massValues, massScalingFactor=1):
    
    metrics = dict()
    metrics['meanTime'] = numpy.mean(timeValues)
    metrics['meanMass'] = numpy.mean(massValues)
    metrics['devTime'] = timeValues - metrics['meanTime']
    metrics['devMass'] = massValues - metrics['meanMass']
    #metrics['devMass'] = (1-metrics['meanMass']/massValues)
    metrics['spreadTime'] = numpy.max(timeValues) - numpy.min(timeValues)
    metrics['spreadMass'] = numpy.max(massValues) - numpy.min(massValues)
    #metrics['spreadMass'] = (1-numpy.min(massValues) / numpy.max(massValues))
    metrics['devEuc'] = numpy.sqrt(numpy.power(metrics['devTime'], 2) +
                                   numpy.power(metrics['devMass']*massScalingFactor, 2)
                                   )
    metrics['meanEuc'] = numpy.mean(metrics['devEuc'])
    metrics['devTime'] = metrics['devTime'].tolist()
    metrics['devMass'] = metrics['devMass'].tolist()
    metrics['devEuc'] = metrics['devEuc'].tolist()
    return metrics
935812
def _interpolate_stream_track_aA(self):
        """"""
        if hasattr(self,'_interpolatedObsTrackAA'):
            return None #Already did this
        #Calculate 1D meanOmega on a fine grid in angle and interpolate
        if not hasattr(self,'_interpolatedThetasTrack'):
            self._interpolate_stream_track()
        dmOs= numpy.array([self.meanOmega(da,oned=True,use_physical=False) 
                          for da in self._interpolatedThetasTrack])
        self._interpTrackAAdmeanOmegaOneD=\
            interpolate.InterpolatedUnivariateSpline(\
            self._interpolatedThetasTrack,dmOs,k=3)
        #Build the interpolated AA
        self._interpolatedObsTrackAA=\
            numpy.empty((len(self._interpolatedThetasTrack),6))
        for ii in range(len(self._interpolatedThetasTrack)):
            self._interpolatedObsTrackAA[ii,:3]=\
                self._progenitor_Omega+dmOs[ii]*self._dsigomeanProgDirection\
                *self._sigMeanSign   
            self._interpolatedObsTrackAA[ii,3:]=\
                self._progenitor_angle+self._interpolatedThetasTrack[ii]\
                *self._dsigomeanProgDirection*self._sigMeanSign
            self._interpolatedObsTrackAA[ii,3:]=\
                numpy.mod(self._interpolatedObsTrackAA[ii,3:],2.*numpy.pi)
        return None
350379
def K(self, X, X2=None):        
    # model :   -a d^2y/dx^2  + b dy/dt + c * y = U
    # kernel Kyy rbf spatiol temporal
    # vyt Y temporal variance  vyx Y spatiol variance   lyt Y temporal lengthscale   lyx Y spatiol lengthscale
    # kernel Kuu doper( doper(Kyy))
    # a   b    c    lyt   lyx    vyx*vyt
        """"""        
        X,slices = X[:,:-1],index_to_slices(X[:,-1])
        if X2 is None:
            X2,slices2 = X,slices
            K = np.zeros((X.shape[0], X.shape[0]))
        else:
            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
            K = np.zeros((X.shape[0], X2.shape[0]))


        tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
        xdist = (X[:,1][:,None] - X2[:,1][None,:])**2

        ttdist = (X[:,0][:,None] - X2[:,0][None,:])
        #rdist = [tdist,xdist]
        #dist = np.abs(X - X2.T)
        vyt = self.variance_Yt
        vyx = self.variance_Yx
        
        lyt=1/(2*self.lengthscale_Yt)
        lyx=1/(2*self.lengthscale_Yx)

        a = self.a ## -a is used in the model, negtive diffusion
        b = self.b
        c = self.c

        kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))

        k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )

        k2 = lambda xdist: ( 4*lyx**2 * (xdist)  - 2*lyx )

        k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )

        k4 = lambda ttdist: 2*lyt*(ttdist)

        for i, s1 in enumerate(slices):
            for j, s2 in enumerate(slices2):
                for ss1 in s1:
                    for ss2 in s2:
                        if i==0 and j==0:
                            K[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        elif i==0 and j==1:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) + b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] )   )
                        elif i==1 and j==1:
                            K[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )* vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                        else:
                            K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) - b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
                            #K[ss1,ss2]= np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) )   )
                            #K[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] )   )
        
        #stop
        return K
1285204
def _create_initial_point(self, Ts, E, entries, weights, Xs):
        
        import lifelines  # kinda hacky but lol

        def transform_ith_param(model, i):
            param = model._fitted_parameters_[i]
            if param <= 0:
                return param
            # technically this is suboptimal for log normal mu, but that's okay.
            return np.log(param)

        name = self.__class__.__name__.replace("AFT", "")
        uni_model = getattr(lifelines, name)()

        if self._censoring_type == CensoringType.RIGHT:
            uni_model.fit_right_censoring(Ts[0], event_observed=E, entry=entries, weights=weights)
        elif self._censoring_type == CensoringType.INTERVAL:
            uni_model.fit_interval_censoring(Ts[0], Ts[1], event_observed=E, entry=entries, weights=weights)
        elif self._censoring_type == CensoringType.LEFT:
            uni_model.fit_left_censoring(Ts[1], event_observed=E, entry=entries, weights=weights)

        # we may use this later in print_summary
        self._ll_null_ = uni_model._log_likelihood

        return np.concatenate(
            [
                # tack on as the intercept
                [0] * (_X.shape[1] - 1) + [transform_ith_param(uni_model, i)]
                for i, _X in enumerate(Xs)
            ]
        )
651930
def funcConvPar(aryDm,
                vecHrf,
                varNumVol):

    
    # In order to avoid an artefact at the end of the time series, we have to
    # concatenate an empty array to both the design matrix and the HRF model
    # before convolution.
    aryDm = np.concatenate((aryDm, np.zeros((aryDm.shape[0], 100))), axis=1)
    vecHrf = np.concatenate((vecHrf, np.zeros((100,))))

    aryDmConv = np.empty((aryDm.shape[0], varNumVol))
    for idx in range(0, aryDm.shape[0]):
        vecDm = aryDm[idx, :]
        # Convolve design matrix with HRF model:
        aryDmConv[idx, :] = np.convolve(vecDm, vecHrf,
                                        mode='full')[:varNumVol]
    return aryDmConv
980034
@Reference(authors = "D. Pelleg, A. Moore", //
      title = "X-means: Extending K-means with Efficient Estimation on the Number of Clusters", //
      booktitle = "Proc. 17th Int. Conf. on Machine Learning (ICML 2000)", //
      url = "http://www.pelleg.org/shared/hp/download/xmeans.ps", //
      bibkey = "DBLP:conf/icml/PellegM00")
  public static <V extends NumberVector> double logLikelihood(Relation<V> relation, Clustering<? extends MeanModel> clustering, NumberVectorDistanceFunction<? super V> distanceFunction) {
    List<? extends Cluster<? extends MeanModel>> clusters = clustering.getAllClusters();
    // number of clusters
    final int m = clusters.size();

    // number of objects in the clustering
    int n = 0;
    // cluster sizes
    int[] n_i = new int[m];
    // total variance
    double d = 0.;
    // variances
    double[] d_i = new double[m];

    // Iterate over clusters:
    Iterator<? extends Cluster<? extends MeanModel>> it = clusters.iterator();
    for(int i = 0; it.hasNext(); ++i) {
      Cluster<? extends MeanModel> cluster = it.next();
      n += n_i[i] = cluster.size();
      d += d_i[i] = varianceOfCluster(cluster, distanceFunction, relation);
    }

    // No remaining variance, if every point is on its own:
    if(n <= m) {
      return Double.NEGATIVE_INFINITY;
    }

    // Total variance (corrected for bias)
    final double logv = FastMath.log(d / (n - m));

    final int dim = RelationUtil.dimensionality(relation);
    // log likelihood of this clustering
    double logLikelihood = 0.;

    // Aggregate
    for(int i = 0; i < m; i++) {
      logLikelihood += n_i[i] * FastMath.log(n_i[i]) // Post. entropy Rn log Rn
          - n_i[i] * .5 * MathUtil.LOGTWOPI // Rn/2 log2pi
          - n_i[i] * dim * .5 * logv // Rn M/2 log sigma^2
          - (d_i[i] - m) * .5; // (Rn-K)/2
    }
    logLikelihood -= n * FastMath.log(n); // Prior entropy, sum_i Rn log R
    return logLikelihood;
  }
63534
def generate_state_data(means, weights):
    
    x_true = np.dot(means, weights)
    sample = np.random.poisson(x_true)
    return sample.astype(float)
72690
function(delta, numPoints, firstTime) {
	let k = 0;
	delta = firstTime ? floor(delta / damp) : delta >> 1;
	delta += floor(delta / numPoints);
	for (/* no initialization */; delta > baseMinusTMin * tMax >> 1; k += base) {
		delta = floor(delta / baseMinusTMin);
	}
	return floor(k + (baseMinusTMin + 1) * delta / (delta + skew));
}
504476
function TimeseriesHistogram( canvas, width, height, left, top ) {
		var graph, data, histogram, edges, axes, annotations, title, text, means;

		// [1] Instantiate a new graph generator and configure:
		graph = xfig.graph( canvas )
			.width( width )
			.height( height )
			.position({
				'left': left,
				'top': top
			})
			.xMin( 0 )
			.xMax( 1 );

		// Create the graph:
		graph.create( 'timeseries-histogram' );

		// Get data:
		d3.json( 'data/timeseries-histogram.data.json', function ( error, json ) {

			// [2] Instantiate a new data generator and configure:
			data = xfig.data( json )
				.accessors( 'x', function ( d ) {
					return d[ 0 ];
				})
				.accessors( 'y', function ( d ) {
					return d[ 1 ];
				});

			// Create edges to define our histogram bins:
			edges = xfig.vector.linspace( -0.025, 1.025, 0.05 );
			
			// Format the data and extract the data to histogram:
			data.format( [ 'x', 'y' ] )
				.extract( function ( d ) {
					return d[ 1 ];
				});

			// Calculate each dataset's mean value and sort the datasets based on their means:
			means = data.mean( function ( d ) {
				return d;
			});
			means = means.map( function ( d, i ) {
				return [ d, i ];
			});
			means.sort( function ( a, b ) {
				return a[0] < b[0] ? -1 : ( a[0] > b[0] ? 1 : 0 );
			});
			means = means.map( function ( d ) {
				return d[ 1 ];
			});

			// Calculate the histogram:
			data.reorder( means )
				.histc( function ( d ) { return d; }, edges );

			// Bind the data instance to the graph:
			graph.data( data )
				.yMin( 0 )
				.yMax( data.size()[0] )
				.zMin( 0 )
				.zMax( data.max( function ( d ) {
					return d[ 1 ];
				}))
				.zRange( ['#ffffff', '#000000'] );

			// [3] Instantiate a new histogram generator and configure:
			histogram = xfig.timeserieshistogram( graph );

			// Create the histogram:
			histogram.create();

			// [4] Instantiate a new axes generator and configure:
			axes = xfig.axes( graph )
				.xLabel( 'percent' )
				.yLabel( 'index' );

			// Create the axes:
			axes.create();

			// [5] Instantiate a new annotations generator and configure:
			annotations = xfig.annotations( graph );

			// Create the annotations element:
			annotations.create();

			// [5.1] Instantiate a new title instance and configure:
			title = annotations.title()
				.top( -30 )
				.left( 0 );

			// Add a (sub)title:
			title.create( 'Subtitle' );

		});

	}
515222
def _forward_step(self, Y, T_X, Var_X, Var_dX, rho_e, sigma2_e, weight):
        
        # We currently only implement diagonal form
        # of covariance matrix for Var_X, Var_dX and T_X, which means
        # each dimension of X is independent and their innovation noise
        # are also independent. Note that log_p_data takes this assumption.
        if Var_X.ndim == 1:
            inv_Var_X = np.diag(1 / Var_X)
            half_log_det_Var_X = np.sum(np.log(Var_X)) / 2.0
            Var_X = np.diag(Var_X)
            # the marginal variance of X
        else:
            half_log_det_Var_X = self._half_log_det(Var_X)
            inv_Var_X = np.linalg.inv(Var_X)
        if Var_dX.ndim == 1:
            inv_Var_dX = np.diag(1 / Var_dX)
            half_log_det_Var_dX = np.sum(np.log(Var_dX)) / 2.0
            Var_dX = np.diag(Var_dX)
            # the marginal variance of Delta X (the change of X from
            # previous time point)
        else:
            inv_Var_dX = np.linalg.inv(Var_dX)
            half_log_det_Var_dX = self._half_log_det(Var_dX)
        if T_X.ndim == 1:
            T_X = np.diag(T_X)
            # Transfer function of X: the expected mean of X at t+1
            # time point is T_x * X
        [n_T, n_V] = np.shape(Y)
        # numbers of time points and voxels
        mu = [None] * n_T
        # posterior mean of X, conditioned on all data up till the current
        # time point
        Gamma_inv = [None] * n_T
        # inverse of poterior Gamma.
        mu_Gamma_inv = [None] * n_T
        # mu * inv(Gamma)
        log_p_data = - np.log(np.pi * 2) * (n_T * n_V) / 2 \
            - half_log_det_Var_X - np.sum(np.log(sigma2_e)) * n_T / 2.0\
            + np.sum(np.log(1 - rho_e**2)) / 2.0 - half_log_det_Var_dX \
            * (n_T - 1)
        # This is the term to be incremented by c_n at each time step.
        # We first add all the fixed terms to it.

        # The following are a few fixed terms.
        Lambda_0 = np.dot(T_X, np.dot(inv_Var_dX, T_X.T)) \
            + np.dot(weight * rho_e**2 / sigma2_e, weight.T)
        H = np.dot(inv_Var_dX, T_X.T) + np.dot(weight * rho_e / sigma2_e,
                                               weight.T)
        Lambda_1 = inv_Var_dX + np.dot(weight / sigma2_e, weight.T)

        Gamma_inv[0] = inv_Var_X + np.dot(
            weight * (1 - rho_e**2) / sigma2_e, weight.T)
        # We might not need this and only use linalg.solve for related terms.
        mu_Gamma_inv[0] = np.dot(
            Y[0, :] * (1 - rho_e**2) / sigma2_e, weight.T)
        mu[0] = np.linalg.solve(Gamma_inv[0], mu_Gamma_inv[0])
        log_p_data -= 0.5 * np.sum(Y[0, :]**2 * (1 - rho_e**2) / sigma2_e)
        # This is the term added for the first time point.

        deltaY = Y[1:, :] - rho_e * Y[:-1, :]
        deltaY_sigma2inv_rho_weightT = np.dot(
            deltaY / sigma2_e * rho_e, weight.T)
        for t in np.arange(1, n_T):
            Gamma_tilde_inv = Lambda_0 + Gamma_inv[t - 1]
            tmp = np.linalg.solve(Gamma_tilde_inv, H.T)
            Gamma_inv[t] = Lambda_1 - np.dot(H, tmp)
            mu_Gamma_inv[t] = np.dot(deltaY[t - 1, :] / sigma2_e, weight.T) \
                + np.dot(mu_Gamma_inv[t - 1]
                         - deltaY_sigma2inv_rho_weightT[t - 1, :], tmp)
            mu[t] = np.linalg.solve(Gamma_inv[t], mu_Gamma_inv[t])
            tmp2 = mu_Gamma_inv[t - 1] - deltaY_sigma2inv_rho_weightT[t - 1, :]
            log_p_data += -self._half_log_det(Gamma_tilde_inv) \
                + np.dot(tmp2, np.linalg.solve(Gamma_tilde_inv, tmp2)) / 2.0
        log_p_data += -self._half_log_det(Gamma_inv[-1]) \
            + np.dot(mu_Gamma_inv[-1], mu[-1]) / 2.0 \
            - np.sum(deltaY**2 / sigma2_e) / 2.0
        return mu, mu_Gamma_inv, Gamma_inv, log_p_data, Lambda_0, \
            Lambda_1, H, deltaY, deltaY_sigma2inv_rho_weightT
125650
def moment2(self):
        """"""
        moment1 = self.moment1
        delays, response = self.delay_response_series
        return statstools.calc_mean_time_deviation(
            delays, response, moment1)
802254
def generate_return_periods(qout_file,
                            return_period_file,
                            num_cpus=multiprocessing.cpu_count(),
                            storm_duration_days=7,
                            method='weibull'):
    
    # get ERA Interim Data Analyzed
    with RAPIDDataset(qout_file) as qout_nc_file:
        print("Setting up Return Periods File ...")
        return_period_nc = Dataset(return_period_file, 'w')

        return_period_nc.createDimension('rivid', qout_nc_file.size_river_id)

        timeSeries_var = \
            return_period_nc.createVariable('rivid', 'i4', ('rivid',))
        timeSeries_var.long_name = (
            'unique identifier for each river reach')

        max_flow_var = \
            return_period_nc.createVariable('max_flow', 'f8', ('rivid',))
        max_flow_var.long_name = 'maximum streamflow'
        max_flow_var.units = 'm3/s'

        if method == 'weibull':

            return_period_20_var = \
                return_period_nc.createVariable('return_period_20',
                                                'f8', ('rivid',))
            return_period_20_var.long_name = '20 year return period flow'
            return_period_20_var.units = 'm3/s'

        if method == 'gumble':

            return_period_100_var = \
                return_period_nc.createVariable('return_period_100',
                                                'f8', ('rivid',))
            return_period_100_var.long_name = '100 year return period flow'
            return_period_100_var.units = 'm3/s'

            return_period_50_var = \
                return_period_nc.createVariable('return_period_50',
                                                'f8', ('rivid',))
            return_period_50_var.long_name = '50 year return period flow'
            return_period_50_var.units = 'm3/s'

            return_period_20_var = \
                return_period_nc.createVariable('return_period_20',
                                                'f8', ('rivid',))
            return_period_20_var.long_name = '20 year return period flow'
            return_period_20_var.units = 'm3/s'

        if method == 'log_pearson':

            return_period_100_var = \
                return_period_nc.createVariable('return_period_100',
                                                'f8', ('rivid',))
            return_period_100_var.long_name = '100 year return period flow'
            return_period_100_var.units = 'm3/s'

            return_period_50_var = \
                return_period_nc.createVariable('return_period_50',
                                                'f8', ('rivid',))
            return_period_50_var.long_name = '50 year return period flow'
            return_period_50_var.units = 'm3/s'

            return_period_25_var = \
                return_period_nc.createVariable('return_period_25',
                                                'f8', ('rivid',))
            return_period_25_var.long_name = '25 year return period flow'
            return_period_25_var.units = 'm3/s'

        return_period_10_var = \
            return_period_nc.createVariable('return_period_10',
                                            'f8', ('rivid',))
        return_period_10_var.long_name = '10 year return period flow'
        return_period_10_var.units = 'm3/s'

        return_period_2_var = \
            return_period_nc.createVariable('return_period_2',
                                            'f8', ('rivid',))
        return_period_2_var.long_name = '2 year return period flow'
        return_period_2_var.units = 'm3/s'

        lat_var = return_period_nc.createVariable('lat', 'f8', ('rivid',),
                                                  fill_value=-9999.0)

        lon_var = return_period_nc.createVariable('lon', 'f8', ('rivid',),
                                                  fill_value=-9999.0)

        add_latlon_metadata(lat_var, lon_var)

        return_period_nc.variables['lat'][:] = \
            qout_nc_file.qout_nc.variables['lat'][:]
        return_period_nc.variables['lon'][:] = \
            qout_nc_file.qout_nc.variables['lon'][:]

        river_id_list = qout_nc_file.get_river_id_array()
        return_period_nc.variables['rivid'][:] = river_id_list

        return_period_nc.return_period_method = method

        return_period_nc.close()

        time_array = qout_nc_file.get_time_array()

    log("Extracting Data and Generating Return Periods ...")
    num_years = int((datetime.utcfromtimestamp(time_array[-1]) -
                     datetime.utcfromtimestamp(time_array[0])).days/365.2425)
    time_steps_per_day = (24 * 3600) / float(
        (datetime.utcfromtimestamp(time_array[1]) -
         datetime.utcfromtimestamp(time_array[0])).total_seconds())
    step = max(1, int(time_steps_per_day * storm_duration_days))

    # generate multiprocessing jobs
    # pylint: disable=no-member
    mp_lock = multiprocessing.Manager().Lock()
    job_combinations = []
    partition_index_list = partition(river_id_list, num_cpus*2)[1]
    for sub_partition_index_list in partition_index_list:
        # pylint: disable=len-as-condition
        if len(sub_partition_index_list) > 0:
            job_combinations.append((qout_file,
                                     return_period_file,
                                     sub_partition_index_list,
                                     step,
                                     num_years,
                                     method,
                                     mp_lock
                                     ))

    pool = multiprocessing.Pool(num_cpus)
    pool.map(generate_single_return_period,
             job_combinations)
    pool.close()
    pool.join()
38063
def _determine_stream_spread_single(sigomatrixEig,
                                    thetasTrack,
                                    sigOmega,
                                    sigAngle,
                                    allinvjacsTrack):
    """"""
    #Estimate the spread in all frequencies and angles
    sigObig2= sigOmega(thetasTrack)**2.
    tsigOdiag= copy.copy(sigomatrixEig[0])
    tsigOdiag[numpy.argmax(tsigOdiag)]= sigObig2
    tsigO= numpy.dot(sigomatrixEig[1],
                     numpy.dot(numpy.diag(tsigOdiag),
                               numpy.linalg.inv(sigomatrixEig[1])))
    #angles
    if hasattr(sigAngle,'__call__'):
        sigangle2= sigAngle(thetasTrack)**2.
    else:
        sigangle2= sigAngle**2.
    tsigadiag= numpy.ones(3)*sigangle2
    tsigadiag[numpy.argmax(tsigOdiag)]= 1.
    tsiga= numpy.dot(sigomatrixEig[1],
                    numpy.dot(numpy.diag(tsigadiag),
                              numpy.linalg.inv(sigomatrixEig[1])))
    #correlations, assume half correlated for now (can be calculated)
    correlations= numpy.diag(0.5*numpy.ones(3))*numpy.sqrt(tsigOdiag*tsigadiag)
    correlations[numpy.argmax(tsigOdiag),numpy.argmax(tsigOdiag)]= 0.
    correlations= numpy.dot(sigomatrixEig[1],
                            numpy.dot(correlations,
                                      numpy.linalg.inv(sigomatrixEig[1])))
    #Now convert
    fullMatrix= numpy.empty((6,6))
    fullMatrix[:3,:3]= tsigO
    fullMatrix[3:,3:]= tsiga
    fullMatrix[3:,:3]= correlations
    fullMatrix[:3,3:]= correlations.T
    return numpy.dot(allinvjacsTrack,numpy.dot(fullMatrix,allinvjacsTrack.T))
587840
def _copy_streamflow_values(self):
        
        log('Creating streamflow variable', 'INFO')
        q_var = self.cf_nc.createVariable(
            self.output_flow_var_name, 'f4', (self.output_id_dim_name, 'time'))
        q_var.long_name = 'Discharge'
        q_var.units = 'm^3/s'
        q_var.coordinates = 'time lat lon z'
        q_var.grid_mapping = 'crs'
        q_var.source = ('Generated by the Routing Application for Parallel '
                        'computatIon of Discharge (RAPID) river routing '
                        'model.')
        q_var.references = 'http://rapid-hub.org/'
        q_var.comment = ('lat, lon, and z values taken at midpoint of river '
                         'reach feature')

        log('Copying streamflow values', 'INFO')
        master_begin_time_step_index = 1
        # to reduce RAM, copy by chunks
        max_2d_dimension = 1000000000  # ~8GB Max
        for raw_nc in self.raw_nc_list:
            max_time_step_size = min(raw_nc.size_time,
                                     max(1, int(float(max_2d_dimension) /
                                                float(raw_nc.size_river_id))))
            raw_nc_begin_time_step_index = 0
            for raw_nc_time_index in \
                    xrange(0, raw_nc.size_time, max_time_step_size):
                time_interval_size = \
                    max(1, min(raw_nc.size_time - raw_nc_time_index,
                               max_time_step_size))

                raw_nc_end_time_step_index = \
                    raw_nc_begin_time_step_index + time_interval_size
                master_end_time_step_index = \
                    master_begin_time_step_index + time_interval_size

                q_var[:,
                      master_begin_time_step_index:master_end_time_step_index]\
                    = raw_nc.get_qout(
                        time_index_start=raw_nc_begin_time_step_index,
                        time_index_end=raw_nc_end_time_step_index)

                master_begin_time_step_index = master_end_time_step_index
                raw_nc_begin_time_step_index = raw_nc_end_time_step_index

        log('Adding initial streamflow values', 'INFO')
        # add initial flow to RAPID output file
        if self.qinit_file and self.rapid_connect_file:
            lookup_table = csv_to_list(self.rapid_connect_file)
            lookup_comids = np.array([int(float(row[0])) for row
                                      in lookup_table])

            init_flow_table = csv_to_list(self.qinit_file)

            for index, comid in enumerate(
                    self.cf_nc.variables[self.output_id_dim_name][:]):
                try:
                    lookup_index = np.where(lookup_comids == comid)[0][0]
                except IndexError:
                    log('COMID %s misssing in rapid_connect file' % comid,
                        'ERROR')
                q_var[index, 0] = float(init_flow_table[lookup_index][0])
        else:
            for index, comid in enumerate(
                    self.cf_nc.variables[self.output_id_dim_name][:]):
                q_var[index, 0] = 0
833806
def _single_multichannel_mp_run(X, Phi, bound, selection_rule, stop_crit,
                                max_iter, verbose=False, pad=0,
                                random_state=None, memory=Memory(None)):
    
    rng = check_random_state(random_state)

    # padding as v stak
    pad = int(pad)
    n_channels = X.shape[0]
    X = np.hstack((np.zeros((n_channels, pad)), X,
                   np.zeros((n_channels, pad))))
    n_samples = X.shape[1]
    n_projs = Phi.doth(X).shape[1]
    err_mse = {}

    # Initialisation
    residual = np.hstack((X.copy(), np.zeros((n_channels,
                                              max(Phi.sizes) / 2))))

    s_rep = np.zeros((n_channels, n_projs))
    X_est = np.zeros((n_channels, n_samples))
    # Main algorithm
    coeffs = np.zeros((n_channels, n_projs))

    it_number = 0
    current_lambda = 1
    for c_idx in range(n_channels):
        err_mse[c_idx] = []
        err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))

    # Decomposition loop: stopping criteria is either SNR or iteration number
    while (current_lambda > bound) & (it_number < max_iter):

        # pick a shift at random : in each size
        rndshifts = {}
        for c_idx in range(n_channels):
            rndshifts[c_idx] = []
        for s_idx, L in enumerate(Phi.sizes):
            shift = rng.randint(low=0, high=L / 4)
            for c_idx in range(n_channels):
                coeffs[c_idx, s_idx * n_samples:(s_idx + 1) * n_samples] = \
                    mdct(residual[c_idx, shift:shift + n_samples], L).ravel()
                rndshifts[c_idx].append(shift)

        # Multichannel mode : we combine projections
        combined = selection_rule(coeffs ** 2)

        # Select a new element
        idx = np.argmax(np.abs(combined))
        # find scale and frequency bin of selected atom
        s_idx = idx // n_samples
        L = Phi.sizes[s_idx]
        F = n_samples // (L // 2)
        frame = (idx - (s_idx * n_samples)) % F
        freq_bin = ((idx - (s_idx * n_samples))) // F

        mdct_wf = memory.cache(mdct_waveform)

        # Update coefficients and residual
        current_lambda_array = np.zeros(n_channels)
        for c_idx in range(n_channels):
            s_rep[c_idx, idx] += coeffs[c_idx, idx]

            # Only one method now : local update via a cached waveform
            pos = (frame * L / 2) - L / 4 + rndshifts[c_idx][s_idx]
            residual[c_idx, pos:pos + L] -= coeffs[c_idx, idx] * \
                mdct_wf(L, freq_bin)

            # also add it to the reconstruction
            X_est[c_idx, pos:pos + L] += coeffs[c_idx, idx] * \
                mdct_wf(L, freq_bin)

            # error computation (err_mse)
            err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))

            current_lambda_array[c_idx] = np.sqrt(
                1. - err_mse[c_idx][-1] / err_mse[c_idx][-2])

        current_lambda = stop_crit(current_lambda_array)

        if verbose:
            print("Iteration %d : Current lambda of %1.4f" % (
                  it_number, current_lambda))
        it_number += 1

    return X_est[:, pad: -pad], err_mse