Python numpy.set_printoptions() Examples

def assert_almost_equal(a, b, rtol=None, atol=None, names=('a', 'b'), equal_nan=False):
    """Test that two numpy arrays are almost equal. Raise exception message if not.

    Parameters
    ----------
    a : np.ndarray
    b : np.ndarray
    threshold : None or float
        The checking threshold. Default threshold will be used if set to ``None``.
    """
    rtol = get_rtol(rtol)
    atol = get_atol(atol)
    if almost_equal(a, b, rtol, atol, equal_nan=equal_nan):
        return
    index, rel = find_max_violation(a, b, rtol, atol)
    np.set_printoptions(threshold=4, suppress=True)
    msg = npt.build_err_msg([a, b],
                            err_msg="Error %f exceeds tolerance rtol=%f, atol=%f. "
                                    " Location of maximum error:%s, a=%f, b=%f"
                            % (rel, rtol, atol, str(index), a[index], b[index]),
                            names=names)
    raise AssertionError(msg) 

def get_hash(self):
        """
        Returns the hash of the molecule.
        """

        m = hashlib.sha1()
        concat = ""

        np.set_printoptions(precision=16)
        for field in self.hash_fields:
            data = getattr(self, field)
            if field == "geometry":
                data = float_prep(data, GEOMETRY_NOISE)
            elif field == "fragment_charges":
                data = float_prep(data, CHARGE_NOISE)
            elif field == "molecular_charge":
                data = float_prep(data, CHARGE_NOISE)
            elif field == "masses":
                data = float_prep(data, MASS_NOISE)

            concat += json.dumps(data, default=lambda x: x.ravel().tolist())

        m.update(concat.encode("utf-8"))
        return m.hexdigest() 

def pytest_configure(config):

    if ASTROPY_HEADER:

        config.option.astropy_header = True

        PYTEST_HEADER_MODULES.pop('h5py', None)
        PYTEST_HEADER_MODULES.pop('Pandas', None)
        PYTEST_HEADER_MODULES['Astropy'] = 'astropy'
        PYTEST_HEADER_MODULES['healpy'] = 'healpy'

        from . import __version__
        packagename = os.path.basename(os.path.dirname(__file__))
        TESTED_VERSIONS[packagename] = __version__

    # Set the Numpy print style to a fixed version to make doctest outputs
    # reproducible.
    try:
        np.set_printoptions(legacy='1.13')
    except TypeError:
        # On older versions of Numpy, the unrecognized 'legacy' option will
        # raise a TypeError.
        pass 

def test_str_repr_legacy(self):
        oldopts = np.get_printoptions()
        np.set_printoptions(legacy='1.13')
        try:
            a = array([0, 1, 2], mask=[False, True, False])
            assert_equal(str(a), '[0 -- 2]')
            assert_equal(repr(a), 'masked_array(data = [0 -- 2],\n'
                                  '             mask = [False  True False],\n'
                                  '       fill_value = 999999)\n')

            a = np.ma.arange(2000)
            a[1:50] = np.ma.masked
            assert_equal(
                repr(a),
                'masked_array(data = [0 -- -- ..., 1997 1998 1999],\n'
                '             mask = [False  True  True ..., False False False],\n'
                '       fill_value = 999999)\n'
            )
        finally:
            np.set_printoptions(**oldopts) 

def test_formatter_reset(self):
        x = np.arange(3)
        np.set_printoptions(formatter={'all':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'int':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        np.set_printoptions(formatter={'all':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'all':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        np.set_printoptions(formatter={'int':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'int_kind':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        x = np.arange(3.)
        np.set_printoptions(formatter={'float':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1.0, 0.0, 1.0])")
        np.set_printoptions(formatter={'float_kind':None})
        assert_equal(repr(x), "array([0., 1., 2.])") 

def test_float_spacing(self):
        x = np.array([1., 2., 3.])
        y = np.array([1., 2., -10.])
        z = np.array([100., 2., -1.])
        w = np.array([-100., 2., 1.])

        assert_equal(repr(x), 'array([1., 2., 3.])')
        assert_equal(repr(y), 'array([  1.,   2., -10.])')
        assert_equal(repr(np.array(y[0])), 'array(1.)')
        assert_equal(repr(np.array(y[-1])), 'array(-10.)')
        assert_equal(repr(z), 'array([100.,   2.,  -1.])')
        assert_equal(repr(w), 'array([-100.,    2.,    1.])')

        assert_equal(repr(np.array([np.nan, np.inf])), 'array([nan, inf])')
        assert_equal(repr(np.array([np.nan, -np.inf])), 'array([ nan, -inf])')

        x = np.array([np.inf, 100000, 1.1234])
        y = np.array([np.inf, 100000, -1.1234])
        z = np.array([np.inf, 1.1234, -1e120])
        np.set_printoptions(precision=2)
        assert_equal(repr(x), 'array([     inf, 1.00e+05, 1.12e+00])')
        assert_equal(repr(y), 'array([      inf,  1.00e+05, -1.12e+00])')
        assert_equal(repr(z), 'array([       inf,  1.12e+000, -1.00e+120])') 

def test_linewidth_str(self):
        a = np.full(18, fill_value=2)
        np.set_printoptions(linewidth=18)
        assert_equal(
            str(a),
            textwrap.dedent("""\
            [2 2 2 2 2 2 2 2
             2 2 2 2 2 2 2 2
             2 2]""")
        )
        np.set_printoptions(linewidth=18, legacy='1.13')
        assert_equal(
            str(a),
            textwrap.dedent("""\
            [2 2 2 2 2 2 2 2 2
             2 2 2 2 2 2 2 2 2]""")
        ) 

def __init__(self):
		# NEST options
		np.set_printoptions(precision=1)
		nest.set_verbosity('M_WARNING')
		nest.ResetKernel()
		nest.SetKernelStatus({"local_num_threads" : 1, "resolution" : p.time_resolution})
		# Create Poisson neurons
		self.spike_generators = nest.Create("poisson_generator", p.resolution[0]*p.resolution[1], params=p.poisson_params)
		self.neuron_pre = nest.Create("parrot_neuron", p.resolution[0]*p.resolution[1])
		# Create motor IAF neurons
		self.neuron_post = nest.Create("iaf_psc_alpha", 2, params=p.iaf_params)
		# Create Output spike detector
		self.spike_detector = nest.Create("spike_detector", 2, params={"withtime": True})
		# Create R-STDP synapses
		self.syn_dict = {"model": "stdp_dopamine_synapse",
						"weight": {"distribution": "uniform", "low": p.w0_min, "high": p.w0_max}}
		self.vt = nest.Create("volume_transmitter")
		nest.SetDefaults("stdp_dopamine_synapse", {"vt": self.vt[0], "tau_c": p.tau_c, "tau_n": p.tau_n, "Wmin": p.w_min, "Wmax": p.w_max, "A_plus": p.A_plus, "A_minus": p.A_minus})
		nest.Connect(self.spike_generators, self.neuron_pre, "one_to_one")
		nest.Connect(self.neuron_pre, self.neuron_post, "all_to_all", syn_spec=self.syn_dict)
		nest.Connect(self.neuron_post, self.spike_detector, "one_to_one")
		# Create connection handles for left and right motor neuron
		self.conn_l = nest.GetConnections(target=[self.neuron_post[0]])
		self.conn_r = nest.GetConnections(target=[self.neuron_post[1]]) 

def test_formatter_reset(self):
        x = np.arange(3)
        np.set_printoptions(formatter={'all':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'int':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        np.set_printoptions(formatter={'all':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'all':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        np.set_printoptions(formatter={'int':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1, 0, 1])")
        np.set_printoptions(formatter={'int_kind':None})
        assert_equal(repr(x), "array([0, 1, 2])")

        x = np.arange(3.)
        np.set_printoptions(formatter={'float':lambda x: str(x-1)})
        assert_equal(repr(x), "array([-1.0, 0.0, 1.0])")
        np.set_printoptions(formatter={'float_kind':None})
        assert_equal(repr(x), "array([ 0.,  1.,  2.])") 

def test_str_repr_legacy(self):
        oldopts = np.get_printoptions()
        np.set_printoptions(legacy='1.13')
        try:
            a = array([0, 1, 2], mask=[False, True, False])
            assert_equal(str(a), '[0 -- 2]')
            assert_equal(repr(a), 'masked_array(data = [0 -- 2],\n'
                                  '             mask = [False  True False],\n'
                                  '       fill_value = 999999)\n')

            a = np.ma.arange(2000)
            a[1:50] = np.ma.masked
            assert_equal(
                repr(a),
                'masked_array(data = [0 -- -- ..., 1997 1998 1999],\n'
                '             mask = [False  True  True ..., False False False],\n'
                '       fill_value = 999999)\n'
            )
        finally:
            np.set_printoptions(**oldopts) 

def shear_test():
    model=FEModel()
    n1=Node(0,0,0)
    n2=Node(0,0,5)
    n3=Node(5,0,5)
    n4=Node(10,0,5)
    n5=Node(10,0,0)
    
    a1=Membrane3(n1,n2,n3,0.25,2e11,0.3,7849)
    a2=Membrane3(n3,n4,n5,0.25,2e11,0.3,7849)
    a3=Membrane3(n1,n3,n5,0.25,2e11,0.3,7849)
    
    model.add_node(n1)
    model.add_node(n2)
    model.add_node(n3)
    model.add_node(n4)
    model.add_node(n5)

    model.add_membrane3(a1)
    model.add_membrane3(a2)
    model.add_membrane3(a3)

    n3.fn=(0,0,-100000,0,0,0)
    n1.dn=n2.dn=n4.dn=n5.dn=[0,0,0,0,0,0]

    model.assemble_KM()
    model.assemble_f()
    model.assemble_boundary()
    res=solve_linear(model)
        
    np.set_printoptions(precision=6,suppress=True)
    print(res)
    print(r"correct answer should be about ???") 

def pseudo_cantilever_test(l=25,h=5):
    """
    This is a cantilever beam with 50x10
    l,h: division on l and h direction
    """
    model=FEModel()
    nodes=[]
    for i in range(h+1):
        for j in range(l+1):
            nodes.append(Node(j*50/l,0,i*10/h))
            model.add_node(nodes[-1])
            
    for i in range(h):
        for j in range(l):
            area1=Membrane3(nodes[i*(l+1)+j],
                           nodes[i*(l+1)+j+1],
                           nodes[(i+1)*(l+1)+j+1],
                           0.25,2e11,0.3,7849)
            area2=Membrane3(nodes[i*(l+1)+j],
                           nodes[(i+1)*(l+1)+j+1],
                           nodes[(i+1)*(l+1)+j],
                           0.25,2e11,0.3,7849)
            if j==0:
                nodes[i*(l+1)+j].dn=[0]*6
                nodes[(i+1)*(l+1)+j].dn=[0]*6

            model.add_membrane3(area1)
            model.add_membrane3(area2)

    nodes[(l+1)*(h+1)-1].fn=(0,0,-100000,0,0,0)
    
    model.assemble_KM()
    model.assemble_f()
    model.assemble_boundary()
    res=solve_linear(model)

    np.set_printoptions(precision=6,suppress=True)
    print(res[(l+1)*(h+1)*6-6:])
    print(r"correct answer should be ???") 

def shear_test4():
    model=FEModel()
    n1=Node(0,0,0)
    n2=Node(0,0,5)
    n3=Node(5,0,5)
    n4=Node(10,0,5)
    n5=Node(10,0,0)
    n6=Node(5,0,0)
    
    a1=Membrane4(n1,n2,n3,n6,0.25,2e11,0.3,7849)
    a2=Membrane4(n3,n4,n5,n6,0.25,2e11,0.3,7849)
    
    model.add_node(n1)
    model.add_node(n2)
    model.add_node(n3)
    model.add_node(n4)
    model.add_node(n5)
    model.add_node(n6)

    model.add_membrane4(a1)
    model.add_membrane4(a2)

    n3.fn=(0,0,-100000,0,0,0)
    n1.dn=n2.dn=n4.dn=n5.dn=[0]*6
    n3.dn=n6.dn=[0,0,None,0,0,0]
    
    model.assemble_KM()
    model.assemble_f()
    model.assemble_boundary()
    res=solve_linear(model)

    np.set_printoptions(precision=6,suppress=True)
    print(res)
    print(r"correct answer should be ???") 

def pseudo_cantilever_test4(l=25,h=5):
    """
    This is a cantilever beam with 50x10
    l,h: division on l and h direction
    """
    model=FEModel()
    nodes=[]
    model=FEModel()
    nodes=[]
    for i in range(h+1):
        for j in range(l+1):
            nodes.append(Node(j*50/l,0,i*10/h))
            model.add_node(nodes[-1])
            
    for i in range(h):
        for j in range(l):
            area=Membrane4(nodes[i*(l+1)+j],
                           nodes[i*(l+1)+j+1],
                           nodes[(i+1)*(l+1)+j+1],
                           nodes[(i+1)*(l+1)+j+1],
                           0.25,2e11,0.3,7849)
            if j==0:
                nodes[i*(l+1)+j].dn=[0]*6
                nodes[(i+1)*(l+1)+j].dn=[0]*6

            model.add_membrane4(area)

    nodes[(l+1)*(h+1)-1].fn=(0,0,-100000,0,0,0)
    
    model.assemble_KM()
    model.assemble_f()
    model.assemble_boundary()
    res=solve_linear(model)

    np.set_printoptions(precision=6,suppress=True)
    print(res[(l+1)*(h+1)*6-6:])
    print(r"correct answer should be ???") 

def position_camera(self, camera_xyz, lookat_xyz, up):
    camera_xyz = np.array(camera_xyz)
    lookat_xyz = np.array(lookat_xyz)
    up = np.array(up)
    lookat_to = lookat_xyz - camera_xyz
    lookat_from = np.array([0, 1., 0.])
    up_from = np.array([0, 0., 1.])
    up_to = up * 1.
    # np.set_printoptions(precision=2, suppress=True)
    # print up_from, lookat_from, up_to, lookat_to
    r = ru.rotate_camera_to_point_at(up_from, lookat_from, up_to, lookat_to)
    R = np.eye(4, dtype=np.float32)
    R[:3,:3] = r

    t = np.eye(4, dtype=np.float32)
    t[:3,3] = -camera_xyz

    view_matrix = np.dot(R.T, t)
    flip_yz = np.eye(4, dtype=np.float32)
    flip_yz[1,1] = 0; flip_yz[2,2] = 0; flip_yz[1,2] = 1; flip_yz[2,1] = -1;
    view_matrix = np.dot(flip_yz, view_matrix)
    view_matrix = view_matrix.T
    # print np.concatenate((R, t, view_matrix), axis=1)
    view_matrix = np.reshape(view_matrix, (-1))
    view_matrix_o = glGetUniformLocation(self.egl_program, 'uViewMatrix')
    glUniformMatrix4fv(view_matrix_o, 1, GL_FALSE, view_matrix)
    return None, None #camera_xyz, q 

def get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None):
    '''Label the occupancies for each orbital for sampled k-points.

    This is a k-point version of scf.hf.SCF.get_occ
    '''
    if mo_energy_kpts is None: mo_energy_kpts = mf.mo_energy

    nkpts = len(mo_energy_kpts)
    nocc = mf.cell.tot_electrons(nkpts) // 2

    mo_energy = np.sort(np.hstack(mo_energy_kpts))
    fermi = mo_energy[nocc-1]
    mo_occ_kpts = []
    for mo_e in mo_energy_kpts:
        mo_occ_kpts.append((mo_e <= fermi).astype(np.double) * 2)

    if nocc < mo_energy.size:
        logger.info(mf, 'HOMO = %.12g  LUMO = %.12g',
                    mo_energy[nocc-1], mo_energy[nocc])
        if mo_energy[nocc-1]+1e-3 > mo_energy[nocc]:
            logger.warn(mf, 'HOMO %.12g == LUMO %.12g',
                        mo_energy[nocc-1], mo_energy[nocc])
    else:
        logger.info(mf, 'HOMO = %.12g', mo_energy[nocc-1])

    if mf.verbose >= logger.DEBUG:
        np.set_printoptions(threshold=len(mo_energy))
        logger.debug(mf, '     k-point                  mo_energy')
        for k,kpt in enumerate(mf.cell.get_scaled_kpts(mf.kpts)):
            logger.debug(mf, '  %2d (%6.3f %6.3f %6.3f)   %s %s',
                         k, kpt[0], kpt[1], kpt[2],
                         mo_energy_kpts[k][mo_occ_kpts[k]> 0],
                         mo_energy_kpts[k][mo_occ_kpts[k]==0])
        np.set_printoptions(threshold=1000)

    return mo_occ_kpts 

def get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None):
    '''Label the occupancies for each orbital for sampled k-points.

    This is a k-point version of scf.hf.SCF.get_occ
    '''
    if mo_energy_kpts is None: mo_energy_kpts = mf.mo_energy

    nkpts = len(mo_energy_kpts)
    nocc = mf.cell.nelectron * nkpts

    mo_energy = np.sort(np.hstack(mo_energy_kpts))
    fermi = mo_energy[nocc-1]
    mo_occ_kpts = []
    for mo_e in mo_energy_kpts:
        mo_occ_kpts.append((mo_e <= fermi).astype(np.double))

    if nocc < mo_energy.size:
        logger.info(mf, 'HOMO = %.12g  LUMO = %.12g',
                    mo_energy[nocc-1], mo_energy[nocc])
        if mo_energy[nocc-1]+1e-3 > mo_energy[nocc]:
            logger.warn(mf, 'HOMO %.12g == LUMO %.12g',
                        mo_energy[nocc-1], mo_energy[nocc])
    else:
        logger.info(mf, 'HOMO = %.12g', mo_energy[nocc-1])

    if mf.verbose >= logger.DEBUG:
        np.set_printoptions(threshold=len(mo_energy))
        logger.debug(mf, '     k-point                  mo_energy')
        for k,kpt in enumerate(mf.cell.get_scaled_kpts(mf.kpts)):
            logger.debug(mf, '  %2d (%6.3f %6.3f %6.3f)   %s %s',
                         k, kpt[0], kpt[1], kpt[2],
                         mo_energy_kpts[k][mo_occ_kpts[k]> 0],
                         mo_energy_kpts[k][mo_occ_kpts[k]==0])
        np.set_printoptions(threshold=1000)

    return mo_occ_kpts 

def get_occ(mf, mo_energy=None, mo_coeff=None):
    '''Label the occupancies for each orbital

    Kwargs:
        mo_energy : 1D ndarray
            Obital energies

        mo_coeff : 2D ndarray
            Obital coefficients

    Examples:

    >>> from pyscf import gto, scf
    >>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
    >>> mf = scf.hf.SCF(mol)
    >>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
    >>> mf.get_occ(energy)
    array([2, 2, 0, 2, 2, 2])
    '''
    if mo_energy is None: mo_energy = mf.mo_energy
    e_idx = numpy.argsort(mo_energy)
    e_sort = mo_energy[e_idx]
    nmo = mo_energy.size
    mo_occ = numpy.zeros(nmo)
    nocc = mf.mol.nelectron // 2
    mo_occ[e_idx[:nocc]] = 2
    if mf.verbose >= logger.INFO and nocc < nmo:
        if e_sort[nocc-1]+1e-3 > e_sort[nocc]:
            logger.warn(mf, 'HOMO %.15g == LUMO %.15g',
                        e_sort[nocc-1], e_sort[nocc])
        else:
            logger.info(mf, '  HOMO = %.15g  LUMO = %.15g',
                        e_sort[nocc-1], e_sort[nocc])

    if mf.verbose >= logger.DEBUG:
        numpy.set_printoptions(threshold=nmo)
        logger.debug(mf, '  mo_energy =\n%s', mo_energy)
        numpy.set_printoptions(threshold=1000)
    return mo_occ 

def get_occ(mf, mo_energy=None, mo_coeff=None):
    if mo_energy is None: mo_energy = mf.mo_energy
    e_idx_a = numpy.argsort(mo_energy[0])
    e_idx_b = numpy.argsort(mo_energy[1])
    e_sort_a = mo_energy[0][e_idx_a]
    e_sort_b = mo_energy[1][e_idx_b]
    nmo = mo_energy[0].size
    n_a, n_b = mf.nelec
    mo_occ = numpy.zeros_like(mo_energy)
    mo_occ[0,e_idx_a[:n_a]] = 1
    mo_occ[1,e_idx_b[:n_b]] = 1
    if mf.verbose >= logger.INFO and n_a < nmo and n_b > 0 and n_b < nmo:
        if e_sort_a[n_a-1]+1e-3 > e_sort_a[n_a]:
            logger.warn(mf, 'alpha nocc = %d  HOMO %.15g >= LUMO %.15g',
                        n_a, e_sort_a[n_a-1], e_sort_a[n_a])
        else:
            logger.info(mf, '  alpha nocc = %d  HOMO = %.15g  LUMO = %.15g',
                        n_a, e_sort_a[n_a-1], e_sort_a[n_a])

        if e_sort_b[n_b-1]+1e-3 > e_sort_b[n_b]:
            logger.warn(mf, 'beta  nocc = %d  HOMO %.15g >= LUMO %.15g',
                        n_b, e_sort_b[n_b-1], e_sort_b[n_b])
        else:
            logger.info(mf, '  beta  nocc = %d  HOMO = %.15g  LUMO = %.15g',
                        n_b, e_sort_b[n_b-1], e_sort_b[n_b])

        if e_sort_a[n_a-1]+1e-3 > e_sort_b[n_b]:
            logger.warn(mf, 'system HOMO %.15g >= system LUMO %.15g',
                        e_sort_b[n_a-1], e_sort_b[n_b])

        numpy.set_printoptions(threshold=nmo)
        logger.debug(mf, '  alpha mo_energy =\n%s', mo_energy[0])
        logger.debug(mf, '  beta  mo_energy =\n%s', mo_energy[1])
        numpy.set_printoptions(threshold=1000)

    if mo_coeff is not None and mf.verbose >= logger.DEBUG:
        ss, s = mf.spin_square((mo_coeff[0][:,mo_occ[0]>0],
                                mo_coeff[1][:,mo_occ[1]>0]), mf.get_ovlp())
        logger.debug(mf, 'multiplicity <S^2> = %.8g  2S+1 = %.8g', ss, s)
    return mo_occ 

def _regular_step(heff, ovlp, xs, lindep, log):
    try:
        e, c = scipy.linalg.eigh(heff[1:,1:], ovlp[1:,1:])
    except scipy.linalg.LinAlgError:
        e, c = lib.safe_eigh(heff[1:,1:], ovlp[1:,1:], lindep)[:2]
    if numpy.any(e < -1e-5):
        log.debug('Negative hessians found %s', e[e<0])

    w, v, seig = lib.safe_eigh(heff, ovlp, lindep)
    if log.verbose >= logger.DEBUG3:
        numpy.set_printoptions(3, linewidth=1000)
        log.debug3('v[0] %s', v[0])
        log.debug3('AH eigs %s', w)
        numpy.set_printoptions(8, linewidth=75)

    #if e[0] < -.1:
    #    sel = 0
    #else:
    # There exists systems that the first eigenvalue of AH is -inf.
    # Dynamically choosing the eigenvectors may be better.
    idx = numpy.where(abs(v[0]) > 0.1)[0]
    sel = idx[0]
    log.debug1('CIAH eigen-sel %s', sel)

    w_t = w[sel]
    xtrial = _dgemv(v[1:,sel]/v[0,sel], xs)
    return xtrial, w_t, v[:,sel], sel, seig 

def polariz_nonin_ave_matelem(self, comega):
    from scipy.sparse import spmatrix 
    """ Computes the non-interacting optical polarizability via the dipole matrix elements."""

    x,y,z = map(spmatrix.toarray, self.dipole_coo())
    i2d = array((x,y,z))
    n = self.mo_occ.shape[-1]
    eemax = max(comega.real)+20.0*max(comega.imag)
    
    p = zeros((len(comega)), dtype=np.complex128) # result to accumulate

    #print(__name__, 'Fermi energy', self.fermi_energy)
    #np.set_printoptions(linewidth=1000)
    for s in range(self.nspin):
      o,e,cc = self.mo_occ[0,s],self.mo_energy[0,s],self.mo_coeff[0,s,:,:,0]
      #print(o[:10])
      #print(e[:10])

      oo1,ee1 = np.subtract.outer(o,o).reshape(n*n), np.subtract.outer(e,e).reshape(n*n)
      idx = unravel_index( np.intersect1d(where(oo1<0.0), where(ee1<eemax)), (n,n))
      ivrt,iocc = array(list(set(idx[0]))), array(list(set(idx[1])))
      voi2d = einsum('nia,ma->nmi', einsum('iab,nb->nia', i2d, cc[ivrt]), cc[iocc])
      t2osc = 2.0/3.0*einsum('voi,voi->vo', voi2d, voi2d)
      t2w =  np.subtract.outer(e[ivrt],e[iocc])
      t2o = -np.subtract.outer(o[ivrt],o[iocc])

      for iw,w in enumerate(comega):
        p[iw] += 0.5*(t2osc*((t2o/(w-t2w))-(t2o/(w+t2w)))).sum()
      
    return p 

def _main():
   # read a file, including orbitals and basis sets, and test
   # if the output orbitals are orthogonal.
   def rmsd(a):
      return np.mean(a.flatten()**2)**.5
   FileName = "benzene.xml"
   #FileName = "/home/cgk/dev/xml-molpro/test1.xml"
   XmlData = ReadMolproXml(FileName,SkipVirtual=True)
   print("Atoms from file [a.u.]:\n{}".format(XmlData.Atoms.MakeXyz(NumFmt="%20.15f",Scale=1/wmme.ToAng)))
   OrbBasis = XmlData.OrbBasis
   #BasisLibs = ["def2-nzvpp-jkfit.libmol"]


   BasisLibs = []
   ic = wmme.FIntegralContext(XmlData.Atoms, XmlData.OrbBasis, FitBasis="univ-JKFIT", BasisLibs=BasisLibs)
   from wmme import mdot
   C = XmlData.Orbs
   S = ic.MakeOverlap()
   print("Orbital matrix shape: {} (loaded from '{}')".format(C.shape, FileName))
   print("Overlap matrix shape: {} (made via WMME)".format(S.shape))
   np.set_printoptions(precision=4,linewidth=10000,edgeitems=3,suppress=False)
   SMo = mdot(C.T, S, C)
   print("Read orbitals:")
   for OrbInfo in XmlData.Orbitals:
      print("{:30s}".format(OrbInfo.Desc))
   print("MO deviation from orthogonality: {:.2e}".format(rmsd(SMo - np.eye(SMo.shape[0]))))

   pass 

def loginstant(self, rho, c_am, v_lm, fmat, jmat, kmat, tnow, it):
        """
        time is logged in atomic units.
        Args:
            rho: complex
                MO density matrix.
            c_am: complex
                Transformation Matrix |AO><MO|
            v_lm: complex
                Transformation Matrix |LAO><MO|
            fmat: complex
                Fock matrix in Lowdin AO basis
            jmat: complex
                Coulomb matrix in AO basis
            kmat: complex
                Exact Exchange in AO basis
            tnow: float
                Current time in propagation in A.U.
            it: int
                Number of iteration of propagation
        Returns:
            tore: str
                |t, dipole(x,y,z), energy|

        """
        np.set_printoptions(precision = 7)
        tore = str(tnow)+" "+str(self.dipole(rho, c_am).real).rstrip("]").lstrip("[")+\
         " " +str(self.energy(transmat(rho,v_lm,-1),fmat, jmat, kmat))

        if it%self.params["StatusEvery"] ==0 or it == self.params["MaxIter"]-1:
            logger.log(self, "t: %f fs    Energy: %f a.u.   Total Density: %f",\
            tnow*FSPERAU,self.energy(transmat(rho,v_lm,-1),fmat, jmat, kmat), \
            2*np.trace(rho))
            logger.log(self, "Dipole moment(X, Y, Z, au): %8.5f, %8.5f, %8.5f",\
             self.dipole(rho, c_am).real[0],self.dipole(rho, c_am).real[1],\
             self.dipole(rho, c_am).real[2])
        return tore 

def run_simnibs(simnibs_struct, cpus=1):
    """Runs a simnnibs problem.

    Parameters:
    --------------
    simnibs_struct: sim_struct.mat.SESSION of str
        SESSION of name of '.mat' file defining the simulation
    cpus: int
        Number of processes to run in parallel (if avaliable)
    """
    np.set_printoptions(precision=4)

    if isinstance(simnibs_struct, str):
        p = read_mat(simnibs_struct)
    else:
        p = simnibs_struct

    out = p.run(cpus=cpus)
    logging.shutdown()
    return out 

def main():
    args = parse_arguments(sys.argv[1:])
    m2m_dir = os.path.abspath(os.path.realpath(os.path.expanduser(args.m2mpath)))
    if not os.path.isdir(m2m_dir):
        raise IOError('Could not find directory: {0}'.format(args.m2mpath))
    if args.out is not None:
        fn_out = os.path.abspath(os.path.realpath(os.path.expanduser(args.out)))
        fn_geo = os.path.splitext(fn_out)[0] + '.geo'
    else:
        fn_out = None
        fn_geo = None

    if args.coords is not None:
        coords = [float(d) for d in args.coords]

    elif args.csv is not None:
        coords = os.path.abspath(os.path.realpath(os.path.expanduser(args.csv)))
        if not os.path.isfile(coords):
            raise IOError('Could not find CSV file: {0}'.format(args.csv))
    else:
        raise argparse.ArgumentTypeError(
            'Plase use either -c or -s')

    if args.coords is not None and args.csv is not None:
        raise argparse.ArgumentError(
            'Please use only -c or -s')

    coords = transformations.warp_coordinates(
        coords, m2m_dir,
        transformation_direction='mni2subject',
        out_name=fn_out,
        transformation_type=args.t,
        out_geo=fn_geo)[1]
    if fn_out is None:
        np.set_printoptions(precision=2,
                            formatter={'float': lambda x: '{0:.2f}'.format(x)})
        print('Transformed coodinates:\n{0}'.format(coords.squeeze())) 

def main():
    args = parse_arguments(sys.argv[1:])
    m2m_dir = os.path.abspath(os.path.realpath(os.path.expanduser(args.m2mpath)))
    if not os.path.isdir(m2m_dir):
        raise IOError('Could not find directory: {0}'.format(args.m2mpath))
    if args.out is not None:
        fn_out = os.path.abspath(os.path.realpath(os.path.expanduser(args.out)))
        fn_geo = os.path.splitext(fn_out)[0] + '.geo'
    else:
        fn_out = None
        fn_geo = None

    if args.coords is not None:
        coords = np.array([float(d) for d in args.coords])

    elif args.csv is not None:
        coords = os.path.abspath(os.path.realpath(os.path.expanduser(args.csv)))
        if not os.path.isfile(coords):
            raise IOError('Could not find CSV file: {0}'.format(args.csv))
    else:
        raise argparse.ArgumentTypeError(
            'Plase use either -c or -s')

    if args.coords is not None and args.csv is not None:
        raise argparse.ArgumentError(
            'Please use only -c or -s')

    coords = transformations.warp_coordinates(
        coords, m2m_dir,
        transformation_direction='subject2mni',
        out_name=fn_out,
        transformation_type=args.t,
        out_geo=fn_geo)[1]
    if fn_out is None:
        np.set_printoptions(precision=2,
                            formatter={'float': lambda x: '{0:.2f}'.format(x)})
        print('Transformed coodinates:\n{0}'.format(coords.squeeze())) 

def display_greedy_policy(self):
        # Display greedy policy:
        #  - rows are s_player
        #  - columns are s_dealer
        print("\nDisplaying greedy policy:")
        np.set_printoptions(precision=3)
        greedy_policy = np.zeros(self.state_dim, dtype=int)
        for s_player in range(self.state_dim[0]):
            for s_dealer in range(self.state_dim[1]):
                greedy_policy[s_player, s_dealer] = np.argmax(self.Q[s_player, s_dealer, :])
        print(greedy_policy)
        print() 

def afterContext(self):
        numpy.set_printoptions(**print_state)

    # Ignore NumPy-specific build files that shouldn't be searched for tests 

def get_printoptions():
    """
    Return the current print options.

    Returns
    -------
    print_opts : dict
        Dictionary of current print options with keys

          - precision : int
          - threshold : int
          - edgeitems : int
          - linewidth : int
          - suppress : bool
          - nanstr : str
          - infstr : str
          - formatter : dict of callables
          - sign : str

        For a full description of these options, see `set_printoptions`.

    See Also
    --------
    set_printoptions, set_string_function

    """
    return _format_options.copy() 

def printoptions(*args, **kwargs):
    """Context manager for setting print options.

    Set print options for the scope of the `with` block, and restore the old
    options at the end. See `set_printoptions` for the full description of
    available options.

    Examples
    --------

    >>> with np.printoptions(precision=2):
    ...     print(np.array([2.0])) / 3
    [0.67]

    The `as`-clause of the `with`-statement gives the current print options:

    >>> with np.printoptions(precision=2) as opts:
    ...      assert_equal(opts, np.get_printoptions())

    See Also
    --------
    set_printoptions, get_printoptions

    """
    opts = np.get_printoptions()
    try:
        np.set_printoptions(*args, **kwargs)
        yield np.get_printoptions()
    finally:
        np.set_printoptions(**opts)