Plots the abundances of H-1, He-4, C-12 and O-16.
198391
def plot4_nolog(self, num):
        
        self.plot_prof_2(num,'H-1',0.,5.)
        self.plot_prof_2(num,'He-4',0.,5.)
        self.plot_prof_2(num,'C-12',0.,5.)
        self.plot_prof_2(num,'O-16',0.,5.)
        pyl.legend(loc=3)
947875
def plot_results(fout_img, goea_results, **kws):
    """"""
    if "{NS}" not in fout_img:
        plt_goea_results(fout_img, goea_results, **kws)
    else:
        # Plot separately by NS: BP, MF, CC
        ns2goea_results = cx.defaultdict(list)
        for rec in goea_results:
            ns2goea_results[rec.NS].append(rec)
        for ns_name, ns_res in ns2goea_results.items():
            fout = fout_img.format(NS=ns_name)
            plt_goea_results(fout, ns_res, **kws)
688558
def map_with_slider(h3: HistogramND, *, show_zero: bool = True, show_values: bool = False, **kwargs) -> dict:
    
    vega = _create_figure(kwargs)

    values_arr = get_data(h3, kwargs.pop("density", None), kwargs.pop("cumulative", None))
    values = values_arr.tolist()
    value_format = get_value_format(kwargs.pop("value_format", None))

    _add_title(h3, vega, kwargs)
    _create_scales(h3, vega, kwargs)
    _create_axes(h3, vega, kwargs)
    _create_cmap_scale(values_arr, vega, kwargs)
    _create_colorbar(vega, kwargs)

    x = h3.get_bin_centers(0)
    y = h3.get_bin_centers(1)
    x1 = h3.get_bin_left_edges(0)
    x2 = h3.get_bin_right_edges(0)
    y1 = h3.get_bin_left_edges(1)
    y2 = h3.get_bin_right_edges(1)

    data = []
    for i in range(h3.shape[0]):
        for j in range(h3.shape[1]):
            for k in range(h3.shape[2]):
                if not show_zero and values[i][j][k] == 0:
                    continue
                item = {
                    "x": float(x[i]),
                    "x1": float(x1[i]),
                    "x2": float(x2[i]),
                    "y": float(y[j]),
                    "y1": float(y1[j]),
                    "y2": float(y2[j]),
                    "k": k,
                    "c": float(values[i][j][k]),
                }
                if show_values:
                    item["label"] = value_format(values[i][j][k])
                data.append(item)

    vega["signals"] = [
        {
            "name": "k", "value": h3.shape[2] // 2,
            "bind": {
                "input": "range", "min": 0, "max": h3.shape[2] - 1,
                "step": 1, "name": (h3.axis_names[2] or "axis2") + " [slice]"
            }
        }
    ]

    vega["data"] = [{
        "name": "table",
        "values": data,
        "transform": [
             {
                 "type": "filter",
                 "expr": "k == datum.k",
             }
        ]
    }]

    vega["marks"] = [
        {
            "type": "rect",
            "from": {"data": "table"},
            "encode": {
                "enter": {
                    "x": {"scale": "xscale", "field": "x1"},
                    "x2": {"scale": "xscale", "field": "x2"},
                    "y": {"scale": "yscale", "field": "y1"},
                    "y2": {"scale": "yscale", "field": "y2"},
                    "fill": {"scale": "color", "field": "c"},
                    "stroke": {"value": 0},
                    # "strokeWidth": {"value": 0},
                    # "fillColor": {"value": "#ffff00"}
                },
                # "update": {
                #     "fillOpacity": {"value": 0.6}
                # },
                # "hover": {
                #     "fillOpacity": {"value": 0.5}
                # }
            }
        }
    ]

    if show_values:
        vega["marks"].append(
            {
                "type": "text",
                "from": {"data": "table"},
                "encode": {
                    "enter": {
                        "align": {"value": "center"},
                        "baseline": {"value": "middle"},
                        "fontSize": {"value": 13},
                        "fontWeight": {"value": "bold"},
                        "text": {"field": "label"},
                        "x": {"scale": "xscale", "field": "x"},
                        "y": {"scale": "yscale", "field": "y"},
                    }
                }
            }
        )

    return vega
214849
def line(h1: Union[Histogram1D, "HistogramCollection"], ax: Axes, *, errors: bool = False, **kwargs):
    """"""


    show_stats = kwargs.pop("show_stats", False)
    show_values = kwargs.pop("show_values", False)
    density = kwargs.pop("density", False)
    cumulative = kwargs.pop("cumulative", False)
    value_format = kwargs.pop("value_format", None)
    text_kwargs = pop_kwargs_with_prefix("text_", kwargs)
    kwargs["label"] = kwargs.get("label", h1.name)

    data = get_data(h1, cumulative=cumulative, density=density)
    _apply_xy_lims(ax, h1, data, kwargs)
    _add_ticks(ax, h1, kwargs)
    _add_labels(ax, h1, kwargs)

    if errors:
        err_data = get_err_data(h1, cumulative=cumulative, density=density)
        ax.errorbar(h1.bin_centers, data, yerr=err_data, fmt=kwargs.pop(
            "fmt", "-"), ecolor=kwargs.pop("ecolor", "black"), **kwargs)
    else:
        ax.plot(h1.bin_centers, data, **kwargs)

    if show_stats:
        _add_stats_box(h1, ax, stats=show_stats)
    if show_values:
        _add_values(ax, h1, data, value_format=value_format, **text_kwargs)
204748
def show_stacking(self):
        """"""
        grp = self.getPseudoBondGroup("pi-Stacking-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        grp.lineType = self.chimera.Dash
        for i, stack in enumerate(self.plcomplex.pistacking):

            m = self.model
            r = m.newResidue("pseudoatoms", " ", 1, " ")
            centroid_prot = m.newAtom("CENTROID", self.chimera.Element("CENTROID"))
            x, y, z = stack.proteinring_center
            centroid_prot.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(centroid_prot)

            centroid_lig = m.newAtom("CENTROID", self.chimera.Element("CENTROID"))
            x, y, z = stack.ligandring_center
            centroid_lig.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(centroid_lig)

            b = grp.newPseudoBond(centroid_lig, centroid_prot)
            b.color = self.colorbyname('forest green')

            self.bs_res_ids += stack.proteinring_atoms
455825
def run(self):
        

        S = self.model.stoichiometry_matrix
        amat = self.model.propensities
        ymat = self.model.species
        n_species = len(ymat)

        # dPdt is matrix of each species differentiated w.r.t. time
        # The code below literally multiplies the stoichiometry matrix to a column vector of propensities
        # from the right (::math::`\frac{dP}{dt} = \mathbf{Sa}`)
        dPdt =  S * amat


        # A Is a matrix of each species (rows) and the derivatives of their stoichiometry matrix rows
        # against each other species
        # Code below computes the matrix A, that is of size `len(ymat) x len(ymat)`, for which each entry
        # ::math::`A_{ik} = \sum_j S_{ij} \frac{\partial a_j}{\partial y_k} = \mathfb{S_i} \frac{\partial \mathbf{a}}{\partial y_k}`
        A = sp.Matrix(len(ymat), len(ymat), lambda i, j: 0)
        for i in range(A.rows):
            for k in range(A.cols):
                A[i, k] = reduce(operator.add, [S[i, j] * sp.diff(amat[j], ymat[k]) for j in range(len(amat))])


        # `diagA` is a matrix that has values sqrt(a[i]) on the diagonal (0 elsewhere)
        diagA = sp.Matrix(len(amat), len(amat), lambda i, j: amat[i] ** sp.Rational(1,2) if i==j else 0)
        # E is stoichiometry matrix times diagA
        E = S * diagA

        variance_terms = []
        cov_matrix = []
        for i in range(len(ymat)):
            row = []
            for j in range(len(ymat)):
                if i <= j:
                    symbol = 'V_{0}_{1}'.format(i, j)
                    variance_terms.append(VarianceTerm(position=(i,j), symbol=symbol))
                else:
                    # Since Vi,j = Vj,i, i.e. covariance are equal, we only record Vi,j but not Vj,i
                    symbol = 'V_{0}_{1}'.format(j, i)
                    variance_terms.append(VarianceTerm(position=(j,i), symbol=symbol))
                row.append(symbol)
            cov_matrix.append(row)

        V = sp.Matrix(cov_matrix)


        # Matrix of variances (diagonal) and covariances of species i and j differentiated wrt time.
        # I.e. if i=j, V_ij is the variance, and if i!=j, V_ij is the covariance between species i and species j
        dVdt = A * V + V * (A.T) + E * (E.T)


        # build ODEProblem object
        rhs_redundant = sp.Matrix([i for i in dPdt] + [i for i in dVdt])


        #generate ODE terms
        n_vectors = [tuple([1 if i==j else 0 for i in range(n_species)]) for j in range(n_species)]
        moment_terms = [Moment(nvec,lhs) for (lhs, nvec) in zip(ymat, n_vectors)]

        ode_description = moment_terms + variance_terms


        non_redundant_idx = []
        ode_terms = []
        # remove repetitive covariances, as Vij = Vji
        for i, cov in enumerate(ode_description):
            if cov in ode_terms:
                continue
            else:
                ode_terms.append(cov)
                non_redundant_idx.append(i)
        rhs = []
        for i in non_redundant_idx:
            rhs.append(rhs_redundant[i])


        out_problem = ODEProblem("LNA", ode_terms, rhs, sp.Matrix(self.model.parameters))


        return out_problem
182376
def _plot_residuals_to_ax(
    data_all, model_ML, ax, e_unit=u.eV, sed=True, errorbar_opts={}
):
    """"""
    if "group" not in data_all.keys():
        data_all["group"] = np.zeros(len(data_all))

    groups = np.unique(data_all["group"])

    MLf_unit, MLsedf = sed_conversion(model_ML[0], model_ML[1].unit, sed)
    MLene = model_ML[0].to(e_unit)
    MLflux = (model_ML[1] * MLsedf).to(MLf_unit)

    ax.axhline(0, color="k", lw=1, ls="--")

    interp = False
    if data_all["energy"].size != MLene.size or not np.allclose(
        data_all["energy"].value, MLene.value
    ):
        interp = True
        from scipy.interpolate import interp1d

        modelfunc = interp1d(MLene.value, MLflux.value, bounds_error=False)

    for g in groups:
        groupidx = np.where(data_all["group"] == g)

        data = data_all[groupidx]

        notul = ~data["ul"]
        df_unit, dsedf = sed_conversion(data["energy"], data["flux"].unit, sed)
        ene = data["energy"].to(e_unit)
        xerr = u.Quantity((data["energy_error_lo"], data["energy_error_hi"]))
        flux = (data["flux"] * dsedf).to(df_unit)
        dflux = (data["flux_error_lo"] + data["flux_error_hi"]) / 2.0
        dflux = (dflux * dsedf).to(df_unit)[notul]

        if interp:
            difference = flux[notul] - modelfunc(ene[notul]) * flux.unit
        else:
            difference = flux[notul] - MLflux[groupidx][notul]

        # wrap around color and marker cycles
        color = color_cycle[int(g) % len(color_cycle)]
        marker = marker_cycle[int(g) % len(marker_cycle)]

        opts = dict(
            zorder=100,
            marker=marker,
            ls="",
            elinewidth=2,
            capsize=0,
            mec=color,
            mew=0.1,
            ms=6,
            color=color,
        )
        opts.update(errorbar_opts)

        ax.errorbar(
            ene[notul].value,
            (difference / dflux).decompose().value,
            yerr=(dflux / dflux).decompose().value,
            xerr=xerr[:, notul].to(e_unit).value,
            **opts
        )

    from matplotlib.ticker import MaxNLocator

    ax.yaxis.set_major_locator(
        MaxNLocator(5, integer="True", prune="upper", symmetric=True)
    )

    ax.set_ylabel(r"$\Delta\sigma$")
    ax.set_xscale("log")
78532
def plot_gos(fout_png, goids, obo_dag, *args, **kws):
    """"""
    engine = kws['engine'] if 'engine' in kws else 'pydot'
    godagsmall = OboToGoDagSmall(goids=goids, obodag=obo_dag).godag
    godagplot = GODagSmallPlot(godagsmall, *args, **kws)
    godagplot.plt(fout_png, engine)
375767
def plot(q: np.ndarray, bodies: List[str], plot_colors: Dict[str, str],
         sim_name: str, fname: Optional[str] = None):
    
    # Get N and number of dims
    N, dims = q.shape

    # Slices for 
    x_slice = slice(0, dims, 3)
    y_slice = slice(1, dims, 3)
    # Convert all distances from meters to astronomical units (AU)
    plot_x = q[:, x_slice] / au2m
    plot_y = q[:, y_slice] / au2m
    
    # Set up chart title and scale
    fig, ax = plt.subplots(figsize=[12,12])
    ax.set_title(f'Inner Planetary Orbits in 2018; Weekly from {sim_name}')
    ax.set_xlabel('x in J2000.0 Frame; Astronomical Units (au)')
    ax.set_ylabel('y in J2000.0 Frame; Astronomical Units (au)')
    # Scale and tick size
    a = 5.0
    da = 1.0
    ticks = np.arange(-a, a+da, da)
    # Set limits and ticks
    ax.set_xlim(-a, a)
    ax.set_ylim(-a, a)
    ax.set_xticks(ticks)
    ax.set_yticks(ticks)

    # Set marker sizes proportional to size of bodies
    radius_earth = radius_tbl['earth']
    markersize_earth = 4.0
    markersize_tbl = {body : cbrt(radius_tbl[body] / radius_earth) * markersize_earth for body in bodies}
    
    # Plot the orbit of each body
    for k, body in enumerate(bodies):
        ax.plot(plot_x[:, k], plot_y[:, k], label=body, color=plot_colors[body], 
                linewidth=0, markersize = markersize_tbl[body], marker='o')
    
    # Legend and grid
    fig.legend(loc=7, bbox_to_anchor=(0.85, 0.5))
    # ax.legend()
    ax.grid()

    # Save plot if a filename was provided
    if fname is not None:
        fig.savefig(fname, bbox_inches='tight')

    # Display plot
    plt.show()
168858
def plot_density(population,  # pylint: disable=too-many-arguments, too-many-locals
                 bins=100, new_fig=True, subplot=111, levels=None, plane='xy',
                 colorlabel='Nodes per unit area', labelfontsize=16,
                 color_map='Reds', no_colorbar=False, threshold=0.01,
                 neurite_type=NeuriteType.basal_dendrite, **kwargs):
    '''
    '''
    fig, ax = common.get_figure(new_fig=new_fig, subplot=subplot)

    H1, xedges1, yedges1 = extract_density(population, plane=plane, bins=bins,
                                           neurite_type=neurite_type)

    mask = H1 < threshold  # mask = H1==0
    H2 = np.ma.masked_array(H1, mask)

    getattr(plt.cm, color_map).set_bad(color='white', alpha=None)

    plots = ax.contourf((xedges1[:-1] + xedges1[1:]) / 2,
                        (yedges1[:-1] + yedges1[1:]) / 2,
                        np.transpose(H2), # / np.max(H2),
                        cmap=getattr(plt.cm, color_map), levels=levels)

    if not no_colorbar:
        cbar = plt.colorbar(plots)
        cbar.ax.set_ylabel(colorlabel, fontsize=labelfontsize)

    kwargs['title'] = kwargs.get('title', '')
    kwargs['xlabel'] = kwargs.get('xlabel', plane[0])
    kwargs['ylabel'] = kwargs.get('ylabel', plane[1])

    return common.plot_style(fig=fig, ax=ax, **kwargs)
1087119
def parse(self):
        ''''''
        for self.line in self.output:
            # Parse general data: charge, multiplicity, coordinates, etc.
            self.index = 0

            if self.line[1:13] == 'Total Charge':
                tokens = self.line.split()
                self.charge = int(tokens[-1])

            if (self.line[1:13] or self.line[0:12]) == 'Multiplicity':
                tokens = self.line.split()
                self.multiplicity = int(tokens[-1])

            if self.line[0:33] == 'CARTESIAN COORDINATES (ANGSTROEM)':
                if not hasattr(self, 'names'):
                    self.names = dict()
                if not hasattr(self, 'coords'):
                    self.coords = dict()
                self.line = self._skip_lines(2)
                names = list()
                coords = list()
                while self.line.strip():
                    tokens = self.line.split()
                    names.append(tokens[0])
                    x = float(tokens[1])
                    y = float(tokens[2])
                    z = float(tokens[3])
                    coords.append((x, y, z))
                    self.line = next(self.output)
                self.names = np.array(names)
                self.coords[self.index] = np.array(coords)

            if self.line[22:50] == 'MULLIKEN POPULATION ANALYSIS':
                if not hasattr(self, 'populations'):
                    self.populations = dict()
                self.line = self._skip_lines(6)
                populations = list()
                while self.line.strip() and 'Sum' not in self.line:
                    tokens = self.line.split()
                    populations.append((float(tokens[-2]), float(tokens[-1])))
                    self.line = next(self.output)
                self.populations['mulliken'][self.index] = np.array(populations) # noqa

            # Parse data from the EPR/NMR module
            if self.line[37:44] == 'EPR/NMR':
                self.eprnmr = dict()

            if self.line[0:19] == 'ELECTRONIC G-MATRIX':
                self.line = self._skip_lines(4)
                self.eprnmr['g']['tensor'] = self._parse_tensor()

            if self.line[0:27] == 'ZERO-FIELD-SPLITTING TENSOR':
                self.line = self._skip_lines(4)
                self.eprnmr['zfs']['tensor'] = self._parse_tensor()

            if self.line[1:8] == 'Nucleus':
                tokens = self.line.split()
                nucleus = int(re.findall(r'\d+', tokens[1])[0])
                while 'Raw HFC' not in self.line:
                    self.line = self._skip_lines(1)
                self.line = self._skip_lines(2)
                self.eprnmr['hfc'][nucleus]['tensor'] = self._parse_tensor()
                self.line = self._skip_lines(1)
                self.eprnmr['hfc'][nucleus]['fc'] = self._parse_components()
                self.eprnmr['hfc'][nucleus]['sd'] = self._parse_components()
                self.line = self._skip_lines(1)
                self.eprnmr['hfc'][nucleus]['orb'] = self._parse_components()
                self.eprnmr['hfc'][nucleus]['dia'] = self._parse_components()

            # Parse data from the MRCI module
            if self.line[36:43] == 'M R C I':
                self.mrci = dict()

            if self.line[1:19] == 'SPIN-SPIN COUPLING':
                self.line = self._skip_lines(4)
                self.mrci['zfs']['ssc']['tensor'] = self._parse_tensor()

            if self.line[1:30] == '2ND ORDER SPIN-ORBIT COUPLING':
                while 'Second' not in self.line:
                    self.line = self._skip_lines(1)
                self.line = self._skip_lines(1)
                self.mrci['zfs']['soc']['second_order']['0']['tensor'] = self._parse_tensor() # noqa
                self.line = self._skip_lines(2)
                self.mrci['zfs']['soc']['second_order']['m']['tensor'] = self._parse_tensor() # noqa
                self.line = self._skip_lines(2)
                self.mrci['zfs']['soc']['second_order']['p']['tensor'] = self._parse_tensor() # noqa

            if self.line[1:42] == 'EFFECTIVE HAMILTONIAN SPIN-ORBIT COUPLING':
                self.line = self._skip_lines(4)
                self.mrci['zfs']['soc']['heff']['tensor'] = self._parse_tensor()
225334
function() {
        var tmp = [], //12
			src = [], //16
			dst = [], //16
			mat = this.toArray(),
			i = 0;

        for (i = 0; i < 4; i++) {
            var i4 = i << 2;
            src[i] = mat[i4];
            src[i + 4] = mat[i4 + 1];
            src[i + 8] = mat[i4 + 2];
            src[i + 12] = mat[i4 + 3];
        }

        // calculate pairs for first 8 elements (cofactors)
        tmp[0] = src[10] * src[15];
        tmp[1] = src[11] * src[14];
        tmp[2] = src[9] * src[15];
        tmp[3] = src[11] * src[13];
        tmp[4] = src[9] * src[14];
        tmp[5] = src[10] * src[13];
        tmp[6] = src[8] * src[15];
        tmp[7] = src[11] * src[12];
        tmp[8] = src[8] * src[14];
        tmp[9] = src[10] * src[12];
        tmp[10] = src[8] * src[13];
        tmp[11] = src[9] * src[12];

        // calculate first 8 elements (cofactors)
        var src0 = src[0],
			src1 = src[1],
			src2 = src[2],
			src3 = src[3],
			src4 = src[4],
			src5 = src[5],
			src6 = src[6],
			src7 = src[7];
		dst[0] = tmp[0] * src5 + tmp[3] * src6 + tmp[4] * src7;
		dst[0] -= tmp[1] * src5 + tmp[2] * src6 + tmp[5] * src7;
		dst[1] = tmp[1] * src4 + tmp[6] * src6 + tmp[9] * src7;
		dst[1] -= tmp[0] * src4 + tmp[7] * src6 + tmp[8] * src7;
		dst[2] = tmp[2] * src4 + tmp[7] * src5 + tmp[10] * src7;
		dst[2] -= tmp[3] * src4 + tmp[6] * src5 + tmp[11] * src7;
		dst[3] = tmp[5] * src4 + tmp[8] * src5 + tmp[11] * src6;
		dst[3] -= tmp[4] * src4 + tmp[9] * src5 + tmp[10] * src6;
		dst[4] = tmp[1] * src1 + tmp[2] * src2 + tmp[5] * src3;
		dst[4] -= tmp[0] * src1 + tmp[3] * src2 + tmp[4] * src3;
		dst[5] = tmp[0] * src0 + tmp[7] * src2 + tmp[8] * src3;
		dst[5] -= tmp[1] * src0 + tmp[6] * src2 + tmp[9] * src3;
		dst[6] = tmp[3] * src0 + tmp[6] * src1 + tmp[11] * src3;
		dst[6] -= tmp[2] * src0 + tmp[7] * src1 + tmp[10] * src3;
		dst[7] = tmp[4] * src0 + tmp[9] * src1 + tmp[10] * src2;
		dst[7] -= tmp[5] * src0 + tmp[8] * src1 + tmp[11] * src2;

        // calculate pairs for second 8 elements (cofactors)
		tmp[0] = src2 * src7;
		tmp[1] = src3 * src6;
		tmp[2] = src1 * src7;
		tmp[3] = src3 * src5;
		tmp[4] = src1 * src6;
		tmp[5] = src2 * src5;
		tmp[6] = src0 * src7;
		tmp[7] = src3 * src4;
		tmp[8] = src0 * src6;
		tmp[9] = src2 * src4;
		tmp[10] = src0 * src5;
		tmp[11] = src1 * src4;

        // calculate second 8 elements (cofactors)
		src0 = src[8];
		src1 = src[9];
		src2 = src[10];
		src3 = src[11];
		src4 = src[12];
		src5 = src[13];
		src6 = src[14];
		src7 = src[15];
		dst[8] = tmp[0] * src5 + tmp[3] * src6 + tmp[4] * src7;
		dst[8] -= tmp[1] * src5 + tmp[2] * src6 + tmp[5] * src7;
		dst[9] = tmp[1] * src4 + tmp[6] * src6 + tmp[9] * src7;
		dst[9] -= tmp[0] * src4 + tmp[7] * src6 + tmp[8] * src7;
		dst[10] = tmp[2] * src4 + tmp[7] * src5 + tmp[10] * src7;
		dst[10] -= tmp[3] * src4 + tmp[6] * src5 + tmp[11] * src7;
		dst[11] = tmp[5] * src4 + tmp[8] * src5 + tmp[11] * src6;
		dst[11] -= tmp[4] * src4 + tmp[9] * src5 + tmp[10] * src6;
		dst[12] = tmp[2] * src2 + tmp[5] * src3 + tmp[1] * src1;
		dst[12] -= tmp[4] * src3 + tmp[0] * src1 + tmp[3] * src2;
		dst[13] = tmp[8] * src3 + tmp[0] * src0 + tmp[7] * src2;
		dst[13] -= tmp[6] * src2 + tmp[9] * src3 + tmp[1] * src0;
		dst[14] = tmp[6] * src1 + tmp[11] * src3 + tmp[3] * src0;
		dst[14] -= tmp[10] * src3 + tmp[2] * src0 + tmp[7] * src1;
		dst[15] = tmp[10] * src2 + tmp[4] * src0 + tmp[9] * src1;
		dst[15] -= tmp[8] * src1 + tmp[11] * src2 + tmp[5] * src0;

		var det = 1.0 / (src[0] * dst[0] + src[1] * dst[1] + src[2] * dst[2] + src[3] * dst[3]);
		for (i = 0, k = 0; i < 4; i++) {
			var m = this.matrix[i];
			for (var j = 0; j < 4; j++) {
				m[j] = dst[k++] * det;
			}
		}
		return this;
    }
708005
def plot_eps_data_hist(self, dfs):
        
        # check if this is a DC inversion
        if 'datum' in dfs[0]:
            dc_inv = True
        else:
            dc_inv = False

        nr_y = len(dfs)
        size_y = 5 / 2.54 * nr_y
        if dc_inv:
            nr_x = 1
        else:
            nr_x = 3
        size_x = 15 / 2.54

        fig, axes = plt.subplots(nr_y, nr_x, figsize=(size_x, size_y))
        axes = np.atleast_2d(axes)

        # plot initial data errors
        df = dfs[0]
        if dc_inv:
            ax = axes[0, 0]
            ax.hist(
                df['datum'] / df['eps_r'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|) / \epsilon_r$')
            ax.set_ylabel(r'count')
        else:
            # complex inversion
            ax = axes[0, 0]
            ax.hist(
                df['-log(|R|)'] / df['eps'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|)$')
            ax.set_ylabel(r'count')

            ax = axes[0, 1]
            ax.hist(
                df['-log(|R|)'] / df['eps_r'],
                100,
            )
            ax.set_xlabel(r'$-log(|R|) / \epsilon_r$')
            ax.set_ylabel(r'count')

            ax = axes[0, 2]
            phase_data = df['-Phase(rad)'] / df['eps_p']
            if not np.all(np.isinf(phase_data) | np.isnan(phase_data)):
                ax.hist(
                    phase_data,
                    100,
                )
            ax.set_xlabel(r'$-\phi[rad] / \epsilon_p$')
            ax.set_ylabel(r'count')

        # iterations
        for it, df in enumerate(dfs[1:]):
            ax = axes[1 + it, 0]
            ax.hist(
                df['psi'],
                100
            )
            rms = np.sqrt(
                1 / df['psi'].shape[0] *
                np.sum(
                    df['psi'] ** 2
                )
            )
            ax.axvline(rms, color='k', linestyle='dashed')
            ax.set_title('iteration: {0}'.format(it))
            ax.set_xlabel('psi')
            ax.set_ylabel(r'count')

            ax = axes[1 + it, 1]
            Rdat = df['Re(d)']
            Rmod = df['Re(f(m))']

            ax.scatter(
                Rdat,
                Rmod,
            )
            ax.set_xlabel(r'$log(R_{data}~[\Omega])$')
            ax.set_ylabel(r'$log(R_{mod}~[\Omega])$')

            ax = axes[1 + it, 2]
            phidat = df['Im(d)']
            phimod = df['Im(f(m))']

            ax.scatter(
                phidat,
                phimod,
            )
            ax.set_xlabel(r'$\phi_{data}~[mrad]$')
            ax.set_ylabel(r'$\phi_{mod}~[mrad]$')

        fig.tight_layout()
        fig.savefig('eps_plot_hist.png', dpi=300)
1072690
def H_donor_count(mol):
    
    mol.require("Valence")
    return sum(1 for _, a in mol.atoms_iter() if a.H_donor)
283272
def plot_eqarea_pars(self, pars, fig):
        
        if pars == {}:
            pass
        elif 'calculation_type' in list(pars.keys()) and pars['calculation_type'] == 'DE-BFP':
            ymin, ymax = fig.get_ylim()
            xmin, xmax = fig.get_xlim()

            D_c, I_c = pmag.circ(pars["specimen_dec"],
                                 pars["specimen_inc"], 90)
            X_c_up, Y_c_up = [], []
            X_c_d, Y_c_d = [], []
            for k in range(len(D_c)):
                XY = pmag.dimap(D_c[k], I_c[k])
                if I_c[k] < 0:
                    X_c_up.append(XY[0])
                    Y_c_up.append(XY[1])
                if I_c[k] > 0:
                    X_c_d.append(XY[0])
                    Y_c_d.append(XY[1])
            fig.plot(X_c_d, Y_c_d, 'b', lw=0.5)
            fig.plot(X_c_up, Y_c_up, 'c', lw=0.5)
            if self.ie_open:
                self.ie.plot(X_c_d, Y_c_d, 'b', lw=0.5)
                self.ie.plot(X_c_up, Y_c_up, 'c', lw=0.5)

            fig.set_xlim(xmin, xmax)
            fig.set_ylim(ymin, ymax)
        # plot best-fit direction
        else:
            if "specimen_dec" in list(pars.keys()) and "specimen_inc" in list(pars.keys()):
                dec = pars["specimen_dec"]
                inc = pars["specimen_inc"]
            elif "dec" in list(pars.keys()) and "inc" in list(pars.keys()):
                dec = pars["dec"]
                inc = pars["inc"]
            else:
                print(("either dec or inc missing from values recived for high level plot, was given %s, aborting" % (
                    str(pars))))
                return
            XY = pmag.dimap(float(dec), float(inc))
            if inc > 0:
                if 'color' in list(pars.keys()):
                    FC = pars['color']
                    EC = pars['color']
                    SIZE = 15*self.GUI_RESOLUTION
                else:
                    FC = 'grey'
                    EC = 'grey'
                    SIZE = 15*self.GUI_RESOLUTION
            else:
                if 'color' in list(pars.keys()):
                    FC = 'white'
                    EC = pars['color']
                    SIZE = 15*self.GUI_RESOLUTION
                else:
                    FC = 'white'
                    EC = 'grey'
                    SIZE = 15*self.GUI_RESOLUTION
            fig.scatter([XY[0]], [XY[1]], marker='o', edgecolor=EC,
                        facecolor=FC, s=SIZE, lw=1, clip_on=False)
            if self.ie_open:
                self.ie.scatter([XY[0]], [XY[1]], marker='o', edgecolor=EC,
                                facecolor=FC, s=SIZE, lw=1, clip_on=False)
1072714
def bbt_simple_plot(self):
        

        # First, sum the different types of alignment counts
        data = OrderedDict()
        cats = OrderedDict()
        for s_name in self.bbt_data:
            data[s_name] = OrderedDict()
            for org in self.bbt_data[s_name]:
                data[s_name][org] = self.bbt_data[s_name][org]['hits'] - self.bbt_data[s_name][org]['shared']
                if org not in cats and org != 'multiMatch' and org != 'noMatch':
                    if org.lower().endswith('.fa'):
                        cname = org[:-3]
                    elif org.lower().endswith('.fasta'):
                        cname = org[:-6]
                    else:
                        cname = org
                    cats[org] = { 'name': cname }

        pconfig = {
            'id': 'biobloom_tools',
            'title': 'BioBloom Tools: Alignment counts per species',
            'ylab': 'Number of hits',
            'hide_zero_cats': False
        }
        cats['multiMatch'] = { 'name': 'Multiple Genomes', 'color': '#820000' }
        cats['noMatch'] = { 'name': 'No Match', 'color': '#cccccc' }

        return bargraph.plot(data, cats, pconfig)
373850
def _adjust_stoichiometry(self):
        
        num_O = len(list(self.particles_by_name('O')))
        num_Si = len(list(self.particles_by_name('Si')))
        n_deletions = num_O - 2*num_Si

        bottom_Os = [atom for atom in self.particles()
                     if atom.name == 'O' and
                        atom.pos[2] < self._O_buffer and
                        len(self.bond_graph.neighbors(atom)) == 1]

        for _ in range(n_deletions):
            O1 = random.choice(bottom_Os)
            bottom_Os.remove(O1)
            self.remove(O1)