numpy.empty

Here are the examples of the python api numpy.empty taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

200 Examples 7

Example 1

Project: fuel
Source File: test_image.py
View license
    def test_format_exceptions(self):
        estream = Random2DRotation(self.example_stream,
                                   which_sources=('source2',))
        bstream = Random2DRotation(self.batch_stream,
                                   which_sources=('source2',))
        assert_raises(ValueError, estream.transform_source_example,
                      numpy.empty((5, 6)), 'source2')
        assert_raises(ValueError, bstream.transform_source_batch,
                      [numpy.empty((7, 6))], 'source2')
        assert_raises(ValueError, bstream.transform_source_batch,
                      [numpy.empty((8, 6))], 'source2')

Example 2

Project: chainer
Source File: transpose_sequence.py
View license
def _transpose(xs, length):
    if length == 0:
        return ()

    xp = cuda.get_array_module(*xs)
    lengths = numpy.empty(length, dtype='i')
    end = length
    for i, x in enumerate(xs):
        lengths[len(x):end] = i
        end = len(x)
    lengths[0:end] = len(xs)

    if xp is numpy:
        dtype = xs[0].dtype
        unit = xs[0].shape[1:]

        outs = tuple([xp.empty((l,) + unit, dtype=dtype) for l in lengths])
        for i, x in enumerate(xs):
            for p, xi in enumerate(x):
                outs[p][i] = xi

    else:
        offsets1 = numpy.empty(len(xs) + 1, dtype='i')
        offsets1[0] = 0
        numpy.cumsum([len(x) for x in xs], out=offsets1[1:])

        offsets2 = numpy.empty(length + 1, dtype='i')
        offsets2[0] = 0
        numpy.cumsum(lengths, dtype='i', out=offsets2[1:])

        x = xp.concatenate(xs, axis=0)
        o = xp.empty_like(x)
        unit = xs[0].size // len(xs[0])
        size = length * len(xs) * unit
        cuda.elementwise(
            'int32 len, int32 unit, raw int32 off1, raw int32 off2, raw T vs',
            'raw T hs',
            '''
            int ind = i / unit;
            int off = i - ind * unit;
            int y = ind / len;
            int x = ind - y * len;
            if (off2[x] + y < off2[x + 1]) {
              hs[(off2[x] + y) * unit + off] = vs[(off1[y] + x) * unit + off];
            }
            ''',
            'transpose_sequence'
        )(length, unit, cuda.to_gpu(offsets1), cuda.to_gpu(offsets2), x, o,
          size=size)
        outs = tuple(xp.split(o, offsets2[1:-1]))

    return outs

Example 3

Project: LV_groundhog
Source File: mainLoop.py
View license
    def roll_vocab_large2small(self):
        # Transfer from large to small parameters
        logger.debug("Called roll_vocab_large2small()")

        self.state['n_sym_source'] = len(self.model.small2large_src)
        logger.debug("n_sym_source=%d" % self.state['n_sym_source'])
        temp = numpy.empty((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32')
        if not self.state['fixed_embeddings']:
            temp_g2 = numpy.empty((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32')
            temp_d2 = numpy.empty((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32')
            if self.state['save_gs']:
                temp_gs = numpy.empty((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32')
        for small in self.model.small2large_src:
            large = self.model.small2large_src[small]
            temp[small] = self.model.large_W_0_enc_approx_embdr[large]
            if not self.state['fixed_embeddings']:
                temp_g2[small] = self.algo.large_W_0_enc_approx_embdr_g2[large]
                temp_d2[small] = self.algo.large_W_0_enc_approx_embdr_d2[large]
                if self.state['save_gs']:
                    temp_gs[small] = self.algo.large_W_0_enc_approx_embdr_gs[large]
        self.model.params[self.model.name2pos['W_0_enc_approx_embdr']].set_value(temp)
        if not self.state['fixed_embeddings']:
            self.algo.gnorm2[self.model.name2pos['W_0_enc_approx_embdr']].set_value(temp_g2)
            self.algo.dnorm2[self.model.name2pos['W_0_enc_approx_embdr']].set_value(temp_d2)
            if self.state['save_gs']:
                self.algo.gs[self.model.name2pos['W_0_enc_approx_embdr']].set_value(temp_gs)
        elif self.state['var_src_len']: # When embeddings are fixed, gnorm2 and dnorm2 must still have the right shape
            self.algo.gnorm2[self.model.name2pos['W_0_enc_approx_embdr']].set_value(numpy.zeros((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32'))
            self.algo.dnorm2[self.model.name2pos['W_0_enc_approx_embdr']].set_value(numpy.zeros((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32'))
            if self.state['save_gs']:
                self.algo.gs[self.model.name2pos['W_0_enc_approx_embdr']].set_value(numpy.zeros((self.state['n_sym_source'], self.state['rank_n_approx']), dtype='float32'))

        temp = self.model.params[self.model.name2pos['W_0_dec_approx_embdr']].get_value()
        if not self.state['fixed_embeddings']:
            temp_g2 = self.algo.gnorm2[self.model.name2pos['W_0_dec_approx_embdr']].get_value()
            temp_d2 = self.algo.dnorm2[self.model.name2pos['W_0_dec_approx_embdr']].get_value()
            if self.state['save_gs']:
                temp_gs = self.algo.gs[self.model.name2pos['W_0_dec_approx_embdr']].get_value()
        for small in self.model.small2large_trgt:
            large = self.model.small2large_trgt[small]
            temp[small] = self.model.large_W_0_dec_approx_embdr[large]
            if not self.state['fixed_embeddings']:
                temp_g2[small] = self.algo.large_W_0_dec_approx_embdr_g2[large]
                temp_d2[small] = self.algo.large_W_0_dec_approx_embdr_d2[large]
                if self.state['save_gs']:
                    temp_gs[small] = self.algo.large_W_0_dec_approx_embdr_gs[large]
        self.model.params[self.model.name2pos['W_0_dec_approx_embdr']].set_value(temp)
        if not self.state['fixed_embeddings']:
            self.algo.gnorm2[self.model.name2pos['W_0_dec_approx_embdr']].set_value(temp_g2)
            self.algo.dnorm2[self.model.name2pos['W_0_dec_approx_embdr']].set_value(temp_d2)
            if self.state['save_gs']:
                self.algo.gs[self.model.name2pos['W_0_dec_approx_embdr']].set_value(temp_gs)

        temp = self.model.params[self.model.name2pos['W2_dec_deep_softmax']].get_value()
        if not self.state['fixed_embeddings']:
            temp_g2 = self.algo.gnorm2[self.model.name2pos['W2_dec_deep_softmax']].get_value()
            temp_d2 = self.algo.dnorm2[self.model.name2pos['W2_dec_deep_softmax']].get_value()
            if self.state['save_gs']:
                temp_gs = self.algo.gs[self.model.name2pos['W2_dec_deep_softmax']].get_value()
        for small in self.model.small2large_trgt:
            large = self.model.small2large_trgt[small]
            temp[:,small] = self.model.large_W2_dec_deep_softmax[:,large]
            if not self.state['fixed_embeddings']:
                temp_g2[:,small] = self.algo.large_W2_dec_deep_softmax_g2[:,large]
                temp_d2[:,small] = self.algo.large_W2_dec_deep_softmax_d2[:,large]
                if self.state['save_gs']:
                    temp_gs[:,small] = self.algo.large_W2_dec_deep_softmax_gs[:,large]
        self.model.params[self.model.name2pos['W2_dec_deep_softmax']].set_value(temp)
        if not self.state['fixed_embeddings']:
            self.algo.gnorm2[self.model.name2pos['W2_dec_deep_softmax']].set_value(temp_g2)
            self.algo.dnorm2[self.model.name2pos['W2_dec_deep_softmax']].set_value(temp_d2)
            if self.state['save_gs']:
                self.algo.gs[self.model.name2pos['W2_dec_deep_softmax']].set_value(temp_gs)

        temp = self.model.params[self.model.name2pos['b_dec_deep_softmax']].get_value()
        if not self.state['fixed_embeddings']:
            temp_g2 = self.algo.gnorm2[self.model.name2pos['b_dec_deep_softmax']].get_value()
            temp_d2 = self.algo.dnorm2[self.model.name2pos['b_dec_deep_softmax']].get_value()
            if self.state['save_gs']:
                temp_gs = self.algo.gs[self.model.name2pos['b_dec_deep_softmax']].get_value()
        for small in self.model.small2large_trgt:
            large = self.model.small2large_trgt[small]
            temp[small] = self.model.large_b_dec_deep_softmax[large]
            if not self.state['fixed_embeddings']:
                temp_g2[small] = self.algo.large_b_dec_deep_softmax_g2[large]
                temp_d2[small] = self.algo.large_b_dec_deep_softmax_d2[large]
                if self.state['save_gs']:
                    temp_gs[small] = self.algo.large_b_dec_deep_softmax_gs[large]
        self.model.params[self.model.name2pos['b_dec_deep_softmax']].set_value(temp)
        if not self.state['fixed_embeddings']:
            self.algo.gnorm2[self.model.name2pos['b_dec_deep_softmax']].set_value(temp_g2)
            self.algo.dnorm2[self.model.name2pos['b_dec_deep_softmax']].set_value(temp_d2)
            if self.state['save_gs']:
                self.algo.gs[self.model.name2pos['b_dec_deep_softmax']].set_value(temp_gs)

Example 4

Project: pyscf
Source File: outcore.py
View license
def general(mol, mo_coeffs, erifile, dataname='eri_mo', tmpdir=None,
            intor='cint2e_sph', aosym='s4', comp=1,
            max_memory=2000, ioblk_size=IOBLK_SIZE, verbose=logger.WARN, compact=True):
    r'''For the given four sets of orbitals, transfer arbitrary spherical AO
    integrals to MO integrals on the fly.

    Args:
        mol : :class:`Mole` object
            AO integrals will be generated in terms of mol._atm, mol._bas, mol._env
        mo_coeffs : 4-item list of ndarray
            Four sets of orbital coefficients, corresponding to the four
            indices of (ij|kl)
        erifile : str or h5py File or h5py Group object
            To store the transformed integrals, in HDF5 format.

    Kwargs
        dataname : str
            The dataset name in the erifile (ref the hierarchy of HDF5 format
            http://www.hdfgroup.org/HDF5/doc1.6/UG/09_Groups.html).  By assigning
            different dataname, the existed integral file can be reused.  If
            the erifile contains the dataname, the new integrals data will
            overwrite the old one.
        tmpdir : str
            The directory where to temporarily store the intermediate data
            (the half-transformed integrals).  By default, it's controlled by
            shell environment variable ``TMPDIR``.  The disk space requirement
            is about  comp*mo_coeffs[0].shape[1]*mo_coeffs[1].shape[1]*nao**2
        intor : str
            Name of the 2-electron integral.  Ref to :func:`getints_by_shell`
            for the complete list of available 2-electron integral names
        aosym : int or str
            Permutation symmetry for the AO integrals

            | 4 or '4' or 's4': 4-fold symmetry (default)
            | '2ij' or 's2ij' : symmetry between i, j in (ij|kl)
            | '2kl' or 's2kl' : symmetry between k, l in (ij|kl)
            | 1 or '1' or 's1': no symmetry
            | 'a4ij' : 4-fold symmetry with anti-symmetry between i, j in (ij|kl) (TODO)
            | 'a4kl' : 4-fold symmetry with anti-symmetry between k, l in (ij|kl) (TODO)
            | 'a2ij' : anti-symmetry between i, j in (ij|kl) (TODO)
            | 'a2kl' : anti-symmetry between k, l in (ij|kl) (TODO)

        comp : int
            Components of the integrals, e.g. cint2e_ip_sph has 3 components.
        max_memory : float or int
            The maximum size of cache to use (in MB), large cache may **not**
            improve performance.
        ioblk_size : float or int
            The block size for IO, large block size may **not** improve performance
        verbose : int
            Print level
        compact : bool
            When compact is True, depending on the four oribital sets, the
            returned MO integrals has (up to 4-fold) permutation symmetry.
            If it's False, the function will abandon any permutation symmetry,
            and return the "plain" MO integrals

    Returns:
        None

    Examples:

    >>> from pyscf import gto
    >>> from pyscf import ao2mo
    >>> import h5py
    >>> def view(h5file, dataname='eri_mo'):
    ...     f5 = h5py.File(h5file)
    ...     print('dataset %s, shape %s' % (str(f5.keys()), str(f5[dataname].shape)))
    ...     f5.close()
    >>> mol = gto.M(atom='O 0 0 0; H 0 1 0; H 0 0 1', basis='sto3g')
    >>> mo1 = numpy.random.random((mol.nao_nr(), 10))
    >>> mo2 = numpy.random.random((mol.nao_nr(), 8))
    >>> mo3 = numpy.random.random((mol.nao_nr(), 6))
    >>> mo4 = numpy.random.random((mol.nao_nr(), 4))
    >>> ao2mo.outcore.general(mol, (mo1,mo2,mo3,mo4), 'oh2.h5')
    >>> view('oh2.h5')
    dataset ['eri_mo'], shape (80, 24)
    >>> ao2mo.outcore.general(mol, (mo1,mo2,mo3,mo3), 'oh2.h5')
    >>> view('oh2.h5')
    dataset ['eri_mo'], shape (80, 21)
    >>> ao2mo.outcore.general(mol, (mo1,mo2,mo3,mo3), 'oh2.h5', compact=False)
    >>> view('oh2.h5')
    dataset ['eri_mo'], shape (80, 36)
    >>> ao2mo.outcore.general(mol, (mo1,mo1,mo2,mo2), 'oh2.h5')
    >>> view('oh2.h5')
    dataset ['eri_mo'], shape (55, 36)
    >>> ao2mo.outcore.general(mol, (mo1,mo1,mo1,mo1), 'oh2.h5', dataname='new')
    >>> view('oh2.h5', 'new')
    dataset ['eri_mo', 'new'], shape (55, 55)
    >>> ao2mo.outcore.general(mol, (mo1,mo1,mo1,mo1), 'oh2.h5', intor='cint2e_ip1_sph', aosym='s1', comp=3)
    >>> view('oh2.h5')
    dataset ['eri_mo', 'new'], shape (3, 100, 100)
    >>> ao2mo.outcore.general(mol, (mo1,mo1,mo1,mo1), 'oh2.h5', intor='cint2e_ip1_sph', aosym='s2kl', comp=3)
    >>> view('oh2.h5')
    dataset ['eri_mo', 'new'], shape (3, 100, 55)
    '''
    time_0pass = (time.clock(), time.time())
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)

    nmoi = mo_coeffs[0].shape[1]
    nmoj = mo_coeffs[1].shape[1]
    nmok = mo_coeffs[2].shape[1]
    nmol = mo_coeffs[3].shape[1]
    nao = mo_coeffs[0].shape[0]
    aosym = _stand_sym_code(aosym)
    if aosym in ('s4', 's2kl'):
        nao_pair = nao * (nao+1) // 2
    else:
        nao_pair = nao * nao

    if (compact and iden_coeffs(mo_coeffs[0], mo_coeffs[1]) and
        aosym in ('s4', 's2ij')):
        nij_pair = nmoi*(nmoi+1) // 2
    else:
        nij_pair = nmoi*nmoj

    klmosym, nkl_pair, mokl, klshape = \
            incore._conc_mos(mo_coeffs[2], mo_coeffs[3],
                             compact and aosym in ('s4', 's2kl'))

#    if nij_pair > nkl_pair:
#        log.warn('low efficiency for AO to MO trans!')

    if isinstance(erifile, str):
        if h5py.is_hdf5(erifile):
            feri = h5py.File(erifile)
            if dataname in feri:
                del(feri[dataname])
        else:
            feri = h5py.File(erifile, 'w')
    else:
        assert(isinstance(erifile, h5py.Group))
        feri = erifile
    if comp == 1:
        chunks = (nmoj,nmol)
        h5d_eri = feri.create_dataset(dataname, (nij_pair,nkl_pair),
                                      'f8', chunks=chunks)
    else:
        chunks = (1,nmoj,nmol)
        h5d_eri = feri.create_dataset(dataname, (comp,nij_pair,nkl_pair),
                                      'f8', chunks=chunks)

    if nij_pair == 0 or nkl_pair == 0:
        if isinstance(erifile, str):
            feri.close()
        return erifile
    log.debug('MO integrals %s are saved in %s/%s', intor, erifile, dataname)
    log.debug('num. MO ints = %.8g, required disk %.8g MB',
              float(nij_pair)*nkl_pair*comp, nij_pair*nkl_pair*comp*8/1e6)

# transform e1
    swapfile = tempfile.NamedTemporaryFile(dir=tmpdir)
    fswap = h5py.File(swapfile.name, 'w')
    half_e1(mol, mo_coeffs, fswap, intor, aosym, comp, max_memory, ioblk_size,
            log, compact)

    time_1pass = log.timer('AO->MO transformation for %s 1 pass'%intor,
                           *time_0pass)

    ioblk_size = max(max_memory*.2, ioblk_size)
    iobuflen = guess_e2bufsize(ioblk_size, nij_pair, max(nao_pair,nkl_pair))[0]
    reading_frame = [numpy.empty((iobuflen,nao_pair)),
                     numpy.empty((iobuflen,nao_pair))]
    def prefetch(icomp, row0, row1, buf):
        if icomp+1 < comp:
            icomp += 1
        else:
            row0, row1 = row1, min(nij_pair, row1+iobuflen)
            icomp = 0
        if row0 < row1:
            _load_from_h5g(fswap['%d'%icomp], row0, row1, buf)
    def async_read(icomp, row0, row1, thread_read):
        buf_current, buf_prefetch = reading_frame
        reading_frame[:] = [buf_prefetch, buf_current]
        if thread_read is None:
            _load_from_h5g(fswap['%d'%icomp], row0, row1, buf_current)
        else:
            thread_read.join()
        thread_read = lib.background_thread(prefetch, icomp, row0, row1, buf_prefetch)
        return buf_current[:row1-row0], thread_read

    def save(icomp, row0, row1, buf):
        if comp == 1:
            h5d_eri[row0:row1] = buf[:row1-row0]
        else:
            h5d_eri[icomp,row0:row1] = buf[:row1-row0]
    def async_write(icomp, row0, row1, buf, thread_io):
        if thread_io is not None:
            thread_io.join()
        thread_io = lib.background_thread(save, icomp, row0, row1, buf)
        return thread_io

    log.debug('step2: kl-pair (ao %d, mo %d), mem %.8g MB, ioblock %.8g MB',
              nao_pair, nkl_pair, iobuflen*nao_pair*8/1e6,
              iobuflen*nkl_pair*8/1e6)

    klaoblks = len(fswap['0'])
    ijmoblks = int(numpy.ceil(float(nij_pair)/iobuflen)) * comp
    ao_loc = mol.ao_loc_nr('cart' in intor)
    ti0 = time_1pass
    bufs1 = numpy.empty((iobuflen,nkl_pair))
    buf_write = numpy.empty_like(bufs1)
    istep = 0
    read_handler = write_handler = None

    for row0, row1 in prange(0, nij_pair, iobuflen):
        nrow = row1 - row0

        for icomp in range(comp):
            istep += 1
            log.debug1('step 2 [%d/%d], [%d,%d:%d], row = %d',
                       istep, ijmoblks, icomp, row0, row1, nrow)

            buf, read_handler = async_read(icomp, row0, row1, read_handler)
            _ao2mo.nr_e2(buf, mokl, klshape, aosym, klmosym,
                         ao_loc=ao_loc, out=bufs1)
            write_handler = async_write(icomp, row0, row1, bufs1, write_handler)
            bufs1, buf_write = buf_write, bufs1  # avoid flushing writing buffer

            ti1 = (time.clock(), time.time())
            log.debug1('step 2 [%d/%d] CPU time: %9.2f, Wall time: %9.2f',
                       istep, ijmoblks, ti1[0]-ti0[0], ti1[1]-ti0[1])
            ti0 = ti1
    write_handler.join()
    fswap.close()
    if isinstance(erifile, str):
        feri.close()

    log.timer('AO->MO transformation for %s 2 pass'%intor, *time_1pass)
    log.timer('AO->MO transformation for %s '%intor, *time_0pass)
    return erifile

Example 5

Project: pyscf
Source File: outcore.py
View license
def half_e1(mol, mo_coeffs, swapfile,
            intor='cint2e_sph', aosym='s4', comp=1,
            max_memory=2000, ioblk_size=IOBLK_SIZE, verbose=logger.WARN, compact=True,
            ao2mopt=None):
    r'''Half transform arbitrary spherical AO integrals to MO integrals
    for the given two sets of orbitals

    Args:
        mol : :class:`Mole` object
            AO integrals will be generated in terms of mol._atm, mol._bas, mol._env
        mo_coeff : ndarray
            Transform (ij|kl) with the same set of orbitals.
        swapfile : str or h5py File or h5py Group object
            To store the transformed integrals, in HDF5 format.  The transformed
            integrals are saved in blocks.

    Kwargs
        intor : str
            Name of the 2-electron integral.  Ref to :func:`getints_by_shell`
            for the complete list of available 2-electron integral names
        aosym : int or str
            Permutation symmetry for the AO integrals

            | 4 or '4' or 's4': 4-fold symmetry (default)
            | '2ij' or 's2ij' : symmetry between i, j in (ij|kl)
            | '2kl' or 's2kl' : symmetry between k, l in (ij|kl)
            | 1 or '1' or 's1': no symmetry
            | 'a4ij' : 4-fold symmetry with anti-symmetry between i, j in (ij|kl) (TODO)
            | 'a4kl' : 4-fold symmetry with anti-symmetry between k, l in (ij|kl) (TODO)
            | 'a2ij' : anti-symmetry between i, j in (ij|kl) (TODO)
            | 'a2kl' : anti-symmetry between k, l in (ij|kl) (TODO)

        comp : int
            Components of the integrals, e.g. cint2e_ip_sph has 3 components.
        verbose : int
            Print level
        max_memory : float or int
            The maximum size of cache to use (in MB), large cache may **not**
            improve performance.
        ioblk_size : float or int
            The block size for IO, large block size may **not** improve performance
        verbose : int
            Print level
        compact : bool
            When compact is True, depending on the four oribital sets, the
            returned MO integrals has (up to 4-fold) permutation symmetry.
            If it's False, the function will abandon any permutation symmetry,
            and return the "plain" MO integrals
        ao2mopt : :class:`AO2MOpt` object
            Precomputed data to improve perfomance

    Returns:
        None

    '''
    time0 = (time.clock(), time.time())
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)

    nao = mo_coeffs[0].shape[0]
    aosym = _stand_sym_code(aosym)
    if aosym in ('s4', 's2ij'):
        nao_pair = nao * (nao+1) // 2
    else:
        nao_pair = nao * nao

    ijmosym, nij_pair, moij, ijshape = \
            incore._conc_mos(mo_coeffs[0], mo_coeffs[1],
                             compact and aosym in ('s4', 's2ij'))

    e1buflen, mem_words, iobuf_words, ioblk_words = \
            guess_e1bufsize(max_memory, ioblk_size, nij_pair, nao_pair, comp)
    ioblk_size = ioblk_words * 8/1e6
# The buffer to hold AO integrals in C code, see line (@)
    aobuflen = max(int((mem_words - 2*comp*e1buflen*nij_pair) // (nao_pair*comp)),
                   IOBUF_ROW_MIN)
    shranges = guess_shell_ranges(mol, (aosym in ('s4', 's2kl')), e1buflen, aobuflen)
    if ao2mopt is None:
        if intor == 'cint2e_sph':
            ao2mopt = _ao2mo.AO2MOpt(mol, intor, 'CVHFnr_schwarz_cond',
                                     'CVHFsetnr_direct_scf')
        else:
            ao2mopt = _ao2mo.AO2MOpt(mol, intor)

    if isinstance(swapfile, str):
        fswap = h5py.File(swapfile, 'w')
    else:
        fswap = swapfile
    for icomp in range(comp):
        g = fswap.create_group(str(icomp)) # for h5py old version

    log.debug('step1: tmpfile %s  %.8g MB', fswap.filename, nij_pair*nao_pair*8/1e6)
    log.debug('step1: (ij,kl) = (%d,%d), mem cache %.8g MB, iobuf %.8g MB',
              nij_pair, nao_pair, mem_words*8/1e6, iobuf_words*8/1e6)
    nstep = len(shranges)
    e1buflen = max([x[2] for x in shranges])

    e2buflen, chunks = guess_e2bufsize(ioblk_size, nij_pair, e1buflen)
    def save(istep, iobuf):
        for icomp in range(comp):
            _transpose_to_h5g(fswap, '%d/%d'%(icomp,istep), iobuf[icomp],
                              e2buflen, None)
    def async_write(istep, iobuf, thread_io):
        if thread_io is not None:
            thread_io.join()
        thread_io = lib.background_thread(save, istep, iobuf)
        return thread_io

    # transform e1
    ti0 = log.timer('Initializing ao2mo.outcore.half_e1', *time0)
    bufs1 = numpy.empty((comp*e1buflen,nao_pair))
    bufs2 = numpy.empty((comp*e1buflen,nij_pair))
    buf_write = numpy.empty_like(bufs2)
    write_handler = None
    for istep,sh_range in enumerate(shranges):
        log.debug1('step 1 [%d/%d], AO [%d:%d], len(buf) = %d', \
                   istep+1, nstep, *(sh_range[:3]))
        buflen = sh_range[2]
        iobuf = numpy.ndarray((comp,buflen,nij_pair), buffer=bufs2)
        nmic = len(sh_range[3])
        p0 = 0
        for imic, aoshs in enumerate(sh_range[3]):
            log.debug2('      fill iobuf micro [%d/%d], AO [%d:%d], len(aobuf) = %d', \
                       imic+1, nmic, *aoshs)
            buf = numpy.ndarray((comp*aoshs[2],nao_pair), buffer=bufs1) # (@)
            _ao2mo.nr_e1fill(intor, aoshs, mol._atm, mol._bas, mol._env,
                             aosym, comp, ao2mopt, out=buf)
            buf = _ao2mo.nr_e1(buf, moij, ijshape, aosym, ijmosym)
            iobuf[:,p0:p0+aoshs[2]] = buf.reshape(comp,aoshs[2],-1)
            p0 += aoshs[2]
        ti0 = log.timer_debug1('gen AO/transform MO [%d/%d]'%(istep+1,nstep), *ti0)

        write_handler = async_write(istep, iobuf, write_handler)
        bufs2, buf_write = buf_write, bufs2  # avoid flushing writing buffer
    write_handler.join()
    bufs1 = bufs2 = None
    if isinstance(swapfile, str):
        fswap.close()
    return swapfile

Example 6

Project: sphere
Source File: permeabilitycalculator.py
View license
    def plotEvolution(self, axis=2, outformat='png'):
        '''
        Plot temporal evolution of parameters on the selected axis.
        Note that the first 5 output files are ignored.
        '''
        skipsteps = 5
        nsteps = self.sim.status() - skipsteps
        self.t_series = numpy.empty(nsteps)
        self.Q_series = numpy.empty((nsteps, 3))
        self.phi_bar_series = numpy.empty(nsteps)
        self.k_series = numpy.empty((nsteps, 3))
        self.K_series = numpy.empty((nsteps, 3))

        print('Reading ' + str(nsteps) + ' output files... '),
        sys.stdout.flush()
        for i in numpy.arange(skipsteps, self.sim.status()):
            self.sim.readstep(i, verbose=False)

            self.t_series[i-skipsteps] = self.sim.time_current[0]

            self.findCrossSectionalFlux()
            self.Q_series[i-skipsteps,:] = self.Q

            self.findMeanPorosity()
            self.phi_bar_series[i-skipsteps] = self.phi_bar

            self.findPermeability()
            self.k_series[i-skipsteps,:] = self.k

            self.findConductivity()
            self.K_series[i-skipsteps,:] = self.K
        print('Done')

        fig = plt.figure()

        plt.subplot(2,2,1)
        plt.xlabel('Time $t$ [s]')
        plt.ylabel('Flux $Q$ [m^3/s]')
        plt.plot(self.t_series, self.Q_series[:,axis])
        #plt.legend()
        plt.grid()

        plt.subplot(2,2,2)
        plt.xlabel('Time $t$ [s]')
        plt.ylabel('Porosity $\phi$ [-]')
        plt.plot(self.t_series, self.phi_bar_series)
        plt.grid()

        plt.subplot(2,2,3)
        plt.xlabel('Time $t$ [s]')
        plt.ylabel('Permeability $k$ [m^2]')
        plt.plot(self.t_series, self.k_series[:,axis])
        plt.grid()

        plt.subplot(2,2,4)
        plt.xlabel('Time $t$ [s]')
        plt.ylabel('Conductivity $K$ [m/s]')
        plt.plot(self.t_series, self.K_series[:,axis])
        plt.grid()

        fig.tight_layout()

        filename = self.sid + '-permeability.' + outformat
        plt.savefig(filename)
        print('Figure saved as "' + filename + '"')

Example 7

Project: Arelle
Source File: UITkTable.py
View license
    def __init__(self, parentWidget, rows, columns, titleRows, titleColumns,
                 tableName=None, browsecmd=None):
        '''
        The initial size of the table (including the header sizes) must be
        supplied at table creation time.
        The contextual menu will have to be created later
        with a 'widget.contextMenu()' command.
        The Tab and Return key are bound to cell navigation.
        '''
        self.data = numpy.empty((rows, columns), dtype=object)
        self.objectIds = numpy.empty((rows, columns), dtype=object)
        self.maxColumnWidths = numpy.empty(columns, dtype=int)
        self.maxRowHeights = numpy.empty(rows, dtype=int)
        self.modifiedCells = dict()
        self.currentCellCoordinates = None
        self.data.fill('')
        self.objectIds.fill('')
        self.maxColumnWidths.fill(0)
        self.maxRowHeights.fill(0)
        self.titleRows = titleRows
        self.titleColumns = titleColumns

        if browsecmd is None:
            if USE_resizeTableCells:
                super(XbrlTable, self).__init__(parentWidget,
                                            rows=rows,
                                            cols=columns,
                                            state='normal',
                                            titlerows=titleRows,
                                            titlecols=titleColumns,
                                            roworigin=0,
                                            colorigin=0,
                                            selectmode='extended',
                                            selecttype='cell',
                                            rowstretch='none',
                                            colstretch='none',
                                            rowheight=self.DEFAULT_ROW_HEIGHT,
                                            colwidth=self.DEFAULT_COLUMN_WIDTH,
                                            flashmode='off',
                                            anchor=NE,
                                            usecommand=1,
                                            background='#d00d00d00',
                                            relief='ridge',
                                            command=self._valueCommand,
                                            takefocus=False,
                                            rowseparator='\n',
                                            wrap=1,
                                            borderwidth=(1, 1, 0, 0),
                                            multiline=1)
            else:
                super(XbrlTable, self).__init__(parentWidget,
                                            rows=rows,
                                            cols=columns,
                                            state='normal',
                                            titlerows=titleRows,
                                            titlecols=titleColumns,
                                            roworigin=0,
                                            colorigin=0,
                                            selectmode='extended',
                                            selecttype='cell',
                                            rowstretch='none',
                                            colstretch='none',
                                            rowheight=-30,
                                            colwidth=15,
                                            flashmode='off',
                                            anchor=NE,
                                            usecommand=1,
                                            background='#d00d00d00',
                                            relief='ridge',
                                            command=self._valueCommand,
                                            takefocus=False,
                                            rowseparator='\n',
                                            wrap=1,
                                            borderwidth=(1, 1, 0, 0),
                                            multiline=1)
        else:
            if USE_resizeTableCells:
                super(XbrlTable, self).__init__(parentWidget,
                                                rows=rows,
                                                cols=columns,
                                                state='normal',
                                                titlerows=titleRows,
                                                titlecols=titleColumns,
                                                roworigin=0,
                                                colorigin=0,
                                                selectmode='extended',
                                                selecttype='cell',
                                                rowstretch='none',
                                                colstretch='none',
                                                rowheight=self.DEFAULT_ROW_HEIGHT,
                                                colwidth=self.DEFAULT_COLUMN_WIDTH,
                                                flashmode='off',
                                                anchor=NE,
                                                usecommand=1,
                                                background='#d00d00d00',
                                                relief='ridge',
                                                command=self._valueCommand,
                                                takefocus=False,
                                                rowseparator='\n',
                                                wrap=1,
                                                multiline=1,
                                                borderwidth=(1, 1, 0, 0),
                                                browsecmd=browsecmd)
            else:
                super(XbrlTable, self).__init__(parentWidget,
                                                rows=rows,
                                                cols=columns,
                                                state='normal',
                                                titlerows=titleRows,
                                                titlecols=titleColumns,
                                                roworigin=0,
                                                colorigin=0,
                                                selectmode='extended',
                                                selecttype='cell',
                                                rowstretch='none',
                                                colstretch='none',
                                                rowheight=-30,
                                                colwidth=15,
                                                flashmode='off',
                                                anchor=NE,
                                                usecommand=1,
                                                background='#d00d00d00',
                                                relief='ridge',
                                                command=self._valueCommand,
                                                takefocus=False,
                                                rowseparator='\n',
                                                wrap=1,
                                                multiline=1,
                                                borderwidth=(1, 1, 0, 0),
                                                browsecmd=browsecmd)

        self.lastMouseCoordinates = None
        self.toolTipShown = False
        self.headerToolTipText = StringVar()
        self.headerToolTip = MyToolTip(self, textvariable=self.headerToolTipText, wraplength=480, follow_mouse=True, state="disabled")

        # Extra key bindings for navigating through the table:
        # Tab: go right
        # Return (and Enter): go down
        self.bind("<Tab>", func=self.cellRight)
        self.bind("<Return>", func=self.cellDown)
        if False:
            self.bind("<Motion>", func=self.mouseMotion)

        # Configure the graphical appearance of the cells by using tags
        self.tag_configure('sel', background = '#b00e00e60',
                           fg='#000000000')
        self.tag_configure('active', background = '#000f70ff0',
                           fg='#000000000')
        self.tag_configure('title', anchor='w', bg='#d00d00d00',
                           fg='#000000000', relief='ridge')
        self.tag_configure(XbrlTable.TG_DISABLED, bg='#d00d00d00',
                           fg='#000000000',
                           relief='flat', state='disabled')

        if rows+columns > 2:
            # The content of the left/top corner cell can already be defined
            topCell = '0,0'
            if titleColumns+titleRows-2 > 0:
                cellSpans = {topCell : '%i,%i'% (titleRows-1, titleColumns-1)}
                self.spans(index=None, **cellSpans)
            self.format_cell(XbrlTable.TG_TOP_LEFT, topCell)
            self.tag_raise(XbrlTable.TG_TOP_LEFT, abovethis='title')
            self.tag_configure(XbrlTable.TG_TOP_LEFT, anchor='ne')
            self.format_cell(XbrlTable.TG_BORDER_ALL, topCell)
            indexValue = {topCell:tableName}
            self.set(**indexValue)

Example 8

Project: Arelle
Source File: UITkTable.py
View license
    def __init__(self, parentWidget, rows, columns, titleRows, titleColumns,
                 tableName=None, browsecmd=None):
        '''
        The initial size of the table (including the header sizes) must be
        supplied at table creation time.
        The contextual menu will have to be created later
        with a 'widget.contextMenu()' command.
        The Tab and Return key are bound to cell navigation.
        '''
        self.data = numpy.empty((rows, columns), dtype=object)
        self.objectIds = numpy.empty((rows, columns), dtype=object)
        self.maxColumnWidths = numpy.empty(columns, dtype=int)
        self.maxRowHeights = numpy.empty(rows, dtype=int)
        self.modifiedCells = dict()
        self.currentCellCoordinates = None
        self.data.fill('')
        self.objectIds.fill('')
        self.maxColumnWidths.fill(0)
        self.maxRowHeights.fill(0)
        self.titleRows = titleRows
        self.titleColumns = titleColumns

        if browsecmd is None:
            if USE_resizeTableCells:
                super(XbrlTable, self).__init__(parentWidget,
                                            rows=rows,
                                            cols=columns,
                                            state='normal',
                                            titlerows=titleRows,
                                            titlecols=titleColumns,
                                            roworigin=0,
                                            colorigin=0,
                                            selectmode='extended',
                                            selecttype='cell',
                                            rowstretch='none',
                                            colstretch='none',
                                            rowheight=self.DEFAULT_ROW_HEIGHT,
                                            colwidth=self.DEFAULT_COLUMN_WIDTH,
                                            flashmode='off',
                                            anchor=NE,
                                            usecommand=1,
                                            background='#d00d00d00',
                                            relief='ridge',
                                            command=self._valueCommand,
                                            takefocus=False,
                                            rowseparator='\n',
                                            wrap=1,
                                            borderwidth=(1, 1, 0, 0),
                                            multiline=1)
            else:
                super(XbrlTable, self).__init__(parentWidget,
                                            rows=rows,
                                            cols=columns,
                                            state='normal',
                                            titlerows=titleRows,
                                            titlecols=titleColumns,
                                            roworigin=0,
                                            colorigin=0,
                                            selectmode='extended',
                                            selecttype='cell',
                                            rowstretch='none',
                                            colstretch='none',
                                            rowheight=-30,
                                            colwidth=15,
                                            flashmode='off',
                                            anchor=NE,
                                            usecommand=1,
                                            background='#d00d00d00',
                                            relief='ridge',
                                            command=self._valueCommand,
                                            takefocus=False,
                                            rowseparator='\n',
                                            wrap=1,
                                            borderwidth=(1, 1, 0, 0),
                                            multiline=1)
        else:
            if USE_resizeTableCells:
                super(XbrlTable, self).__init__(parentWidget,
                                                rows=rows,
                                                cols=columns,
                                                state='normal',
                                                titlerows=titleRows,
                                                titlecols=titleColumns,
                                                roworigin=0,
                                                colorigin=0,
                                                selectmode='extended',
                                                selecttype='cell',
                                                rowstretch='none',
                                                colstretch='none',
                                                rowheight=self.DEFAULT_ROW_HEIGHT,
                                                colwidth=self.DEFAULT_COLUMN_WIDTH,
                                                flashmode='off',
                                                anchor=NE,
                                                usecommand=1,
                                                background='#d00d00d00',
                                                relief='ridge',
                                                command=self._valueCommand,
                                                takefocus=False,
                                                rowseparator='\n',
                                                wrap=1,
                                                multiline=1,
                                                borderwidth=(1, 1, 0, 0),
                                                browsecmd=browsecmd)
            else:
                super(XbrlTable, self).__init__(parentWidget,
                                                rows=rows,
                                                cols=columns,
                                                state='normal',
                                                titlerows=titleRows,
                                                titlecols=titleColumns,
                                                roworigin=0,
                                                colorigin=0,
                                                selectmode='extended',
                                                selecttype='cell',
                                                rowstretch='none',
                                                colstretch='none',
                                                rowheight=-30,
                                                colwidth=15,
                                                flashmode='off',
                                                anchor=NE,
                                                usecommand=1,
                                                background='#d00d00d00',
                                                relief='ridge',
                                                command=self._valueCommand,
                                                takefocus=False,
                                                rowseparator='\n',
                                                wrap=1,
                                                multiline=1,
                                                borderwidth=(1, 1, 0, 0),
                                                browsecmd=browsecmd)

        self.lastMouseCoordinates = None
        self.toolTipShown = False
        self.headerToolTipText = StringVar()
        self.headerToolTip = MyToolTip(self, textvariable=self.headerToolTipText, wraplength=480, follow_mouse=True, state="disabled")

        # Extra key bindings for navigating through the table:
        # Tab: go right
        # Return (and Enter): go down
        self.bind("<Tab>", func=self.cellRight)
        self.bind("<Return>", func=self.cellDown)
        if False:
            self.bind("<Motion>", func=self.mouseMotion)

        # Configure the graphical appearance of the cells by using tags
        self.tag_configure('sel', background = '#b00e00e60',
                           fg='#000000000')
        self.tag_configure('active', background = '#000f70ff0',
                           fg='#000000000')
        self.tag_configure('title', anchor='w', bg='#d00d00d00',
                           fg='#000000000', relief='ridge')
        self.tag_configure(XbrlTable.TG_DISABLED, bg='#d00d00d00',
                           fg='#000000000',
                           relief='flat', state='disabled')

        if rows+columns > 2:
            # The content of the left/top corner cell can already be defined
            topCell = '0,0'
            if titleColumns+titleRows-2 > 0:
                cellSpans = {topCell : '%i,%i'% (titleRows-1, titleColumns-1)}
                self.spans(index=None, **cellSpans)
            self.format_cell(XbrlTable.TG_TOP_LEFT, topCell)
            self.tag_raise(XbrlTable.TG_TOP_LEFT, abovethis='title')
            self.tag_configure(XbrlTable.TG_TOP_LEFT, anchor='ne')
            self.format_cell(XbrlTable.TG_BORDER_ALL, topCell)
            indexValue = {topCell:tableName}
            self.set(**indexValue)

Example 9

Project: pyunicorn
Source File: cross_recurrence_plot.py
View license
    def __init__(self, x, y, metric="supremum", normalize=False,
                 sparse_rqa=False, silence_level=0, **kwds):
        """
        Initialize an instance of CrossRecurrencePlot.

        .. note::
           For a cross recurrence plot, time series x and y generally do
           **not** need to have the same length!

        Either recurrence threshold ``threshold`` or recurrence rate
        ``recurrence_rate`` have to be given as keyword arguments.

        Embedding is only supported for scalar time series. If embedding
        dimension ``dim`` and delay ``tau`` are **both** given as keyword
        arguments, embedding is applied. Multidimensional time series are
        processed as is by default. The same delay embedding is applied to both
        time series x and y.

        :type x: 2D array (time, dimension)
        :arg x: One of the time series to be analyzed, can be scalar or
            multi-dimensional.
        :type y: 2D array (time, dimension)
        :arg y: One of the time series to be analyzed, can be scalar or
            multi-dimensional.
        :arg str metric: The metric for measuring distances in phase space
            ("manhattan", "euclidean", "supremum").
        :arg bool normalize: Decide whether to normalize both time series to
            zero mean and unit standard deviation.
        :arg number silence_level: Inverse level of verbosity of the object.
        :arg number threshold: The recurrence threshold keyword for generating
            the cross recurrence plot using a fixed threshold.
        :arg number recurrence_rate: The recurrence rate keyword for generating
            the cross recurrence plot using a fixed recurrence rate.
        :arg number dim: The embedding dimension.
        :arg number tau: The embedding delay.
        """
        threshold = kwds.get("threshold")
        recurrence_rate = kwds.get("recurrence_rate")

        RecurrencePlot.__init__(
            self, np.empty((2, 0)), metric=metric, normalize=normalize,
            sparse_rqa=sparse_rqa, threshold=threshold,
            recurrence_rate=recurrence_rate, silence_level=silence_level)

        self.CR = None
        """The cross recurrence matrix."""
        self.N = 0
        """The length of the embedded time series x."""
        self.M = 0
        """The length of the embedded time series y."""

        #  Store time series
        self.x = x.copy().astype("float32")
        """The time series x."""
        self.y = y.copy().astype("float32")
        """The time series y."""

        #  Reshape time series
        self.x.shape = (self.x.shape[0], -1)
        self.y.shape = (self.y.shape[0], -1)

        #  Normalize time series
        if normalize:
            self.normalize_time_series(self.x)
            self.normalize_time_series(self.y)

        #  Get embedding dimension and delay from **kwds
        dim = kwds.get("dim")
        tau = kwds.get("tau")

        if dim is not None and tau is not None:
            #  Embed the time series
            self.x_embedded = self.embed_time_series(self.x, dim, tau)
            """The embedded time series x."""
            self.y_embedded = self.embed_time_series(self.y, dim, tau)
            """The embedded time series y."""
        else:
            self.x_embedded = self.x
            self.y_embedded = self.y

        #  construct recurrence plot accordingly to threshold / recurrence rate
        if threshold is not None:
            #  Calculate the recurrence matrix R using the radius of
            #  neighborhood threshold
            CrossRecurrencePlot.set_fixed_threshold(self, threshold)
        elif recurrence_rate is not None:
            #  Calculate the recurrence matrix R using a fixed recurrence rate
            CrossRecurrencePlot.\
                set_fixed_recurrence_rate(self, recurrence_rate)
        else:
            raise NameError(
                "Please give either threshold or recurrence_rate " +
                "to construct the cross recurrence plot!")

Example 10

Project: pyunicorn
Source File: joint_recurrence_plot.py
View license
    def __init__(self, x, y, metric=("supremum", "supremum"),
                 normalize=False, lag=0, silence_level=0, **kwds):
        """
        Initialize an instance of JointRecurrencePlot.

        .. note::
           For a joint recurrence plot, time series x and y need to have the
           same length!

        Either recurrence thresholds ``threshold``/``threshold_std`` or
        recurrence rates ``recurrence_rate`` have to be given as keyword
        arguments.

        Embedding is only supported for scalar time series. If embedding
        dimension ``dim`` and delay ``tau`` are **both** given as keyword
        arguments, embedding is applied. Multidimensional time series are
        processed as is by default.

        :type x: 2D array (time, dimension)
        :arg x: Time series x to be analyzed (scalar/multi-dimensional).
        :type y: 2D array (time, dimension)
        :arg y: Time series y to be analyzed (scalar/multi-dimensional).
        :type metric: tuple of string
        :arg metric: The metric for measuring distances in phase space
            ("manhattan", "euclidean", "supremum"). Give separately for each
            time series.
        :type normalize: tuple of bool
        :arg normalize: Decide whether to normalize the time series to zero
            mean and unit standard deviation. Give separately for each time
            series.
        :arg number lag: To create a delayed version of the JRP.
        :arg number silence_level: Inverse level of verbosity of the object.
        :type threshold: tuple of number
        :keyword threshold: The recurrence threshold keyword for generating the
            recurrence plot using a fixed threshold.  Give separately for each
            time series.
        :type threshold_std: tuple of number
        :keyword threshold_std: The recurrence threshold keyword for generating
            the recurrence plot using a fixed threshold in units of the time
            series' STD. Give separately for each time series.
        :type recurrence_rate: tuple of number
        :keyword recurrence_rate: The recurrence rate keyword for generating
            the recurrence plot using a fixed recurrence rate. Give separately
            for each time series.
        :type dim: tuple of number
        :keyword dim: Embedding dimension. Give separately for each time
            series.
        :type tau: tuple of number
        :keyword tau: Embedding delay. Give separately for each time series.
        """
        threshold = kwds.get("threshold")
        threshold_std = kwds.get("threshold_std")
        recurrence_rate = kwds.get("recurrence_rate")

        RecurrencePlot.__init__(
            self, np.empty((2, 0)), metric=metric[0], normalize=normalize,
            threshold=threshold[0] if threshold else 0,
            recurrence_rate=recurrence_rate, silence_level=silence_level)

        self.JR = None
        """The joint recurrence matrix."""
        self.N = 0
        """The length of both embedded time series x and y."""

        #  Check for consistency: x and y need to have the same length
        if x.shape[0] == y.shape[0]:

            #  Store time series
            self.x = x.copy().astype("float32")
            """The time series x."""
            self.y = y.copy().astype("float32")
            """The time series y."""

            #  Reshape time series
            self.x.shape = (self.x.shape[0], -1)
            self.y.shape = (self.y.shape[0], -1)

            #  Normalize time series
            if normalize:
                self.normalize_time_series(self.x)
                self.normalize_time_series(self.y)

            #  Store lag
            self.lag = lag
            """Used to create a delayed JRP."""

            #  Get embedding dimension and delay from **kwds
            dim = kwds.get("dim")
            tau = kwds.get("tau")

            if dim is not None and tau is not None:
                #  Embed the time series
                self.x_embedded = \
                    self.embed_time_series(self.x, dim[0], tau[0])
                """The embedded time series x."""
                self.y_embedded = \
                    self.embed_time_series(self.y, dim[1], tau[1])
                """The embedded time series y."""
            else:
                self.x_embedded = self.x
                self.y_embedded = self.y

            #  Prune embedded time series to same length
            #  (number of state vectors)
            min_N = min(self.x_embedded.shape[0], self.y_embedded.shape[0])
            self.x_embedded = self.x_embedded[:min_N, :]
            self.y_embedded = self.y_embedded[:min_N, :]

            #  construct recurrence plot accordingly to
            #  threshold / recurrence rate
            if np.abs(lag) > x.shape[0]:
                #  Lag must be smaller than size of recurrence plot
                raise ValueError(
                    "Delay value (lag) must not exceed length of time series!")
            elif threshold is not None:
                #  Calculate the recurrence matrix R using the radius of
                #  neighborhood threshold
                JointRecurrencePlot.set_fixed_threshold(self, threshold)
            elif threshold_std is not None:
                #  Calculate the recurrence matrix R using the radius of
                #  neighborhood threshold in units of the time series' STD
                JointRecurrencePlot.\
                    set_fixed_threshold_std(self, threshold_std)
            elif recurrence_rate is not None:
                #  Calculate the recurrence matrix R using a fixed recurrence
                #  rate
                JointRecurrencePlot.\
                    set_fixed_recurrence_rate(self, recurrence_rate)
            else:
                raise NameError(
                    "Please give either threshold or recurrence_rate " +
                    "to construct the joint recurrence plot!")

            #  No treatment of missing values yet!
            self.missing_values = False

        else:
            raise ValueError(
                "Both time series x and y need to have the same length!")

Example 11

Project: pyunicorn
Source File: surrogates.py
View license
    def twin_surrogates(self, original_data, dimension, delay, threshold,
                        min_dist=7):
        """
        Return surrogates using the twin surrogate method.

        Scalar twin surrogates are created by isolating the first component
        (dimension) of the twin surrogate trajectories.

        Twin surrogates share linear and nonlinear properties with the original
        time series, since they correspond to realizations of trajectories of
        the same dynamical systems with different initial conditions.

        References: [Thiel2006]_ [*], [Marwan2007]_.

        The twin lists of all time series are cached to facilitate a faster
        generation of several surrogates for each time series. Hence,
        :meth:`clear_cache` has to be called before generating twin surrogates
        from a different set of time series!

        :type original_data: 2D array [index, time]
        :arg original_data: The original time series.
        :arg int dimension: The embedding dimension.
        :arg int delay: The embedding delay.
        :arg float threshold: The recurrence threshold.
        :arg number min_dist: The minimum temporal distance for twins.
        :rtype: 2D array [index, time]
        :return: the twin surrogates.
        """
        #  The algorithm proceeds in several steps:
        #  1. Embed the original_data time series, using time delay embedding
        #     for simplicity. Use the same dimension and time delay delay for
        #     all time series for simplicity. Determine delay using time
        #     delayed mutual information and d using false nearest neighbors
        #     methods.
        #  2. Use the algorithm proposed in [*] to find twins
        #  3. Reconstruct one-dimensional twin surrogate time series
        (N, n_time) = original_data.shape

        #  Make sure that twins are calculated only once
        if self._twins_cached:
            twins = self._twins
        else:
            embedding = self.embed_time_series_array(original_data,
                                                     dimension, delay)
            twins = self.twins(embedding, threshold, min_dist)
            self._twins = twins
            self._twins_cached = True

        surrogates = np.empty(original_data.shape)

        code = r"""
        int i, j, k, new_k, n_twins, rand;

        //  Initialize random number generator
        srand48(time(0));

        for (i = 0; i < N; i++) {
            //  Get the twin list for time series i
            py::list twins_i = PyList_GetItem(twins, i);

            //  Randomly choose a starting point in the original_data
            //  trajectory.
            k = floor(drand48() * n_time);

            j = 0;

            while (j < n_time) {
                surrogates(i,j) = original_data(i,k);

                //  Get the list of twins of sample k in the original_data
                //  time series.
                py::list twins_ik = PyList_GetItem(twins_i,k);

                //  Get the number of twins of k
                n_twins = PyList_Size(twins_ik);

                //  If k has no twins, go to the next sample k+1. If k has
                //  twins at m, choose among m+1 and k+1 with equal probability
                if (n_twins == 0)
                    k++;
                else {
                    //  Generate a random integer between 0 and n_twins
                    rand = floor(drand48() * (n_twins + 1));

                    //  If rand = n_twins go to sample k+1, otherwise jump
                    //  to the future of one of the twins.
                    if (rand == n_twins)
                        k++;
                    else {
                        k = twins_ik[rand];
                        k++;
                    }

                }

                //  If the new k >= n_time, choose a new random starting point
                //  in the original_data time series.
                if (k >= n_time) {
                    do {
                        new_k = floor(drand48() * n_time);
                    }
                    while (k == new_k);

                    k = new_k;
                }

                j++;
            }
        }
        """
        weave_inline(locals(), code,
                     ['N', 'n_time', 'original_data', 'twins', 'surrogates'])
        return surrogates

Example 12

Project: OpenPNM
Source File: __Cubic__.py
View license
    def __init__(self, shape=None, template=None, spacing=[1, 1, 1],
                 connectivity=6, **kwargs):
        super().__init__(**kwargs)

        if shape is not None:
            arr = np.atleast_3d(np.empty(shape))
        elif template is not None:
            arr = sp.array(template, ndmin=3, dtype=bool)
        else:
            arr = np.atleast_3d(np.empty([1, 1, 1]))

        # Store original network shape
        self._shape = sp.shape(arr)
        # Store network spacing
        self._spacing = sp.ones(3)*sp.array(spacing, ndmin=1)

        points = np.array([i for i, v in np.ndenumerate(arr)], dtype=float)
        points += 0.5
        points *= spacing

        I = np.arange(arr.size).reshape(arr.shape)

        face_joints = [
            (I[:, :, :-1], I[:, :, 1:]),
            (I[:, :-1], I[:, 1:]),
            (I[:-1], I[1:]),
        ]

        corner_joints = [
            (I[:-1, :-1, :-1], I[1:, 1:, 1:]),
            (I[:-1, :-1, 1:], I[1:, 1:, :-1]),
            (I[:-1, 1:, :-1], I[1:, :-1, 1:]),
            (I[1:, :-1, :-1], I[:-1, 1:, 1:]),
        ]

        edge_joints = [
            (I[:, :-1, :-1], I[:, 1:, 1:]),
            (I[:, :-1, 1:], I[:, 1:, :-1]),
            (I[:-1, :, :-1], I[1:, :, 1:]),
            (I[1:, :, :-1], I[:-1, :, 1:]),
            (I[1:, 1:, :], I[:-1, :-1, :]),
            (I[1:, :-1, :], I[:-1, 1:, :]),
        ]

        if connectivity == 6:
            joints = face_joints
        elif connectivity == 8:
            joints = corner_joints
        elif connectivity == 12:
            joints = edge_joints
        elif connectivity == 14:
            joints = face_joints + corner_joints
        elif connectivity == 18:
            joints = face_joints + edge_joints
        elif connectivity == 20:
            joints = edge_joints + corner_joints
        elif connectivity == 26:
            joints = face_joints + corner_joints + edge_joints
        else:
            raise Exception('Invalid connectivity receieved. Must be 6, 8, 12, 14, '
                            '18, 20 or 26')

        I = np.arange(arr.size).reshape(arr.shape)
        tails, heads = [], []
        for T, H in joints:
            tails.extend(T.flat)
            heads.extend(H.flat)
        pairs = np.vstack([tails, heads]).T

        self['pore.coords'] = points
        self['throat.conns'] = pairs
        self['pore.all'] = np.ones(len(self['pore.coords']), dtype=bool)
        self['throat.all'] = np.ones(len(self['throat.conns']), dtype=bool)
        self['pore.index'] = sp.arange(0, len(self['pore.coords']))

        self._label_surfaces()

        # If an image was sent as 'template', then trim network to image shape
        if template is not None:
            self.trim(~arr.flatten())

Example 13

View license
def auc_metric(solution, prediction, task=BINARY_CLASSIFICATION):
    """
    Normarlized Area under ROC curve (AUC).

    Return Gini index = 2*AUC-1 for  binary classification problems.
    Should work for a vector of binary 0/1 (or -1/1)"solution" and any discriminant values
    for the predictions. If solution and prediction are not vectors, the AUC
    of the columns of the matrices are computed and averaged (with no weight).
    The same for all classification problems (in fact it treats well only the
    binary and multilabel classification problems).
    :param solution:
    :param prediction:
    :param task:
    :return:
    """
    if task == BINARY_CLASSIFICATION:
        if len(solution.shape) == 1:
            # Solution won't be touched - no copy
            solution = solution.reshape((-1, 1))
        elif len(solution.shape) == 2:
            if solution.shape[1] > 1:
                raise ValueError('Solution array must only contain one class '
                                 'label, but contains %d' % solution.shape[1])
        else:
            raise ValueError('Solution.shape %s' % solution.shape)
        solution = solution.copy()

        if len(prediction.shape) == 2:
            if prediction.shape[1] > 2:
                raise ValueError('A prediction array with probability values '
                                 'for %d classes is not a binary '
                                 'classification problem' % prediction.shape[1])
            elif prediction.shape[1] == 2:
                # Prediction will be copied into a new binary array - no copy
                prediction = prediction[:, 1].reshape((-1, 1))
        else:
            raise ValueError('Invalid prediction shape %s' % prediction.shape)

    elif task == MULTICLASS_CLASSIFICATION:
        if len(solution.shape) == 1:
            solution = create_multiclass_solution(solution, prediction)
        elif len(solution.shape) == 2:
            if solution.shape[1] > 1:
                raise ValueError('Solution array must only contain one class '
                                 'label, but contains %d' % solution.shape[1])
            else:
                solution = create_multiclass_solution(solution.reshape((-1, 1)),
                                                      prediction)
        else:
            raise ValueError('Solution.shape %s' % solution.shape)
    elif task == MULTILABEL_CLASSIFICATION:
        solution = solution.copy()
    else:
        raise NotImplementedError('auc_metric does not support task type %s'
                                  % task)

    solution, prediction = normalize_array(solution, prediction.copy())

    label_num = solution.shape[1]
    auc = np.empty(label_num)
    for k in range(label_num):
        r_ = scipy.stats.rankdata(prediction[:, k])
        s_ = solution[:, k]
        if sum(s_) == 0:
            print(
                'WARNING: no positive class example in class {}'.format(k + 1))
        npos = np.sum(s_ == 1)
        nneg = np.sum(s_ < 1)
        auc[k] = (np.sum(r_[s_ == 1]) - npos * (npos + 1) / 2) / (nneg * npos)
    auc[~np.isfinite(auc)] = 0
    return 2 * np.mean(auc) - 1

Example 14

Project: pupil
Source File: square_marker_detect.py
View license
    def track_in_frame(self,gray_img,grid_size,min_marker_perimeter=40,aperture=11,visualize=False):
        observed_markers = self.detect_in_frame(gray_img,grid_size,min_marker_perimeter,aperture,visualize)

        # flatten observed marker into state entry format
        map_fn = self.make_raw_to_flat_map(gray_img.T.shape)
        observed_markers = map(map_fn, observed_markers)
        distances = np.empty((len(self.state), len(observed_markers)))
        for n_i, hist_m in enumerate(self.state):
            for n_ip1, new_m in enumerate(observed_markers):
                distances[n_i, n_ip1] = self.distance(hist_m, new_m)

        hist_to_match = np.ones(len(self.state)).astype(bool)
        observ_to_match =np.ones(len(observed_markers)).astype(bool)

        while hist_to_match.any() and observ_to_match.any():
            match_dist = np.min(distances)
            if match_dist > self.max_match_dist:
                break # do not match markers that are to distant to each other
            match_idx = np.argmin(distances)
            matched_hist_idx, matched_observ_idx = np.unravel_index(match_idx, distances.shape)

            self._merge_marker(
                self.state[matched_hist_idx],
                observed_markers[matched_observ_idx])

            # remove rows and columns
            distances[matched_hist_idx, :] = 2
            distances[:, matched_observ_idx] = 2
            hist_to_match[matched_hist_idx] = False
            observ_to_match[matched_observ_idx] = False

        prev_pts = []
        unmatched_history = filter(
            lambda (_,unmatched): unmatched,
            izip(self.state,hist_to_match))
        for hist_marker, unmatched in unmatched_history:
            # use optical flow to track unmatched markers...
            norm_verts = hist_marker[self.vert_slc].reshape((4,2))
            verts = (norm_verts * gray_img.T.shape).astype(np.float32)
            prev_pts.append(verts)

        if self.prev_img is not None and prev_pts:
            prev_pts = np.vstack(prev_pts)
            new_pts, flow_found, err = cv2.calcOpticalFlowPyrLK(
                self.prev_img, gray_img, prev_pts,
                minEigThreshold=0.01,**lk_params)
            for marker_idx in xrange(flow_found.shape[0]/4):
                hist_marker = unmatched_history[marker_idx][0]
                hist_marker[self.loc_conf_idx] -= self.unmatched_penalty
                m_slc = slice(marker_idx*4,marker_idx*4+4)
                if flow_found[m_slc].sum() >= 4 :
                    found, _  = np.where(flow_found[m_slc])
                    # calculate differences
                    old_verts = prev_pts[m_slc][found,:]
                    new_verts = new_pts[m_slc][found,:]
                    vert_difs = new_verts - old_verts

                    # calc mean dif
                    mean_dif = vert_difs.mean(axis=0)
                    # take 3 closest difs
                    mean_dif_dist = np.linalg.norm(mean_dif - vert_difs,axis=1)
                    closest_mean_dif = np.argsort(mean_dif_dist)[:-1]
                    # recalc mean dif
                    mean_dif = vert_difs[closest_mean_dif].mean(axis=0)

                    m_id = bin_list_to_int(np.round(hist_marker[self.id_slc]).astype(int))
                    # apply mean dif and normalize
                    proj_verts = (prev_pts[m_slc] + mean_dif) / gray_img.T.shape
                    hist_marker[self.vert_slc] = proj_verts.flatten()
                else:
                    # penalize again, if optical flow did not work
                    hist_marker[self.loc_conf_idx] -= 10*self.unmatched_penalty

        for observ_marker, unmatched in zip(observed_markers, observ_to_match):
            # add unmatched oberservations
            if unmatched: self._append_marker(observ_marker)

        # purge low confidence markers
        for m_idx, m_hist in reversedEnumerate(self.state):
            if (self.marker_id_confidence(m_hist) < self.id_purge_threshold or
                m_hist[self.loc_conf_idx] < self.loc_purge_threshold):
                del self.state[m_idx]

        tracked_markers = self.extract_markers()
        for tracked_m in tracked_markers:
            norm_verts = tracked_m['norm_verts']
            # cv2.getPerspectiveTransform needs np.float32
            tracked_m['verts'] = (norm_verts * gray_img.T.shape).astype(np.float32)
            tracked_m['centroid'] = np.mean(tracked_m['verts'], axis=0).reshape((2,))
            tracked_m['perimeter'] = cv2.arcLength(tracked_m['verts'],closed=True)
            # What is this for?
            # tracked_m['img'] = gray_img

        self.prev_img = gray_img

        return tracked_markers

Example 15

Project: pycalphad
Source File: calculate.py
View license
def _compute_phase_values(phase_obj, components, variables, statevar_dict,
                          points, func, output, maximum_internal_dof, broadcast=True, fake_points=False,
                          largest_energy=None):
    """
    Calculate output values for a particular phase.

    Parameters
    ----------
    phase_obj : Phase
        Phase object from a thermodynamic database.
    components : list
        Names of components to consider in the calculation.
    variables : list
        Names of variables in the phase's internal degrees of freedom.
    statevar_dict : OrderedDict {str -> float or sequence}
        Mapping of state variables to desired values. This will broadcast if necessary.
    points : ndarray
        Inputs to 'func', except state variables. Columns should be in 'variables' order.
    func : callable
        Function of state variables and 'variables'.
        See 'make_callable' docstring for details.
    output : string
        Desired name of the output result in the Dataset.
    maximum_internal_dof : int
        Largest number of internal degrees of freedom of any phase. This is used
        to guarantee different phase's Datasets can be concatenated.
    broadcast : bool
        If True, broadcast state variables against each other to create a grid.
        If False, assume state variables are given as equal-length lists (or single-valued).
    fake_points : bool, optional (Default: False)
        If True, the first few points of the output surface will be fictitious
        points used to define an equilibrium hyperplane guaranteed to be above
        all the other points. This is used for convex hull computations.

    Returns
    -------
    Dataset of the output attribute as a function of state variables

    Examples
    --------
    None yet.
    """
    if broadcast:
        # Broadcast compositions and state variables along orthogonal axes
        # This lets us eliminate an expensive Python loop
        statevars = np.meshgrid(*itertools.chain(statevar_dict.values(),
                                                     [np.empty(points.shape[-2])]),
                                    sparse=True, indexing='ij')[:-1]
        points = broadcast_to(points, tuple(len(np.atleast_1d(x)) for x in statevar_dict.values()) + points.shape[-2:])
    else:
        statevars = list(np.atleast_1d(x) for x in statevar_dict.values())
        statevars_ = []
        for statevar in statevars:
            if (len(statevar) != len(points)) and (len(statevar) == 1):
                statevar = np.repeat(statevar, len(points))
            if (len(statevar) != len(points)) and (len(statevar) != 1):
                raise ValueError('Length of state variable list and number of given points must be equal when '
                                 'broadcast=False.')
            statevars_.append(statevar)
        statevars = statevars_
    # func may only have support for vectorization along a single axis (no broadcasting)
    # we need to force broadcasting and flatten the result before calling
    bc_statevars = [np.ascontiguousarray(broadcast_to(x, points.shape[:-1]).reshape(-1)) for x in statevars]
    pts = points.reshape(-1, points.shape[-1]).T

    phase_output = func(*itertools.chain(bc_statevars, pts))
    if isinstance(phase_output, (float, int)):
        phase_output = broadcast_to(phase_output, points.shape[:-1])
    phase_output = np.asarray(phase_output, dtype=np.float)
    phase_output.shape = points.shape[:-1]
    if fake_points:
        phase_output = np.concatenate((broadcast_to(largest_energy, points.shape[:-2] + (len(components),)), phase_output), axis=-1)
        phase_names = np.concatenate((broadcast_to('_FAKE_', points.shape[:-2] + (len(components),)),
                                      np.full(points.shape[:-1], phase_obj.name, dtype='U' + str(len(phase_obj.name)))), axis=-1)
    else:
        phase_names = np.full(points.shape[:-1], phase_obj.name, dtype='U'+str(len(phase_obj.name)))

    # Map the internal degrees of freedom to global coordinates
    # Normalize site ratios by the sum of site ratios times a factor
    # related to the site fraction of vacancies
    site_ratio_normalization = np.zeros(points.shape[:-1])
    for idx, sublattice in enumerate(phase_obj.constituents):
        vacancy_column = np.ones(points.shape[:-1])
        if 'VA' in set(sublattice):
            var_idx = variables.index(v.SiteFraction(phase_obj.name, idx, 'VA'))
            vacancy_column -= points[..., :, var_idx]
        site_ratio_normalization += phase_obj.sublattices[idx] * vacancy_column

    phase_compositions = np.empty(points.shape[:-1] + (len(components),))
    for col, comp in enumerate(components):
        avector = [float(vxx.species == comp) * \
            phase_obj.sublattices[vxx.sublattice_index] for vxx in variables]
        phase_compositions[..., :, col] = np.divide(np.dot(points[..., :, :], avector),
                                               site_ratio_normalization)
    if fake_points:
        phase_compositions = np.concatenate((np.broadcast_to(np.eye(len(components)), points.shape[:-2] + (len(components), len(components))), phase_compositions), axis=-2)

    coordinate_dict = {'component': components}
    # Resize 'points' so it has the same number of columns as the maximum
    # number of internal degrees of freedom of any phase in the calculation.
    # We do this so that everything is aligned for concat.
    # Waste of memory? Yes, but the alternatives are unclear.
    if fake_points:
        expanded_points = np.full(points.shape[:-2] + (len(components)+points.shape[-2], maximum_internal_dof), np.nan)
        expanded_points[..., len(components):, :points.shape[-1]] = points
    else:
        expanded_points = np.full(points.shape[:-1] + (maximum_internal_dof,), np.nan)
        expanded_points[..., :points.shape[-1]] = points
    if broadcast:
        coordinate_dict.update({key: np.atleast_1d(value) for key, value in statevar_dict.items()})
        output_columns = [str(x) for x in statevar_dict.keys()] + ['points']
    else:
        output_columns = ['points']
    data_arrays = {'X': (output_columns + ['component'], phase_compositions),
                   'Phase': (output_columns, phase_names),
                   'Y': (output_columns + ['internal_dof'], expanded_points),
                   output: (['dim_'+str(i) for i in range(len(phase_output.shape) - len(output_columns))] + output_columns, phase_output)
                   }
    if not broadcast:
        # Add state variables as data variables rather than as coordinates
        for sym, vals in zip(statevar_dict.keys(), statevars):
            data_arrays.update({sym: (output_columns, vals)})

    return Dataset(data_arrays, coords=coordinate_dict)

Example 16

Project: pycalphad
Source File: lower_convex_hull.py
View license
def old_lower_convex_hull(global_grid, result_array, verbose=False):
    """
    Find the simplices on the lower convex hull satisfying the specified
    conditions in the result array.

    Parameters
    ----------
    global_grid : Dataset
        A sample of the energy surface of the system.
    result_array : Dataset
        This object will be modified!
        Coordinates correspond to conditions axes.
    verbose : bool
        Display details to stdout. Useful for debugging.

    Returns
    -------
    None. Results are written to result_array.

    Notes
    -----
    This routine will not check if any simplex is degenerate.
    Degenerate simplices will manifest with duplicate or NaN indices.

    Examples
    --------
    None yet.
    """
    indep_conds = sorted([x for x in sorted(result_array.coords.keys()) if x in ['T', 'P']])
    comp_conds = sorted([x for x in sorted(result_array.coords.keys()) if x.startswith('X_')])
    pot_conds = sorted([x for x in sorted(result_array.coords.keys()) if x.startswith('MU_')])
    # force conditions to have particular ordering
    conditions = indep_conds + pot_conds + comp_conds
    trial_shape = (len(result_array.coords['component']),)
    trial_points = None
    _initialize_array(global_grid, result_array)

    # Enforce ordering of shape if this is the first iteration
    #if result_array.attrs['hull_iterations'] == 1:
    #    result_array['points'] = result_array['points'].transpose(*(conditions + ['vertex']))
    #    result_array['GM'] = result_array['GM'].transpose(*conditions)
    #    result_array['NP'] = result_array['NP'].transpose(*(conditions + ['vertex']))

    # Determine starting combinations of chemical potentials and compositions
    # TODO: Check Gibbs phase rule compliance

    if len(pot_conds) > 0:
        raise NotImplementedError('Chemical potential conditions are not yet supported')

    # FIRST CASE: Only composition conditions specified
    # We only need to compute the dependent composition value directly
    # Initialize trial points as lowest energy point in the system
    if (len(comp_conds) > 0) and (len(pot_conds) == 0):
        trial_points = np.empty(result_array['GM'].T.shape)
        trial_points.fill(np.inf)
        trial_points[...] = global_grid['GM'].argmin(dim='points').values.T
        trial_points = trial_points.T
        comp_values = cartesian([result_array.coords[cond] for cond in comp_conds])
        # Insert dependent composition value
        # TODO: Handle W(comp) as well as X(comp) here
        specified_components = set([x[2:] for x in comp_conds])
        dependent_component = set(result_array.coords['component'].values) - specified_components
        dependent_component = list(dependent_component)
        if len(dependent_component) != 1:
            raise ValueError('Number of dependent components is different from one')
        insert_idx = sorted(result_array.coords['component'].values).index(dependent_component[0])
        comp_values = np.concatenate((comp_values[..., :insert_idx],
                                      1 - np.sum(comp_values, keepdims=True, axis=-1),
                                      comp_values[..., insert_idx:]),
                                     axis=-1)
        # Prevent compositions near an edge from going negative
        comp_values[np.nonzero(comp_values < MIN_SITE_FRACTION)] = MIN_SITE_FRACTION*10
        # TODO: Assumes N=1
        comp_values /= comp_values.sum(axis=-1, keepdims=True)
        #print(comp_values)

    # SECOND CASE: Only chemical potential conditions specified
    # TODO: Implementation of chemical potential

    # THIRD CASE: Mixture of composition and chemical potential conditions
    # TODO: Implementation of mixed conditions

    if trial_points is None:
        raise ValueError('Invalid conditions')

    # factored out via profiling
    result_array_GM_values = result_array.GM.values
    result_array_points_values = result_array.points.values
    result_array_MU_values = result_array.MU.values
    result_array_NP_values = result_array.NP.values
    global_grid_GM_values = global_grid.GM.values
    global_grid_X_values = global_grid.X.values

    df_shape = result_array_GM_values.shape + (len(global_grid.points),)
    driving_forces = np.zeros(df_shape, dtype=np.float)
    trial_simplices = np.empty(result_array_points_values.shape +
                               (result_array_points_values.shape[-1],), dtype=np.int)
    max_iterations = 200
    iterations = 0
    while iterations < max_iterations:
        iterations += 1
        # Initialize trial simplices with values from best guess simplices
        trial_simplices[..., :, :] = result_array_points_values[..., np.newaxis, :]
        # Trial simplices will be the current simplex with each vertex
        #     replaced by the trial point
        # Exactly one of those simplices will contain a given test point,
        #     excepting edge cases
        trial_simplices.T[np.diag_indices(trial_shape[0])] = trial_points.T
        #print('trial_simplices.shape', trial_simplices.shape)
        #print('global_grid.X.values.shape', global_grid.X.values.shape)
        flat_statevar_indices = np.unravel_index(np.arange(np.multiply.reduce(result_array_MU_values.shape)),
                                                 result_array_MU_values.shape)[:len(indep_conds)]
        #print('flat_statevar_indices', flat_statevar_indices)
        trial_matrix = global_grid_X_values[np.index_exp[flat_statevar_indices +
                                                         (trial_simplices.reshape(-1, trial_simplices.shape[-1]).T,)]]
        trial_matrix = np.rollaxis(trial_matrix, 0, -1)
        #print('trial_matrix', trial_matrix)
        # Partially ravel the array to make indexing operations easier
        trial_matrix.shape = (-1,) + trial_matrix.shape[-2:]

        # We have to filter out degenerate simplices before
        #     phase fraction computation
        # This is because even one degenerate simplex causes the entire tensor
        #     to be singular
        nondegenerate_indices = np.all(np.linalg.svd(trial_matrix,
                                                     compute_uv=False) > 1e-09,
                                       axis=-1, keepdims=True)
        #('NONDEGENERATE INDICES', nondegenerate_indices)
        # Determine how many trial simplices remain for each target point.
        # In principle this would always be one simplex per point, but once
        # some target values reach equilibrium, trial_points starts
        # to contain points already on our best guess simplex.
        # This causes trial_simplices to create degenerate simplices.
        # We can safely filter them out since those target values are
        # already at equilibrium.
        #sum_array = np.sum(nondegenerate_indices, axis=-1, dtype=np.int)
        #index_array = np.repeat(np.arange(trial_matrix.shape[0], dtype=np.int),
        #                        sum_array)
        index_array = np.arange(trial_matrix.shape[0], dtype=np.int)
        comp_shape = trial_simplices.shape[:len(indep_conds)+len(pot_conds)] + \
                     (comp_values.shape[0], trial_simplices.shape[-2])

        comp_indices = np.unravel_index(index_array, comp_shape)[len(indep_conds)+len(pot_conds)]
        fractions = np.full(result_array_points_values.shape +
                            (result_array_points_values.shape[-1],), -1.)
        fractions[np.unravel_index(index_array, fractions.shape[:-1])] = \
            stacked_lstsq(np.swapaxes(trial_matrix[index_array], -2, -1),
                          comp_values[comp_indices])
        fractions /= fractions.sum(axis=-1, keepdims=True)
        #print('fractions', fractions)
        # A simplex only contains a point if its barycentric coordinates
        # (phase fractions) are non-negative.
        bounding_indices = np.all(fractions >= -MIN_SITE_FRACTION*100, axis=-1)
        #print('BOUNDING INDICES', bounding_indices)
        if ~np.any(bounding_indices):
            raise ValueError('Desired composition is not inside any candidate simplex. This is a bug.')
        multiple_success_trials = np.sum(bounding_indices, axis=-1, dtype=np.int, keepdims=False) != 1
        #print('MULTIPLE SUCCESS TRIALS SHAPE', np.nonzero(multiple_success_trials))
        if np.any(multiple_success_trials):
            saved_trial = np.zeros_like(multiple_success_trials, dtype=np.int)
            # Case of only degenerate simplices (zero bounding)
            # Choose trial with "least negative" fraction
            zero_success_indices = np.logical_and(~nondegenerate_indices.reshape(bounding_indices.shape),
                                                  multiple_success_trials[..., np.newaxis])
            saved_trial[np.nonzero(zero_success_indices.any(axis=-1))] = \
                np.argmax(fractions[np.nonzero(zero_success_indices.any(axis=-1))].min(axis=-1), axis=-1)
            # Case of multiple bounding non-degenerate simplices
            # Choose the first one. This addresses gh-28.
            multiple_bounding_indices = \
                np.logical_and(np.logical_and(bounding_indices, nondegenerate_indices.reshape(bounding_indices.shape)),
                               multiple_success_trials[..., np.newaxis])
            #print('MULTIPLE SUCCESS TRIALS.shape', multiple_success_trials.shape)
            #print('BOUNDING INDICES.shape', bounding_indices.shape)
            #print('MULTIPLE_BOUNDING_INDICES.shape', multiple_bounding_indices.shape)
            saved_trial[np.nonzero(multiple_bounding_indices.any(axis=-1))] = \
                np.argmax(multiple_bounding_indices[np.nonzero(multiple_bounding_indices.any(axis=-1))], axis=-1)
            #print('SAVED TRIAL.shape', saved_trial.shape)
            #print('BOUNDING INDICES BEFORE', bounding_indices)
            bounding_indices[np.nonzero(multiple_success_trials)] = False
            #print('BOUNDING INDICES FALSE', bounding_indices)
            bounding_indices[np.nonzero(multiple_success_trials) + \
                             (saved_trial[np.nonzero(multiple_success_trials)],)] = True
            #print('BOUNDING INDICES AFTER', bounding_indices)
        fractions.shape = (-1, fractions.shape[-1])
        bounding_indices.shape = (-1,)
        index_array = np.arange(trial_matrix.shape[0], dtype=np.int)[bounding_indices]

        raveled_simplices = trial_simplices.reshape((-1,) + trial_simplices.shape[-1:])
        candidate_simplices = raveled_simplices[index_array, :]
        #print('candidate_simplices', candidate_simplices)

        # We need to convert the flat index arrays into multi-index tuples.
        # These tuples will tell us which state variable combinations are relevant
        # for the calculation. We can drop the last dimension, 'trial'.
        #print('trial_simplices.shape[:-1]', trial_simplices.shape[:-1])
        statevar_indices = np.unravel_index(index_array, trial_simplices.shape[:-1]
                                            )[:len(indep_conds)+len(pot_conds)]
        aligned_energies = global_grid_GM_values[statevar_indices + (candidate_simplices.T,)].T
        statevar_indices = tuple(x[..., np.newaxis] for x in statevar_indices)
        #print('statevar_indices', statevar_indices)
        aligned_compositions = global_grid_X_values[np.index_exp[statevar_indices + (candidate_simplices,)]]
        #print('aligned_compositions', aligned_compositions)
        #print('aligned_energies', aligned_energies)
        candidate_potentials = stacked_lstsq(aligned_compositions.astype(np.float, copy=False),
                                             aligned_energies.astype(np.float, copy=False))
        #print('candidate_potentials', candidate_potentials)
        #logger.debug('candidate_simplices: %s', candidate_simplices)
        comp_indices = np.unravel_index(index_array, comp_shape)[len(indep_conds)+len(pot_conds)]
        #print('comp_values[comp_indices]', comp_values[comp_indices])
        candidate_energies = np.multiply(candidate_potentials,
                                         comp_values[comp_indices]).sum(axis=-1)
        #print('candidate_energies', candidate_energies)

        # Generate a matrix of energies comparing our calculations for this iteration
        # to each other.
        # 'conditions' axis followed by a 'trial' axis
        # Empty values are filled in with infinity
        comparison_matrix = np.empty([int(trial_matrix.shape[0] / trial_shape[0]),
                                      trial_shape[0]])
        if comparison_matrix.shape[0] != aligned_compositions.shape[0]:
            raise ValueError('Arrays have become misaligned. This is a bug. Try perturbing your composition conditions '
                             'by a small amount (1e-4). If you would like, you can report this issue to the development'
                             ' team and they will fix it for future versions.')
        comparison_matrix.fill(np.inf)
        comparison_matrix[np.divide(index_array, trial_shape[0]).astype(np.int),
                          np.mod(index_array, trial_shape[0])] = candidate_energies
        #print('comparison_matrix', comparison_matrix)

        # If a condition point is all infinities, it means we did not calculate it
        # We should filter those out from any comparisons
        calculated_indices = ~np.all(comparison_matrix == np.inf, axis=-1)
        # Extract indices for trials with the lowest energy for each target point
        lowest_energy_indices = np.argmin(comparison_matrix[calculated_indices],
                                          axis=-1)
        # Filter conditions down to only those calculated this iteration
        calculated_conditions_indices = np.arange(comparison_matrix.shape[0])[calculated_indices]

        #print('comparison_matrix[calculated_conditions_indices,lowest_energy_indices]',comparison_matrix[calculated_conditions_indices,
        #                                    lowest_energy_indices])
        # This has to be greater-than-or-equal because, in the case where
        # the desired condition is right on top of a simplex vertex (gh-28), there
        # will be no change in energy changing a "_FAKE_" vertex to a real one.
        is_lower_energy = comparison_matrix[calculated_conditions_indices,
                                            lowest_energy_indices] <= \
            result_array_GM_values.flat[calculated_conditions_indices]
        #print('is_lower_energy', is_lower_energy)

        # These are the conditions we will update this iteration
        final_indices = calculated_conditions_indices[is_lower_energy]
        #print('final_indices', final_indices)
        # Convert to multi-index form so we can index the result array
        final_multi_indices = np.unravel_index(final_indices,
                                               result_array_GM_values.shape)

        updated_potentials = candidate_potentials[is_lower_energy]
        result_array_points_values[final_multi_indices] = candidate_simplices[is_lower_energy]
        result_array_GM_values[final_multi_indices] = candidate_energies[is_lower_energy]
        result_array_MU_values[final_multi_indices] = updated_potentials
        result_array_NP_values[final_multi_indices] = \
            fractions[np.nonzero(bounding_indices)][is_lower_energy]
        #print('result_array.GM.values', result_array.GM.values)

        # By profiling, it's faster to recompute all driving forces in-place
        # versus doing fancy indexing to only update "changed" driving forces
        # This avoids the creation of expensive temporaries
        np.einsum('...i,...i',
                  result_array_MU_values[..., np.newaxis, :],
                  global_grid_X_values[np.index_exp[...] + ((np.newaxis,) * len(comp_conds)) + np.index_exp[:, :]],
                  out=driving_forces)
        np.subtract(driving_forces,
                    global_grid_GM_values[np.index_exp[...] + ((np.newaxis,) * len(comp_conds)) + np.index_exp[:]],
                    out=driving_forces)


        # Update trial points to choose points with largest remaining driving force
        trial_points = np.argmax(driving_forces, axis=-1)
        #print('trial_points', trial_points)
        #logger.debug('trial_points: %s', trial_points)

        # If all driving force (within some tolerance) is consumed, we found equilibrium
        if np.all(driving_forces <= DRIVING_FORCE_TOLERANCE):
            return result_array
    if verbose:
        print('Max hull iterations exceeded. Remaining driving force: ', driving_forces.max())
    return result_array

Example 17

Project: pydy
Source File: ode_function_generators.py
View license
    def __init__(self, right_hand_side, coordinates, speeds, constants,
                 mass_matrix=None, coordinate_derivatives=None,
                 specifieds=None, linear_sys_solver='numpy',
                 constants_arg_type=None, specifieds_arg_type=None):
        """Generates a numerical function which can evaluate the right hand
        side of the first order ordinary differential equations from a
        system described by one of the following three symbolic forms:

            [1] x' = F(x, t, r, p)

            [2] M(x, p) x' = F(x, t, r, p)

            [3] M(q, p) u' = F(q, u, t, r, p)
                q' = G(q, u, t, r, p)

        where

            x : states, i.e. [q, u]
            t : time
            r : specified (exogenous) inputs
            p : constants
            q : generalized coordinates
            u : generalized speeds
            M : mass matrix (full or minimum)
            F : right hand side (full or minimum)
            G : right hand side of the kinematical differential equations

        The generated function is of the form F(x, t, p) or F(x, t, r, p)
        depending on whether the system has specified inputs or not.

        Parameters
        ==========
        right_hand_side : SymPy Matrix, shape(n, 1)
            A column vector containing the symbolic expressions for the
            right hand side of the ordinary differential equations. If the
            right hand side has been solved for symbolically then only F is
            required, see form [1]; if not then the mass matrix must also be
            supplied, see forms [2, 3].
        coordinates : sequence of SymPy Functions
            The generalized coordinates. These must be ordered in the same
            order as the rows in M, F, and/or G and be functions of time.
        speeds : sequence of SymPy Functions
            The generalized speeds. These must be ordered in the same order
            as the rows in M, F, and/or G and be functions of time.
        constants : sequence of SymPy Symbols
            All of the constants present in the equations of motion. The
            order does not matter.
        mass_matrix : sympy.Matrix, shape(n, n), optional
            This can be either the "full" mass matrix as in [2] or the
            "minimal" mass matrix as in [3]. The rows and columns must be
            ordered to match the order of the coordinates and speeds. In the
            case of the full mass matrix, the speeds should always be
            ordered before the speeds, i.e. x = [q, u].
        coordinate_derivatives : sympy.Matrix, shape(m, 1), optional
            If the "minimal" mass matrix, form [3], is supplied, then this
            column vector represents the right hand side of the kinematical
            differential equations.
        specifieds : sequence of SymPy Functions
            The specified exogenous inputs to the system. These should be
            functions of time and the order does not matter.
        linear_sys_solver : string or function
            Specify either `numpy` or `scipy` to use the linear solvers
            provided in each package or supply a function that solves a
            linear system Ax=b with the call signature x = solve(A, b). For
            example, if you need to use custom kwargs for the SciPy solver,
            pass in a lambda function that wraps the solver and sets them.
        constants_arg_type : string
            The generated function accepts two different types of arguments
            for the numerical values of the constants: either a ndarray of
            the constants values in the correct order or a dictionary
            mapping the constants symbols to the numerical values. If None,
            this is determined inside of the generated function and can
            cause a significant slow down for performance critical code. If
            you know apriori what arg types you need to support choose
            either ``array`` or ``dictionary``. Note that ``array`` is
            faster than ``dictionary``.
        specifieds_arg_type : string
            The generated function accepts three different types of
            arguments for the numerical values of the specifieds: either a
            ndarray of the specifieds values in the correct order, a
            function that generates the correctly ordered ndarray, or a
            dictionary mapping the specifieds symbols or tuples of thereof
            to floats, ndarrays, or functions. If None, this is determined
            inside of the generated function and can cause a significant
            slow down for performance critical code. If you know apriori
            what arg types you want to support choose either ``array``,
            ``function``, or ``dictionary``. The speed of each, from fast to
            slow, are ``array``, ``function``, ``dictionary``, None.

        Notes
        =====
        The generated function still supports the pre-0.3.0 extra argument
        style, i.e. args = {'constants': ..., 'specified': ...}, but only if
        ``constants_arg_type`` and ``specifieds_arg_type`` are both set to
        None. This functionality is deprecated and will be removed in 0.4.0,
        so it's best to adjust your code to support the new argument types.
        See the docstring for the generated function for more info on the
        new style of arguments.

        """

        self.right_hand_side = right_hand_side
        self.coordinates = coordinates
        self.speeds = speeds
        self.constants = constants
        self.mass_matrix = mass_matrix
        self.coordinate_derivatives = coordinate_derivatives
        self.specifieds = specifieds
        self.linear_sys_solver = linear_sys_solver
        self.constants_arg_type = constants_arg_type
        self.specifieds_arg_type = specifieds_arg_type

        # As the order of the constants and specifieds arguments if not
        # important, allow Sets to be used as input. However, the order must be
        # maintained and converted to a Sequence.
        if constants is not None and not isinstance(constants, Sequence):
            self.constants = tuple(constants)
        if specifieds is not None and not isinstance(specifieds, Sequence):
            self.specifieds = tuple(specifieds)

        self.system_type = self._deduce_system_type(
            mass_matrix=mass_matrix,
            coordinate_derivatives=coordinate_derivatives)

        self.num_coordinates = len(coordinates)
        self.num_speeds = len(speeds)
        self.num_states = self.num_coordinates + self.num_speeds
        self.num_constants = len(constants)

        if self.specifieds is None:
            self.num_specifieds = 0
            self.specifieds_arg_type = None
        else:
            self.num_specifieds = len(specifieds)

        # These are pre-allocated storage for the numerical values used in
        # some of the rhs() evaluations.
        self._constants_values = np.empty(self.num_constants)
        self._specifieds_values = np.empty(self.num_specifieds)

        self._check_system_consitency()

Example 18

Project: pymc3
Source File: advi_minibatch.py
View license
def advi_minibatch(vars=None, start=None, model=None, n=5000, n_mcsamples=1,
                   minibatch_RVs=None, minibatch_tensors=None,
                   minibatches=None, local_RVs=None, observed_RVs=None,
                   encoder_params=[], total_size=None, optimizer=None,
                   learning_rate=.001, epsilon=.1, random_seed=None,
                   verbose=1):
    """Perform mini-batch ADVI.

    This function implements a mini-batch ADVI with the meanfield
    approximation. Autoencoding variational inference is also supported.

    The log probability terms for mini-batches, corresponding to RVs in
    minibatch_RVs, are scaled to (total_size) / (the number of samples in each
    mini-batch), where total_size is an argument for the total data size.

    minibatch_tensors is a list of tensors (can be shared variables) to which
    mini-batch samples are set during the optimization. In most cases, these
    tensors are observations for RVs in the model.

    local_RVs and observed_RVs are used for autoencoding variational Bayes.
    Both of these RVs are associated with each of given samples.
    The difference is that local_RVs are unkown and their posterior
    distributions are approximated.

    local_RVs are Ordered dict, whose keys and values are RVs and a tuple of
    two objects. The first is the theano expression of variational parameters
    (mean and log of std) of the approximate posterior, which are encoded from
    given samples by an arbitrary deterministic function, e.g., MLP. The other
    one is a scaling constant to be multiplied to the log probability term
    corresponding to the RV.

    observed_RVs are also Ordered dict with RVs as the keys, but whose values
    are only the scaling constant as in local_RVs. In this case, total_size is
    ignored.

    If local_RVs is None (thus not using autoencoder), the following two
    settings are equivalent:

    - observed_RVs=OrderedDict([(rv, total_size / minibatch_size)])
    - minibatch_RVs=[rv], total_size=total_size

    where minibatch_size is minibatch_tensors[0].shape[0].

    The variational parameters and the parameters of the autoencoder are
    simultaneously optimized with given optimizer, which is a function that
    returns a dictionary of parameter updates as provided to Theano function.
    See the docstring of pymc3.variational.advi().

    Parameters
    ----------
    vars : object
        List of random variables. If None, variational posteriors (normal
        distribution) are fit for all RVs in the given model.
    start : Dict or None
        Initial values of parameters (variational means).
    model : Model
        Probabilistic model.
    n : int
        Number of interations updating parameters.
    n_mcsamples : int
        Number of Monte Carlo samples to approximate ELBO.
    minibatch_RVs : list of ObservedRVs
        Random variables in the model for which mini-batch tensors are set.
        When this argument is given, both of arguments local_RVs and 
        observed_RVs must be None.
    minibatch_tensors : list of (tensors or shared variables)
        Tensors used to create ObservedRVs in minibatch_RVs.
    minibatches : generator of list
        Generates a set of minibatches when calling next().
        The length of the returned list must be the same with the number of
        random variables in `minibatch_tensors`.
    total_size : int
        Total size of training samples. This is used to appropriately scale the
        log likelihood terms corresponding to mini-batches in ELBO.
    local_RVs : Ordered dict
        Include encoded variational parameters and a scaling constant for
        the corresponding RV. See the above description.
    observed_RVs : Ordered dict
        Include a scaling constant for the corresponding RV. See the above
        description
    encoder_params : list of theano shared variables
        Parameters of encoder.
    optimizer : (loss, list of shared variables) -> dict or OrderedDict
        A function that returns parameter updates given loss and shared
        variables of parameters. If :code:`None` (default), a default
        Adagrad optimizer is used with parameters :code:`learning_rate`
        and :code:`epsilon` below.
    learning_rate: float
        Base learning rate for adagrad. This parameter is ignored when
        an optimizer is given.
    epsilon : float
        Offset in denominator of the scale of learning rate in Adagrad.
        This parameter is ignored when an optimizer is given.
    random_seed : int
        Seed to initialize random state.

    Returns
    -------
    ADVIFit
        Named tuple, which includes 'means', 'stds', and 'elbo_vals'.
    """
    theano.config.compute_test_value = 'ignore'

    model = modelcontext(model)
    vars = inputvars(vars if vars is not None else model.vars)
    start = start if start is not None else model.test_point
    check_discrete_rvs(vars)
    _check_minibatches(minibatch_tensors, minibatches)

    # Prepare optimizer
    if optimizer is None:
        optimizer = adagrad_optimizer(learning_rate, epsilon)

    # For backward compatibility in how input arguments are given
    local_RVs, observed_RVs = _get_rvss(minibatch_RVs, local_RVs, observed_RVs,
                                        minibatch_tensors, total_size)

    # Replace local_RVs with transformed variables
    ds = model.deterministics

    def get_transformed(v):
        if v in ds:
            return v.transformed
        return v
    local_RVs = OrderedDict(
        [(get_transformed(v), (uw, s)) for v, (uw, s) in local_RVs.items()]
    )

    # Get global variables
    global_RVs = list(set(vars) - set(list(local_RVs) + list(observed_RVs)))

    # Ordering for concatenation of random variables
    global_order = ArrayOrdering([v for v in global_RVs])
    local_order = ArrayOrdering([v for v in local_RVs])

    # ELBO wrt variational parameters
    inarray_g, uw_g, replace_g = _join_global_RVs(global_RVs, global_order)
    inarray_l, uw_l, replace_l = _join_local_RVs(local_RVs, local_order)
    logpt = _make_logpt(global_RVs, local_RVs, observed_RVs, model)
    replace = replace_g
    replace.update(replace_l)
    logp = theano.clone(logpt, replace, strict=False)
    elbo = _elbo_t(logp, uw_g, uw_l, inarray_g, inarray_l,
                   n_mcsamples, random_seed)
    del logpt

    # Replacements tensors of variational parameters in the graph
    replaces = dict()

    # Variational parameters for global RVs
    if 0 < len(global_RVs):
        uw_global_shared, bij = _init_uw_global_shared(start, global_RVs,
                                                       global_order)
        replaces.update({uw_g: uw_global_shared})

    # Variational parameters for local RVs, encoded from samples in
    # mini-batches
    if 0 < len(local_RVs):
        uws = [uw for _, (uw, _) in local_RVs.items()]
        uw_local_encoded = tt.concatenate([uw[0].ravel() for uw in uws] +
                                          [uw[1].ravel() for uw in uws])
        replaces.update({uw_l: uw_local_encoded})

    # Replace tensors of variational parameters in ELBO
    elbo = theano.clone(elbo, OrderedDict(replaces), strict=False)

    # Replace input shared variables with tensors
    def is_shared(t):
        return isinstance(t, theano.compile.sharedvalue.SharedVariable)
    tensors = [(t.type() if is_shared(t) else t) for t in minibatch_tensors]
    updates = OrderedDict(
        {t: t_ for t, t_ in zip(minibatch_tensors, tensors) if is_shared(t)}
    )
    elbo = theano.clone(elbo, updates, strict=False)

    # Create parameter update function used in the training loop
    params = encoder_params
    if 0 < len(global_RVs):
        params += [uw_global_shared]
    updates = OrderedDict(optimizer(loss=-1 * elbo, param=params))
    f = theano.function(tensors, elbo, updates=updates)

    # Optimization loop
    elbos = np.empty(n)
    for i in range(n):
        e = f(*next(minibatches))
        elbos[i] = e
        if verbose and not i % (n // 10):
            if not i:
                print('Iteration {0} [{1}%]: ELBO = {2}'.format(
                    i, 100 * i // n, e.round(2)))
            else:
                avg_elbo = elbos[i - n // 10:i].mean()
                print('Iteration {0} [{1}%]: Average ELBO = {2}'.format(
                    i, 100 * i // n, avg_elbo.round(2)))

    if verbose:
        print('Finished [100%]: ELBO = {}'.format(elbos[-1].round(2)))

    # Variational parameters of global RVs
    if 0 < len(global_RVs):
        l = int(uw_global_shared.get_value(borrow=True).size / 2)
        u = bij.rmap(uw_global_shared.get_value(borrow=True)[:l])
        w = bij.rmap(uw_global_shared.get_value(borrow=True)[l:])
        # w is in log space
        for var in w.keys():
            w[var] = np.exp(w[var])
    else:
        u = dict()
        w = dict()

    return ADVIFit(u, w, elbos)

Example 19

Project: pyqtgraph
Source File: GLSurfacePlotItem.py
View license
    def setData(self, x=None, y=None, z=None, colors=None):
        """
        Update the data in this surface plot. 
        
        ==============  =====================================================================
        **Arguments:**
        x,y             1D arrays of values specifying the x,y positions of vertexes in the
                        grid. If these are omitted, then the values will be assumed to be
                        integers.
        z               2D array of height values for each grid vertex.
        colors          (width, height, 4) array of vertex colors.
        ==============  =====================================================================
        
        All arguments are optional.
        
        Note that if vertex positions are updated, the normal vectors for each triangle must 
        be recomputed. This is somewhat expensive if the surface was initialized with smooth=False
        and very expensive if smooth=True. For faster performance, initialize with 
        computeNormals=False and use per-vertex colors or a normal-independent shader program.
        """
        if x is not None:
            if self._x is None or len(x) != len(self._x):
                self._vertexes = None
            self._x = x
        
        if y is not None:
            if self._y is None or len(y) != len(self._y):
                self._vertexes = None
            self._y = y
        
        if z is not None:
            #if self._x is None:
                #self._x = np.arange(z.shape[0])
                #self._vertexes = None
            #if self._y is None:
                #self._y = np.arange(z.shape[1])
                #self._vertexes = None
                
            if self._x is not None and z.shape[0] != len(self._x):
                raise Exception('Z values must have shape (len(x), len(y))')
            if self._y is not None and z.shape[1] != len(self._y):
                raise Exception('Z values must have shape (len(x), len(y))')
            self._z = z
            if self._vertexes is not None and self._z.shape != self._vertexes.shape[:2]:
                self._vertexes = None
        
        if colors is not None:
            self._colors = colors
            self._meshdata.setVertexColors(colors)
        
        if self._z is None:
            return
        
        updateMesh = False
        newVertexes = False
        
        ## Generate vertex and face array
        if self._vertexes is None:
            newVertexes = True
            self._vertexes = np.empty((self._z.shape[0], self._z.shape[1], 3), dtype=float)
            self.generateFaces()
            self._meshdata.setFaces(self._faces)
            updateMesh = True
        
        ## Copy x, y, z data into vertex array
        if newVertexes or x is not None:
            if x is None:
                if self._x is None:
                    x = np.arange(self._z.shape[0])
                else:
                    x = self._x
            self._vertexes[:, :, 0] = x.reshape(len(x), 1)
            updateMesh = True
        
        if newVertexes or y is not None:
            if y is None:
                if self._y is None:
                    y = np.arange(self._z.shape[1])
                else:
                    y = self._y
            self._vertexes[:, :, 1] = y.reshape(1, len(y))
            updateMesh = True
        
        if newVertexes or z is not None:
            self._vertexes[...,2] = self._z
            updateMesh = True
        
        ## Update MeshData
        if updateMesh:
            self._meshdata.setVertexes(self._vertexes.reshape(self._vertexes.shape[0]*self._vertexes.shape[1], 3))
            self.meshDataChanged()

Example 20

Project: python-control
Source File: frdata.py
View license
    def __init__(self, *args, **kwargs):
        """Construct an FRD object

        The default constructor is FRD(d, w), where w is an iterable of
        frequency points, and d is the matching frequency data.

        If d is a single list, 1d array, or tuple, a SISO system description
        is assumed. d can also be

        To call the copy constructor, call FRD(sys), where sys is a
        FRD object.

        To construct frequency response data for an existing LTI
        object, other than an FRD, call FRD(sys, omega)

        """
        smooth = kwargs.get('smooth', False)

        if len(args) == 2:
            if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
                # not an FRD, but still a system, second argument should be
                # the frequency range
                otherlti = args[0]
                self.omega = array(args[1], dtype=float)
                self.omega.sort()
                numfreq = len(self.omega)

                # calculate frequency response at my points
                self.fresp = empty(
                    (otherlti.outputs, otherlti.inputs, numfreq),
                    dtype=complex)
                for k, w in enumerate(self.omega):
                    self.fresp[:, :, k] = otherlti.evalfr(w)

            else:
                # The user provided a response and a freq vector
                self.fresp = array(args[0], dtype=complex)
                if len(self.fresp.shape) == 1:
                    self.fresp = self.fresp.reshape(1, 1, len(args[0]))
                self.omega = array(args[1], dtype=float)
                if len(self.fresp.shape) != 3 or \
                        self.fresp.shape[-1] != self.omega.shape[-1] or \
                        len(self.omega.shape) != 1:
                    raise TypeError(
                        "The frequency data constructor needs a 1-d or 3-d"
                        " response data array and a matching frequency vector"
                        " size")

        elif len(args) == 1:
            # Use the copy constructor.
            if not isinstance(args[0], FRD):
                raise TypeError(
                    "The one-argument constructor can only take in"
                    " an FRD object.  Received %s." % type(args[0]))
            self.omega = args[0].omega
            self.fresp = args[0].fresp
        else:
            raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))

        # create interpolation functions
        if smooth:
            self.ifunc = empty((self.fresp.shape[0], self.fresp.shape[1]),
                               dtype=tuple)
            for i in range(self.fresp.shape[0]):
                for j in range(self.fresp.shape[1]):
                    self.ifunc[i,j],u = splprep(
                        u=self.omega, x=[real(self.fresp[i, j, :]),
                                         imag(self.fresp[i, j, :])],
                        w=1.0/(absolute(self.fresp[i, j, :])+0.001), s=0.0)
        else:
            self.ifunc = None
        LTI.__init__(self, self.fresp.shape[1], self.fresp.shape[0])

Example 21

Project: mpop
Source File: HRWimage.py
View license
def HRWimage( HRW_data, obj_area, hrw_channels=None, min_correlation=None, cloud_type=None, style='barbs', \
                        barb_length=None, color_mode='channel', pivot='middle', legend=True, legend_loc=3):

    """Create a PIL image from high resolution wind data
       required input:
         HRW data [HRW_class]:      HRW_class instant containing the data, see mpop/satin/nwcsaf_hrw_hdf.py
         obj_area [area_class]:     instant of area class, returned by area_def
       optional input
       hrw_channels [string array]: giving the channels that are used derive the HRW vectors
                                    e.g. hrw_channels['HRV','WV_073']
       min_correlation [int]:       minimum correlation of tracking, if below arrow is not shown
       cloud_type [int array]:      cloud types of the wind vectors, e.g. cloud_type=[8,10,11]
       style [string]:              different styles of plotting
                                    style='barbs' or style='5min_displacement' or style='15min_displacement'
       color_mode [string]:         choose color of the wind symbols, possible choises:
                                    color_mode='channel'      -> one color per SEVIRI channel used to derive HRW
                                    color_mode='pressure'     -> colorcoded cloud top pressure
                                    color_mode='temperature'  -> colorcoded cloud top temperature
                                    color_mode='cloud_type'   -> NWC-SAF cloud types
                                    color_mode='correlation'  80 ... 100
                                    color_mode='conf_nwp'      70 ... 100
                                    color_mode='conf_no_nwp'   70 ... 100
       pivot [string]:              position of the barb, e.g. pivot='middle' == center of barb at origin
       legend [True or False] :     show legend or not 
       legend_loc [string or int]:  location of the legend 
                                    upper right    1
                                    upper left     2
                                    lower left     3
                                    lower right    4
                                    right          5
                                    center left    6
                                    center right   7
                                    lower center   8
                                    upper center   9
                                    center         10
                                    best
    """

    #if min_correlation != None:
    #    print "    filter for min_correlation = ", min_correlation
    #    inds = where(HRW_data.correlation > min_correlation)
    #    HRW_data = HRW_data.subset(inds)
     
    print "... create HRWimage, color_mode = ", color_mode

    # get a empty figure with transparent background, no axis and no margins outside the diagram
    fig, ax = prepare_figure(obj_area)

    # define arrow properties 
    head_width  = 0.006 * min(obj_area.x_size,obj_area.x_size)
    head_length = 2 * head_width
    m_per_s_to_knots = 1.944

    #barb_length = 0.008 * min(obj_area.x_size,obj_area.x_size)
    
    if barb_length == None:
        n_winds = len(HRW_data.wind_id)
        if n_winds < 300:
            barb_length = 5.68
        elif n_winds < 500:
            barb_length = 5.43
        elif n_winds < 700:
            barb_length = 5.18
        elif n_winds < 900:
            barb_length = 4.68
        else:          
            barb_length = 4.00
    print "barb_length", barb_length

    if color_mode == 'channel':
        classes = ('HRV',          'VIS008 ', 'WV_062 ',   'WV_073 ',   'IR_120 ')
        colors   = ['mediumorchid', 'red',     'limegreen', 'darkgreen', 'darkturquoise']
    elif color_mode == 'pressure':
        classes = ['<200hPa',  '200-300hPa','300-400hPa','400-500hPa','500-600hPa','600-700hPa','700-800hPa', '800-900hPa','>900hPa']
        colors   = ['darksalmon', 'red'     ,'darkorange','yellow'    ,'lime'      ,'seagreen',  'deepskyblue','blue',      'mediumorchid']
        classes = tuple(['CTP '+cl for cl in classes])
    elif color_mode == 'cloud_type' or color_mode == 'cloudtype':
        classes=['non-processed','cloud free land', 'cloud free sea', 'land snow', 'sea ice',\
                 'very low cum.', 'very low', 'low cum.', 'low', 'med cum.', 'med', 'high cum.', 'high', 'very high cum.', 'very high', \
                 'sem. thin', 'sem. med.', 'sem. thick', 'sem. above', 'broken', 'undefined']

        colors = empty( (len(classes),3), dtype=int )
        colors[ 0,:] = [100, 100, 100]
        colors[ 1,:] = [  0, 120,   0]
        colors[ 2,:] = [  0,   0,   0]
        colors[ 3,:] = [250, 190, 250]
        colors[ 4,:] = [220, 160, 220]
        colors[ 5,:] = [255, 150,   0]
        colors[ 6,:] = [255, 100,   0]
        colors[ 7,:] = [255, 220,   0]
        colors[ 8,:] = [255, 180,   0]
        colors[ 9,:] = [255, 255, 140]
        colors[10,:] = [240, 240,   0]
        colors[11,:] = [250, 240, 200]
        colors[12,:] = [215, 215, 150]
        colors[13,:] = [255, 255, 255]
        colors[14,:] = [230, 230, 230]
        colors[15,:] = [  0,  80, 215]
        colors[16,:] = [  0, 180, 230]
        colors[17,:] = [  0, 240, 240]
        colors[18,:] = [ 90, 200, 160]
        colors[19,:] = [200,   0, 200]
        colors[20,:] = [ 95,  60,  30]
        colors = colors/255.
    elif color_mode in ['correlation','conf_nwp','conf_no_nwp']:
        classes  = ['<70',    '<75',     '<80',    '<85',   '<90',  '<95',   '>95' ]
        colors   = ['indigo', 'darkred', 'red','darkorange','gold', 'lime', 'green']
        classes = tuple([color_mode+' '+cl for cl in classes])
    else:
          print "*** Error in HRW_streamplot (mpop/imageo/HRWimage.py)"
          print "    unknown color_mode"
          quit()

    for wid in range(len(HRW_data.wind_id)):

        if color_mode == 'channel':

            if HRW_data.channel[wid].find('HRV') != -1: # HRV
                barbcolor = colors[0]
            elif HRW_data.channel[wid].find('VIS008') != -1: #  0.8 micro m
                barbcolor = colors[1]
            elif HRW_data.channel[wid].find('WV_062') != -1: #  6.2 micro m
                barbcolor = colors[2]
            elif HRW_data.channel[wid].find('WV_073') != -1: #  7.3 micro m
                barbcolor = colors[3]
            elif HRW_data.channel[wid].find('IR_120') != -1: # 12.0 micro m
                barbcolor = colors[4]

        elif color_mode == 'pressure':

            if HRW_data.pressure[wid] < 20000:
                barbcolor = colors[0]
            elif HRW_data.pressure[wid] < 30000:
                barbcolor = colors[1]
            elif HRW_data.pressure[wid] < 40000:
                barbcolor = colors[2]
            elif HRW_data.pressure[wid] < 50000:
                barbcolor = colors[3] 
            elif HRW_data.pressure[wid] < 60000:
                barbcolor = colors[4]
            elif HRW_data.pressure[wid] < 70000:
                barbcolor = colors[5]
            elif HRW_data.pressure[wid] < 80000:
                barbcolor = colors[6]
            elif HRW_data.pressure[wid] < 90000:
                barbcolor = colors[7]
            else:
                barbcolor = colors[8]

        elif color_mode == 'cloud_type' or color_mode == 'cloudtype':

            barbcolor = list(colors[HRW_data.cloud_type[wid], :])

        elif color_mode in ['correlation','conf_nwp','conf_no_nwp']:
            if color_mode == 'correlation':
                cdata = HRW_data.correlation
            elif color_mode == 'conf_nwp':
                cdata = HRW_data.conf_nwp
            elif color_mode == 'conf_no_nwp':
                cdata = HRW_data.conf_no_nwp

            if cdata[wid] < 70:
                barbcolor = colors[0]
            elif cdata[wid] < 75:
                barbcolor = colors[1]
            elif cdata[wid] < 80:
                barbcolor = colors[2]
            elif cdata[wid] < 85:
                barbcolor = colors[3] 
            elif cdata[wid] < 90:
                barbcolor = colors[4]
            elif cdata[wid] < 95:
                barbcolor = colors[5]
            else:
                barbcolor = colors[6]
        else:
              print "*** Error in HRW_streamplot (mpop/imageo/HRWimage.py)"
              print "    unknown color_mode"
              quit()
      
        x0, y0 = obj_area.get_xy_from_lonlat( HRW_data.lon[wid], HRW_data.lat[wid], outside_error=False) #, return_int=True

        u = HRW_data.wind_speed[wid] * -1 * sin(radians(HRW_data.wind_direction[wid])) 
        v = HRW_data.wind_speed[wid] * -1 * cos(radians(HRW_data.wind_direction[wid]))

        #print '%6s %3d %10.7f %10.7f %7.2f %7.1f %8.1f %10s' % (HRW_data.channel[wid], HRW_data.wind_id[wid], \
        #                                                        HRW_data.lon[wid], HRW_data.lat[wid], \
        #                                                        HRW_data.wind_speed[wid]*m_per_s_to_knots, \
        #                                                        HRW_data.wind_direction[wid], HRW_data.pressure[wid], barbcolor)


        if style == 'barbs':
            u = HRW_data.wind_speed[wid] * -1 * sin(radians(HRW_data.wind_direction[wid])) * m_per_s_to_knots
            v = HRW_data.wind_speed[wid] * -1 * cos(radians(HRW_data.wind_direction[wid])) * m_per_s_to_knots
            ax.barbs(x0, obj_area.y_size - y0, u * m_per_s_to_knots, v * m_per_s_to_knots, length = barb_length, pivot='middle', barbcolor=barbcolor)

        elif style == '5min_displacement' or style == '15min_displacement':
            if style == '5min_displacement':
                t_in_s =  5*60
            elif style == '15min_displacement':
                t_in_s = 15*60
            dx = u * t_in_s / obj_area.pixel_size_x
            dy = v * t_in_s / obj_area.pixel_size_y
            ax.arrow(x0, y0, dx, dy, head_width = head_width, head_length = head_length, fc=barbcolor, ec=barbcolor)

    if legend:

        rcParams['legend.handlelength'] = 0
        rcParams['legend.numpoints'] = 1

        # create blank rectangle
        rec = Rectangle((0, 0), 0, 0, fc="w", fill=False, edgecolor='none', linewidth=0)

        ##  *fontsize*: [size in points | 'xx-small' | 'x-small' | 'small' |
        ##              'medium' | 'large' | 'x-large' | 'xx-large']


        alpha=1.0
        bbox={'facecolor':'white', 'alpha':alpha, 'pad':10}

        print "... add legend: color is a function of ",  color_mode
        
        recs = empty( len(classes), dtype=object)
        recs[:] = rec 

        #if color_mode == 'pressure':
        #    recs = [rec, rec, rec, rec, rec, rec, rec, rec, rec]
        #if color_mode == 'channel':
        #    recs = [rec, rec, rec, rec, rec]
        #if color_mode in ['correlation','conf_nwp','conf_no_nwp']:
        #    recs = [rec, rec, rec, rec, rec, rec, rec]

        size=12
        if color_mode=='cloud_type':
            size=10

        leg = ax.legend(recs, classes, loc=legend_loc, prop={'size':size})

        for color,text in zip(colors,leg.get_texts()):
            text.set_color(color)


    return fig2img ( fig )

Example 22

Project: mpop
Source File: HRWimage.py
View license
def HRW_2dfield( HRW_data, obj_area, interpol_method=None, hrw_channels=None, min_correlation=None, level=''):

    print "... calculate 2d wind field (HRW_2dfield)"

    if min_correlation != None:
        print "    filter for min_correlation = ", min_correlation
        inds = where(HRW_data.correlation > min_correlation)
        HRW_data.subset(inds)

    xx, yy = obj_area.get_xy_from_lonlat( HRW_data.lon, HRW_data.lat, outside_error=False, return_int=False) #, return_int=True

    yy = obj_area.y_size - yy  

    uu = - HRW_data.wind_speed * sin(radians(HRW_data.wind_direction))
    vv = - HRW_data.wind_speed * cos(radians(HRW_data.wind_direction))

    # get rid of all vectors outside of the field 
    index = nonzero(xx)
    xx = xx[index]
    yy = yy[index]
    uu = uu[index]
    vv = vv[index]

    points = transpose(append([xx], [yy], axis=0))
    #print type(uu), uu.shape
    #print type(points), points.shape
    #print points[0], yy[0], xx[0]
    #print uu[0]

    nx = obj_area.x_size
    ny = obj_area.y_size

    x2 = arange(nx)
    y2 = (ny-1) - arange(ny)

    grid_x, grid_y = meshgrid(x2, y2)
    
    if interpol_method == None:
        # we need at least 2 winds to interpolate 
        if uu.size < 4:  
            print "*** Warning, not wnough wind data available, n_winds = ", uu.size
            fake = empty(grid_x.shape)
            fake[:,:] = nan
            HRW_data.interpol_method = None
            return fake, fake

        elif uu.size < 50:
            interpol_method = "RBF"

        else:
            interpol_method = "linear + nearest"
     
    #interpol_method = "nearest"
    #interpol_method = "cubic + nearest" # might cause unrealistic overshoots
    #interpol_method = "kriging"
    #interpol_method = "..."

    print "... min windspeed (org data): ", HRW_data.wind_speed.min()
    print "... max windspeed (org data): ", HRW_data.wind_speed.max()

    for i_iteration in [0,1]:

        if interpol_method == "nearest":

            print '... fill with nearest neighbour'
            # griddata, see http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
            grid_u1x = griddata(points, uu, (grid_x, grid_y), method='nearest')
            grid_v1x = griddata(points, vv, (grid_x, grid_y), method='nearest')

        elif interpol_method == "RBF":

            print '... inter- and extrapolation using radial basis functions'
            # https://www.youtube.com/watch?v=_cJLVhdj0j4
            print "... start Rbf"
            from scipy.interpolate import Rbf
            # rbfu = Rbf(xx, yy, uu, epsilon=0.1) #
            rbfu = Rbf(xx, yy, uu, epsilon=0.2)
            grid_u1x = rbfu(grid_x, grid_y)
            rbfv = Rbf(xx, yy, vv, epsilon=0.1) #
            grid_v1x = rbfv(grid_x, grid_y)
            print "... finish Rbf"
            # !very! slow for a large number of observations 

        elif interpol_method == "linear + nearest" or interpol_method == "cubic + nearest":

            if interpol_method == "linear + nearest":
                print '... calculate linear interpolation'
                # griddata, see http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html
                grid_u1 = griddata(points, uu, (grid_x, grid_y), method='linear')
                grid_v1 = griddata(points, vv, (grid_x, grid_y), method='linear')
            elif interpol_method == "cubic + nearest":
                # smoother, but can cause unrealistic overshoots
                print '... calculate cubic interpolation'
                grid_u1 = griddata(points, uu, (grid_x, grid_y), method='cubic')
                grid_v1 = griddata(points, vv, (grid_x, grid_y), method='cubic')
            else:
                print "*** Error in mpop/imageo/HRWimage.py"
                print "    unknown interpolation method: ", interpol_method
                quit()

            if 1==1:
                # use faster function to extrapolate with closest neighbour
                print "... fill outside area with closest value"
                grid_u1x = fill_with_closest_pixel(grid_u1, invalid=None) 
                grid_v1x = fill_with_closest_pixel(grid_v1, invalid=None) 
            else:
                # use griddata to extrapolate with closest neighbour
                points2 = transpose(append([grid_x.flatten()], [grid_y.flatten()], axis=0))
                print type(grid_x.flatten()), grid_x.flatten().shape
                print type(points2), points2.shape
                mask = ~isnan(grid_v1.flatten())
                inds = where(mask)[0]
                grid_u1x = griddata(points2[inds], grid_u1.flatten()[inds], (grid_x, grid_y), method='nearest')
                grid_v1x = griddata(points2[inds], grid_v1.flatten()[inds], (grid_x, grid_y), method='nearest')

            if 1==0:
                # add othermost points as additional data
                y_add = [0,    0, ny-1, ny-1]
                x_add = [0, nx-1,    0, nx-1]
                for (i,j) in zip(x_add,y_add):
                    uu = append(uu, grid_u0[i,j])
                    vv = append(vv, grid_v0[i,j])
                xx = append(xx, x_add)
                yy = append(yy, y_add)
                points = transpose(append([yy], [xx], axis=0))

                print 'calc extent1'
                grid_u1e = griddata(points, uu, (grid_x, grid_y), method='linear')
                grid_v1e = griddata(points, vv, (grid_x, grid_y), method='linear')

        else:
            print "*** Error in mpop/imageo/HRWimage.py"
            print "    unknown interpol_method", interpol_method
            quit()

        ##http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
        ##http://stackoverflow.com/questions/3526514/problem-with-2d-interpolation-in-scipy-non-rectangular-grid
        #print "SmoothBivariateSpline:"
        #from scipy.interpolate import SmoothBivariateSpline
        #fitu = SmoothBivariateSpline( xx, yy, uu, s=1000) # , kx=3, ky=3, s = smooth * z2sum m 
        #from numpy import empty 
        #grid_u_SBP = empty(grid_x.shape)
        #for i in range(0,nx-1):       # starting upper right going down
        #    for j in range(0,ny-1):   # starting lower right going right
        #        #print i,j
        #        grid_u_SBP[j,i] = fitu(j,i)

        #grid_u_SBP = np.array([k.predict([x,y]) for x,y in zip(np.ravel(grid_x), np.ravel(grid_y))])
        #grid_u_SBP = grid_u_SBP.reshape(grid_x.shape)

        ##print x2
        ##print y2
        #grid_u_SBP = fitu(x2,y2)
        ##print "grid_u_SBP.shape", grid_u_SBP.shape
        ###print grid_u_SBP
        #print "End SmoothBivariateSpline:"

        #print "bisplrep:"
        #from scipy import interpolate
        #tck = interpolate.bisplrep(xx, yy, uu)
        #grid_u_BSR = interpolate.bisplev(grid_x[:,0], grid_y[0,:], tck)
        #print grid_u_BSR.shape
        #print "bisplrep"
        #print "grid_v1x.shape", grid_v1x.shape

        extent=(0,nx,0,ny)
        origin='lower'
        origin='upper'
        origin=None 

        # show different stages of 2d inter- and extra-polation 
        if 1==0:
            print 'make matplotlib.pyplot'
            import matplotlib.pyplot as plt
            vmin=-10
            vmax=10
            fig = plt.figure()
            plt.subplot(221)
            plt.title('u '+interpol_method)
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_u1x, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.subplot(222)
            plt.title('v '+interpol_method)
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_v1x, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()

            # standard calculation for comparison 
            print '... calculate linear interpolation'
            grid_u1 = griddata(points, uu, (grid_x, grid_y), method='linear')
            grid_v1 = griddata(points, vv, (grid_x, grid_y), method='linear')
            grid_u1xx = fill_with_closest_pixel(grid_u1, invalid=None) 
            grid_v1xx = fill_with_closest_pixel(grid_v1, invalid=None) 

            plt.subplot(223)
            plt.title('U Linear+Nearest')
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            plt.imshow(grid_u1xx, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.subplot(224)
            plt.title('V Linear+Nearest') 
            plt.plot(points[:,0], ny-1-points[:,1], 'k.', ms=1)
            #plt.title('Cubic')
            plt.imshow(grid_v1xx, origin=origin, vmin=vmin, vmax=vmax) #, extent=extent
            #plt.colorbar()
            plt.gcf().set_size_inches(6, 6)
            #plt.show()  # does not work with AGG
            tmpfile="test_hrw"+level+".png"
            fig.savefig(tmpfile)
            print "display "+tmpfile+" &"


        if grid_u1x.min() < -150 or grid_v1x.min() < -150 or grid_u1x.max() > 150 or grid_v1x.max() > 150:
            print "*** Warning, numerical instability detected, interpolation method: ", interpol_method
            print "    min u windspeed (u 2dimensional): ", grid_u1x.min()
            print "    min v windspeed (v 2dimensional): ", grid_v1x.min()
            print "    max u windspeed (u 2dimensional): ", grid_u1x.max()
            print "    max v windspeed (v 2dimensional): ", grid_v1x.max()
            interpol_method = "glinear + nearest"
            print "... try another interpolation method: ", interpol_method
        else:
            # (hopefully) numerical stable interpolation, exit the interpolation loop
            break 

    HRW_data.interpol_method = interpol_method

    return grid_u1x, grid_v1x

Example 23

Project: python-qinfer
Source File: clustering.py
View license
def particle_clusters(
        particle_locations, particle_weights=None,
        eps=0.5, min_particles=5, metric='euclidean',
        weighted=False, w_pow=0.5,
        quiet=True
    ):
    """
    Yields an iterator onto tuples ``(cluster_label, cluster_particles)``,
    where ``cluster_label`` is an `int` identifying the cluster (or ``NOISE``
    for the particles lying outside of all clusters), and where
    ``cluster_particles`` is an array of ``dtype`` `bool` specifying the indices
    of all particles in that cluster. That is, particle ``i`` is in the cluster
    if ``cluster_particles[i] == True``.
    """
    
    
    if weighted == True and particle_weights is None:
        raise ValueError("Weights must be specified for weighted clustering.")
        
    # Allocate new arrays to hold the weights and locations.        
    new_weights = np.empty(particle_weights.shape)
    new_locs    = np.empty(particle_locations.shape)
    
    # Calculate and possibly reweight the metric.
    if weighted:
        M = sklearn.metrics.pairwise.pairwise_distances(particle_locations, metric=metric)
        M = metrics.weighted_pairwise_distances(M, particle_weights, w_pow=w_pow)
    
        # Create and run a SciKit-Learn DBSCAN clusterer.
        clusterer = sklearn.cluster.DBSCAN(
            min_samples=min_particles,
            eps=eps,
            metric='precomputed'
        )
        cluster_labels = clusterer.fit_predict(M)
    else:
        clusterer = sklearn.cluster.DBSCAN(
            min_samples=min_particles,
            eps=eps,
            metric=metric
        )
        cluster_labels = clusterer.fit_predict(particle_locations)
    
    # Find out how many clusters were identified.
    # Cluster counting logic from:
    # [http://scikit-learn.org/stable/auto_examples/cluster/plot_dbscan.html].
    is_noise = -1 in cluster_labels
    n_clusters = len(set(cluster_labels)) - (1 if is_noise else 0)
    
    # If more than 10% of the particles were labeled as NOISE,
    # warn.
    n_noise = np.sum(cluster_labels == -1)
    if n_noise / particle_weights.shape[0] >= 0.1:
        warnings.warn("More than 10% of the particles were classified as NOISE. Consider increasing the neighborhood size ``eps``.", ResamplerWarning)
    
    # Print debugging info.
    if not quiet:
        print("[Clustering] DBSCAN identified {} cluster{}. "\
              "{} particles identified as NOISE.".format(
                  n_clusters, "s" if n_clusters > 1 else "", n_noise
              ))
        
    # Loop over clusters, calling the secondary resampler for each.
    # The loop should include -1 if noise was found.
    for idx_cluster in range(-1 if is_noise else 0, n_clusters):
        # Grab a boolean array identifying the particles in a  particular
        # cluster.
        this_cluster = cluster_labels == idx_cluster
        
        yield idx_cluster, this_cluster

Example 24

Project: QuantEcon.py
Source File: ddp.py
View license
def backward_induction(ddp, T, v_term=None):
    r"""
    Solve by backward induction a :math:`T`-period finite horizon
    discrete dynamic program with stationary reward and transition
    probability functions :math:`r` and :math:`q` and discount factor
    :math:`\beta \in [0, 1]`.

    The optimal value functions :math:`v^*_0, \ldots, v^*_T` and policy
    functions :math:`\sigma^*_0, \ldots, \sigma^*_{T-1}` are obtained by
    :math:`v^*_T = v_T`, and

    .. math::

        v^*_{t-1}(s) = \max_{a \in A(s)} r(s, a) +
            \beta \sum_{s' \in S} q(s'|s, a) v^*_t(s')
            \quad (s \in S)

    and

    .. math::

        \sigma^*_{t-1}(s) \in \operatorname*{arg\,max}_{a \in A(s)}
            r(s, a) + \beta \sum_{s' \in S} q(s'|s, a) v^*_t(s')
            \quad (s \in S)

    for :math:`t = T, \ldots, 1`, where the terminal value function
    :math:`v_T` is exogenously given.

    Parameters
    ----------
    ddp : DiscreteDP
        DiscreteDP instance storing reward array `R`, transition
        probability array `Q`, and discount factor `beta`.

    T : scalar(int)
        Number of decision periods.

    v_term : array_like(float, ndim=1), optional(default=None)
        Terminal value function, of length equal to n (the number of
        states). If None, it defaults to the vector of zeros.

    Returns
    -------
    vs : ndarray(float, ndim=2)
        Array of shape (T+1, n) where `vs[t]` contains the optimal
        value function at period `t = 0, ..., T`.

    sigmas : ndarray(int, ndim=2)
        Array of shape (T, n) where `sigmas[t]` contains the optimal
        policy function at period `t = 0, ..., T-1`.

    """
    n = ddp.num_states
    vs = np.empty((T+1, n))
    sigmas = np.empty((T, n), dtype=int)

    if v_term is None:
        v_term = np.zeros(n)
    vs[T, :] = v_term

    for t in range(T, 0, -1):
        ddp.bellman_operator(vs[t, :], Tv=vs[t-1, :], sigma=sigmas[t-1, :])

    return vs, sigmas

Example 25

Project: zipline
Source File: yahoo.py
View license
def yahoo_equities(symbols, start=None, end=None):
    """Create a data bundle ingest function from a set of symbols loaded from
    yahoo.

    Parameters
    ----------
    symbols : iterable[str]
        The ticker symbols to load data for.
    start : datetime, optional
        The start date to query for. By default this pulls the full history
        for the calendar.
    end : datetime, optional
        The end date to query for. By default this pulls the full history
        for the calendar.

    Returns
    -------
    ingest : callable
        The bundle ingest function for the given set of symbols.

    Examples
    --------
    This code should be added to ~/.zipline/extension.py

    .. code-block:: python

       from zipline.data.bundles import yahoo_equities, register

       symbols = (
           'AAPL',
           'IBM',
           'MSFT',
       )
       register('my_bundle', yahoo_equities(symbols))

    Notes
    -----
    The sids for each symbol will be the index into the symbols sequence.
    """
    # strict this in memory so that we can reiterate over it
    symbols = tuple(symbols)

    def ingest(environ,
               asset_db_writer,
               minute_bar_writer,  # unused
               daily_bar_writer,
               adjustment_writer,
               calendar,
               start_session,
               end_session,
               cache,
               show_progress,
               output_dir,
               # pass these as defaults to make them 'nonlocal' in py2
               start=start,
               end=end):
        if start is None:
            start = start_session
        if end is None:
            end = None

        metadata = pd.DataFrame(np.empty(len(symbols), dtype=[
            ('start_date', 'datetime64[ns]'),
            ('end_date', 'datetime64[ns]'),
            ('auto_close_date', 'datetime64[ns]'),
            ('symbol', 'object'),
        ]))

        def _pricing_iter():
            sid = 0
            with maybe_show_progress(
                    symbols,
                    show_progress,
                    label='Downloading Yahoo pricing data: ') as it, \
                    requests.Session() as session:
                for symbol in it:
                    path = _cachpath(symbol, 'ohlcv')
                    try:
                        df = cache[path]
                    except KeyError:
                        df = cache[path] = DataReader(
                            symbol,
                            'yahoo',
                            start,
                            end,
                            session=session,
                        ).sort_index()

                    # the start date is the date of the first trade and
                    # the end date is the date of the last trade
                    start_date = df.index[0]
                    end_date = df.index[-1]
                    # The auto_close date is the day after the last trade.
                    ac_date = end_date + pd.Timedelta(days=1)
                    metadata.iloc[sid] = start_date, end_date, ac_date, symbol

                    df.rename(
                        columns={
                            'Open': 'open',
                            'High': 'high',
                            'Low': 'low',
                            'Close': 'close',
                            'Volume': 'volume',
                        },
                        inplace=True,
                    )
                    yield sid, df
                    sid += 1

        daily_bar_writer.write(_pricing_iter(), show_progress=show_progress)

        symbol_map = pd.Series(metadata.symbol.index, metadata.symbol)

        # Hardcode the exchange to "YAHOO" for all assets and (elsewhere)
        # register "YAHOO" to resolve to the NYSE calendar, because these are
        # all equities and thus can use the NYSE calendar.
        metadata['exchange'] = "YAHOO"
        asset_db_writer.write(equities=metadata)

        adjustments = []
        with maybe_show_progress(
                symbols,
                show_progress,
                label='Downloading Yahoo adjustment data: ') as it, \
                requests.Session() as session:
            for symbol in it:
                path = _cachpath(symbol, 'adjustment')
                try:
                    df = cache[path]
                except KeyError:
                    df = cache[path] = DataReader(
                        symbol,
                        'yahoo-actions',
                        start,
                        end,
                        session=session,
                    ).sort_index()

                df['sid'] = symbol_map[symbol]
                adjustments.append(df)

        adj_df = pd.concat(adjustments)
        adj_df.index.name = 'date'
        adj_df.reset_index(inplace=True)

        splits = adj_df[adj_df.action == 'SPLIT']
        splits = splits.rename(
            columns={'value': 'ratio', 'date': 'effective_date'},
        )
        splits.drop('action', axis=1, inplace=True)

        dividends = adj_df[adj_df.action == 'DIVIDEND']
        dividends = dividends.rename(
            columns={'value': 'amount', 'date': 'ex_date'},
        )
        dividends.drop('action', axis=1, inplace=True)
        # we do not have this data in the yahoo dataset
        dividends['record_date'] = pd.NaT
        dividends['declared_date'] = pd.NaT
        dividends['pay_date'] = pd.NaT

        adjustment_writer.write(splits=splits, dividends=dividends)

    return ingest

Example 26

Project: sharedmem
Source File: sort.py
View license
def argsort(data, out=None, chunksize=None, 
        baseargsort=None, 
        argmerge=None, np=None):
    """
     parallel argsort, like numpy.argsort

     use sizeof(intp) * len(data) as scratch space

     use baseargsort for serial sort 
         ind = baseargsort(data)

     use argmerge to merge
         def argmerge(data, A, B, out):
             ensure data[out] is sorted
             and out[:] = A join B

     TODO: shall try to use the inplace merge mentioned in 
            http://keithschwarz.com/interesting/code/?dir=inplace-merge.
    """
    if baseargsort is None:
        baseargsort = lambda x:x.argsort()

    if argmerge is None:
        argmerge = default_argmerge

    if chunksize is None:
        chunksize = 1024 * 1024 * 16

    if out is None:
        arg1 = numpy.empty(len(data), dtype='intp')
        out = arg1
    else:
        assert out.dtype == numpy.dtype('intp')
        assert len(out) == len(data)
        arg1 = out

    if np is None:
        np = sharedmem.cpu_count()

    if np <= 1 or len(data) < chunksize: 
        out[:] = baseargsort(data)
        return out

    CHK = [slice(i, i + chunksize) for i in range(0, len(data), chunksize)]
    DUMMY = slice(len(data), len(data))
    if len(CHK) % 2: CHK.append(DUMMY)
    with sharedmem.TPool() as pool:
        def work(i):
            C = CHK[i]
            start, stop, step = C.indices(len(data))
            arg1[C] = baseargsort(data[C])
            arg1[C] += start
        pool.map(work, range(len(CHK)))
  
    arg2 = numpy.empty_like(arg1)
  
    flip = 0
    while len(CHK) > 1:
        with sharedmem.TPool() as pool:
            def work(i):
                C1 = CHK[i]
                C2 = CHK[i+1]
                start1, stop1, step1 = C1.indices(len(data))
                start2, stop2, step2 = C2.indices(len(data))
        #        print 'argmerge', start1, stop1, start2, stop2
                assert start2 == stop1
                argmerge(data, arg1[C1], arg1[C2], arg2[start1:stop2])
                return slice(start1, stop2)
            CHK = pool.map(work, range(0, len(CHK), 2))
            arg1, arg2 = arg2, arg1
            flip = flip + 1
        if len(CHK) == 1: break
        if len(CHK) % 2: CHK.append(DUMMY)
    if flip % 2 != 0:
        # only even flips out ends up pointing to arg2 and needs to be
        # copied
        out[:] = arg1
    return out

Example 27

Project: pygame_sdl2
Source File: pixelcopy_test.py
View license
    def test_surface_to_array_2d(self):
        try:
            from numpy import empty, dtype
        except ImportError:
            return

        palette = self.test_palette
        alpha_color = (0, 0, 0, 128)

        dst_dims = self.surf_size
        destinations = [empty(dst_dims, t) for t in self.dst_types]
        if (pygame.get_sdl_byteorder() == pygame.LIL_ENDIAN):
            swapped_dst = empty(dst_dims, dtype('>u4'))
        else:
            swapped_dst = empty(dst_dims, dtype('<u4'))

        for surf in self.sources:
            src_bytesize = surf.get_bytesize()
            for dst in destinations:
                if dst.itemsize < src_bytesize:
                    self.assertRaises(ValueError, surface_to_array, dst, surf)
                    continue
                dst[...] = 0
                self.assertFalse(surf.get_locked())
                surface_to_array(dst, surf)
                self.assertFalse(surf.get_locked())
                for posn, i in self.test_points:
                    sp = unsigned32(surf.get_at_mapped(posn))
                    dp = dst[posn]
                    self.assertEqual(dp, sp,
                                     "%s != %s: flags: %i"
                                     ", bpp: %i, dtype: %s,  posn: %s" %
                                     (dp, sp,
                                      surf.get_flags(), surf.get_bitsize(),
                                      dst.dtype,
                                      posn))

                if surf.get_masks()[3]:
                    posn = (2, 1)
                    surf.set_at(posn, alpha_color)
                    surface_to_array(dst, surf)
                    sp = unsigned32(surf.get_at_mapped(posn))
                    dp = dst[posn]
                    self.assertEqual(dp, sp, "%s != %s: bpp: %i" %
                                     (dp, sp, surf.get_bitsize()))

            swapped_dst[...] = 0
            self.assertFalse(surf.get_locked())
            surface_to_array(swapped_dst, surf)
            self.assertFalse(surf.get_locked())
            for posn, i in self.test_points:
                sp = unsigned32(surf.get_at_mapped(posn))
                dp = swapped_dst[posn]
                self.assertEqual(dp, sp,
                                 "%s != %s: flags: %i"
                                 ", bpp: %i, dtype: %s,  posn: %s" %
                                 (dp, sp,
                                  surf.get_flags(), surf.get_bitsize(),
                                  dst.dtype,
                                  posn))

            if surf.get_masks()[3]:
                posn = (2, 1)
                surf.set_at(posn, alpha_color)
                self.assertFalse(surf.get_locked())
                surface_to_array(swapped_dst, surf)
                self.assertFalse(surf.get_locked())
                sp = unsigned32(surf.get_at_mapped(posn))
                dp = swapped_dst[posn]
                self.assertEqual(dp, sp, "%s != %s: bpp: %i" %
                                     (dp, sp, surf.get_bitsize()))

Example 28

View license
    def section_coordinates(self, coordinate_list):
        """
        :param coordinate_list: [[x, z], [x,z]]
        :return: create an array that represents the width over the cross section.
        """

        # Closes the cross section if necessary
        if coordinate_list[0] != coordinate_list[-1]:
            coordinate_list.append(coordinate_list[0])

        self.coordinate_list = []
        for i in range(len(coordinate_list)):
            self.coordinate_list.append(Point(coordinate_list[i][0], coordinate_list[i][1]))

        # Check the input
        self.valid_cross_section()

        # The empty array that will be filled with width values
        self.width_array = np.empty(self.nValues)
        count = 0

        # determine the max z-value in the polygon and set this to the height
        self.height = 0
        for i in coordinate_list:
            if i[1] > self.height:
                self.height = i[1]

        self.create_height_array()

        # interpolate between the given points.
        for section_height in self.height_array:
            for val in range(len(coordinate_list)):
                if coordinate_list[val][1] >= section_height:
                    left_higher_coordinate = coordinate_list[val]
                    left_lower_coordinate = coordinate_list[val - 1]
                    break

            for val in range(1, len(coordinate_list)):
                if coordinate_list[-val][1] >= section_height:
                    right_higher_coordinate = coordinate_list[-val]
                    right_lower_coordinate = coordinate_list[-val + 1]
                    break

            diff_height_left = left_higher_coordinate[1] - left_lower_coordinate[1]
            diff_width_left = left_higher_coordinate[0] - left_lower_coordinate[0]
            diff_height_right = right_higher_coordinate[1] - right_lower_coordinate[1]
            diff_width_right = right_higher_coordinate[0] - right_lower_coordinate[0]

            if diff_height_left and diff_height_right > 1e-6:

                fraction = (section_height - left_lower_coordinate[1]) / diff_height_left
                x_left = diff_width_left * fraction + left_lower_coordinate[0]

                fraction = (section_height - right_lower_coordinate[1]) / diff_height_right
                x_right = diff_width_right * fraction + right_lower_coordinate[0]
                self.width_array[count] = (x_left + x_right)
                self.xValuesPerSection[count] = np.array([x_left, x_right])
            count += 1

        # set the height of one section
        self.section_height = self.height / self.nValues

Example 29

Project: attention-lvcsr
Source File: neighbours.py
View license
    def perform(self, node, inp, out_):
        ten4, neib_shape, neib_step = inp
        z, = out_
        # GpuImages2Neibs should not run this perform in DebugMode
        if type(self) != Images2Neibs:
            raise theano.gof.utils.MethodNotDefined()

        def CEIL_INTDIV(a, b):
            if a % b:
                return (a // b) + 1
            else:
                return a // b

        grid_c = -1  # number of patch in height
        grid_d = -1  # number of patch in width
        assert ten4.ndim == 4
        assert neib_shape.ndim == 1
        assert neib_shape.shape[0] == 2
        assert neib_step.ndim == 1
        assert neib_step.shape[0] == 2
        c, d = neib_shape
        step_x, step_y = neib_step
        mode = self.mode

        if mode == "wrap_centered":
            if (c % 2 != 1) or (d % 2 != 1):
                raise TypeError(
                    "Images2Neibs:"
                    " in mode wrap_centered need patch with odd shapes")

            if (ten4.shape[2] < c) or (ten4.shape[3] < d):
                raise TypeError(
                    "Images2Neibs: in wrap_centered mode, don't support"
                    " image shapes smaller then the patch shapes:"
                    " neib_shape=(%d,%d), ten4[2:]=[%d,%d]" %
                    (c, d, ten4.shape[2], ten4.shape[3]))
            grid_c = CEIL_INTDIV(ten4.shape[2], step_x)
            grid_d = CEIL_INTDIV(ten4.shape[3], step_y)

        elif mode == "valid":
            if (ten4.shape[2] < c) or (((ten4.shape[2] - c) % step_x) != 0):
                raise TypeError(
                    "neib_shape[0]=%d, neib_step[0]=%d and"
                    " ten4.shape[2]=%d not consistent" %
                    (c, step_x, ten4.shape[2]))
            if (ten4.shape[3] < d) or (((ten4.shape[3] - d) % step_y) != 0):
                raise TypeError(
                    "neib_shape[1]=%d, neib_step[1]=%d and"
                    " ten4.shape[3]=%d not consistent" %
                    (d, step_y, ten4.shape[3]))
            # number of patch in height
            grid_c = 1 + ((ten4.shape[2] - c) // step_x)
            # number of patch in width
            grid_d = 1 + ((ten4.shape[3] - d) // step_y)
        elif mode == "ignore_borders":
            # number of patch in height
            grid_c = 1 + ((ten4.shape[2] - c) // step_x)
            # number of patch in width
            grid_d = 1 + ((ten4.shape[3] - d) // step_y)
        else:
            raise TypeError("Images2Neibs: unknow mode '%s'" % mode)

        z_dim0 = grid_c * grid_d * ten4.shape[1] * ten4.shape[0]
        z_dim1 = c * d
        z[0] = numpy.empty((z_dim0, z_dim1), dtype=node.outputs[0].dtype)

        nb_batch = ten4.shape[0]
        nb_stack = ten4.shape[1]
        height = ten4.shape[2]
        width = ten4.shape[3]

        wrap_centered_idx_shift_x = c // 2
        wrap_centered_idx_shift_y = d // 2
        for n in range(nb_batch):
            for s in range(nb_stack):
                # loop over the number of patch in height
                for a in range(grid_c):
                    # loop over the number of patch in width
                    for b in range(grid_d):
                        z_row = b + grid_d * (a + grid_c * (s + nb_stack * n))
                        for i in range(c):
                            ten4_2 = i + a * step_x
                            if mode == "wrap_centered":
                                ten4_2 -= wrap_centered_idx_shift_x
                                if ten4_2 < 0:
                                    ten4_2 += height
                                elif ten4_2 >= height:
                                    ten4_2 -= height
                            for j in range(d):
                                ten4_3 = j + b * step_y
                                if mode == "wrap_centered":
                                    ten4_3 -= wrap_centered_idx_shift_y
                                    if ten4_3 < 0:
                                        ten4_3 += width
                                    elif ten4_3 >= width:
                                        ten4_3 -= width
                                z_col = j + d * i

                                z[0][z_row, z_col] = ten4[n, s, ten4_2, ten4_3]

Example 30

Project: lfd
Source File: tps_registration_parallel.py
View license
def tps_rpm_bij_grid(x_knd, y_lmd, n_iter = 20, reg_init = .1, reg_final = .001, rad_init = .1, rad_final = .005, rot_reg = 1e-3, 
                     plotting = False, plot_cb = None, x_weights = None, y_weights = None, outlierprior = .1, outlierfrac = 2e-1, 
                     parallel = False, ppservers=(), partition_step=None):
    """
    If parallel=True the computation is split among the cores of the local machine and the cores of the cluster
    So, if parallel=True and ppservers=(), the computation is split only among the local cores
    If parallel=True and partition_step is especified, the computation is split in blocks of partition_step x partition_step
    """
    #TODO warn if default parameters are different?
    if len(x_knd.shape) == 2:
        x_knd.resize((1,)+x_knd.shape)
    if len(y_lmd.shape) == 2:
        y_lmd.resize((1,)+y_lmd.shape)
    k = x_knd.shape[0]
    l = y_lmd.shape[0]
    tps_tups = np.empty((k,l), dtype=object)
     
    if not parallel:
        for i in range(k):
            for j in range(l):
                tps_tups[i,j] = tps_rpm_bij(x_knd[i], y_lmd[j], n_iter=n_iter, reg_init=reg_init, reg_final=reg_final, 
                                            rad_init=rad_init, rad_final=rad_final, rot_reg=rot_reg, 
                                            plotting=plotting, plot_cb=plot_cb, x_weights=x_weights, y_weights=y_weights, 
                                            outlierprior=outlierprior, outlierfrac=outlierfrac)
    else:
        job_server = pp.Server(ppservers=ppservers)
        #TODO check if servers on nodes are running
        print "Starting pp with", job_server.get_ncpus(), "local workers"
 
        # make sure to change the order of optional arguments if the function's signature is changed
        # parallel=False
        opt_arg_vals = (n_iter, reg_init, reg_final, rad_init, rad_final, rot_reg, plotting, plot_cb, x_weights, y_weights, outlierprior, outlierfrac, False, ppservers)
        dep_funcs = ()
        modules = ("import numpy as np",
                   "import scipy.spatial.distance as ssd")
        myglobals = {'tps_rpm_bij':registration.tps_rpm_bij, 
                 'loglinspace':registration.loglinspace,
                 'Transformation':Transformation,
                 'ThinPlateSpline':ThinPlateSpline,
                 'tps_eval':tps.tps_eval,
                 'tps_kernel_matrix2':tps.tps_kernel_matrix2,
                 'tps_apply_kernel':tps.tps_apply_kernel,
                 'balance_matrix3':registration.balance_matrix3,
                 'fit_ThinPlateSpline':registration.fit_ThinPlateSpline,
                 'tps_fit3':tps.tps_fit3,
                 'tps_kernel_matrix':tps.tps_kernel_matrix,
                 'solve_eqp1':tps.solve_eqp1,
                 'tps_cost':tps.tps_cost
                 }
        if partition_step:
            partitions = [(i,min(i+partition_step,k),j,min(j+partition_step,l)) for i in range(0,k,partition_step) for j in range(0,l,partition_step)]
        else:
            partitions = [(i,i+1,0,l) for i in range(k)]
        jobs = [job_server.submit(tps_rpm_bij_grid, (x_knd[i_start:i_end], y_lmd[j_start:j_end])+opt_arg_vals, dep_funcs, modules, globals=myglobals) 
                for (i_start,i_end,j_start,j_end) in partitions]
        for ((i_start,i_end,j_start,j_end),job) in zip(partitions, jobs):
            tps_tups[i_start:i_end,j_start:j_end] = job()
        job_server.print_stats()
     
    return tps_tups

Example 31

Project: wikidump
Source File: segment.py
View license
def create_dataset( num_doc
                  , distribution
                  , num_seg 
                  , text_transformer
                  , exclude = []
                  , seed = None
                  ):
  """
  Create a multilingual dataset from wikipedia data based on segments
  of monolingual documents.

  @param num_doc: number of documents in intended dataset
  @param distribution: Mapping from class name to proportion
  @param segments: Sequence of segments and their relative sizes
  @param exclude: Docids to be excluded from the final dataset
  """
  logger = logging.getLogger('wikidump.segment.create_dataset')
  logger.info('Creating dataset')
  docmap = {}
  classmap = {}
  # Set up a parser for each class
  parser = {}
  classlabels = []
  class_probs = numpy.empty( len(distribution), dtype=float) 
  dump_indices = [] 
  denominator = 0 
  logger.info('Setting up parsers')
  for i,classlabel in enumerate(distribution):
    parser[classlabel] = page_parser(classlabel)
    classlabels.append(classlabel)
    denominator += distribution[classlabel] 
    class_probs[i] = denominator 
    dumpsize = dumpSize(classlabel)
    indices = range(dumpsize)
    random.shuffle(indices)
    dump_indices.append(indices)
    logger.debug('Done for %s, %d indices', classlabel, dumpsize)
  logger.info('Parsers ready')

  # Normalize to a CDF
  class_probs /= denominator

  random.seed(seed)

  logger.info('Obtaining documents')
  for doc_num in range(num_doc):
    doc = ""
    classes = set() 
    segids = []
    for seg_num in range(num_seg):
      seg_class_index = numpy.searchsorted(class_probs, random.random())
      seg_class = classlabels[seg_class_index] 
      classes.add(seg_class)
      doc_index = dump_indices[seg_class_index].pop(0)
      segid = "_".join((seg_class, str(doc_index)))
      while segid in exclude:
        logger.debug('Ignored %s: in exclude', segid)
        try:
          doc_index = dump_indices[seg_class_index].pop(0)
        except IndexError:
          raise ValueError, "No more documents in class %s available" % classlabel 
        segid = "_".join((seg_class, str(doc_index)))
      segids.append(segid)
      content = text_transformer(parser[seg_class].get_entry(doc_index))
      seg_size = len(content)/num_seg
      seg_start = seg_size * seg_num
      seg_end = seg_start + seg_size
      doc += content[seg_start:seg_end] 
    docid = '-'.join(segids) 
    docmap[docid] = doc
    classmap[docid] = list(classes)
    logger.debug('Index: %d ID: %s', doc_num, docid)
  return docmap, classmap

Example 32

Project: SALib
Source File: dgsm.py
View license
def analyze(problem, X, Y, num_resamples=1000,
            conf_level=0.95, print_to_console=False):
    """Calculates Derivative-based Global Sensitivity Measure on model outputs.
    
    Returns a dictionary with keys 'vi', 'vi_std', 'dgsm', and 'dgsm_conf',
    where each entry is a list of size D (the number of parameters) containing
    the indices in the same order as the parameter file.
          
    Parameters
    ----------
    problem : dict
        The problem definition
    X : numpy.matrix
        The NumPy matrix containing the model inputs
    Y : numpy.array
        The NumPy array containing the model outputs
    num_resamples : int
        The number of resamples used to compute the confidence
        intervals (default 1000)
    conf_level : float
        The confidence interval level (default 0.95)
    print_to_console : bool
        Print results directly to console (default False)
        
    References
    ----------
    .. [1] Sobol, I. M. and S. Kucherenko (2009). "Derivative based global
           sensitivity measures and their link with global sensitivity
           indices." Mathematics and Computers in Simulation, 79(10):3009-3017,
           doi:10.1016/j.matcom.2009.01.023.
    """

    D = problem['num_vars']

    if Y.size % (D + 1) == 0:
        N = int(Y.size / (D + 1))
    else:
        raise RuntimeError("Incorrect number of samples in model output file.")

    if not 0 < conf_level < 1:
        raise RuntimeError("Confidence level must be between 0-1.")

    base = np.empty(N)
    X_base = np.empty((N, D))
    perturbed = np.empty((N, D))
    X_perturbed = np.empty((N, D))
    step = D + 1

    base = Y[0:Y.size:step]
    X_base = X[0:Y.size:step, :]
    for j in range(D):
        perturbed[:, j] = Y[(j + 1):Y.size:step]
        X_perturbed[:, j] = X[(j + 1):Y.size:step, j]

    # First order (+conf.) and Total order (+conf.)
    keys = ('vi', 'vi_std', 'dgsm', 'dgsm_conf')
    S = dict((k, np.empty(D)) for k in keys)
    if print_to_console:
        print("Parameter %s %s %s %s" % keys)

    for j in range(D):
        S['vi'][j], S['vi_std'][j] = calc_vi(
            base, perturbed[:, j], X_perturbed[:, j] - X_base[:, j])
        S['dgsm'][j], S['dgsm_conf'][j] = calc_dgsm(base, perturbed[:, j], X_perturbed[
                                                    :, j] - X_base[:, j], problem['bounds'][j], num_resamples, conf_level)

        if print_to_console:
            print("%s %f %f %f %f" % (
                problem['names'][j], S['vi'][j], S['vi_std'][j], S['dgsm'][j], S['dgsm_conf'][j]))

    return S

Example 33

Project: SALib
Source File: fast.py
View license
def analyze(problem, Y, M=4, print_to_console=False):
    """Performs the Fourier Amplitude Sensitivity Test (FAST) on model outputs.
    
    Returns a dictionary with keys 'S1' and 'ST', where each entry is a list of
    size D (the number of parameters) containing the indices in the same order
    as the parameter file.
    
    Parameters
    ----------
    problem : dict
        The problem definition
    Y : numpy.array
        A NumPy array containing the model outputs
    M : int
        The interference parameter, i.e., the number of harmonics to sum in
        the Fourier series decomposition (default 4)
    print_to_console : bool
        Print results directly to console (default False)
        
    References
    ----------
    .. [1] Cukier, R. I., C. M. Fortuin, K. E. Shuler, A. G. Petschek, and J. H.
           Schaibly (1973).  "Study of the sensitivity of coupled reaction
           systems to uncertainties in rate coefficients."  J. Chem. Phys.,
           59(8):3873-3878, doi:10.1063/1.1680571.
           
    .. [2] Saltelli, A., S. Tarantola, and K. P.-S. Chan (1999).  "A
          Quantitative Model-Independent Method for Global Sensitivity
          Analysis of Model Output."  Technometrics, 41(1):39-56,
          doi:10.1080/00401706.1999.10485594.
          
    Examples
    --------
    >>> X = fast_sampler.sample(problem, 1000)
    >>> Y = Ishigami.evaluate(X)
    >>> Si = fast.analyze(problem, Y, print_to_console=False)
    """

    D = problem['num_vars']

    if Y.size % (D) == 0:
        N = int(Y.size / D)
    else:
        print("""
            Error: Number of samples in model output file must be a multiple of D, 
            where D is the number of parameters in your parameter file.
          """)
        exit()

    # Recreate the vector omega used in the sampling
    omega = np.empty([D])
    omega[0] = math.floor((N - 1) / (2 * M))
    m = math.floor(omega[0] / (2 * M))

    if m >= (D - 1):
        omega[1:] = np.floor(np.linspace(1, m, D - 1))
    else:
        omega[1:] = np.arange(D - 1) % m + 1

    # Calculate and Output the First and Total Order Values
    if print_to_console:
        print("Parameter First Total")
    Si = dict((k, [None] * D) for k in ['S1', 'ST'])
    for i in range(D):
        l = np.arange(i * N, (i + 1) * N)
        Si['S1'][i] = compute_first_order(Y[l], N, M, omega[0])
        Si['ST'][i] = compute_total_order(Y[l], N, omega[0])
        if print_to_console:
            print("%s %f %f" % 
                  (problem['names'][i], Si['S1'][i], Si['ST'][i]))
    return Si

Example 34

Project: SALib
Source File: saltelli.py
View license
def sample(problem, N, calc_second_order=True):
    """Generates model inputs using Saltelli's extension of the Sobol sequence.

    Returns a NumPy matrix containing the model inputs using Saltelli's sampling
    scheme.  Saltelli's scheme extends the Sobol sequence in a way to reduce
    the error rates in the resulting sensitivity index calculations.  If
    calc_second_order is False, the resulting matrix has N * (D + 2)
    rows, where D is the number of parameters.  If calc_second_order is True,
    the resulting matrix has N * (2D + 2) rows.  These model inputs are
    intended to be used with :func:`SALib.analyze.sobol.analyze`.

    Parameters
    ----------
    problem : dict
        The problem definition
    N : int
        The number of samples to generate
    calc_second_order : bool
        Calculate second-order sensitivities (default True)
    """
    D = problem['num_vars']
    groups = problem.get('groups')

    if not groups:
        Dg = problem['num_vars']
    else:
        Dg = len(set(groups))
        G, group_names = compute_groups_matrix(groups, D)

    # How many values of the Sobol sequence to skip
    skip_values = 1000

    # Create base sequence - could be any type of sampling
    base_sequence = sobol_sequence.sample(N + skip_values, 2 * D)

    if calc_second_order:
        saltelli_sequence = np.empty([(2 * Dg + 2) * N, D])
    else:
        saltelli_sequence = np.empty([(Dg + 2) * N, D])
    index = 0

    for i in range(skip_values, N + skip_values):

        # Copy matrix "A"
        for j in range(D):
            saltelli_sequence[index, j] = base_sequence[i, j]

        index += 1

        # Cross-sample elements of "B" into "A"
        for k in range(Dg):
            for j in range(D):
                if (not groups and j == k) or (groups and group_names[k] == groups[j]):
                    saltelli_sequence[index, j] = base_sequence[i, j + D]
                else:
                    saltelli_sequence[index, j] = base_sequence[i, j]

            index += 1

        # Cross-sample elements of "A" into "B"
        # Only needed if you're doing second-order indices (true by default)
        if calc_second_order:
            for k in range(Dg):
                for j in range(D):
                    if (not groups and j == k) or (groups and group_names[k] == groups[j]):
                        saltelli_sequence[index, j] = base_sequence[i, j]
                    else:
                        saltelli_sequence[index, j] = base_sequence[i, j + D]

                index += 1

        # Copy matrix "B"
        for j in range(D):
            saltelli_sequence[index, j] = base_sequence[i, j + D]

        index += 1
    if not problem.get('dists'):
        # scaling values out of 0-1 range with uniform distributions
        scale_samples(saltelli_sequence, problem['bounds'])
        return saltelli_sequence
    else:
        # scaling values to other distributions based on inverse CDFs
        scaled_saltelli = nonuniform_scale_samples(saltelli_sequence, problem['bounds'], problem['dists'])
        return scaled_saltelli

Example 35

Project: sci-wms
Source File: ugrid.py
View license
    def getmap(self, layer, request):
        time_index, time_value = self.nearest_time(layer, request.GET['time'])
        wgs84_bbox = request.GET['wgs84_bbox']

        with self.dataset() as nc:
            data_obj = nc.variables[layer.access_name]
            data_location = data_obj.location
            mesh_name = data_obj.mesh

            ug = UGrid.from_ncfile(self.topology_file, mesh_name=mesh_name)
            coords = np.empty(0)
            if data_location == 'node':
                coords = ug.nodes
            elif data_location == 'face':
                coords = ug.face_coordinates
            elif data_location == 'edge':
                coords = ug.edge_coordinates

            lon = coords[:, 0]
            lat = coords[:, 1]

            # Calculate any vector padding if we need to
            padding = None
            vector_step = request.GET['vectorstep']
            if request.GET['image_type'] == 'vectors':
                padding_factor = calc_safety_factor(request.GET['vectorscale'])
                padding = calc_lon_lat_padding(lon, lat, padding_factor) * vector_step

            # Calculate the boolean spatial mask to slice with
            bool_spatial_idx = data_handler.ugrid_lat_lon_subset_idx(lon, lat,
                                                                     bbox=wgs84_bbox.bbox,
                                                                     padding=padding)

            # Randomize vectors to subset if we need to
            if request.GET['image_type'] == 'vectors' and vector_step > 1:
                num_vec = int(bool_spatial_idx.size / vector_step)
                step = int(bool_spatial_idx.size / num_vec)
                bool_spatial_idx[np.where(bool_spatial_idx==True)][0::step] = False

            # If no triangles intersect the field of view, return a transparent tile
            if not np.any(bool_spatial_idx):
                logger.warning("No triangles in field of view, returning empty tile.")
                return self.empty_response(layer, request)

            if isinstance(layer, Layer):
                if (len(data_obj.shape) == 3):
                    z_index, z_value = self.nearest_z(layer, request.GET['elevation'])
                    data = data_obj[time_index, z_index, :]
                elif (len(data_obj.shape) == 2):
                    data = data_obj[time_index, :]
                elif len(data_obj.shape) == 1:
                    data = data_obj[:]
                else:
                    logger.debug("Dimension Mismatch: data_obj.shape == {0} and time = {1}".format(data_obj.shape, time_value))
                    return self.empty_response(layer, request)

                if request.GET['image_type'] in ['pcolor', 'contours', 'filledcontours']:
                    # Avoid triangles with nan values
                    bool_spatial_idx[np.isnan(data)] = False

                    # Get the faces to plot
                    faces = ug.faces[:]
                    face_idx = data_handler.face_idx_from_node_idx(faces, bool_spatial_idx)
                    faces_subset = faces[face_idx]
                    tri_subset = Tri.Triangulation(lon, lat, triangles=faces_subset)

                    if request.GET['image_type'] == 'pcolor':
                        return mpl_handler.tripcolor_response(tri_subset, data, request, data_location=data_location)
                    else:
                        return mpl_handler.tricontouring_response(tri_subset, data, request)
                elif request.GET['image_type'] in ['filledhatches', 'hatches']:
                    raise NotImplementedError('matplotlib does not support hatching on triangular grids... sorry!')
                else:
                    raise NotImplementedError('Image type "{}" is not supported.'.format(request.GET['image_type']))

            elif isinstance(layer, VirtualLayer):
                # Data needs to be [var1,var2] where var are 1D (nodes only, elevation and time already handled)
                data = []
                for l in layer.layers:
                    data_obj = nc.variables[l.var_name]
                    if (len(data_obj.shape) == 3):
                        z_index, z_value = self.nearest_z(layer, request.GET['elevation'])
                        data.append(data_obj[time_index, z_index, bool_spatial_idx])
                    elif (len(data_obj.shape) == 2):
                        data.append(data_obj[time_index, bool_spatial_idx])
                    elif len(data_obj.shape) == 1:
                        data.append(data_obj[bool_spatial_idx])
                    else:
                        logger.debug("Dimension Mismatch: data_obj.shape == {0} and time = {1}".format(data_obj.shape, time_value))
                        return self.empty_response(layer, request)

                if request.GET['image_type'] == 'vectors':
                    return mpl_handler.quiver_response(lon[bool_spatial_idx],
                                                       lat[bool_spatial_idx],
                                                       data[0],
                                                       data[1],
                                                       request)
                else:
                    raise NotImplementedError('Image type "{}" is not supported.'.format(request.GET['image_type']))

Example 36

Project: sci-wms
Source File: ugrid_tides.py
View license
    def get_tidal_vectors(self, layer, time, bbox, vector_scale=None, vector_step=None):

        vector_scale = vector_scale or 1
        vector_step = vector_step or 1

        with netCDF4.Dataset(self.topology_file) as nc:
            data_obj = nc.variables[layer.access_name]
            data_location = getattr(data_obj, 'location', 'node')
            mesh_name = data_obj.mesh

            ug = UGrid.from_nc_dataset(nc, mesh_name=mesh_name)
            coords = np.empty(0)
            if data_location == 'node':
                coords = ug.nodes
            elif data_location == 'face':
                coords = ug.face_coordinates
            elif data_location == 'edge':
                coords = ug.edge_coordinates

            lon = coords[:, 0]
            lat = coords[:, 1]

            padding_factor = calc_safety_factor(vector_scale)
            padding = calc_lon_lat_padding(lon, lat, padding_factor) * vector_step
            spatial_idx = data_handler.ugrid_lat_lon_subset_idx(lon, lat,
                                                                bbox=bbox,
                                                                padding=padding)

            tnames = nc.get_variables_by_attributes(standard_name='tide_constituent')[0]
            tfreqs = nc.get_variables_by_attributes(standard_name='tide_frequency')[0]

            from utide import _ut_constants_fname
            from utide.utilities import loadmatbunch
            con_info = loadmatbunch(_ut_constants_fname)['const']

            # Get names from the utide constant file
            utide_const_names = [ e.strip() for e in con_info['name'].tolist() ]

            # netCDF4-python is returning ugly arrays of bytes...
            names = []
            for n in tnames[:]:
                z = ''.join([ x.decode('utf-8') for x in n.tolist() if x ]).strip()
                names.append(z)

            if 'STEADY' in names:
                names[names.index('STEADY')] = 'Z0'
            extract_names = list(set(utide_const_names).intersection(set(names)))

            ntides = data_obj.shape[data_obj.dimensions.index('ntides')]
            extract_mask = np.zeros(shape=(ntides,), dtype=bool)
            for n in extract_names:
                extract_mask[names.index(n)] = True

            if not spatial_idx.any() or not extract_mask.any():
                e = np.ma.empty(0)
                return e, e, e, e

            ua = nc.variables['u'][extract_mask, spatial_idx]
            va = nc.variables['v'][extract_mask, spatial_idx]
            up = nc.variables['u_phase'][extract_mask, spatial_idx]
            vp = nc.variables['v_phase'][extract_mask, spatial_idx]
            freqs = tfreqs[extract_mask]

            omega = freqs * 3600  # Convert from radians/s to radians/hour.

            from utide.harmonics import FUV
            from matplotlib.dates import date2num
            v, u, f = FUV(t=np.array([date2num(time) + 366.1667]),
                          tref=np.array([0]),
                          lind=np.array([ utide_const_names.index(x) for x in extract_names ]),
                          lat=55,  # Reference latitude for 3rd order satellites (degrees) (55 is fine always)
                          ngflgs=[0, 0, 0, 0])  # [NodsatLint NodsatNone GwchLint GwchNone]

            s = calendar.timegm(time.timetuple()) / 60 / 60.
            v, u, f = map(np.squeeze, (v, u, f))
            v = v * 2 * np.pi  # Convert phase in radians.
            u = u * 2 * np.pi  # Convert phase in radians.

            U = (f * ua.T * np.cos(v + s * omega + u - up.T * np.pi / 180)).sum(axis=1)
            V = (f * va.T * np.cos(v + s * omega + u - vp.T * np.pi / 180)).sum(axis=1)

            return U, V, lon[spatial_idx], lat[spatial_idx]

Example 37

Project: scikit-beam
Source File: binned_statistic.py
View license
    def __call__(self, values):
        """
        Parameters
        ----------
        values : array_like
            The values on which the statistic will be computed.  This must be
            the same shape as `sample` in the constructor.

        Returns
        -------
        statistic : array
            The values of the selected statistic in each bin.
        """

        self.result = np.empty(self.nbin.prod(), float)
        if self.statistic == 'mean':
            self.result.fill(np.nan)
            flatsum = np.bincount(self.xy, values)
            a = self.flatcount.nonzero()
            self.result[a] = flatsum[a] / self.flatcount[a]
        elif self.statistic == 'std':
            self.result.fill(0)
            flatsum = np.bincount(self.xy, values)
            flatsum2 = np.bincount(self.xy, values ** 2)
            a = self.flatcount.nonzero()
            self.result[a] = np.sqrt(flatsum2[a] / self.flatcount[a] -
                                     (flatsum[a] / self.flatcount[a]) ** 2)
        elif self.statistic == 'count':
            self.result.fill(0)
            a = np.arange(len(self.flatcount))
            self.result[a] = self.flatcount
        elif self.statistic == 'sum':
            self.result.fill(0)
            flatsum = np.bincount(self.xy, values)
            a = np.arange(len(flatsum))
            self.result[a] = flatsum
        elif self.statistic == 'median':
            self.result.fill(np.nan)
            for i in np.unique(self.xy):
                self.result[i] = np.median(values[self.xy == i])
        elif callable(self.statistic):
            with warnings.catch_warnings():
                # Numpy generates a warnings for mean/std/... with empty list
                warnings.filterwarnings('ignore', category=RuntimeWarning)
                old = np.seterr(invalid='ignore')
                try:
                    null = self.statistic([])
                except:
                    null = np.nan
                np.seterr(**old)
            self.result.fill(null)
            for i in np.unique(self.xy):
                self.result[i] = self.statistic(values[self.xy == i])

        # Shape into a proper matrix
        self.result = self.result.reshape(np.sort(self.nbin))
        ni = np.copy(self.ni)
        for i in np.arange(self.nbin.size):
            j = ni.argsort()[i]
            self.result = self.result.swapaxes(i, j)
            ni[i], ni[j] = ni[j], ni[i]

        # Remove outliers (indices 0 and -1 for each dimension).
        core = self.D * [slice(1, -1)]
        self.result = self.result[core]

        if (self.result.shape != self.nbin - 2).any():
            raise RuntimeError('Internal Shape Error')

        return self.result

Example 38

View license
def compare_ssim(X, Y, win_size=None, gradient=False,
                 dynamic_range=None, multichannel=False,
                 gaussian_weights=False, full=False, **kwargs):
    """Compute the mean structural similarity index between two images.

    Parameters
    ----------
    X, Y : ndarray
        Image.  Any dimensionality.
    win_size : int or None
        The side-length of the sliding window used in comparison.  Must be an
        odd value.  If `gaussian_weights` is True, this is ignored and the
        window size will depend on `sigma`.
    gradient : bool, optional
        If True, also return the gradient.
    dynamic_range : int, optional
        The dynamic range of the input image (distance between minimum and
        maximum possible values).  By default, this is estimated from the image
        data-type.
    multichannel : bool, optional
        If True, treat the last dimension of the array as channels. Similarity
        calculations are done independently for each channel then averaged.
    gaussian_weights : bool, optional
        If True, each patch has its mean and variance spatially weighted by a
        normalized Gaussian kernel of width sigma=1.5.
    full : bool, optional
        If True, return the full structural similarity image instead of the
        mean value.

    Other Parameters
    ----------------
    use_sample_covariance : bool
        if True, normalize covariances by N-1 rather than, N where N is the
        number of pixels within the sliding window.
    K1 : float
        algorithm parameter, K1 (small constant, see [1]_)
    K2 : float
        algorithm parameter, K2 (small constant, see [1]_)
    sigma : float
        sigma for the Gaussian when `gaussian_weights` is True.

    Returns
    -------
    mssim : float
        The mean structural similarity over the image.
    grad : ndarray
        The gradient of the structural similarity index between X and Y [2]_.
        This is only returned if `gradient` is set to True.
    S : ndarray
        The full SSIM image.  This is only returned if `full` is set to True.

    Notes
    -----
    To match the implementation of Wang et. al. [1]_, set `gaussian_weights`
    to True, `sigma` to 1.5, and `use_sample_covariance` to False.

    References
    ----------
    .. [1] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P.
       (2004). Image quality assessment: From error visibility to
       structural similarity. IEEE Transactions on Image Processing,
       13, 600-612.
       https://ece.uwaterloo.ca/~z70wang/publications/ssim.pdf,
       DOI:10.1.1.11.2477

    .. [2] Avanaki, A. N. (2009). Exact global histogram specification
       optimized for structural similarity. Optical Review, 16, 613-621.
       http://arxiv.org/abs/0901.0065,
       DOI:10.1007/s10043-009-0119-z

    """
    if not X.dtype == Y.dtype:
        raise ValueError('Input images must have the same dtype.')

    if not X.shape == Y.shape:
        raise ValueError('Input images must have the same dimensions.')

    if multichannel:
        # loop over channels
        args = dict(win_size=win_size,
                    gradient=gradient,
                    dynamic_range=dynamic_range,
                    multichannel=False,
                    gaussian_weights=gaussian_weights,
                    full=full)
        args.update(kwargs)
        nch = X.shape[-1]
        mssim = np.empty(nch)
        if gradient:
            G = np.empty(X.shape)
        if full:
            S = np.empty(X.shape)
        for ch in range(nch):
            ch_result = compare_ssim(X[..., ch], Y[..., ch], **args)
            if gradient and full:
                mssim[..., ch], G[..., ch], S[..., ch] = ch_result
            elif gradient:
                mssim[..., ch], G[..., ch] = ch_result
            elif full:
                mssim[..., ch], S[..., ch] = ch_result
            else:
                mssim[..., ch] = ch_result
        mssim = mssim.mean()
        if gradient and full:
            return mssim, G, S
        elif gradient:
            return mssim, G
        elif full:
            return mssim, S
        else:
            return mssim

    K1 = kwargs.pop('K1', 0.01)
    K2 = kwargs.pop('K2', 0.03)
    sigma = kwargs.pop('sigma', 1.5)
    if K1 < 0:
        raise ValueError("K1 must be positive")
    if K2 < 0:
        raise ValueError("K2 must be positive")
    if sigma < 0:
        raise ValueError("sigma must be positive")
    use_sample_covariance = kwargs.pop('use_sample_covariance', True)

    if win_size is None:
        if gaussian_weights:
            win_size = 11  # 11 to match Wang et. al. 2004
        else:
            win_size = 7   # backwards compatibility

    if np.any((np.asarray(X.shape) - win_size) < 0):
        raise ValueError(
            "win_size exceeds image extent.  If the input is a multichannel "
            "(color) image, set multichannel=True.")

    if not (win_size % 2 == 1):
        raise ValueError('Window size must be odd.')

    if dynamic_range is None:
        dmin, dmax = dtype_range[X.dtype.type]
        dynamic_range = dmax - dmin

    ndim = X.ndim

    if gaussian_weights:
        # sigma = 1.5 to approximately match filter in Wang et. al. 2004
        # this ends up giving a 13-tap rather than 11-tap Gaussian
        filter_func = gaussian_filter
        filter_args = {'sigma': sigma}

    else:
        filter_func = uniform_filter
        filter_args = {'size': win_size}

    # ndimage filters need floating point data
    X = X.astype(np.float64)
    Y = Y.astype(np.float64)

    NP = win_size ** ndim

    # filter has already normalized by NP
    if use_sample_covariance:
        cov_norm = NP / (NP - 1)  # sample covariance
    else:
        cov_norm = 1.0  # population covariance to match Wang et. al. 2004

    # compute (weighted) means
    ux = filter_func(X, **filter_args)
    uy = filter_func(Y, **filter_args)

    # compute (weighted) variances and covariances
    uxx = filter_func(X * X, **filter_args)
    uyy = filter_func(Y * Y, **filter_args)
    uxy = filter_func(X * Y, **filter_args)
    vx = cov_norm * (uxx - ux * ux)
    vy = cov_norm * (uyy - uy * uy)
    vxy = cov_norm * (uxy - ux * uy)

    R = dynamic_range
    C1 = (K1 * R) ** 2
    C2 = (K2 * R) ** 2

    A1, A2, B1, B2 = ((2 * ux * uy + C1,
                       2 * vxy + C2,
                       ux ** 2 + uy ** 2 + C1,
                       vx + vy + C2))
    D = B1 * B2
    S = (A1 * A2) / D

    # to avoid edge effects will ignore filter radius strip around edges
    pad = (win_size - 1) // 2

    # compute (weighted) mean of ssim
    mssim = crop(S, pad).mean()

    if gradient:
        # The following is Eqs. 7-8 of Avanaki 2009.
        grad = filter_func(A1 / D, **filter_args) * X
        grad += filter_func(-S / B2, **filter_args) * Y
        grad += filter_func((ux * (A2 - A1) - uy * (B2 - B1) * S) / D,
                            **filter_args)
        grad *= (2 / X.size)

        if full:
            return mssim, grad, S
        else:
            return mssim, grad
    else:
        if full:
            return mssim, S
        else:
            return mssim

Example 39

Project: scikit-image
Source File: _geometric.py
View license
def warp_coords(coord_map, shape, dtype=np.float64):
    """Build the source coordinates for the output of a 2-D image warp.

    Parameters
    ----------
    coord_map : callable like GeometricTransform.inverse
        Return input coordinates for given output coordinates.
        Coordinates are in the shape (P, 2), where P is the number
        of coordinates and each element is a ``(row, col)`` pair.
    shape : tuple
        Shape of output image ``(rows, cols[, bands])``.
    dtype : np.dtype or string
        dtype for return value (sane choices: float32 or float64).

    Returns
    -------
    coords : (ndim, rows, cols[, bands]) array of dtype `dtype`
            Coordinates for `scipy.ndimage.map_coordinates`, that will yield
            an image of shape (orows, ocols, bands) by drawing from source
            points according to the `coord_transform_fn`.

    Notes
    -----

    This is a lower-level routine that produces the source coordinates for 2-D
    images used by `warp()`.

    It is provided separately from `warp` to give additional flexibility to
    users who would like, for example, to re-use a particular coordinate
    mapping, to use specific dtypes at various points along the the
    image-warping process, or to implement different post-processing logic
    than `warp` performs after the call to `ndi.map_coordinates`.


    Examples
    --------
    Produce a coordinate map that shifts an image up and to the right:

    >>> from skimage import data
    >>> from scipy.ndimage import map_coordinates
    >>>
    >>> def shift_up10_left20(xy):
    ...     return xy - np.array([-20, 10])[None, :]
    >>>
    >>> image = data.astronaut().astype(np.float32)
    >>> coords = warp_coords(shift_up10_left20, image.shape)
    >>> warped_image = map_coordinates(image, coords)

    """
    shape = safe_as_int(shape)
    rows, cols = shape[0], shape[1]
    coords_shape = [len(shape), rows, cols]
    if len(shape) == 3:
        coords_shape.append(shape[2])
    coords = np.empty(coords_shape, dtype=dtype)

    # Reshape grid coordinates into a (P, 2) array of (row, col) pairs
    tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T

    # Map each (row, col) pair to the source image according to
    # the user-provided mapping
    tf_coords = coord_map(tf_coords)

    # Reshape back to a (2, M, N) coordinate grid
    tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)

    # Place the y-coordinate mapping
    _stackcopy(coords[1, ...], tf_coords[0, ...])

    # Place the x-coordinate mapping
    _stackcopy(coords[0, ...], tf_coords[1, ...])

    if len(shape) == 3:
        coords[2, ...] = range(shape[2])

    return coords

Example 40

Project: scikit-learn
Source File: bench_plot_omp_lars.py
View license
def compute_bench(samples_range, features_range):

    it = 0

    results = dict()
    lars = np.empty((len(features_range), len(samples_range)))
    lars_gram = lars.copy()
    omp = lars.copy()
    omp_gram = lars.copy()

    max_it = len(samples_range) * len(features_range)
    for i_s, n_samples in enumerate(samples_range):
        for i_f, n_features in enumerate(features_range):
            it += 1
            n_informative = n_features / 10
            print('====================')
            print('Iteration %03d of %03d' % (it, max_it))
            print('====================')
            # dataset_kwargs = {
            #     'n_train_samples': n_samples,
            #     'n_test_samples': 2,
            #     'n_features': n_features,
            #     'n_informative': n_informative,
            #     'effective_rank': min(n_samples, n_features) / 10,
            #     #'effective_rank': None,
            #     'bias': 0.0,
            # }
            dataset_kwargs = {
                'n_samples': 1,
                'n_components': n_features,
                'n_features': n_samples,
                'n_nonzero_coefs': n_informative,
                'random_state': 0
            }
            print("n_samples: %d" % n_samples)
            print("n_features: %d" % n_features)
            y, X, _ = make_sparse_coded_signal(**dataset_kwargs)
            X = np.asfortranarray(X)

            gc.collect()
            print("benchmarking lars_path (with Gram):", end='')
            sys.stdout.flush()
            tstart = time()
            G = np.dot(X.T, X)  # precomputed Gram matrix
            Xy = np.dot(X.T, y)
            lars_path(X, y, Xy=Xy, Gram=G, max_iter=n_informative)
            delta = time() - tstart
            print("%0.3fs" % delta)
            lars_gram[i_f, i_s] = delta

            gc.collect()
            print("benchmarking lars_path (without Gram):", end='')
            sys.stdout.flush()
            tstart = time()
            lars_path(X, y, Gram=None, max_iter=n_informative)
            delta = time() - tstart
            print("%0.3fs" % delta)
            lars[i_f, i_s] = delta

            gc.collect()
            print("benchmarking orthogonal_mp (with Gram):", end='')
            sys.stdout.flush()
            tstart = time()
            orthogonal_mp(X, y, precompute=True,
                          n_nonzero_coefs=n_informative)
            delta = time() - tstart
            print("%0.3fs" % delta)
            omp_gram[i_f, i_s] = delta

            gc.collect()
            print("benchmarking orthogonal_mp (without Gram):", end='')
            sys.stdout.flush()
            tstart = time()
            orthogonal_mp(X, y, precompute=False,
                          n_nonzero_coefs=n_informative)
            delta = time() - tstart
            print("%0.3fs" % delta)
            omp[i_f, i_s] = delta

    results['time(LARS) / time(OMP)\n (w/ Gram)'] = (lars_gram / omp_gram)
    results['time(LARS) / time(OMP)\n (w/o Gram)'] = (lars / omp)
    return results

Example 41

Project: scikit-learn
Source File: mutual_info_.py
View license
def _estimate_mi(X, y, discrete_features='auto', discrete_target=False,
                 n_neighbors=3, copy=True, random_state=None):
    """Estimate mutual information between the features and the target.

    Parameters
    ----------
    X : array_like or sparse matrix, shape (n_samples, n_features)
        Feature matrix.

    y : array_like, shape (n_samples,)
        Target vector.

    discrete_features : {'auto', bool, array_like}, default 'auto'
        If bool, then determines whether to consider all features discrete
        or continuous. If array, then it should be either a boolean mask
        with shape (n_features,) or array with indices of discrete features.
        If 'auto', it is assigned to False for dense `X` and to True for
        sparse `X`.

    discrete_target : bool, default False
        Whether to consider `y` as a discrete variable.

    n_neighbors : int, default 3
        Number of neighbors to use for MI estimation for continuous variables,
        see [1]_ and [2]_. Higher values reduce variance of the estimation, but
        could introduce a bias.

    copy : bool, default True
        Whether to make a copy of the given data. If set to False, the initial
        data will be overwritten.

    random_state : int seed, RandomState instance or None, default None
        The seed of the pseudo random number generator for adding small noise
        to continuous variables in order to remove repeated values.

    Returns
    -------
    mi : ndarray, shape (n_features,)
        Estimated mutual information between each feature and the target.
        A negative value will be replaced by 0.

    References
    ----------
    .. [1] A. Kraskov, H. Stogbauer and P. Grassberger, "Estimating mutual
           information". Phys. Rev. E 69, 2004.
    .. [2] B. C. Ross "Mutual Information between Discrete and Continuous
           Data Sets". PLoS ONE 9(2), 2014.
    """
    X, y = check_X_y(X, y, accept_sparse='csc', y_numeric=not discrete_target)
    n_samples, n_features = X.shape

    if discrete_features == 'auto':
        discrete_features = issparse(X)

    if isinstance(discrete_features, bool):
        discrete_mask = np.empty(n_features, dtype=bool)
        discrete_mask.fill(discrete_features)
    else:
        discrete_features = np.asarray(discrete_features)
        if discrete_features.dtype != 'bool':
            discrete_mask = np.zeros(n_features, dtype=bool)
            discrete_mask[discrete_features] = True
        else:
            discrete_mask = discrete_features

    continuous_mask = ~discrete_mask
    if np.any(continuous_mask) and issparse(X):
        raise ValueError("Sparse matrix `X` can't have continuous features.")

    rng = check_random_state(random_state)
    if np.any(continuous_mask):
        if copy:
            X = X.copy()

        if not discrete_target:
            X[:, continuous_mask] = scale(X[:, continuous_mask],
                                          with_mean=False, copy=False)

        # Add small noise to continuous features as advised in Kraskov et. al.
        X = X.astype(float)
        means = np.maximum(1, np.mean(np.abs(X[:, continuous_mask]), axis=0))
        X[:, continuous_mask] += 1e-10 * means * rng.randn(
                n_samples, np.sum(continuous_mask))

    if not discrete_target:
        y = scale(y, with_mean=False)
        y += 1e-10 * np.maximum(1, np.mean(np.abs(y))) * rng.randn(n_samples)

    mi = [_compute_mi(x, y, discrete_feature, discrete_target) for
          x, discrete_feature in moves.zip(_iterate_columns(X), discrete_mask)]

    return np.array(mi)

Example 42

Project: scipy
Source File: _fitpack_impl.py
View license
def splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None,
           full_output=0, per=0, quiet=1):
    """
    Find the B-spline representation of 1-D curve.

    Given the set of data points ``(x[i], y[i])`` determine a smooth spline
    approximation of degree k on the interval ``xb <= x <= xe``.

    Parameters
    ----------
    x, y : array_like
        The data points defining a curve y = f(x).
    w : array_like, optional
        Strictly positive rank-1 array of weights the same length as x and y.
        The weights are used in computing the weighted least-squares spline
        fit. If the errors in the y values have standard-deviation given by the
        vector d, then w should be 1/d. Default is ones(len(x)).
    xb, xe : float, optional
        The interval to fit.  If None, these default to x[0] and x[-1]
        respectively.
    k : int, optional
        The order of the spline fit. It is recommended to use cubic splines.
        Even order splines should be avoided especially with small s values.
        1 <= k <= 5
    task : {1, 0, -1}, optional
        If task==0 find t and c for a given smoothing factor, s.

        If task==1 find t and c for another value of the smoothing factor, s.
        There must have been a previous call with task=0 or task=1 for the same
        set of data (t will be stored an used internally)

        If task=-1 find the weighted least square spline for a given set of
        knots, t. These should be interior knots as knots on the ends will be
        added automatically.
    s : float, optional
        A smoothing condition. The amount of smoothness is determined by
        satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where g(x)
        is the smoothed interpolation of (x,y). The user can use s to control
        the tradeoff between closeness and smoothness of fit. Larger s means
        more smoothing while smaller values of s indicate less smoothing.
        Recommended values of s depend on the weights, w. If the weights
        represent the inverse of the standard-deviation of y, then a good s
        value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is
        the number of datapoints in x, y, and w. default : s=m-sqrt(2*m) if
        weights are supplied. s = 0.0 (interpolating) if no weights are
        supplied.
    t : array_like, optional
        The knots needed for task=-1. If given then task is automatically set
        to -1.
    full_output : bool, optional
        If non-zero, then return optional outputs.
    per : bool, optional
        If non-zero, data points are considered periodic with period x[m-1] -
        x[0] and a smooth periodic spline approximation is returned. Values of
        y[m-1] and w[m-1] are not used.
    quiet : bool, optional
        Non-zero to suppress messages.
        This parameter is deprecated; use standard Python warning filters
        instead.

    Returns
    -------
    tck : tuple
        (t,c,k) a tuple containing the vector of knots, the B-spline
        coefficients, and the degree of the spline.
    fp : array, optional
        The weighted sum of squared residuals of the spline approximation.
    ier : int, optional
        An integer flag about splrep success. Success is indicated if ier<=0.
        If ier in [1,2,3] an error occurred but was not raised. Otherwise an
        error is raised.
    msg : str, optional
        A message corresponding to the integer flag, ier.

    Notes
    -----
    See splev for evaluation of the spline and its derivatives.

    The user is responsible for assuring that the values of *x* are unique.
    Otherwise, *splrep* will not return sensible results.

    See Also
    --------
    UnivariateSpline, BivariateSpline
    splprep, splev, sproot, spalde, splint
    bisplrep, bisplev

    Notes
    -----
    See splev for evaluation of the spline and its derivatives. Uses the
    FORTRAN routine curfit from FITPACK.

    If provided, knots `t` must satisfy the Schoenberg-Whitney conditions,
    i.e., there must be a subset of data points ``x[j]`` such that
    ``t[j] < x[j] < t[j+k+1]``, for ``j=0, 1,...,n-k-2``.

    References
    ----------
    Based on algorithms described in [1]_, [2]_, [3]_, and [4]_:

    .. [1] P. Dierckx, "An algorithm for smoothing, differentiation and
       integration of experimental data using spline functions",
       J.Comp.Appl.Maths 1 (1975) 165-184.
    .. [2] P. Dierckx, "A fast algorithm for smoothing data on a rectangular
       grid while using spline functions", SIAM J.Numer.Anal. 19 (1982)
       1286-1304.
    .. [3] P. Dierckx, "An improved algorithm for curve fitting with spline
       functions", report tw54, Dept. Computer Science,K.U. Leuven, 1981.
    .. [4] P. Dierckx, "Curve and surface fitting with splines", Monographs on
       Numerical Analysis, Oxford University Press, 1993.

    Examples
    --------

    >>> import matplotlib.pyplot as plt
    >>> from scipy.interpolate import splev, splrep
    >>> x = np.linspace(0, 10, 10)
    >>> y = np.sin(x)
    >>> tck = splrep(x, y)
    >>> x2 = np.linspace(0, 10, 200)
    >>> y2 = splev(x2, tck)
    >>> plt.plot(x, y, 'o', x2, y2)
    >>> plt.show()

    """
    if task <= 0:
        _curfit_cache = {}
    x, y = map(atleast_1d, [x, y])
    m = len(x)
    if w is None:
        w = ones(m, float)
        if s is None:
            s = 0.0
    else:
        w = atleast_1d(w)
        if s is None:
            s = m - sqrt(2*m)
    if not len(w) == m:
        raise TypeError('len(w)=%d is not equal to m=%d' % (len(w), m))
    if (m != len(y)) or (m != len(w)):
        raise TypeError('Lengths of the first three arguments (x,y,w) must '
                        'be equal')
    if not (1 <= k <= 5):
        raise TypeError('Given degree of the spline (k=%d) is not supported. '
                        '(1<=k<=5)' % k)
    if m <= k:
        raise TypeError('m > k must hold')
    if xb is None:
        xb = x[0]
    if xe is None:
        xe = x[-1]
    if not (-1 <= task <= 1):
        raise TypeError('task must be -1, 0 or 1')
    if t is not None:
        task = -1
    if task == -1:
        if t is None:
            raise TypeError('Knots must be given for task=-1')
        numknots = len(t)
        _curfit_cache['t'] = empty((numknots + 2*k + 2,), float)
        _curfit_cache['t'][k+1:-k-1] = t
        nest = len(_curfit_cache['t'])
    elif task == 0:
        if per:
            nest = max(m + 2*k, 2*k + 3)
        else:
            nest = max(m + k + 1, 2*k + 3)
        t = empty((nest,), float)
        _curfit_cache['t'] = t
    if task <= 0:
        if per:
            _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(8 + 5*k),), float)
        else:
            _curfit_cache['wrk'] = empty((m*(k + 1) + nest*(7 + 3*k),), float)
        _curfit_cache['iwrk'] = empty((nest,), intc)
    try:
        t = _curfit_cache['t']
        wrk = _curfit_cache['wrk']
        iwrk = _curfit_cache['iwrk']
    except KeyError:
        raise TypeError("must call with task=1 only after"
                        " call with task=0,-1")
    if not per:
        n, c, fp, ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk,
                                        xb, xe, k, s)
    else:
        n, c, fp, ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s)
    tck = (t[:n], c[:n], k)
    if ier <= 0 and not quiet:
        _mess = (_iermess[ier][0] + "\tk=%d n=%d m=%d fp=%f s=%f" %
                 (k, len(t), m, fp, s))
        warnings.warn(RuntimeWarning(_mess))
    if ier > 0 and not full_output:
        if ier in [1, 2, 3]:
            warnings.warn(RuntimeWarning(_iermess[ier][0]))
        else:
            try:
                raise _iermess[ier][1](_iermess[ier][0])
            except KeyError:
                raise _iermess['unknown'][1](_iermess['unknown'][0])
    if full_output:
        try:
            return tck, fp, ier, _iermess[ier][0]
        except KeyError:
            return tck, fp, ier, _iermess['unknown'][0]
    else:
        return tck

Example 43

Project: scipy
Source File: special_matrices.py
View license
def pascal(n, kind='symmetric', exact=True):
    """
    Returns the n x n Pascal matrix.

    The Pascal matrix is a matrix containing the binomial coefficients as
    its elements.

    Parameters
    ----------
    n : int
        The size of the matrix to create; that is, the result is an n x n
        matrix.
    kind : str, optional
        Must be one of 'symmetric', 'lower', or 'upper'.
        Default is 'symmetric'.
    exact : bool, optional
        If `exact` is True, the result is either an array of type
        numpy.uint64 (if n < 35) or an object array of Python long integers.
        If `exact` is False, the coefficients in the matrix are computed using
        `scipy.special.comb` with `exact=False`.  The result will be a floating
        point array, and the values in the array will not be the exact
        coefficients, but this version is much faster than `exact=True`.

    Returns
    -------
    p : (n, n) ndarray
        The Pascal matrix.

    See Also
    --------
    invpascal

    Notes
    -----
    See http://en.wikipedia.org/wiki/Pascal_matrix for more information
    about Pascal matrices.

    .. versionadded:: 0.11.0

    Examples
    --------
    >>> from scipy.linalg import pascal
    >>> pascal(4)
    array([[ 1,  1,  1,  1],
           [ 1,  2,  3,  4],
           [ 1,  3,  6, 10],
           [ 1,  4, 10, 20]], dtype=uint64)
    >>> pascal(4, kind='lower')
    array([[1, 0, 0, 0],
           [1, 1, 0, 0],
           [1, 2, 1, 0],
           [1, 3, 3, 1]], dtype=uint64)
    >>> pascal(50)[-1, -1]
    25477612258980856902730428600L
    >>> from scipy.special import comb
    >>> comb(98, 49, exact=True)
    25477612258980856902730428600L

    """

    from scipy.special import comb
    if kind not in ['symmetric', 'lower', 'upper']:
        raise ValueError("kind must be 'symmetric', 'lower', or 'upper'")

    if exact:
        if n >= 35:
            L_n = np.empty((n, n), dtype=object)
            L_n.fill(0)
        else:
            L_n = np.zeros((n, n), dtype=np.uint64)
        for i in range(n):
            for j in range(i + 1):
                L_n[i, j] = comb(i, j, exact=True)
    else:
        L_n = comb(*np.ogrid[:n, :n])

    if kind is 'lower':
        p = L_n
    elif kind is 'upper':
        p = L_n.T
    else:
        p = np.dot(L_n, L_n.T)

    return p

Example 44

Project: scipy
Source File: _solvers.py
View license
def solve_continuous_are(a, b, q, r, e=None, s=None, balanced=True):
    """
    Solves the continuous-time algebraic Riccati equation (CARE).

    The CARE is defined as

    .. math::

          X A + A^H X - X B R^{-1} B^H X + Q = 0

    The limitations for a solution to exist are :

        * All eigenvalues of :math:`A` on the right half plane, should be
          controllable.

        * The associated hamiltonian pencil (See Notes), should have
          eigenvalues sufficiently away from the imaginary axis.

    Moreover, if ``e`` or ``s`` is not precisely ``None``, then the
    generalized version of CARE

    .. math::

          E^HXA + A^HXE - (E^HXB + S) R^{-1} (B^HXE + S^H) + Q = 0

    is solved. When omitted, ``e`` is assumed to be the identity and ``s``
    is assumed to be the zero matrix with sizes compatible with ``a`` and
    ``b`` respectively.

    Parameters
    ----------
    a : (M, M) array_like
        Square matrix
    b : (M, N) array_like
        Input
    q : (M, M) array_like
        Input
    r : (N, N) array_like
        Nonsingular square matrix
    e : (M, M) array_like, optional
        Nonsingular square matrix
    s : (M, N) array_like, optional
        Input
    balanced : bool, optional
        The boolean that indicates whether a balancing step is performed
        on the data. The default is set to True.

    Returns
    -------
    x : (M, M) ndarray
        Solution to the continuous-time algebraic Riccati equation.

    Raises
    ------
    LinAlgError
        For cases where the stable subspace of the pencil could not be
        isolated. See Notes section and the references for details.

    See Also
    --------
    solve_discrete_are : Solves the discrete-time algebraic Riccati equation

    Notes
    -----
    The equation is solved by forming the extended hamiltonian matrix pencil,
    as described in [1]_, :math:`H - \lambda J` given by the block matrices ::

        [ A    0    B ]             [ E   0    0 ]
        [-Q  -A^H  -S ] - \lambda * [ 0  E^H   0 ]
        [ S^H B^H   R ]             [ 0   0    0 ]

    and using a QZ decomposition method.

    In this algorithm, the fail conditions are linked to the symmetry
    of the product :math:`U_2 U_1^{-1}` and condition number of
    :math:`U_1`. Here, :math:`U` is the 2m-by-m matrix that holds the
    eigenvectors spanning the stable subspace with 2m rows and partitioned
    into two m-row matrices. See [1]_ and [2]_ for more details.

    In order to improve the QZ decomposition accuracy, the pencil goes
    through a balancing step where the sum of absolute values of
    :math:`H` and :math:`J` entries (after removing the diagonal entries of
    the sum) is balanced following the recipe given in [3]_.

    .. versionadded:: 0.11.0

    References
    ----------
    .. [1]  P. van Dooren , "A Generalized Eigenvalue Approach For Solving
       Riccati Equations.", SIAM Journal on Scientific and Statistical
       Computing, Vol.2(2), DOI: 10.1137/0902010

    .. [2] A.J. Laub, "A Schur Method for Solving Algebraic Riccati
       Equations.", Massachusetts Institute of Technology. Laboratory for
       Information and Decision Systems. LIDS-R ; 859. Available online :
       http://hdl.handle.net/1721.1/1301

    .. [3] P. Benner, "Symplectic Balancing of Hamiltonian Matrices", 2001,
       SIAM J. Sci. Comput., 2001, Vol.22(5), DOI: 10.1137/S1064827500367993

    """

    # Validate input arguments
    a, b, q, r, e, s, m, n, r_or_c, gen_are = _are_validate_args(
                                                     a, b, q, r, e, s, 'care')

    H = np.empty((2*m+n, 2*m+n), dtype=r_or_c)
    H[:m, :m] = a
    H[:m, m:2*m] = 0.
    H[:m, 2*m:] = b
    H[m:2*m, :m] = -q
    H[m:2*m, m:2*m] = -a.conj().T
    H[m:2*m, 2*m:] = 0. if s is None else -s
    H[2*m:, :m] = 0. if s is None else s.conj().T
    H[2*m:, m:2*m] = b.conj().T
    H[2*m:, 2*m:] = r

    if gen_are:
        J = block_diag(e, e.conj().T, np.zeros_like(r, dtype=r_or_c))
    else:
        J = block_diag(np.eye(2*m), np.zeros_like(r, dtype=r_or_c))

    if balanced:
        # xGEBAL does not remove the diagonals before scaling. Also
        # to avoid destroying the Symplectic structure, we follow Ref.3
        M = np.abs(H) + np.abs(J)
        M[np.diag_indices_from(M)] = 0.
        _, (sca, _) = matrix_balance(M, separate=1, permute=0)
        # do we need to bother?
        if not np.allclose(sca, np.ones_like(sca)):
            # Now impose diag(D,inv(D)) from Benner where D is
            # square root of s_i/s_(n+i) for i=0,....
            sca = np.log2(sca)
            # NOTE: Py3 uses "Bankers Rounding: round to the nearest even" !!
            s = np.round((sca[m:2*m] - sca[:m])/2)
            sca = 2 ** np.r_[s, -s, sca[2*m:]]
            # Elementwise multiplication via broadcasting.
            elwisescale = sca[:, None] * np.reciprocal(sca)
            H *= elwisescale
            J *= elwisescale

    # Deflate the pencil to 2m x 2m ala Ref.1, eq.(55)
    q, r = qr(H[:, -n:])
    H = q[:, n:].conj().T.dot(H[:, :2*m])
    J = q[:2*m, n:].conj().T.dot(J[:2*m, :2*m])

    # Decide on which output type is needed for QZ
    out_str = 'real' if r_or_c == float else 'complex'

    _, _, _, _, _, u = ordqz(H, J, sort='lhp', overwrite_a=True,
                             overwrite_b=True, check_finite=False,
                             output=out_str)

    # Get the relevant parts of the stable subspace basis
    if gen_are:
        u, _ = qr(np.vstack((e.dot(u[:m, :m]), u[m:, :m])))
    u00 = u[:m, :m]
    u10 = u[m:, :m]

    # Solve via back-substituion after checking the condition of u00
    up, ul, uu = lu(u00)
    if 1/cond(uu) < np.spacing(1.):
        raise LinAlgError('Failed to find a finite solution.')

    # Exploit the triangular structure
    x = solve_triangular(ul.conj().T,
                         solve_triangular(uu.conj().T,
                                          u10.conj().T,
                                          lower=True),
                         unit_diagonal=True,
                         ).conj().T.dot(up.conj().T)
    if balanced:
        x *= sca[:m, None] * sca[:m]

    # Check the deviation from symmetry for success
    u_sym = u00.conj().T.dot(u10)
    n_u_sym = norm(u_sym, 1)
    u_sym = u_sym - u_sym.conj().T
    sym_threshold = np.max([np.spacing(1000.), n_u_sym])

    if norm(u_sym, 1) > sym_threshold:
        raise LinAlgError('The associated Hamiltonian pencil has eigenvalues '
                          'too close to the imaginary axis')

    return (x + x.conj().T)/2

Example 45

Project: scipy
Source File: slsqp.py
View license
def _minimize_slsqp(func, x0, args=(), jac=None, bounds=None,
                    constraints=(),
                    maxiter=100, ftol=1.0E-6, iprint=1, disp=False,
                    eps=_epsilon, callback=None,
                    **unknown_options):
    """
    Minimize a scalar function of one or more variables using Sequential
    Least SQuares Programming (SLSQP).

    Options
    -------
    ftol : float
        Precision goal for the value of f in the stopping criterion.
    eps : float
        Step size used for numerical approximation of the jacobian.
    disp : bool
        Set to True to print convergence messages. If False,
        `verbosity` is ignored and set to 0.
    maxiter : int
        Maximum number of iterations.

    """
    _check_unknown_options(unknown_options)
    fprime = jac
    iter = maxiter
    acc = ftol
    epsilon = eps

    if not disp:
        iprint = 0

    # Constraints are triaged per type into a dictionnary of tuples
    if isinstance(constraints, dict):
        constraints = (constraints, )

    cons = {'eq': (), 'ineq': ()}
    for ic, con in enumerate(constraints):
        # check type
        try:
            ctype = con['type'].lower()
        except KeyError:
            raise KeyError('Constraint %d has no type defined.' % ic)
        except TypeError:
            raise TypeError('Constraints must be defined using a '
                            'dictionary.')
        except AttributeError:
            raise TypeError("Constraint's type must be a string.")
        else:
            if ctype not in ['eq', 'ineq']:
                raise ValueError("Unknown constraint type '%s'." % con['type'])

        # check function
        if 'fun' not in con:
            raise ValueError('Constraint %d has no function defined.' % ic)

        # check jacobian
        cjac = con.get('jac')
        if cjac is None:
            # approximate jacobian function.  The factory function is needed
            # to keep a reference to `fun`, see gh-4240.
            def cjac_factory(fun):
                def cjac(x, *args):
                    return approx_jacobian(x, fun, epsilon, *args)
                return cjac
            cjac = cjac_factory(con['fun'])

        # update constraints' dictionary
        cons[ctype] += ({'fun': con['fun'],
                         'jac': cjac,
                         'args': con.get('args', ())}, )

    exit_modes = {-1: "Gradient evaluation required (g & a)",
                    0: "Optimization terminated successfully.",
                    1: "Function evaluation required (f & c)",
                    2: "More equality constraints than independent variables",
                    3: "More than 3*n iterations in LSQ subproblem",
                    4: "Inequality constraints incompatible",
                    5: "Singular matrix E in LSQ subproblem",
                    6: "Singular matrix C in LSQ subproblem",
                    7: "Rank-deficient equality constraint subproblem HFTI",
                    8: "Positive directional derivative for linesearch",
                    9: "Iteration limit exceeded"}

    # Wrap func
    feval, func = wrap_function(func, args)

    # Wrap fprime, if provided, or approx_jacobian if not
    if fprime:
        geval, fprime = wrap_function(fprime, args)
    else:
        geval, fprime = wrap_function(approx_jacobian, (func, epsilon))

    # Transform x0 into an array.
    x = asfarray(x0).flatten()

    # Set the parameters that SLSQP will need
    # meq, mieq: number of equality and inequality constraints
    meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['eq']]))
    mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args'])) for c in cons['ineq']]))
    # m = The total number of constraints
    m = meq + mieq
    # la = The number of constraints, or 1 if there are no constraints
    la = array([1, m]).max()
    # n = The number of independent variables
    n = len(x)

    # Define the workspaces for SLSQP
    n1 = n + 1
    mineq = m - meq + n1 + n1
    len_w = (3*n1+m)*(n1+1)+(n1-meq+1)*(mineq+2) + 2*mineq+(n1+mineq)*(n1-meq) \
            + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1
    len_jw = mineq
    w = zeros(len_w)
    jw = zeros(len_jw)

    # Decompose bounds into xl and xu
    if bounds is None or len(bounds) == 0:
        xl = np.empty(n, dtype=float)
        xu = np.empty(n, dtype=float)
        xl.fill(np.nan)
        xu.fill(np.nan)
    else:
        bnds = array(bounds, float)
        if bnds.shape[0] != n:
            raise IndexError('SLSQP Error: the length of bounds is not '
                             'compatible with that of x0.')

        bnderr = where(bnds[:, 0] > bnds[:, 1])[0]
        if bnderr.any():
            raise ValueError('SLSQP Error: lb > ub in bounds %s.' %
                             ', '.join(str(b) for b in bnderr))
        xl, xu = bnds[:, 0], bnds[:, 1]

        # Mark infinite bounds with nans; the Fortran code understands this
        infbnd = ~isfinite(bnds)
        xl[infbnd[:, 0]] = np.nan
        xu[infbnd[:, 1]] = np.nan

    # Initialize the iteration counter and the mode value
    mode = array(0,int)
    acc = array(acc,float)
    majiter = array(iter,int)
    majiter_prev = 0

    # Print the header if iprint >= 2
    if iprint >= 2:
        print("%5s %5s %16s %16s" % ("NIT","FC","OBJFUN","GNORM"))

    while 1:

        if mode == 0 or mode == 1:  # objective and constraint evaluation requird

            # Compute objective function
            try:
                fx = float(np.asarray(func(x)))
            except:
                raise ValueError("Objective function must return a scalar")
            # Compute the constraints
            if cons['eq']:
                c_eq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['eq']])
            else:
                c_eq = zeros(0)
            if cons['ineq']:
                c_ieq = concatenate([atleast_1d(con['fun'](x, *con['args']))
                                     for con in cons['ineq']])
            else:
                c_ieq = zeros(0)

            # Now combine c_eq and c_ieq into a single matrix
            c = concatenate((c_eq, c_ieq))

        if mode == 0 or mode == -1:  # gradient evaluation required

            # Compute the derivatives of the objective function
            # For some reason SLSQP wants g dimensioned to n+1
            g = append(fprime(x),0.0)

            # Compute the normals of the constraints
            if cons['eq']:
                a_eq = vstack([con['jac'](x, *con['args'])
                               for con in cons['eq']])
            else:  # no equality constraint
                a_eq = zeros((meq, n))

            if cons['ineq']:
                a_ieq = vstack([con['jac'](x, *con['args'])
                                for con in cons['ineq']])
            else:  # no inequality constraint
                a_ieq = zeros((mieq, n))

            # Now combine a_eq and a_ieq into a single a matrix
            if m == 0:  # no constraints
                a = zeros((la, n))
            else:
                a = vstack((a_eq, a_ieq))
            a = concatenate((a,zeros([la,1])),1)

        # Call SLSQP
        slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)

        # call callback if major iteration has incremented
        if callback is not None and majiter > majiter_prev:
            callback(x)

        # Print the status of the current iterate if iprint > 2 and the
        # major iteration has incremented
        if iprint >= 2 and majiter > majiter_prev:
            print("%5i %5i % 16.6E % 16.6E" % (majiter,feval[0],
                                               fx,linalg.norm(g)))

        # If exit mode is not -1 or 1, slsqp has completed
        if abs(mode) != 1:
            break

        majiter_prev = int(majiter)

    # Optimization loop complete.  Print status if requested
    if iprint >= 1:
        print(exit_modes[int(mode)] + "    (Exit mode " + str(mode) + ')')
        print("            Current function value:", fx)
        print("            Iterations:", majiter)
        print("            Function evaluations:", feval[0])
        print("            Gradient evaluations:", geval[0])

    return OptimizeResult(x=x, fun=fx, jac=g[:-1], nit=int(majiter),
                          nfev=feval[0], njev=geval[0], status=int(mode),
                          message=exit_modes[int(mode)], success=(mode == 0))

Example 46

Project: iris
Source File: trajectory.py
View license
def interpolate(cube, sample_points, method=None):
    """
    Extract a sub-cube at the given n-dimensional points.

    Args:

    * cube
        The source Cube.

    * sample_points
        A sequence of coordinate (name) - values pairs.

    Kwargs:

    * method
        Request "linear" interpolation (default) or "nearest" neighbour.
        Only nearest neighbour is available when specifying multi-dimensional coordinates.


    For example::
    
        sample_points = [('latitude', [45, 45, 45]), ('longitude', [-60, -50, -40])]
        interpolated_cube = interpolate(cube, sample_points)

    """
    if method not in [None, "linear", "nearest"]:
        raise ValueError("Unhandled interpolation specified : %s" % method)

    # Convert any coordinate names to coords
    points = []
    for coord, values in sample_points:
        if isinstance(coord, six.string_types):
            coord = cube.coord(coord)
        points.append((coord, values))
    sample_points = points

    # Do all value sequences have the same number of values?
    coord, values = sample_points[0]
    trajectory_size = len(values)
    for coord, values in sample_points[1:]:
        if len(values) != trajectory_size:
            raise ValueError('Lengths of coordinate values are inconsistent.')

    # Which dimensions are we squishing into the last dimension?
    squish_my_dims = set()
    for coord, values in sample_points:
        dims = cube.coord_dims(coord)
        for dim in dims:
            squish_my_dims.add(dim)

    # Derive the new cube's shape by filtering out all the dimensions we're about to sample,
    # and then adding a new dimension to accommodate all the sample points.
    remaining = [(dim, size) for dim, size in enumerate(cube.shape) if dim not in squish_my_dims]
    new_data_shape = [size for dim, size in remaining]
    new_data_shape.append(trajectory_size)

    # Start with empty data and then fill in the "column" of values for each trajectory point.
    new_cube = iris.cube.Cube(np.empty(new_data_shape))
    new_cube.metadata = cube.metadata

    # Derive the mapping from the non-trajectory source dimensions to their
    # corresponding destination dimensions.
    remaining_dims = [dim for dim, size in remaining]
    dimension_remap = {dim: i for i, dim in enumerate(remaining_dims)}

    # Record a mapping from old coordinate IDs to new coordinates,
    # for subsequent use in creating updated aux_factories.
    coord_mapping = {}

    # Create all the non-squished coords
    for coord in cube.dim_coords:
        src_dims = cube.coord_dims(coord)
        if squish_my_dims.isdisjoint(src_dims):
            dest_dims = [dimension_remap[dim] for dim in src_dims]
            new_coord = coord.copy()
            new_cube.add_dim_coord(new_coord, dest_dims)
            coord_mapping[id(coord)] = new_coord
    for coord in cube.aux_coords:
        src_dims = cube.coord_dims(coord)
        if squish_my_dims.isdisjoint(src_dims):
            dest_dims = [dimension_remap[dim] for dim in src_dims]
            new_coord = coord.copy()
            new_cube.add_aux_coord(new_coord, dest_dims)
            coord_mapping[id(coord)] = new_coord

    # Create all the squished (non derived) coords, not filled in yet.
    trajectory_dim = len(remaining_dims)
    for coord in cube.dim_coords + cube.aux_coords:
        src_dims = cube.coord_dims(coord)
        if not squish_my_dims.isdisjoint(src_dims):
            points = np.array([coord.points.flatten()[0]] * trajectory_size)
            new_coord = iris.coords.AuxCoord(points,
                                             standard_name=coord.standard_name,
                                             long_name=coord.long_name,
                                             units=coord.units,
                                             bounds=None,
                                             attributes=coord.attributes,
                                             coord_system=coord.coord_system)
            new_cube.add_aux_coord(new_coord, trajectory_dim)
            coord_mapping[id(coord)] = new_coord

    for factory in cube.aux_factories:
        new_cube.add_aux_factory(factory.updated(coord_mapping))

    # Are the given coords all 1-dimensional? (can we do linear interp?)
    for coord, values in sample_points:
        if coord.ndim > 1:
            if method == "linear":
                raise iris.exceptions.CoordinateMultiDimError("Cannot currently perform linear interpolation for multi-dimensional coordinates.")
            method = "nearest"
            break

    # Use a cache with _nearest_neighbour_indices_ndcoords()
    cache = {}

    for i in range(trajectory_size):
        point = [(coord, values[i]) for coord, values in sample_points]

        if method in ["linear", None]:
            column = linear_regrid(cube, point)
            new_cube.data[..., i] = column.data
        elif method == "nearest":
            column_index = _nearest_neighbour_indices_ndcoords(cube, point, cache=cache)
            column = cube[column_index]
            new_cube.data[..., i] = column.data

        # Fill in the empty squashed (non derived) coords.
        for column_coord in column.dim_coords + column.aux_coords:
            src_dims = cube.coord_dims(column_coord)
            if not squish_my_dims.isdisjoint(src_dims):
                if len(column_coord.points) != 1:
                    raise Exception("Expected to find exactly one point. Found %d" % len(column_coord.points))
                new_cube.coord(column_coord.name()).points[i] = column_coord.points[0]

    return new_cube

Example 47

Project: iris
Source File: _interpolate_private.py
View license
def regrid(source_cube, grid_cube, mode='bilinear', **kwargs):
    """
    Returns a new cube with values derived from the source_cube on the horizontal grid specified
    by the grid_cube.

    Fundamental input requirements:
        1) Both cubes must have a CoordSystem.
        2) The source 'x' and 'y' coordinates must not share data dimensions with any other coordinates.

    In addition, the algorithm currently used requires:
        3) Both CS instances must be compatible:
            i.e. of the same type, with the same attribute values, and with compatible coordinates.
        4) No new data dimensions can be created.
        5) Source cube coordinates to map to a single dimension.

    Args:

    * source_cube:
        An instance of :class:`iris.cube.Cube` which supplies the source data and metadata.
    * grid_cube:
        An instance of :class:`iris.cube.Cube` which supplies the horizontal grid definition.

    Kwargs:

    * mode (string):
        Regridding interpolation algorithm to be applied, which may be one of the following:

            * 'bilinear' for bi-linear interpolation (default), see :func:`iris.analysis.interpolate.linear`.
            * 'nearest' for nearest neighbour interpolation.

    Returns:
        A new :class:`iris.cube.Cube` instance.

    .. note::

        The masked status of values are currently ignored.  See :func:\
`~iris.experimental.regrid.regrid_bilinear_rectilinear_src_and_grid`
        for regrid support with mask awareness.

    .. deprecated:: 1.10

        Please use :meth:`iris.cube.Cube.regrid` instead, with an appropriate
        regridding scheme:

        *   For mode='bilinear', simply use the :class:`~iris.analysis.Linear`
            scheme.

        *   For mode='nearest', use the :class:`~iris.analysis.Nearest` scheme,
            with extrapolation_mode='extrapolate', but be aware of the
            following possible differences:

            *   Any missing result points, i.e. those which match source points
                which are masked or NaN, are returned as as NaN values by this
                routine.  The 'Nearest' scheme, however, represents missing
                results as masked points in a masked array.
                *Which* points are missing is unchanged.

            *   Longitude wrapping for this routine is controlled by the
                'circular' property of the x coordinate.
                The 'Nearest' scheme, however, *always* wraps any coords with
                modular units, such as (correctly formed) longitudes.
                Thus, behaviour can be different if "x_coord.circular" is
                False :  In that case, if the original non-longitude-wrapped
                operation is required, it can be replicated by converting all
                X and Y coordinates' units to '1' and removing their coordinate
                systems.

    """
    if mode == 'bilinear':
        scheme = iris.analysis.Linear(**kwargs)
        return source_cube.regrid(grid_cube, scheme)

    # Condition 1
    source_cs = source_cube.coord_system(iris.coord_systems.CoordSystem)
    grid_cs = grid_cube.coord_system(iris.coord_systems.CoordSystem)
    if (source_cs is None) != (grid_cs is None):
        raise ValueError("The source and grid cubes must both have a CoordSystem or both have None.")

    # Condition 2: We can only have one x coordinate and one y coordinate with the source CoordSystem, and those coordinates
    # must be the only ones occupying their respective dimension
    source_x = source_cube.coord(axis='x', coord_system=source_cs)
    source_y = source_cube.coord(axis='y', coord_system=source_cs)

    source_x_dims = source_cube.coord_dims(source_x)
    source_y_dims = source_cube.coord_dims(source_y)

    source_x_dim = None
    if source_x_dims:
        if len(source_x_dims) > 1:
            raise ValueError('The source x coordinate may not describe more than one data dimension.')
        source_x_dim = source_x_dims[0]
        dim_sharers = ', '.join([coord.name() for coord in source_cube.coords(contains_dimension=source_x_dim) if coord is not source_x])
        if dim_sharers:
            raise ValueError('No coordinates may share a dimension (dimension %s) with the x '
                             'coordinate, but (%s) do.' % (source_x_dim, dim_sharers))

    source_y_dim = None
    if source_y_dims:
        if len(source_y_dims) > 1:
            raise ValueError('The source y coordinate may not describe more than one data dimension.')
        source_y_dim = source_y_dims[0]
        dim_sharers = ', '.join([coord.name() for coord in source_cube.coords(contains_dimension=source_y_dim) if coord is not source_y])
        if dim_sharers:
            raise ValueError('No coordinates may share a dimension (dimension %s) with the y '
                             'coordinate, but (%s) do.' % (source_y_dim, dim_sharers))

    if source_x_dim is not None and source_y_dim == source_x_dim:
        raise ValueError('The source x and y coords may not describe the same data dimension.')


    # Condition 3
    # Check for compatible horizontal CSs. Currently that means they're exactly the same except for the coordinate
    # values.
    # The same kind of CS ...
    compatible = (source_cs == grid_cs)
    if compatible:
        grid_x = grid_cube.coord(axis='x', coord_system=grid_cs)
        grid_y = grid_cube.coord(axis='y', coord_system=grid_cs)
        compatible = source_x.is_compatible(grid_x) and \
            source_y.is_compatible(grid_y)
    if not compatible:
        raise ValueError("The new grid must be defined on the same coordinate system, and have the same coordinate "
                         "metadata, as the source.")

    # Condition 4
    if grid_cube.coord_dims(grid_x) and not source_x_dims or \
            grid_cube.coord_dims(grid_y) and not source_y_dims:
        raise ValueError("The new grid must not require additional data dimensions.")

    x_coord = grid_x.copy()
    y_coord = grid_y.copy()


    #
    # Adjust the data array to match the new grid.
    #

    # get the new shape of the data
    new_shape = list(source_cube.shape)
    if source_x_dims:
        new_shape[source_x_dims[0]] = grid_x.shape[0]
    if source_y_dims:
        new_shape[source_y_dims[0]] = grid_y.shape[0]

    new_data = np.empty(new_shape, dtype=source_cube.data.dtype)

    # Prepare the index pattern which will be used to insert a single "column" of data.
    # NB. A "column" is a slice constrained to a single XY point, which therefore extends over *all* the other axes.
    # For an XYZ cube this means a column only extends over Z and corresponds to the normal definition of "column".
    indices = [slice(None, None)] * new_data.ndim

    if mode == 'bilinear':
        # Perform bilinear interpolation, passing through any keywords.
        points_dict = [(source_x, list(x_coord.points)), (source_y, list(y_coord.points))]
        new_data = linear(source_cube, points_dict, **kwargs).data
    else:
        # Perform nearest neighbour interpolation on each column in turn.
        for iy, y in enumerate(y_coord.points):
            for ix, x in enumerate(x_coord.points):
                column_pos = [(source_x,  x), (source_y, y)]
                column_data = extract_nearest_neighbour(source_cube, column_pos).data
                if source_y_dim is not None:
                    indices[source_y_dim] = iy
                if source_x_dim is not None:
                    indices[source_x_dim] = ix
                new_data[tuple(indices)] = column_data

    # Special case to make 0-dimensional results take the same form as NumPy
    if new_data.shape == ():
        new_data = new_data.flat[0]

    # Start with just the metadata and the re-sampled data...
    new_cube = iris.cube.Cube(new_data)
    new_cube.metadata = source_cube.metadata

    # ... and then copy across all the unaffected coordinates.

    # Record a mapping from old coordinate IDs to new coordinates,
    # for subsequent use in creating updated aux_factories.
    coord_mapping = {}

    def copy_coords(source_coords, add_method):
        for coord in source_coords:
            if coord is source_x or coord is source_y:
                continue
            dims = source_cube.coord_dims(coord)
            new_coord = coord.copy()
            add_method(new_coord, dims)
            coord_mapping[id(coord)] = new_coord

    copy_coords(source_cube.dim_coords, new_cube.add_dim_coord)
    copy_coords(source_cube.aux_coords, new_cube.add_aux_coord)

    for factory in source_cube.aux_factories:
        new_cube.add_aux_factory(factory.updated(coord_mapping))

    # Add the new coords
    if source_x in source_cube.dim_coords:
        new_cube.add_dim_coord(x_coord, source_x_dim)
    else:
        new_cube.add_aux_coord(x_coord, source_x_dims)

    if source_y in source_cube.dim_coords:
        new_cube.add_dim_coord(y_coord, source_y_dim)
    else:
        new_cube.add_aux_coord(y_coord, source_y_dims)

    return new_cube

Example 48

Project: sfepy
Source File: dof_info.py
View license
    def map_equations(self, bcs, field, ts, functions, problem=None,
                      warn=False):
        """
        Create the mapping of active DOFs from/to all DOFs.

        Parameters
        ----------
        bcs : Conditions instance
            The Dirichlet or periodic boundary conditions (single
            condition instances). The dof names in the conditions must
            already be canonized.
        field : Field instance
            The field of the variable holding the DOFs.
        ts : TimeStepper instance
            The time stepper.
        functions : Functions instance
            The registered functions.
        problem : Problem instance, optional
            The problem that can be passed to user functions as a context.
        warn : bool, optional
            If True, warn about BC on non-existent nodes.

        Returns
        -------
        active_bcs : set
            The set of boundary conditions active in the current time.

        Notes
        -----
        - Periodic bc: master and slave DOFs must belong to the same
          field (variables can differ, though).
        """
        if bcs is None:
            self.val_ebc = nm.empty((0,), dtype=field.dtype)
            self._init_empty()
            return set()

        eq_ebc = nm.zeros((self.var_di.n_dof,), dtype=nm.int32)
        val_ebc = nm.zeros((self.var_di.n_dof,), dtype=field.dtype)
        master_slave = nm.zeros((self.var_di.n_dof,), dtype=nm.int32)
        chains = []

        active_bcs = set()
        for bc in bcs:
            # Skip conditions that are not active in the current time.
            if not is_active_bc(bc, ts=ts, functions=functions):
                continue

            active_bcs.add(bc.key)

            if isinstance(bc, EssentialBC):
                ntype = 'EBC'
                region = bc.region

            else:
                ntype = 'EPBC'
                region = bc.regions[0]

            if warn:
                clean_msg = ('warning: ignoring nonexistent %s node (%s) in '
                             % (ntype, self.var_di.var_name))
            else:
                clean_msg = None

            # Get master region nodes.
            master_nod_list = field.get_dofs_in_region(region)
            if len(master_nod_list) == 0:
                continue

            if ntype == 'EBC': # EBC.
                dofs, val = bc.dofs
                ##
                # Evaluate EBC values.
                fun = get_condition_value(val, functions, 'EBC', bc.name)
                if isinstance(fun, Function):
                    aux = fun
                    fun = lambda coors: aux(ts, coors,
                                            bc=bc, problem=problem)

                nods, vv = field.set_dofs(fun, region, len(dofs), clean_msg)

                eq = expand_nodes_to_equations(nods, dofs, self.dof_names)
                # Duplicates removed here...
                eq_ebc[eq] = 1
                if vv is not None: val_ebc[eq] = vv

            else: # EPBC.
                region = bc.regions[1]
                slave_nod_list = field.get_dofs_in_region(region)

                nmaster = nm.unique(nm.hstack(master_nod_list))
                # Treat fields not covering the whole domain.
                if nmaster[0] == -1:
                    nmaster = nmaster[1:]

                nslave = nm.unique(nm.hstack(slave_nod_list))
                # Treat fields not covering the whole domain.
                if nslave[0] == -1:
                    nslave = nslave[1:]

                ## print nmaster + 1
                ## print nslave + 1
                if nmaster.shape != nslave.shape:
                    msg = 'EPBC list lengths do not match!\n(%s,\n %s)' %\
                          (nmaster, nslave)
                    raise ValueError(msg)

                if (nmaster.shape[0] == 0) and (nslave.shape[0] == 0):
                    continue

                mcoor = field.get_coor(nmaster)
                scoor = field.get_coor(nslave)

                fun = get_condition_value(bc.match, functions, 'EPBC', bc.name)
                if isinstance(fun, Function):
                    i1, i2 = fun(mcoor, scoor)

                else:
                    i1, i2 = fun

                ## print nm.c_[mcoor[i1], scoor[i2]]
                ## print nm.c_[nmaster[i1], nslave[i2]] + 1

                meq = expand_nodes_to_equations(nmaster[i1], bc.dofs[0],
                                                self.dof_names)
                seq = expand_nodes_to_equations(nslave[i2], bc.dofs[1],
                                                self.dof_names)

                m_assigned = nm.where(master_slave[meq] != 0)[0]
                s_assigned = nm.where(master_slave[seq] != 0)[0]
                if m_assigned.size or s_assigned.size: # Chain EPBC.
                    aux = master_slave[meq[m_assigned]]
                    sgn = nm.sign(aux)
                    om_chain = zip(meq[m_assigned], (aux - sgn) * sgn)
                    chains.extend(om_chain)

                    aux = master_slave[seq[s_assigned]]
                    sgn = nm.sign(aux)
                    os_chain = zip(seq[s_assigned], (aux - sgn) * sgn)
                    chains.extend(os_chain)

                    m_chain = zip(meq[m_assigned], seq[m_assigned])
                    chains.extend(m_chain)

                    msd = nm.setdiff1d(s_assigned, m_assigned)
                    s_chain = zip(meq[msd], seq[msd])
                    chains.extend(s_chain)

                    msa = nm.union1d(m_assigned, s_assigned)
                    ii = nm.setdiff1d(nm.arange(meq.size), msa)
                    master_slave[meq[ii]] = seq[ii] + 1
                    master_slave[seq[ii]] = - meq[ii] - 1

                else:
                    master_slave[meq] = seq + 1
                    master_slave[seq] = - meq - 1

        chains = group_chains(chains)
        resolve_chains(master_slave, chains)

        ii = nm.argwhere(eq_ebc == 1)
        self.eq_ebc = nm.atleast_1d(ii.squeeze())
        self.val_ebc = nm.atleast_1d(val_ebc[ii].squeeze())
        self.master = nm.argwhere(master_slave > 0).squeeze()
        self.slave = master_slave[self.master] - 1

        assert_((self.eq_ebc.shape == self.val_ebc.shape))
        self.eq[self.eq_ebc] = -2
        self.eq[self.master] = -1
        self.eqi = nm.compress(self.eq >= 0, self.eq)
        self.eq[self.eqi] = nm.arange(self.eqi.shape[0], dtype=nm.int32)
        self.eq[self.master] = self.eq[self.slave]
        self.n_eq = self.eqi.shape[0]
        self.n_ebc = self.eq_ebc.shape[0]
        self.n_epbc = self.master.shape[0]

        return active_bcs

Example 49

Project: sfepy
Source File: fields.py
View license
    def evaluate_at(self, coors, source_vals, mode='val', strategy='general',
                    close_limit=0.1, get_cells_fun=None, cache=None,
                    ret_cells=False, ret_status=False, ret_ref_coors=False,
                    verbose=False):
        """
        Evaluate source DOF values corresponding to the field in the given
        coordinates using the field interpolation.

        Parameters
        ----------
        coors : array, shape ``(n_coor, dim)``
            The coordinates the source values should be interpolated into.
        source_vals : array, shape ``(n_nod, n_components)``
            The source DOF values corresponding to the field.
        mode : {'val', 'grad'}, optional
            The evaluation mode: the field value (default) or the field value
            gradient.
        strategy : {'general', 'convex'}, optional
            The strategy for finding the elements that contain the
            coordinates. For convex meshes, the 'convex' strategy might be
            faster than the 'general' one.
        close_limit : float, optional
            The maximum limit distance of a point from the closest
            element allowed for extrapolation.
        get_cells_fun : callable, optional
            If given, a function with signature ``get_cells_fun(coors, cmesh,
            **kwargs)`` returning cells and offsets that potentially contain
            points with the coordinates `coors`. Applicable only when
            `strategy` is 'general'. When not given,
            :func:`get_potential_cells()
            <sfepy.discrete.common.global_interp.get_potential_cells>` is used.
        cache : Struct, optional
            To speed up a sequence of evaluations, the field mesh and other
            data can be cached. Optionally, the cache can also contain the
            reference element coordinates as `cache.ref_coors`, `cache.cells`
            and `cache.status`, if the evaluation occurs in the same
            coordinates repeatedly. In that case the mesh related data are
            ignored. See :func:`Field.get_evaluate_cache()
            <sfepy.discrete.fem.fields_base.FEField.get_evaluate_cache()>`.
        ret_ref_coors : bool, optional
            If True, return also the found reference element coordinates.
        ret_status : bool, optional
            If True, return also the enclosing cell status for each point.
        ret_cells : bool, optional
            If True, return also the cell indices the coordinates are in.
        verbose : bool
            If False, reduce verbosity.

        Returns
        -------
        vals : array
            The interpolated values with shape ``(n_coor, n_components)`` or
            gradients with shape ``(n_coor, n_components, dim)`` according to
            the `mode`. If `ret_status` is False, the values where the status
            is greater than one are set to ``numpy.nan``.
        ref_coors : array
            The found reference element coordinates, if `ret_ref_coors` is True.
        cells : array
            The cell indices, if `ret_ref_coors` or `ret_cells` or `ret_status`
            are True.
        status : array
            The status, if `ret_ref_coors` or `ret_status` are True, with the
            following meaning: 0 is success, 1 is extrapolation within
            `close_limit`, 2 is extrapolation outside `close_limit`, 3 is
            failure, 4 is failure due to non-convergence of the Newton
            iteration in tensor product cells. If close_limit is 0, then for
            the 'general' strategy the status 5 indicates points outside of the
            field domain that had no potential cells.
        """
        from sfepy.discrete.common.global_interp import get_ref_coors
        from sfepy.discrete.common.extmods.crefcoors import evaluate_in_rc

        output('evaluating in %d points...' % coors.shape[0], verbose=verbose)

        ref_coors, cells, status = get_ref_coors(self, coors,
                                                 strategy=strategy,
                                                 close_limit=close_limit,
                                                 get_cells_fun=get_cells_fun,
                                                 cache=cache,
                                                 verbose=verbose)

        tt = time.clock()

        # Interpolate to the reference coordinates.
        if mode == 'val':
            vals = nm.empty((coors.shape[0], source_vals.shape[1], 1),
                            dtype=source_vals.dtype)
            cmode = 0

        elif mode == 'grad':
            vals = nm.empty((coors.shape[0], source_vals.shape[1],
                             coors.shape[1]),
                            dtype=source_vals.dtype)
            cmode = 1

        ctx = self.create_basis_context()

        evaluate_in_rc(vals, ref_coors, cells, status, source_vals,
                       self.get_econn('volume', self.region), cmode, ctx)
        output('interpolation: %f s' % (time.clock()-tt),verbose=verbose)

        output('...done',verbose=verbose)

        if mode == 'val':
            vals.shape = (coors.shape[0], source_vals.shape[1])

        if not ret_status:
            ii = nm.where(status > 1)[0]
            vals[ii] = nm.nan

        if ret_ref_coors:
            return vals, ref_coors, cells, status

        elif ret_status:
            return vals, cells, status

        elif ret_cells:
            return vals, cells

        else:
            return vals

Example 50

Project: SHARPpy
Source File: thermo.py
View license
    def drawConvectiveIndices(self, qp):
        '''
        This handles the drawing of the parcel indices.
        
        Parameters
        ----------
        qp: QtGui.QPainter object
        
        '''
        ## initialize a white pen with thickness 2 and a solid line
        pen = QtGui.QPen(QtCore.Qt.white, 1, QtCore.Qt.SolidLine)
        qp.setPen(pen)
        qp.setFont(self.label_font)
        ## make the initial x pixel coordinate relative to the frame
        ## width.
        x1 = self.brx / 8
        y1 = self.ylast + self.tpad
        ## get the indices rounded to the nearest int, conver to strings
        ## Start with the surface based parcel.
        capes = np.empty(4, dtype="S10") # Only allow 4 parcels to be displayed
        cins = np.empty(4, dtype="S10")
        lcls = np.empty(4, dtype="S10")
        lis = np.empty(4, dtype="S10")
        lfcs = np.empty(4, dtype="S10")
        els = np.empty(4, dtype="S10")
        texts = self.pcl_types
        for i, key in enumerate(texts):
            parcel = self.parcels[key]
            capes[i] = tab.utils.INT2STR(parcel.bplus) # Append the CAPE value
            cins[i] = tab.utils.INT2STR(parcel.bminus)
            lcls[i] = tab.utils.INT2STR(parcel.lclhght)
            lis[i] = tab.utils.INT2STR(parcel.li5)
            lfcs[i] = tab.utils.INT2STR(parcel.lfchght)
            els[i] = tab.utils.INT2STR(parcel.elhght)
        ## Now that we have all the data, time to plot the text in their
        ## respective columns.
        
        ## PCL type
        pcl_index = 0
        for i, text in enumerate(texts):
            self.bounds[i,0] = y1
            if text == self.pcl_types[self.skewt_pcl]:
                pen = QtGui.QPen(QtCore.Qt.cyan, 1, QtCore.Qt.SolidLine)
                qp.setPen(pen)
                pcl_index = i
            else:
                pen = QtGui.QPen(QtCore.Qt.white, 1, QtCore.Qt.SolidLine)
                qp.setPen(pen)
            rect = QtCore.QRect(0, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace
            self.bounds[i,1] = y1
        ## CAPE
        y1 = self.ylast + self.tpad
        for i, text in enumerate(capes):
            try:
                if pcl_index == i and int(text) >= 4000:
                    color = QtCore.Qt.magenta
                elif pcl_index == i and int(text) >= 3000:
                    color=QtCore.Qt.red
                elif pcl_index == i and int(text) >= 2000:
                    color=QtCore.Qt.yellow
                else:
                    color=QtCore.Qt.white
            except:
                color=QtCore.Qt.white
            pen = QtGui.QPen(color, 1, QtCore.Qt.SolidLine)
            qp.setPen(pen)
            rect = QtCore.QRect(x1*1, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace
        ## CINH
        y1 = self.ylast + self.tpad
        for i, text in enumerate(cins):
            try:
                if int(capes[i]) > 0 and pcl_index == i and int(text) >= -50:
                    color = QtCore.Qt.green
                elif int(capes[i]) > 0 and pcl_index == i and int(text) >= -100:
                    color=QtGui.QColor('#996600')
                elif int(capes[i]) > 0 and pcl_index == i and int(text) < -100:
                    color=QtGui.QColor('#993333')
                else:
                    color=QtCore.Qt.white
            except:
                color=QtCore.Qt.white
            pen = QtGui.QPen(color, 1, QtCore.Qt.SolidLine)
            qp.setPen(pen)
            rect = QtCore.QRect(x1*2, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace

        ## LCL
        y1 = self.ylast + self.tpad
        for i, text in enumerate(lcls):
            try:
                if int(text) < 1000 and pcl_index == i and texts[i] == "ML":
                    color = QtCore.Qt.green
                elif int(text) < 1500 and pcl_index == i and texts[i] == "ML":
                    color=QtGui.QColor('#996600')
                elif int(text) < 2000 and pcl_index == i and texts[i] == "ML":
                    color=QtGui.QColor('#993333')
                else:
                    color=QtCore.Qt.white
            except:
                color=QtCore.Qt.white
            pen = QtGui.QPen(color, 1, QtCore.Qt.SolidLine)
            qp.setPen(pen)
            rect = QtCore.QRect(x1*3, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace

        pen = QtGui.QPen(QtCore.Qt.white, 1, QtCore.Qt.SolidLine)
        qp.setPen(pen)
        ## LI
        y1 = self.ylast + self.tpad
        for text in lis:
            rect = QtCore.QRect(x1*4, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace
        ## LFC
        y1 = self.ylast + self.tpad
        for text in lfcs:
            rect = QtCore.QRect(x1*5, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace
        ## EL
        y1 = self.ylast + self.tpad
        for text in els:
            rect = QtCore.QRect(x1*6, y1, x1*2, self.label_height)
            qp.drawText(rect, QtCore.Qt.TextDontClip | QtCore.Qt.AlignCenter, text)
            vspace = self.label_height
            if platform.system() == "Windows":
                vspace += self.label_metrics.descent()
            y1 += vspace
            self.ylast = y1
        qp.drawLine(0, y1+2, self.brx, y1+2)
        color=QtGui.QColor('#996633')
        pen = QtGui.QPen(color, .5, QtCore.Qt.SolidLine)
        qp.setPen(pen)
        qp.drawLine(2, self.bounds[self.skewt_pcl,0]-1, x1*6 + x1*2 - 1, self.bounds[self.skewt_pcl,0]-1)
        qp.drawLine(2, self.bounds[self.skewt_pcl,0]-1, 2, self.bounds[self.skewt_pcl,1])
        qp.drawLine(2, self.bounds[self.skewt_pcl,1], x1*6 + x1*2 - 1, self.bounds[self.skewt_pcl,1])
        qp.drawLine(x1*6 + x1*2 -1 , self.bounds[self.skewt_pcl,0]-1, x1*6 + x1*2 -1, self.bounds[self.skewt_pcl,1])