Return a Graph of a quantile-quantile (QQ) plot and confidence band
320950
def qqgraph(h1, h2, quantiles=None):
    
    if quantiles is None:
        quantiles = max(min(len(h1), len(h2)) / 2, 1)
    nq = quantiles
    # position where to compute the quantiles in [0, 1]
    xq = array('d', [0.] * nq)
    # array to contain the quantiles
    yq1 = array('d', [0.] * nq)
    # array to contain the quantiles
    yq2 = array('d', [0.] * nq)

    for i in range(nq):
        xq[i] = float(i + 1) / nq

    h1.GetQuantiles(nq, yq1, xq)
    h2.GetQuantiles(nq, yq2, xq)

    xq_plus = array('d', [0.] * nq)
    xq_minus = array('d', [0.] * nq)
    yq2_plus = array('d', [0.] * nq)
    yq2_minus = array('d', [0.] * nq)

    """
    KS_cv: KS critical value

               1.36
    KS_cv = -----------
             sqrt( N )

    Where 1.36 is for alpha = 0.05 (confidence level 1-5%=95%, about 2 sigma)

    For 1 sigma (alpha=0.32, CL=68%), the value in the nominator is 0.9561,
    it is gotten by GetCriticalValue(1, 1 - 0.68).

    Notes
    -----

    * For 1-sample KS test (data and theoretic), N should be n

    * For 2-sample KS test (2 data set), N should be sqrt(m*n/(m+n))!
      Here is the case m or n (size of samples) should be effective size
      for a histogram

    * Critical value here is valid for only for sample size >= 80 (some
      references say 35) which means, for example, for a unweighted histogram,
      it must have more than 80 (or 35) entries filled and then confidence
      band is reliable.

    """

    esum1 = effective_sample_size(h1)
    esum2 = effective_sample_size(h2)

    # one sigma band
    KS_cv = (critical_value(1, 1 - 0.68) /
             sqrt((esum1 * esum2) / (esum1 + esum2)))

    for i in range(nq):
        # upper limit
        xq_plus[i] = float(xq[i] + KS_cv)
        # lower limit
        xq_minus[i] = float(xq[i] - KS_cv)

    h2.GetQuantiles(nq, yq2_plus, xq_plus)
    h2.GetQuantiles(nq, yq2_minus, xq_minus)

    yq2_err_plus = array('d', [0.] * nq)
    yq2_err_minus = array('d', [0.] * nq)
    for i in range(nq):
        yq2_err_plus[i] = yq2_plus[i] - yq2[i]
        yq2_err_minus[i] = yq2[i] - yq2_minus[i]

    # forget the last point, so number of points: (nq - 1)
    gr = Graph(nq - 1)
    for i in range(nq - 1):
        gr[i] = (yq1[i], yq2[i])
        # confidence level band
        gr.SetPointEYlow(i, yq2_err_minus[i])
        gr.SetPointEYhigh(i, yq2_err_plus[i])

    return gr
743166
def plot_result(x_p, y_p, y_p_e, smoothed_data, smoothed_data_diff, filename=None):
    ''' 
    '''
    logging.info('Plot results')
    plt.close()

    p1 = plt.errorbar(x_p * analysis_configuration['vcal_calibration'], y_p, yerr=y_p_e, fmt='o')  # plot data with error bars
    p2, = plt.plot(x_p * analysis_configuration['vcal_calibration'], smoothed_data, '-r')  # plot smoothed data
    factor = np.amax(y_p) / np.amin(smoothed_data_diff) * 1.1
    p3, = plt.plot(x_p * analysis_configuration['vcal_calibration'], factor * smoothed_data_diff, '-', lw=2)  # plot differentiated data
    mpv_index = np.argmax(-analysis_utils.smooth_differentiation(x_p, y_p, weigths=1 / y_p_e, order=3, smoothness=analysis_configuration['smoothness'], derivation=1))
    p4, = plt.plot([x_p[mpv_index] * analysis_configuration['vcal_calibration'], x_p[mpv_index] * analysis_configuration['vcal_calibration']], [0, factor * smoothed_data_diff[mpv_index]], 'k-', lw=2)
    text = 'MPV ' + str(int(x_p[mpv_index] * analysis_configuration['vcal_calibration'])) + ' e'
    plt.text(1.01 * x_p[mpv_index] * analysis_configuration['vcal_calibration'], -10. * smoothed_data_diff[mpv_index], text, ha='left')
    plt.legend([p1, p2, p3, p4], ['data', 'smoothed spline', 'spline differentiation', text], prop={'size': 12}, loc=0)
    plt.title('\'Single hit cluster\'-occupancy for different pixel thresholds')
    plt.xlabel('Pixel threshold [e]')
    plt.ylabel('Single hit cluster occupancy [a.u.]')
    plt.ylim(0, np.amax(y_p) * 1.15)
    if filename is None:
        plt.show()
    else:
        filename.savefig(plt.gcf())
    return smoothed_data_diff
729502
def param_describe(params, quant=95, axis=0):
    """"""
    par = np.mean(params, axis=axis)
    lo, up = perc(quant)
    p_up = np.percentile(params, up, axis=axis)
    p_lo = np.percentile(params, lo, axis=axis)
    return par, p_lo, p_up
1325658
def plot_fit(self,intervals=False,**kwargs):
        
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))
        series_type = kwargs.get('series_type','Smoothed')

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            date_index = copy.deepcopy(self.index)
            date_index = date_index[:self.y.shape[0]+1]

            if series_type == 'Smoothed':
                mu, V = self.smoothed_state(self.data,self.latent_variables.get_z_values())
            elif series_type == 'Filtered':
                mu, V, _, _, _ = self._model(self.data,self.latent_variables.get_z_values())
            else:
                mu, V = self.smoothed_state(self.data,self.latent_variables.get_z_values())

            # Create smoothed/filtered aggregate series
            _, Z, _, _, _ = self._ss_matrices(self.latent_variables.get_z_values())
            smoothed_series = np.zeros(self.y.shape[0])

            for t in range(0,self.y.shape[0]):
                smoothed_series[t] = np.dot(Z[t],mu[:,t])

            plt.figure(figsize=figsize) 
            
            plt.subplot(self.z_no+1, 1, 1)
            plt.title(self.y_name + " Raw and " + series_type)  
            plt.plot(date_index,self.data,label='Data')
            plt.plot(date_index,smoothed_series,label=series_type,c='black')
            plt.legend(loc=2)

            for coef in range(0,self.z_no-1):
                V_coef = V[0][coef][:-1]    
                plt.subplot(self.z_no+1, 1, 2+coef)
                plt.title("Beta " + self.X_names[coef]) 

                if intervals == True:
                    alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                    plt.fill_between(date_index[5:], mu[coef,0:mu.shape[1]-1][5:] + 1.98*np.sqrt(V_coef[5:]), mu[coef,0:mu.shape[1]-1][5:] - 1.98*np.sqrt(V_coef[5:]), alpha=0.15,label='95% C.I.') 
                plt.plot(date_index,mu[coef,0:mu.shape[1]-1],label='Data')
                plt.legend(loc=2)               
            
            plt.subplot(self.z_no+1, 1, self.z_no+1)
            plt.title("Measurement Error")
            plt.plot(date_index,self.data-smoothed_series,label='Irregular')
            plt.legend(loc=2)   

            plt.show()
1364870
def plot_fit(self,intervals=True,**kwargs):
        
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            date_index = copy.deepcopy(self.index)
            date_index = date_index[self.integ:self.data_original.shape[0]+1]
            t_params = self.transform_z()

            if self.model_name2 == 'Skewt':
                states_upper_95 = (self.states + 1.98*np.sqrt(self.states_ses)) + ((t_params[-3] - (1.0/t_params[-3]))*t_params[-2]*gas.SkewtScore.tv_variate_exp(t_params[-1]))
                states_lower_95 = (self.states - 1.98*np.sqrt(self.states_ses)) + ((t_params[-3] - (1.0/t_params[-3]))*t_params[-2]*gas.SkewtScore.tv_variate_exp(t_params[-1])) 
                mean_states = self.states + ((t_params[-3] - (1.0/t_params[-3]))*t_params[-2]*gas.SkewtScore.tv_variate_exp(t_params[-1]))
            else:
                states_upper_95 = (self.states + 1.98*np.sqrt(self.states_ses)) 
                states_lower_95 = (self.states - 1.98*np.sqrt(self.states_ses))
                mean_states = self.states

            plt.figure(figsize=figsize) 
            
            plt.subplot(2, 1, 1)
            plt.title(self.data_name + " Raw and Smoothed") 

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index, self.link(states_lower_95), self.link(states_upper_95), alpha=0.15,label='95% C.I.')   

            plt.plot(date_index, self.data, label='Data')
            plt.plot(date_index, self.link(mean_states), label='Smoothed', c='black')          
            plt.legend(loc=2)
            
            plt.subplot(2, 1, 2)
            plt.title(self.data_name + " Local Level")  

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index, self.link(states_lower_95), self.link(states_upper_95), alpha=0.15,label='95% C.I.')   

            plt.plot(date_index, self.link(mean_states), label='Smoothed State')
            plt.legend(loc=2)
            plt.show()
56268
def sequence_quality_plot (self):
        

        data = dict()
        for s_name in self.fastqc_data:
            try:
                data[s_name] = {self.avg_bp_from_range(d['base']): d['mean'] for d in self.fastqc_data[s_name]['per_base_sequence_quality']}
            except KeyError:
                pass
        if len(data) == 0:
            log.debug('sequence_quality not found in FastQC reports')
            return None

        pconfig = {
            'id': 'fastqc_per_base_sequence_quality_plot',
            'title': 'FastQC: Mean Quality Scores',
            'ylab': 'Phred Score',
            'xlab': 'Position (bp)',
            'ymin': 0,
            'xDecimals': False,
            'tt_label': '<b>Base {point.x}</b>: {point.y:.2f}',
            'colors': self.get_status_cols('per_base_sequence_quality'),
            'yPlotBands': [
                {'from': 28, 'to': 100, 'color': '#c3e6c3'},
                {'from': 20, 'to': 28, 'color': '#e6dcc3'},
                {'from': 0, 'to': 20, 'color': '#e6c3c3'},
            ]
        }
        self.add_section (
            name = 'Sequence Quality Histograms',
            anchor = 'fastqc_per_base_sequence_quality',
            description = 'The mean quality value across each base position in the read.',
            helptext = '''
            To enable multiple samples to be plotted on the same graph, only the mean quality
            scores are plotted (unlike the box plots seen in FastQC reports).

            Taken from the [FastQC help](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/2%20Per%20Base%20Sequence%20Quality.html):

            _The y-axis on the graph shows the quality scores. The higher the score, the better
            the base call. The background of the graph divides the y axis into very good quality
            calls (green), calls of reasonable quality (orange), and calls of poor quality (red).
            The quality of calls on most platforms will degrade as the run progresses, so it is
            common to see base calls falling into the orange area towards the end of a read._
            ''',
            plot = linegraph.plot(data, pconfig)
        )
588645
def isa_from_graph(graph: nx.Graph, oneq_type='Xhalves', twoq_type='CZ') -> ISA:
    
    all_qubits = list(range(max(graph.nodes) + 1))
    qubits = [Qubit(i, type=oneq_type, dead=i not in graph.nodes) for i in all_qubits]
    edges = [Edge(sorted((a, b)), type=twoq_type, dead=False) for a, b in graph.edges]
    return ISA(qubits, edges)
256484
def freqpoly_plot(data):
        """"""
        rel_data = OrderedDict()
        for key, val in data.items():
            tot = sum(val.values(), 0)
            rel_data[key] = {k: v / tot for k, v in val.items()}
        fplotconfig = {
            'data_labels': [
                {'name': 'Absolute', 'ylab': 'Frequency', 'xlab': 'Merged Read Length'},
                {'name': 'Relative', 'ylab': 'Relative Frequency', 'xlab': 'Merged Read Length'}
                ],
            'id': 'flash_freqpoly_plot', 'title': 'FLASh: Frequency of merged read lengths',
            'colors': dict(zip(data.keys(), MultiqcModule.get_colors(len(data))))
            }
        return linegraph.plot([data, rel_data], fplotconfig)
629569
def plot_q_flows(network):
    
    cmap_line = plt.cm.jet
    
    q_flows_max = abs(network.lines_t.q0.abs().max()/(network.lines.s_nom))
    
    ll = network.plot(line_colors = q_flows_max, line_cmap = cmap_line)
    boundaries = [min(q_flows_max), max(q_flows_max)]
    v = np.linspace(boundaries[0], boundaries[1], 101)

    cb = plt.colorbar(ll[1], boundaries=v,
                      ticks=v[0:101:10])

    cb.set_clim(vmin=boundaries[0], vmax=boundaries[1])
823922
def scale_errors(cov_axes, confidence_level=0.95):
    
    dof = len(cov_axes)
    x2t = chi2.ppf(confidence_level,dof)
    return N.sqrt(x2t*cov_axes)
63928
def _qqplot_bar(M=1000000, alphaLevel = 0.05, distr = 'log10'):
	""""""
	mRange=10**(sp.arange(sp.log10(0.5),sp.log10(M-0.5)+0.1,0.1));#should be exp or 10**?
	numPts=len(mRange);
	betaalphaLevel=sp.zeros(numPts);#down in the plot
	betaOneMinusalphaLevel=sp.zeros(numPts);#up in the plot
	betaInvHalf=sp.zeros(numPts);
	for n in range(numPts):
	   m=mRange[n]; #numPLessThanThresh=m;
	   betaInvHalf[n]=st.beta.ppf(0.5,m,M-m);
	   betaalphaLevel[n]=st.beta.ppf(alphaLevel,m,M-m);
	   betaOneMinusalphaLevel[n]=st.beta.ppf(1-alphaLevel,m,M-m);
	betaDown=betaInvHalf-betaalphaLevel;
	betaUp=betaOneMinusalphaLevel-betaInvHalf;

	theoreticalPvals=mRange/M;
	return betaUp, betaDown, theoreticalPvals
218621
def median_hilow(series, confidence_interval=0.95):
    
    tail = (1 - confidence_interval) / 2
    return pd.DataFrame({'y': [np.median(series)],
                         'ymin': np.percentile(series, 100 * tail),
                         'ymax': np.percentile(series, 100 * (1 - tail))})
527754
def summary(self):
        
        ci = 1 - self.alpha
        lower_upper_bounds = self._compute_confidence_bounds_of_parameters()
        df = pd.DataFrame(index=self._fitted_parameter_names)
        df["coef"] = self._fitted_parameters_
        df["se(coef)"] = self._compute_standard_errors().loc["se"]
        df["lower %g" % ci] = lower_upper_bounds["lower-bound"]
        df["upper %g" % ci] = lower_upper_bounds["upper-bound"]
        df["p"] = self._compute_p_values()
        with np.errstate(invalid="ignore", divide="ignore"):
            df["-log2(p)"] = -np.log2(df["p"])
        return df
317137
def plot_fit(self, intervals=True, **kwargs):
        
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))
        series_type = kwargs.get('series_type','Smoothed')

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            date_index = copy.deepcopy(self.index)
            date_index = date_index[self.integ:self.data_original.shape[0]+1]

            if series_type == 'Smoothed':
                mu, V = self.smoothed_state(self.data,self.latent_variables.get_z_values())
            elif series_type == 'Filtered':
                mu, V, _, _, _ = self._model(self.data,self.latent_variables.get_z_values())
            else:
                mu, V = self.smoothed_state(self.data,self.latent_variables.get_z_values())

            mu = mu[0][:-1]
            V = V.ravel()

            plt.figure(figsize=figsize) 
            
            plt.subplot(3, 1, 1)
            plt.title(self.data_name + " Raw and " + series_type)   

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index[2:], mu[2:] + 1.98*np.sqrt(V[:-1][2:]), mu[2:] - 1.98*np.sqrt(V[:-1][2:]), alpha=0.15,label='95% C.I.') 

            plt.plot(date_index,self.data,label='Data')
            plt.plot(date_index,mu,label=series_type,c='black')
            plt.legend(loc=2)
            
            plt.subplot(3, 1, 2)
            plt.title(self.data_name + " Local Level")  

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index[2:], mu[2:] + 1.98*np.sqrt(V[:-1][2:]), mu[2:] - 1.98*np.sqrt(V[:-1][2:]), alpha=0.15,label='95% C.I.') 

            plt.plot(date_index,mu,label='Local Level')
            plt.legend(loc=2)
            plt.subplot(3, 1, 3)
            plt.title("Measurement Noise")  
            plt.plot(date_index,self.data-mu)
            plt.show()
121090
def plot_fit(self,intervals=True,**kwargs):
        
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))
        series_type = kwargs.get('series_type','Smoothed')

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            date_index = copy.deepcopy(self.index)
            date_index = date_index[self.integ:self.data_original.shape[0]+1]

            if series_type == 'Smoothed':
                mu, V= self.smoothed_state(self.data,self.latent_variables.get_z_values())
            elif series_type == 'Filtered':
                mu, V, _, _, _ = self._model(self.data,self.latent_variables.get_z_values())
            else:
                mu, V = self.smoothed_state(self.data,self.latent_variables.get_z_values())

            mu0 = mu[0][:-1]
            mu1 = mu[1][:-1]
            Vlev = V[0][0][:-1]
            Vtrend = V[0][1][:-1]   

            plt.figure(figsize=figsize) 
            
            plt.subplot(2, 2, 1)
            plt.title(self.data_name + " Raw and " + series_type)   

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index[2:], mu0[2:] + 1.98*np.sqrt(Vlev[2:]), mu0[2:] - 1.98*np.sqrt(Vlev[2:]), alpha=0.15,label='95% C.I.')   

            plt.plot(date_index,self.data,label='Data')
            plt.plot(date_index,mu0,label=series_type,c='black')
            plt.legend(loc=2)
            
            plt.subplot(2, 2, 2)

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index[2:], mu0[2:] + 1.98*np.sqrt(Vlev[2:]), mu0[2:] - 1.98*np.sqrt(Vlev[2:]), alpha=0.15,label='95% C.I.')   

            plt.title(self.data_name + " Local Level")  
            plt.plot(date_index,mu0,label='Local Level')
            plt.legend(loc=2)
            
            plt.subplot(2, 2, 3)

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index[2:], mu1[2:] + 1.98*np.sqrt(Vtrend[2:]), mu1[2:] - 1.98*np.sqrt(Vtrend[2:]), alpha=0.15,label='95% C.I.')   

            plt.title(self.data_name + " Trend")    
            plt.plot(date_index,mu1,label='Stochastic Trend')
            plt.legend(loc=2)

            plt.subplot(2, 2, 4)
            plt.title("Measurement Noise")  
            plt.plot(date_index[1:self.data.shape[0]],self.data[1:self.data.shape[0]]-mu0[1:self.data.shape[0]],label='Irregular term')
            plt.show()
3903
def plot_fit(self,intervals=True,**kwargs):
        
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))

        if self.latent_variables.estimated is False:
            raise Exception("No latent variables estimated!")
        else:
            date_index = copy.deepcopy(self.index)
            date_index = date_index[self.integ:self.data_original.shape[0]+1]

            states_0_upper_95 = self.states[0] + 1.98*np.sqrt(self.states_var[0])
            states_0_lower_95 = self.states[0] - 1.98*np.sqrt(self.states_var[0])
            states_1_upper_95 = self.states[1] + 1.98*np.sqrt(self.states_var[1])
            states_1_lower_95 = self.states[1] - 1.98*np.sqrt(self.states_var[1])

            plt.figure(figsize=figsize) 
            
            plt.subplot(2, 2, 1)
            plt.title(self.data_name + " Raw and Smoothed") 

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index, self.link(states_0_lower_95), self.link(states_0_upper_95), alpha=0.15,label='95% C.I.')   

            plt.plot(date_index,self.data,label='Data')
            plt.plot(date_index,self.link(self.states[0]),label="Smoothed",c='black')
            plt.legend(loc=2)
            
            plt.subplot(2, 2, 2)
            plt.title(self.data_name + " Local Level")  

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index, self.link(states_0_lower_95), self.link(states_0_upper_95), alpha=0.15,label='95% C.I.')   

            plt.plot(date_index,self.link(self.states[0]),label='Smoothed State')
            plt.legend(loc=2)
            
            plt.subplot(2, 2, 3)
            plt.title(self.data_name + " Trend")    

            if intervals == True:
                alpha =[0.15*i/float(100) for i in range(50,12,-2)]
                plt.fill_between(date_index, states_1_lower_95, states_1_upper_95, alpha=0.15,label='95% C.I.') 

            plt.plot(date_index,self.states[1],label='Smoothed State')
            plt.legend(loc=2)
            
            plt.show()
682125
def min_freq(self):
        
        i = np.unravel_index(np.argmin(self.bands), self.bands.shape)

        return self.qpoints[i[1]], self.bands[i]