numpy.dot

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

171 Examples 7

Example 1

Project: apogee Source File: cannon.py
def _polyfit_coeffs(spec,specerr,scatter,labelA,return_cov=False):
    """For a given scatter, return the best-fit coefficients"""
    Y= spec/(specerr**2.+scatter**2.)
    ATY= numpy.dot(labelA.T,Y)
    CiA= labelA*numpy.tile(1./(specerr**2.+scatter**2.),(labelA.shape[1],1)).T
    ATCiA= numpy.dot(labelA.T,CiA)
    ATCiAinv= linalg.inv(ATCiA)
    if return_cov:
        return (numpy.dot(ATCiAinv,ATY),ATCiAinv)
    else:
        return numpy.dot(ATCiAinv,ATY)

Example 2

Project: FreeCAD_drawing_dimensioning Source File: dimensionSvgConstructor.py
def arrowHeadSVG( tipPos, d, L1, L2, W, clr='blue'):
    d2 = numpy.dot( [[ 0, -1],[ 1, 0]], d) #same as rotate2D( d, pi/2 )
    R = numpy.array( [d, d2]).transpose()
    p2 = numpy.dot( R, [ L1,    W/2.0 ]) + tipPos
    p3 = numpy.dot( R, [ L1+L2, 0     ]) + tipPos
    p4 = numpy.dot( R, [ L1,   -W/2.0 ]) + tipPos
    return '<polygon points="%f,%f %f,%f %f,%f %f,%f" style="fill:%s;stroke:%s;stroke-width:0" />' % (tipPos[0], tipPos[1], p2[0], p2[1], p3[0], p3[1], p4[0], p4[1], clr, clr)

Example 3

Project: universalSmashSystem Source File: collisionBox.py
    def ejectionDirections(self, _other, _dx=0, _dy=0):
        test_rect = self.current_ecb.rect.copy()
        test_rect.x += _dx
        test_rect.y += _dy

        distances = directionalDisplacements(test_rect, _other)

        working_list = filter(lambda e: numpy.dot(e[0], e[1]) >= 0, distances)
        reference_list = copy.deepcopy(working_list)
        for element in reference_list:
            working_list = filter(lambda k: numpy.dot(k[0], element[1]) != numpy.dot(element[0], element[1]) or numpy.allclose(k[0], element[0]), working_list)
        return working_list

Example 4

Project: pyscf Source File: icmpspt.py
def icmpspt(mc, pttype="NEVPT2", energyE0=0.0, rdmM=0, frozen=0, PTM=1000, PTincore=False, fciExtraLine=[], have3RDM=False, root=0, nroots=1, verbose=None, AAAVsplit=1):

    #remove the -1 state
    import os
    os.system("rm %s/node0/Rotation*.state-1.tmp"%(mc.fcisolver.scratchDirectory))
    os.system("rm %s/node0/wave*.-1.tmp"%(mc.fcisolver.scratchDirectory))

#    if type(mc.fcisolver) is not dmrgci.DMRGCI:
#        if (mc.fcisolver.fcibase_class is not dmrgci.DMRGCI):
#            print "this works with dmrgscf and not regular mcscf"
#            exit(0)

    if (pttype != "NEVPT2" and AAAVsplit != 1):
        print "AAAVsplit only works with CASSCF natural orbitals and NEVPT2"
        exit(0)

    mc.fcisolver.startM = 100
    mc.fcisolver.maxM = max(rdmM,501)
    mc.fcisolver.clearSchedule()
    mc.fcisolver.restart = False

    if (not have3RDM):
        mc.fcisolver.has_threepdm = False

        #we will redo the calculations so lets get ride of -1 states
        import os
        os.system("rm %s/node0/Rotation-*.state-1.tmp"%(mc.fcisolver.scratchDirectory))
        os.system("rm %s/node0/wave-*.-1.tmp"%(mc.fcisolver.scratchDirectory))
        os.system("rm %s/node0/RestartReorder.dat_1"%(mc.fcisolver.scratchDirectory))
    else:
        mc.fcisolver.has_threepdm = True

    mc.fcisolver.generate_schedule()
    mc.fcisolver.extraline = []
    if (PTincore):
        mc.fcisolver.extraline.append('do_npdm_in_core')
    mc.fcisolver.extraline += fciExtraLine


    if (len(mc.fcisolver.orbsym) == 0 and mc.fcisolver.mol.symmetry):
        mcscf.casci_symm.label_symmetry_(mc, mc.mo_coeff)
    ericas = mc.get_h2cas()
    h1e = reduce(numpy.dot, (mc.mo_coeff.T, mc.get_hcore(), mc.mo_coeff))
    dmcore = numpy.dot(mc.mo_coeff[:,:mc.ncore], mc.mo_coeff[:,:mc.ncore].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore = reduce(numpy.dot, (mc.mo_coeff.T, vj-vk*0.5, mc.mo_coeff))
    h1effcas = h1e+vhfcore

    dmrgci.writeIntegralFile(mc.fcisolver, h1effcas[mc.ncore:mc.ncore+mc.ncas, mc.ncore:mc.ncore+mc.ncas], ericas, mc.ncas, mc.nelecas)

    dm1eff = numpy.zeros(shape=(mc.ncas, mc.ncas)) #this is the state average density which is needed in NEVPT2
 
    #loop over all states besides the current root
    if (pttype == "NEVPT2" and nroots>1):
        stateIter = range(nroots)
        stateIter.remove(root)
        for istate in stateIter:
            dm3 = mc.fcisolver.make_rdm3(state=istate, norb=mc.ncas, nelec=mc.nelecas, dt=float_precision)    
            nelec = mc.nelecas[0]+mc.nelecas[1]
            dm2 = numpy.einsum('ijklmk', dm3)/(nelec-2)
            dm1 = numpy.einsum('ijkj', dm2)/(nelec-1)
            dm1eff += dm1

    #now add the contributaion due to the current root
    dm3 = mc.fcisolver.make_rdm3(state=root, norb=mc.ncas, nelec=mc.nelecas, dt=float_precision)    
    nelec = mc.nelecas[0]+mc.nelecas[1]
    dm2 = numpy.einsum('ijklmk', dm3)/(nelec-2)
    dm1 = numpy.einsum('ijkj', dm2)/(nelec-1)
    dm1eff += dm1
    dm1eff = dm1eff/(1.0*nroots)
    import os
    os.system("mkdir int")    
    numpy.save("int/E3",dm3)
    numpy.save("int/E3B.npy", dm3.transpose(0,3,1,4,2,5))
    numpy.save("int/E3C.npy", dm3.transpose(5,0,2,4,1,3))
    del dm3

    #backup the restartreorder file to -1. this is because responseaaav and responseaaac both overwrite this file
    #this means that when we want to restart a calculation after lets say responseaaav didnt finish, the new calculaitons
    #will use the restartreorder file that was written by the incomplete responseaaav run instead of the original dmrg run.
    reorderf1 = "%s/node0/RestartReorder.dat_1"%(mc.fcisolver.scratchDirectory)
    reorderf = "%s/node0/RestartReorder.dat"%(mc.fcisolver.scratchDirectory)
    import os.path
    reorder1present = os.path.isfile(reorderf1) 
    if (reorder1present):
        from subprocess import check_call
        output = check_call("cp %s %s"%(reorderf1, reorderf), shell=True)
    else :
        from subprocess import check_call
        check_call("cp %s %s"%(reorderf, reorderf1), shell=True)
    reorder = numpy.loadtxt("%s/node0/RestartReorder.dat"%(mc.fcisolver.scratchDirectory))


    if (pttype == "NEVPT2") :
        norbs, energyE0 = writeNevpt2Integrals(mc, dm1, dm2, dm1eff, AAAVsplit, frozen)
        sys.stdout.flush()
        print "wrote the integrals to disk"

        for k in range(AAAVsplit):
            writeDMRGConfFile(mc.nelecas[0], mc.nelecas[1], mc.ncore, mc.ncas,  norbs,
                              mc.fcisolver, PTM, "AAAV", mc.fcisolver.memory, mc.fcisolver.num_thrds, reorder, fciExtraLine, aaavsplit=AAAVsplit, aaavIter=k, root=root, name = "NEVPT2")
        writeDMRGConfFile(mc.nelecas[0], mc.nelecas[1], mc.ncore-frozen, mc.ncas,  norbs-frozen,
                          mc.fcisolver, PTM, "AAAC", mc.fcisolver.memory, mc.fcisolver.num_thrds, reorder, fciExtraLine,root=root, name = "NEVPT2")
        sys.stdout.flush()

        totalE = 0.0;
        totalE += executeNEVPT(nelec, mc.ncore, mc.ncas, frozen, mc.fcisolver.memory)# executeMRLCC(nelec, mc.ncore, mc.ncas)

        try:
            for k in range(AAAVsplit):
                outfile, infile = "responseNEVPT2_aaav%d.out"%(k), "responseNEVPT2_aaav%d.conf"%(k)
                output = check_call("%s  %s  %s > %s"%(mc.fcisolver.mpiprefix, mc.fcisolver.executable, infile, outfile), shell=True)
                file1 = open("%s/node0/dmrg.e"%(mc.fcisolver.scratchDirectory),"rb")
                import struct
                energy = struct.unpack('d', file1.read(8))[0]
                file1.close()
                totalE += energy
                print "perturber AAAV%i --  %18.9e"%(k, energy)
                sys.stdout.flush()

            if (mc.ncore-frozen != 0):
                outfile, infile = "responseNEVPT2_aaac.out", "responseNEVPT2_aaac.conf"
                output = check_call("%s  %s  %s > %s"%(mc.fcisolver.mpiprefix, mc.fcisolver.executable, infile, outfile), shell=True)
                file1 = open("%s/node0/dmrg.e"%(mc.fcisolver.scratchDirectory),"rb")
                energy = struct.unpack('d', file1.read(8))[0]
                file1.close()
                totalE += energy
                print "perturber AAAC --  %18.9e"%(energy)

        except ValueError:
            print(output)

        from subprocess import check_call
        return totalE
    else :
        #this is a bad way to do it, the problem is
        #that pyscf works with double precision and
        #
        #energyE0 = writeMRLCCIntegrals(mc, dm1, dm2)
        #sys.stdout.flush()
        energyE0, norbs = writeNumpyforMRLCC(mc, dm1, dm2, frozen) 
        sys.stdout.flush()
        writeDMRGConfFile(mc.nelecas[0], mc.nelecas[1], mc.ncore, mc.ncas,  norbs,
                          mc.fcisolver, PTM, "AAAV", mc.fcisolver.memory, mc.fcisolver.num_thrds, reorder, fciExtraLine, root=root, name="MRLCC")
        writeDMRGConfFile(mc.nelecas[0], mc.nelecas[1], mc.ncore-frozen, mc.ncas,  norbs-frozen,
                          mc.fcisolver, PTM, "AAAC", mc.fcisolver.memory, mc.fcisolver.num_thrds, reorder, fciExtraLine, root=root, name="MRLCC")
        totalE = 0.0
        totalE +=  executeMRLCC(nelec, mc.ncore, mc.ncas, frozen, mc.fcisolver.memory)
        from subprocess import check_call
        try:
            outfile, infile = "responseMRLCC_aaav0.out", "responseMRLCC_aaav0.conf"
            output = check_call("%s  %s  %s > %s"%(mc.fcisolver.mpiprefix, mc.fcisolver.executable, infile, outfile), shell=True)
            file1 = open("%s/node0/dmrg.e"%(mc.fcisolver.scratchDirectory),"rb")
            import struct
            energy = struct.unpack('d', file1.read(8))[0]
            file1.close()
            totalE += energy
            print "perturber AAAV --  %18.9e"%(energy)
        except ValueError:
            print "perturber AAAV -- NA"
            #exit()

        try:
            if (mc.ncore-frozen != 0):
                outfile, infile = "responseMRLCC_aaac.out", "responseMRLCC_aaac.conf"
                output = check_call("%s  %s  %s > %s"%(mc.fcisolver.mpiprefix, mc.fcisolver.executable, infile, outfile), shell=True)
                file1 = open("%s/node0/dmrg.e"%(mc.fcisolver.scratchDirectory),"rb")
                energy = struct.unpack('d', file1.read(8))[0]
                file1.close()
                totalE += energy
                print "perturber AAAC --  %18.9e"%(energy)
        except ValueError:
            print "perturber AAAC -- NA"

        print "total PT  -- %18.9e"%(totalE)
        return totalE

Example 5

Project: pyscf Source File: newton_ah.py
def newton(mf):
    import pyscf.scf
    from pyscf.mcscf import mc1step_symm
    if mf.__class__.__doc__ is None:
        doc = ''
    else:
        doc = mf.__class__.__doc__
    class RHF(mf.__class__):
        __doc__ = doc + \
        '''
        Attributes for Newton solver:
            max_cycle_inner : int
                AH iterations within eacy macro iterations. Default is 10
            max_stepsize : int
                The step size for orbital rotation.  Small step is prefered.  Default is 0.05.
            canonicalization : bool
                To control whether to canonicalize the orbitals optimized by
                Newton solver.  Default is True.
        '''
        def __init__(self):
            self._scf = mf
            self.max_cycle_inner = 10
            self.max_stepsize = .05
            self.canonicalization = True

            self.ah_start_tol = 5.
            self.ah_start_cycle = 1
            self.ah_level_shift = 0
            self.ah_conv_tol = 1e-12
            self.ah_lindep = 1e-14
            self.ah_max_cycle = 30
            self.ah_grad_trust_region = 2.5
# * Classic AH can be simulated by setting
#               max_cycle_micro_inner = 1
#               ah_start_tol = 1e-7
#               max_stepsize = 1.5
#               ah_grad_trust_region = 1e6
            self.ah_decay_rate = .8
            self.keyframe_interval = 5
            self.keyframe_interval_rate = 1.
            self_keys = set(self.__dict__.keys())

            self.__dict__.update(mf.__dict__)
            self._keys = self_keys.union(mf._keys)

        def dump_flags(self):
            log = logger.Logger(self.stdout, self.verbose)
            log.info('\n')
            self._scf.dump_flags()
            log.info('cuem**** %s Newton solver flags ********', self._scf.__class__)
            log.info('SCF tol = %g', self.conv_tol)
            log.info('conv_tol_grad = %s',    self.conv_tol_grad)
            log.info('max. SCF cycles = %d', self.max_cycle)
            log.info('direct_scf = %s', self._scf.direct_scf)
            if self._scf.direct_scf:
                log.info('direct_scf_tol = %g', self._scf.direct_scf_tol)
            if self.chkfile:
                log.info('chkfile to save SCF result = %s', self.chkfile)
            log.info('max_cycle_inner = %d',  self.max_cycle_inner)
            log.info('max_stepsize = %g', self.max_stepsize)
            log.info('ah_start_tol = %g',     self.ah_start_tol)
            log.info('ah_level_shift = %g',   self.ah_level_shift)
            log.info('ah_conv_tol = %g',      self.ah_conv_tol)
            log.info('ah_lindep = %g',        self.ah_lindep)
            log.info('ah_start_cycle = %d',   self.ah_start_cycle)
            log.info('ah_max_cycle = %d',     self.ah_max_cycle)
            log.info('ah_grad_trust_region = %g', self.ah_grad_trust_region)
            log.info('keyframe_interval = %d', self.keyframe_interval)
            log.info('keyframe_interval_rate = %g', self.keyframe_interval_rate)
            log.info('augmented hessian decay rate = %g', self.ah_decay_rate)
            log.info('max_memory %d MB (current use %d MB)',
                     self.max_memory, lib.current_memory()[0])

        def get_fock(self, h1e, s1e, vhf, dm, cycle=-1, adiis=None,
                     diis_start_cycle=None, level_shift_factor=None,
                     damp_factor=None):
            return h1e + vhf

        def kernel(self, mo_coeff=None, mo_occ=None):
            if mo_coeff is None:
                mo_coeff = self.mo_coeff
            else:  # save initial guess because some methods may access them
                self.mo_coeff = mo_coeff
            if mo_occ is None:
                mo_occ = self.mo_occ
            else:
                self.mo_occ = mo_occ
            cput0 = (time.clock(), time.time())

            self.build(self.mol)
            self.dump_flags()
            self.converged, self.e_tot, \
                    self.mo_energy, self.mo_coeff, self.mo_occ = \
                    kernel(self, mo_coeff, mo_occ, conv_tol=self.conv_tol,
                           conv_tol_grad=self.conv_tol_grad,
                           max_cycle=self.max_cycle,
                           callback=self.callback, verbose=self.verbose)

            logger.timer(self, 'Second order SCF', *cput0)
            self._finalize()
            return self.e_tot

        def from_dm(self, dm):
            '''Transform density matrix to the initial guess'''
            if isinstance(dm, str):
                sys.stderr.write('Newton solver reads density matrix from chkfile %s\n' % dm)
                dm = self.from_chk(dm, True)
            mol = self._scf.mol
            h1e = self._scf.get_hcore(mol)
            s1e = self._scf.get_ovlp(mol)
            vhf = self._scf.get_veff(mol, dm)
            fock = self.get_fock(h1e, s1e, vhf, dm, 0, None)
            mo_energy, mo_coeff = self.eig(fock, s1e)
            mo_occ = self.get_occ(mo_energy, mo_coeff)
            return mo_coeff, mo_occ

        def gen_g_hop(self, mo_coeff, mo_occ, fock_ao=None, h1e=None):
            mol = self._scf.mol
            if mol.symmetry:
# label the symmetry every time calling gen_g_hop because kernel function may
# change mo_coeff ordering after calling self.eig
                self._orbsym = symm.label_orb_symm(mol, mol.irrep_id,
                        mol.symm_orb, mo_coeff, s=self.get_ovlp())
            return gen_g_hop_rhf(self, mo_coeff, mo_occ, fock_ao)

        def update_rotate_matrix(self, dx, mo_occ, u0=1):
            dr = pyscf.scf.hf.unpack_uniq_var(dx, mo_occ)
            if self._scf.mol.symmetry:
                dr = mc1step_symm._symmetrize(dr, self._orbsym, None)
            return numpy.dot(u0, expmat(dr))

        def rotate_mo(self, mo_coeff, u, log=None):
            mo = numpy.dot(mo_coeff, u)
            if log is not None and log.verbose >= logger.DEBUG:
                s = reduce(numpy.dot, (mo[:,self.mo_occ>0].T, self.get_ovlp(),
                                       self.mo_coeff[:,self.mo_occ>0]))
                log.debug('Overlap to initial guess, SVD = %s',
                          _effective_svd(s, 1e-5))
                log.debug('Overlap to last step, SVD = %s',
                          _effective_svd(u, 1e-5))
            return mo

    if isinstance(mf, pyscf.scf.rohf.ROHF):
        class ROHF(RHF):
            __doc__ = RHF.__doc__
            def gen_g_hop(self, mo_coeff, mo_occ, fock_ao=None, h1e=None):
                mol = self._scf.mol
                if mol.symmetry:
# label the symmetry every time calling gen_g_hop because kernel function may
# change mo_coeff ordering after calling self.eig
                    self._orbsym = symm.label_orb_symm(mol, mol.irrep_id,
                            mol.symm_orb, mo_coeff, s=self.get_ovlp())
                return gen_g_hop_rohf(self, mo_coeff, mo_occ, fock_ao)

            def get_fock(self, h1e, s1e, vhf, dm, cycle=-1, adiis=None,
                         diis_start_cycle=None, level_shift_factor=None,
                         damp_factor=None):
                fock = h1e + vhf
                self._focka_ao = self._scf._focka_ao = fock[0]  # needed by ._scf.eig
                self._fockb_ao = self._scf._fockb_ao = fock[1]  # needed by ._scf.eig
                self._dm_ao = dm  # needed by .eig
                return fock

            def eig(self, fock, s1e):
                f = (self._focka_ao, self._fockb_ao)
                f = pyscf.scf.rohf.get_roothaan_fock(f, self._dm_ao, s1e)
                return self._scf.eig(f, s1e)
                #fc = numpy.dot(fock[0], mo_coeff)
                #mo_energy = numpy.einsum('pk,pk->k', mo_coeff, fc)
                #return mo_energy
        return ROHF()

    elif isinstance(mf, pyscf.scf.uhf.UHF):
        class UHF(RHF):
            __doc__ = RHF.__doc__
            def gen_g_hop(self, mo_coeff, mo_occ, fock_ao=None, h1e=None):
                mol = self._scf.mol
                if mol.symmetry:
                    ovlp_ao = self.get_ovlp()
                    self._orbsym =(symm.label_orb_symm(mol, mol.irrep_id,
                                            mol.symm_orb, mo_coeff[0], s=ovlp_ao),
                                   symm.label_orb_symm(mol, mol.irrep_id,
                                            mol.symm_orb, mo_coeff[1], s=ovlp_ao))
                return gen_g_hop_uhf(self, mo_coeff, mo_occ, fock_ao)

            def update_rotate_matrix(self, dx, mo_occ, u0=1):
                occidxa = mo_occ[0] > 0
                occidxb = mo_occ[1] > 0
                viridxa = ~occidxa
                viridxb = ~occidxb
                nmo = len(occidxa)
                nocca = occidxa.sum()
                nvira = nmo - nocca

                dra = numpy.zeros((nmo,nmo))
                drb = numpy.zeros((nmo,nmo))
                dra[viridxa[:,None]&occidxa] = dx[:nocca*nvira]
                drb[viridxb[:,None]&occidxb] = dx[nocca*nvira:]
                dra = dra - dra.T
                drb = drb - drb.T

                if self._scf.mol.symmetry:
                    dra = mc1step_symm._symmetrize(dra, self._orbsym[0], None)
                    drb = mc1step_symm._symmetrize(drb, self._orbsym[1], None)
                if isinstance(u0, int) and u0 == 1:
                    return numpy.asarray((expmat(dra), expmat(drb)))
                else:
                    return numpy.asarray((numpy.dot(u0[0], expmat(dra)),
                                          numpy.dot(u0[1], expmat(drb))))

            def rotate_mo(self, mo_coeff, u, log=None):
                return numpy.asarray((numpy.dot(mo_coeff[0], u[0]),
                                      numpy.dot(mo_coeff[1], u[1])))
        return UHF()

    elif isinstance(mf, pyscf.scf.dhf.UHF):
        raise RuntimeError('Not support Dirac-HF')

    else:
        return RHF()

Example 6

Project: pyscf Source File: icmpspt.py
def writeNevpt2Integrals(mc, dm1, dm2, dm1eff, aaavsplit, frozen):
    eris = _ERIS(mc, mc.mo_coeff)
    nocc = mc.ncore + mc.ncas
    ncore = mc.ncore
    nact = mc.ncas
    norbs = eris['ppaa'].shape[0]
    nvirt = norbs-ncore-nact
    eris_sp={}
    #eris_sp['cvcv']= numpy.zeros(shape=(1,1)) #numpy.zeros(shape=(mc.ncore, norbs-mc.ncore-mc.ncas, mc.ncore, norbs-mc.ncas-mc.ncore))
    eris_sp['h1eff']= 1.*eris['h1eff'] #numpy.zeros(shape=(norbs, norbs))
    int1 = reduce(numpy.dot, (mc.mo_coeff.T, mc.get_hcore(), mc.mo_coeff)) 
    #eris_sp['h1eff'] = makeheff(int1, eris['papa'], eris['ppaa'], dm1, ncore, nvirt)

    eris_sp['h1eff'][:mc.ncore,:mc.ncore] += numpy.einsum('abcd,cd', eris['ppaa'][:mc.ncore, :mc.ncore,:,:], dm1eff)
    eris_sp['h1eff'][:mc.ncore,:mc.ncore] -= numpy.einsum('abcd,bd', eris['papa'][:mc.ncore, :, :mc.ncore,:], dm1eff)*0.5
    eris_sp['h1eff'][mc.ncas+mc.ncore:,mc.ncas+mc.ncore:] += numpy.einsum('abcd,cd', eris['ppaa'][mc.ncas+mc.ncore:,mc.ncas+mc.ncore:,:,:], dm1eff)
    eris_sp['h1eff'][mc.ncas+mc.ncore:,mc.ncas+mc.ncore:] -= numpy.einsum('abcd,bd', eris['papa'][mc.ncas+mc.ncore:,:,mc.ncas+mc.ncore:,:], dm1eff)*0.5


    offdiagonal = 0.0
    #zero out off diagonal core
    for k in range(mc.ncore):
        for l in range(mc.ncore):
            if(k != l):
                offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))
    #zero out off diagonal virtuals
    for k in range(mc.ncore+mc.ncas, norbs):
        for l in range(mc.ncore+mc.ncas,norbs):
            if(k != l):
                offdiagonal = max(abs(offdiagonal), abs(eris_sp['h1eff'][k,l] ))

    if (abs(offdiagonal) > 1e-6):
        print "WARNING cuem* Have to use natorual orbitals from casscf ****"
        print "offdiagonal elements ", offdiagonal


    eriscvcv = eris['cvcv']
    if (not isinstance(eris['cvcv'], type(eris_sp['h1eff']))):
        #eriscvcv = h5py.File(eris['cvcv'].name,'r')["eri_mo"]
        eriscvcv = lib.chkfile.load(eris['cvcv'].name, "eri_mo")#h5py.File(eris['cvcv'].name,'r')["eri_mo"]
    eris_sp['cvcv'] = eriscvcv.reshape(mc.ncore, norbs-mc.ncore-mc.ncas, mc.ncore, norbs-mc.ncore-mc.ncas)


    import os
    os.system("mkdir int")    

    numpy.save("int/W:caac", numpy.asfortranarray(eris['papa'][frozen:mc.ncore, :, frozen:mc.ncore, :].transpose(0,3,1,2)))
    numpy.save("int/W:aeca", numpy.asfortranarray(eris['papa'][frozen:mc.ncore, :, mc.ncore+mc.ncas:, :].transpose(1,2,0,3)))
    numpy.save("int/W:ccaa", numpy.asfortranarray(eris['papa'][frozen:mc.ncore, :, frozen:mc.ncore, :].transpose(0,2,1,3)))
    numpy.save("int/W:eeaa", numpy.asfortranarray(eris['papa'][mc.ncore+mc.ncas:, :, mc.ncore+mc.ncas:, :].transpose(0,2,1,3)))
    numpy.save("int/W:caca", numpy.asfortranarray(eris['ppaa'][frozen:mc.ncore, frozen:mc.ncore, :, :].transpose(0,2,1,3))) 
    numpy.save("int/W:eaca", numpy.asfortranarray(eris['ppaa'][mc.ncore+mc.ncas:, frozen:mc.ncore, :, :].transpose(0,2,1,3))) 
    numpy.save("int/W:eecc", numpy.asfortranarray(eris_sp['cvcv'][frozen:,:,frozen:,:].transpose(1,3,0,2))) 
    numpy.save("int/W:ccae", numpy.asfortranarray(eris['pacv'][frozen:mc.ncore,:,frozen:,:].transpose(0,2,1,3))) 

    numpy.save("int/W:aaaa", numpy.asfortranarray(eris['ppaa'][mc.ncore:mc.ncore+mc.ncas, mc.ncore:mc.ncore+mc.ncas, :, :].transpose(0,2,1,3)))
    numpy.save("int/W:eeca", numpy.asfortranarray(eris['pacv'][mc.ncore+mc.ncas:, :, frozen:, :].transpose(3,0,2,1)))
    numpy.save("int/int1eff", numpy.asfortranarray(eris_sp['h1eff'][frozen:,frozen:]))
    numpy.save("int/E1.npy", numpy.asfortranarray(dm1))
    numpy.save("int/E2.npy", numpy.asfortranarray(dm2))
    #numpy.save("int/E3",dm3)
    #numpy.save("int/E3B.npy", dm3.transpose(0,3,1,4,2,5))
    #numpy.save("int/E3C.npy", dm3.transpose(5,0,2,4,1,3))

    nc = mc.ncore+mc.ncas
    energyE0 = 1.0*numpy.einsum('ij,ij', eris_sp['h1eff'][ncore:nc, ncore:nc], dm1)
    energyE0 += 0.5*numpy.einsum('ijkl,ijkl', eris['ppaa'][mc.ncore:mc.ncore+mc.ncas, mc.ncore:mc.ncore+mc.ncas, :, :].transpose(0,2,1,3), dm2)
    mo = mc.mo_coeff
    dmcore = numpy.dot(mo[:,:ncore], mo[:,:ncore].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    energyE0 += numpy.einsum('ij,ji', dmcore, mc.get_hcore()) \
                  + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5
    energyE0 += mc.mol.energy_nuc()
    print "Energy = ", energyE0

    from pyscf import symm
    mol = mc.mol
    orbsymout=[]
    orbsym = []
    if (mol.symmetry):
        orbsym = symm.label_orb_symm(mc.mol, mc.mol.irrep_id,
                                     mc.mol.symm_orb, mc.mo_coeff, s=mc._scf.get_ovlp())
    if mol.symmetry and orbsym:
        if mol.groupname.lower() == 'dooh':
            orbsymout = [dmrg_sym.IRREP_MAP['D2h'][i % 10] for i in orbsym]
        elif mol.groupname.lower() == 'coov':
            orbsymout = [dmrg_sym.IRREP_MAP['C2v'][i % 10] for i in orbsym]
        else:
            orbsymout = [dmrg_sym.IRREP_MAP[mol.groupname][i] for i in orbsym]
    else:
        orbsymout = []


    ncore = mc.ncore
    mo = mc.mo_coeff

    dmcore = numpy.dot(mo[:,:ncore], mo[:,:ncore].T)*2
    vj, vk = mc._scf.get_jk(mc.mol, dmcore)
    vhfcore = reduce(numpy.dot, (mo.T, vj-vk*0.5, mo))
    energy_core = numpy.einsum('ij,ji', dmcore, mc.get_hcore()) \
                  + numpy.einsum('ij,ji', dmcore, vj-0.5*vk) * .5
    print energy_core
    virtOrbs = range(ncore+mc.ncas, eris_sp['h1eff'].shape[0])
    chunks = len(virtOrbs)/aaavsplit
    virtRange = [virtOrbs[i:i+chunks] for i in xrange(0, len(virtOrbs), chunks)]
    for K in range(aaavsplit):
        currentOrbs = range(ncore, mc.ncas+ncore)+virtRange[K]
        fout = open('FCIDUMP_aaav%d'%(K),'w')
        #tools.fcidump.write_head(fout, eris_sp['h1eff'].shape[0]-ncore, mol.nelectron-2*ncore, orbsym= orbsymout[ncore:])

        tools.fcidump.write_head(fout, mc.ncas+len(virtRange[K]), mol.nelectron-2*ncore, orbsym= (orbsymout[ncore:ncore+mc.ncas]+orbsymout[virtRange[K][0]:virtRange[K][-1]+1]) )
        ij = ncore*(ncore+1)/2
        for i in range(len(currentOrbs)):
            for j in range(ncore, mc.ncas+ncore):
                for k in range(mc.ncas):
                    for l in range(k+1):
                        I = currentOrbs[i]
                        if abs(eris['ppaa'][I,j,k,l]) > 1.e-8 :
                            fout.write(' %17.9e %4d %4d %4d %4d\n' \
                                           % (eris['ppaa'][I,j,k,l], i+1, j+1-ncore, k+1, l+1))
    
        h1eff = numpy.zeros(shape=(mc.ncas+len(virtRange[K]), mc.ncas+len(virtRange[K])))
        h1eff[:mc.ncas, :mc.ncas] = eris_sp['h1eff'][ncore:ncore+mc.ncas,ncore:ncore+mc.ncas]
        h1eff[mc.ncas:, mc.ncas:] = eris_sp['h1eff'][virtRange[K][0]:virtRange[K][-1]+1, virtRange[K][0]:virtRange[K][-1]+1]
        h1eff[:mc.ncas, mc.ncas:] = eris_sp['h1eff'][ncore:ncore+mc.ncas, virtRange[K][0]:virtRange[K][-1]+1]
        h1eff[mc.ncas:, :mc.ncas] = eris_sp['h1eff'][virtRange[K][0]:virtRange[K][-1]+1, ncore:ncore+mc.ncas]

        tools.fcidump.write_hcore(fout, h1eff, mc.ncas+len(virtRange[K]), tol=1e-8)
        #tools.fcidump.write_hcore(fout, eris_sp['h1eff'][virtRange[K][0]:virtRange[K][-1]+1, virtRange[K][0]:virtRange[K][-1]+1], len(virtRange[K]), tol=1e-8)
        fout.write(' %17.9e  0  0  0  0\n' %( mol.energy_nuc()+energy_core-energyE0))
        fout.close()

    nc = ncore+mc.ncas
    fout = open('FCIDUMP_aaac','w')
    tools.fcidump.write_head(fout, nc-frozen, mol.nelectron-2*frozen, orbsym= orbsymout[frozen:nc])
    for i in range(frozen,nc):
        for j in range(ncore, nc):
            for k in range(ncore, nc):
                for l in range(ncore,k+1):
                    if abs(eris['ppaa'][i,j,k-ncore,l-ncore]) > 1.e-8 :
                        fout.write(' %17.9e %4d %4d %4d %4d\n' \
                                   % (eris['ppaa'][i,j,k-ncore,l-ncore], i+1-frozen, j+1-frozen, k+1-frozen, l+1-frozen))

    dmrge = energyE0-mol.energy_nuc()-energy_core

    ecore_aaac = 0.0;
    for i in range(frozen,ncore):
        ecore_aaac += 2.0*eris_sp['h1eff'][i,i]
    tools.fcidump.write_hcore(fout, eris_sp['h1eff'][frozen:nc,frozen:nc], nc-frozen, tol=1e-8)
    fout.write(' %17.9e  0  0  0  0\n' %( -dmrge-ecore_aaac))
    fout.close()


    return norbs, energyE0

Example 7

Project: pyscf Source File: icosmo.py
def cosmo_for_casci(mc, cosmo):
    oldCAS = mc.__class__
    cosmo.initialization(cosmo.mol)
    if cosmo.dm is not None:
        cosmo._dm_guess = cosmo.dm
        vcosmo = cosmo.cosmo_fock(cosmo.dm)

    class CAS(oldCAS):
        def __init__(self):
            self.__dict__.update(mc.__dict__)

        def get_hcore(self, mol=None):
            if cosmo.dm is not None:
                v1 = vcosmo
            else:
                if cosmo._dm_guess is None:  # Initial guess
                    na = self.ncore + self.nelecas[0]
                    nb = self.ncore + self.nelecas[1]
                    #log.Initial('Initial DM: na,nb,nelec=',na,nb,na+nb)
                    dm =(numpy.dot(self.mo_coeff[:,:na], self.mo_coeff[:,:na].T)
                       + numpy.dot(self.mo_coeff[:,:nb], self.mo_coeff[:,:nb].T))
                else:
                    dm = cosmo._dm_guess
                v1 = cosmo.cosmo_fock(dm)
            cosmo._v = v1
            return self._scf.get_hcore(mol) + v1

        def kernel(self, mo_coeff=None, ci0=None):
            if mo_coeff is None:
                mo_coeff = self.mo_coeff
            else:
                self.mo_coeff = mo_coeff
            if ci0 is None:
                ci0 = self.ci

            if self.verbose >= lib.logger.WARN:
                self.check_sanity()
            self.dump_flags()

            if self.mol.symmetry:
                mcscf.casci_symm.label_symmetry_(self, self.mo_coeff)
            log = lib.logger.Logger(self.stdout, self.verbose)

            def casci_iter(mo_coeff, ci, cycle):
                # casci.kernel call get_hcore, which initialized cosmo._v
                e_tot, e_cas, fcivec = mcscf.casci.kernel(self, mo_coeff,
                                                          ci0=ci0, verbose=log)
                if isinstance(self.e_cas, (float, numpy.number)):
                    casdm1 = self.fcisolver.make_rdm1(fcivec, self.ncas, self.nelecas)
                else:
                    log.debug('COSMO responses to DM of state %d', cosmo.casci_state_id)
                    c = fcivec[cosmo.casci_state_id]
                    casdm1 = self.fcisolver.make_rdm1(c, self.ncas, self.nelecas)

                mocore = mo_coeff[:,:self.ncore]
                mocas = mo_coeff[:,self.ncore:self.ncore+self.ncas]
                cosmo._dm_guess = reduce(numpy.dot, (mocas, casdm1, mocas.T))
                cosmo._dm_guess += numpy.dot(mocore, mocore.T) * 2
                edup = numpy.einsum('ij,ij', cosmo._v, cosmo._dm_guess)
                # Substract <VP> to get E0, then add Ediel
                e_tot = e_tot - edup + cosmo.ediel

                log.debug('COSMO E_diel = %.15g', cosmo.ediel)

                if (log.verbose >= lib.logger.INFO and
                    hasattr(self.fcisolver, 'spin_square')):
                    ss = self.fcisolver.spin_square(fcivec, self.ncas, self.nelecas)
                    if isinstance(e_cas, (float, numpy.number)):
                        log.info('COSMO_cycle %d CASCI E = %.15g  E(CI) = %.15g  S^2 = %.7f',
                                 cycle, e_tot, e_cas, ss[0])
                    else:
                        for i, e in enumerate(e_cas):
                            log.info('COSMO_cycle %d CASCI root %d  E = %.15g  '
                                     'E(CI) = %.15g  S^2 = %.7f',
                                     cycle, i, e_tot[i], e, ss[0][i])
                else:
                    if isinstance(e_cas, (float, numpy.number)):
                        log.note('COSMO_cycle %d CASCI E = %.15g  E(CI) = %.15g',
                                 cycle, e_tot, e_cas)
                    else:
                        for i, e in enumerate(e_cas):
                            log.note('COSMO_cycle %d CASCI root %d  E = %.15g  E(CI) = %.15g',
                                     cycle, i, e_tot[i], e)
                return e_tot, e_cas, fcivec

            if cosmo.dm is not None:
                self.e_tot, self.e_cas, self.ci = casci_iter(mo_coeff, ci0, 0)
            else:
                e_tot = 0
                for icycle in range(cosmo.casci_max_cycle):
                    self.e_tot, self.e_cas, self.ci = casci_iter(mo_coeff, ci0,
                                                                 icycle)
                    if numpy.all(abs(self.e_tot-e_tot) < cosmo.casci_conv_tol):
                        log.debug('    delta E(CAS) = %s', self.e_tot-e_tot)
                        break
                    e_tot, ci0 = self.e_tot, self.ci

            self._finalize()
            return self.e_tot, self.e_cas, self.ci, self.mo_coeff, self.mo_energy

        def _finalize(self):
            if isinstance(self.e_tot, (float, numpy.number)):
                cosmo.e_tot = self.e_tot
                self.e_tot = cosmo.cosmo_occ(self.make_rdm1())
            else:
                cosmo.e_tot = self.e_tot[cosmo.casci_state_id]
                c = self.ci[cosmo.casci_state_id]
# Assuming the COSMO correction e_cosmo are the same for multiple roots.
                e_cosmo = (cosmo.cosmo_occ(self.make_rdm1(self.mo_coeff, c)) -
                           self.e_tot[cosmo.casci_state_id])
                self.e_tot += e_cosmo

            if self.canonicalization:
                log = lib.logger.Logger(self.stdout, self.verbose)
                if isinstance(self.e_cas, (float, numpy.number)):
                    self.canonicalize_(self.mo_coeff, self.ci,
                                       cas_natorb=self.natorb, verbose=log)
                else:
                    self.canonicalize_(self.mo_coeff, self.ci[0],
                                       cas_natorb=self.natorb, verbose=log)
            return oldCAS._finalize(self)

    return CAS()

Example 8

Project: deap Source File: cma.py
Function: update
    def update(self, population):
        """Update the current covariance matrix strategy from the
        *population*.

        :param population: A list of individuals from which to update the
                           parameters.
        """
        population.sort(key=lambda ind: ind.fitness, reverse=True)

        old_centroid = self.centroid
        self.centroid = numpy.dot(self.weights, population[0:self.mu])

        c_diff = self.centroid - old_centroid

        # Cumulation : update evolution path
        self.ps = (1 - self.cs) * self.ps \
            + sqrt(self.cs * (2 - self.cs) * self.mueff) / self.sigma \
            * numpy.dot(self.B, (1. / self.diagD)
                        * numpy.dot(self.B.T, c_diff))

        hsig = float((numpy.linalg.norm(self.ps) /
                      sqrt(1. - (1. - self.cs) ** (2. * (self.update_count + 1.))) / self.chiN
                      < (1.4 + 2. / (self.dim + 1.))))

        self.update_count += 1

        self.pc = (1 - self.cc) * self.pc + hsig \
            * sqrt(self.cc * (2 - self.cc) * self.mueff) / self.sigma \
            * c_diff

        # Update covariance matrix
        artmp = population[0:self.mu] - old_centroid
        self.C = (1 - self.ccov1 - self.ccovmu + (1 - hsig)
                  * self.ccov1 * self.cc * (2 - self.cc)) * self.C \
            + self.ccov1 * numpy.outer(self.pc, self.pc) \
            + self.ccovmu * numpy.dot((self.weights * artmp.T), artmp) \
            / self.sigma ** 2

        self.sigma *= numpy.exp((numpy.linalg.norm(self.ps) / self.chiN - 1.)
                                * self.cs / self.damps)

        self.diagD, self.B = numpy.linalg.eigh(self.C)
        indx = numpy.argsort(self.diagD)

        self.cond = self.diagD[indx[-1]] / self.diagD[indx[0]]

        self.diagD = self.diagD[indx] ** 0.5
        self.B = self.B[:, indx]
        self.BD = self.B * self.diagD

Example 9

Project: pyscf Source File: davidson.py
def dgeev(abop, x0, precond, type=1, tol=1e-14, max_cycle=50, max_space=12,
          lindep=1e-14, max_memory=2000, dot=numpy.dot, callback=None,
          nroots=1, lessio=False, verbose=logger.WARN):
    '''Davidson diagonalization method to solve  A c = e B c.

    Args:
        abop : function([x]) => ([array_like_x], [array_like_x])
            abop applies two matrix vector multiplications and returns tuple (Ax, Bx)
        x0 : 1D array
            Initial guess
        precond : function(dx, e, x0) => array_like_dx
            Preconditioner to generate new trial vector.
            The argument dx is a residual vector ``a*x0-e*x0``; e is the current
            eigenvalue; x0 is the current eigenvector.

    Kwargs:
        tol : float
            Convergence tolerance.
        max_cycle : int
            max number of iterations.
        max_space : int
            space size to hold trial vectors.
        lindep : float
            Linear dependency threshold.  The function is terminated when the
            smallest eigenvalue of the metric of the trial vectors is lower
            than this threshold.
        max_memory : int or float
            Allowed memory in MB.
        dot : function(x, y) => scalar
            Inner product
        callback : function(envs_dict) => None
            callback function takes one dict as the argument which is
            generated by the builtin function :func:`locals`, so that the
            callback function can access all local variables in the current
            envrionment.
        nroots : int
            Number of eigenvalues to be computed.  When nroots > 1, it affects
            the shape of the return value
        lessio : bool
            How to compute a*x0 for current eigenvector x0.  There are two
            ways to compute a*x0.  One is to assemble the existed a*x.  The
            other is to call aop(x0).  The default is the first method which
            needs more IO and less computational cost.  When IO is slow, the
            second method can be considered.

    Returns:
        e : list of floats
            The lowest :attr:`nroots` eigenvalues.
        c : list of 1D arrays
            The lowest :attr:`nroots` eigenvectors.
    '''
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(sys.stdout, verbose)

    def qr(xs):
        #q, r = numpy.linalg.qr(numpy.asarray(xs).T)
        #q = [qi/numpy_helper.norm(qi)
        #     for i, qi in enumerate(q.T) if r[i,i] > 1e-7]
        qs = [xs[0]/numpy_helper.norm(xs[0])]
        for i in range(1, len(xs)):
            xi = xs[i].copy()
            for j in range(len(qs)):
                xi -= qs[j] * numpy.dot(qs[j], xi)
            norm = numpy_helper.norm(xi)
            if norm > 1e-7:
                qs.append(xi/norm)
        return qs

    toloose = numpy.sqrt(tol) * 1e-2

    if isinstance(x0, numpy.ndarray) and x0.ndim == 1:
        x0 = [x0]
    max_cycle = min(max_cycle, x0[0].size)
    max_space = max_space + nroots * 2
    # max_space*3 for holding ax, bx and xs, nroots*3 for holding axt, bxt and xt
    _incore = max_memory*1e6/x0[0].nbytes > max_space*3+nroots*3
    heff = numpy.empty((max_space,max_space), dtype=x0[0].dtype)
    seff = numpy.empty((max_space,max_space), dtype=x0[0].dtype)
    fresh_start = True

    for icyc in range(max_cycle):
        if fresh_start:
            if _incore:
                xs = []
                ax = []
                bx = []
            else:
                xs = linalg_helper._Xlist()
                ax = linalg_helper._Xlist()
                bx = linalg_helper._Xlist()
            space = 0
# Orthogonalize xt space because the basis of subspace xs must be orthogonal
# but the eigenvectors x0 are very likely non-orthogonal when A is non-Hermitian.
            xt, x0 = qr(x0), None
            e = numpy.zeros(nroots)
            fresh_start = False
        elif len(xt) > 1:
            xt = qr(xt)
            xt = xt[:40]  # 40 trial vectors at most

        axt, bxt = abop(xt)
        if type > 1:
            axt = abop(bxt)[0]
        for k, xi in enumerate(xt):
            xs.append(xt[k])
            ax.append(axt[k])
            bx.append(bxt[k])
        rnow = len(xt)
        head, space = space, space+rnow

        if type == 1:
            for i in range(space):
                if head <= i < head+rnow:
                    for k in range(i-head+1):
                        heff[head+k,i] = dot(xt[k].conj(), axt[i-head])
                        heff[i,head+k] = heff[head+k,i].conj()
                        seff[head+k,i] = dot(xt[k].conj(), bxt[i-head])
                        seff[i,head+k] = seff[head+k,i].conj()
                else:
                    for k in range(rnow):
                        heff[head+k,i] = dot(xt[k].conj(), ax[i])
                        heff[i,head+k] = heff[head+k,i].conj()
                        seff[head+k,i] = dot(xt[k].conj(), bx[i])
                        seff[i,head+k] = seff[head+k,i].conj()
        else:
            for i in range(space):
                if head <= i < head+rnow:
                    for k in range(i-head+1):
                        heff[head+k,i] = dot(bxt[k].conj(), axt[i-head])
                        heff[i,head+k] = heff[head+k,i].conj()
                        seff[head+k,i] = dot(xt[k].conj(), bxt[i-head])
                        seff[i,head+k] = seff[head+k,i].conj()
                else:
                    for k in range(rnow):
                        heff[head+k,i] = dot(bxt[k].conj(), ax[i])
                        heff[i,head+k] = heff[head+k,i].conj()
                        seff[head+k,i] = dot(xt[k].conj(), bx[i])
                        seff[i,head+k] = seff[head+k,i].conj()

        w, v = scipy.linalg.eigh(heff[:space,:space], seff[:space,:space])
        if space < nroots or e.size != nroots:
            de = w[:nroots]
        else:
            de = w[:nroots] - e
        e = w[:nroots]

        x0 = []
        ax0 = []
        bx0 = []
        if lessio and not _incore:
            for k, ek in enumerate(e):
                x0.append(xs[space-1] * v[space-1,k])
            for i in reversed(range(space-1)):
                xsi = xs[i]
                for k, ek in enumerate(e):
                    x0[k] += v[i,k] * xsi
            ax0, bx0 = abop(x0)
            if type > 1:
                ax0 = abop(bx0)[0]
        else:
            for k, ek in enumerate(e):
                x0 .append(xs[space-1] * v[space-1,k])
                ax0.append(ax[space-1] * v[space-1,k])
                bx0.append(bx[space-1] * v[space-1,k])
            for i in reversed(range(space-1)):
                xsi = xs[i]
                axi = ax[i]
                bxi = bx[i]
                for k, ek in enumerate(e):
                    x0 [k] += v[i,k] * xsi
                    ax0[k] += v[i,k] * axi
                    bx0[k] += v[i,k] * bxi

        ide = numpy.argmax(abs(de))
        if abs(de[ide]) < tol:
            log.debug('converge %d %d  e= %s  max|de|= %4.3g',
                      icyc, space, e, de[ide])
            break

        dx_norm = []
        xt = []
        for k, ek in enumerate(e):
            if type == 1:
                dxtmp = ax0[k] - ek * bx0[k]
            else:
                dxtmp = ax0[k] - ek * x0[k]
            xt.append(dxtmp)
            dx_norm.append(numpy_helper.norm(dxtmp))

        if max(dx_norm) < toloose:
            log.debug('converge %d %d  |r|= %4.3g  e= %s  max|de|= %4.3g',
                      icyc, space, max(dx_norm), e, de[ide])
            break

        # remove subspace linear dependency
        for k, ek in enumerate(e):
            if dx_norm[k] > toloose:
                xt[k] = precond(xt[k], e[0], x0[k])
                xt[k] *= 1/numpy_helper.norm(xt[k])
            else:
                xt[k] = None
        xt = [xi for xi in xt if xi is not None]
        for i in range(space):
            for xi in xt:
                xsi = xs[i]
                xi -= xsi * numpy.dot(xi, xsi)
        norm_min = 1
        for i,xi in enumerate(xt):
            norm = numpy_helper.norm(xi)
            if norm > toloose:
                xt[i] *= 1/norm
                norm_min = min(norm_min, norm)
            else:
                xt[i] = None
        xt = [xi for xi in xt if xi is not None]
        if len(xt) == 0:
            log.debug('Linear dependency in trial subspace')
            break
        log.debug('davidson %d %d  |r|= %4.3g  e= %s  max|de|= %4.3g  lindep= %4.3g',
                  icyc, space, max(dx_norm), e, de[ide], norm)

        fresh_start = fresh_start or (space+len(xt) > max_space)

        if callable(callback):
            callback(locals())

    if type == 3:
        for k in range(nroots):
            x0[k] = abop(x0[k])[1]

    if nroots == 1:
        return e[0], x0[0]
    else:
        return e, x0

Example 10

Project: pyscf Source File: davidson.py
def eig(aop, x0, precond, tol=1e-14, max_cycle=50, max_space=12,
        lindep=1e-14, max_memory=2000, dot=numpy.dot, callback=None,
        nroots=1, lessio=False, left=False, pick=pickeig,
        verbose=logger.WARN):
    ''' A X = X w
    '''
    assert(not left)
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(sys.stdout, verbose)

    def qr(xs):
        #q, r = numpy.linalg.qr(numpy.asarray(xs).T)
        #q = [qi/numpy_helper.norm(qi)
        #     for i, qi in enumerate(q.T) if r[i,i] > 1e-7]
        qs = [xs[0]/numpy_helper.norm(xs[0])]
        for i in range(1, len(xs)):
            xi = xs[i].copy()
            for j in range(len(qs)):
                xi -= qs[j] * numpy.dot(qs[j], xi)
            norm = numpy_helper.norm(xi)
            if norm > 1e-7:
                qs.append(xi/norm)
        return qs

    toloose = numpy.sqrt(tol) * 1e-2

    if isinstance(x0, numpy.ndarray) and x0.ndim == 1:
        x0 = [x0]
    max_cycle = min(max_cycle, x0[0].size)
    max_space = max_space + nroots * 2
    # max_space*2 for holding ax and xs, nroots*2 for holding axt and xt
    _incore = max_memory*1e6/x0[0].nbytes > max_space*2+nroots*2
    heff = numpy.empty((max_space,max_space), dtype=x0[0].dtype)
    fresh_start = True

    for icyc in range(max_cycle):
        if fresh_start:
            if _incore:
                xs = []
                ax = []
            else:
                xs = linalg_helper._Xlist()
                ax = linalg_helper._Xlist()
            space = 0
# Orthogonalize xt space because the basis of subspace xs must be orthogonal
# but the eigenvectors x0 are very likely non-orthogonal when A is non-Hermitian.
            xt, x0 = qr(x0), None
            e = numpy.zeros(nroots)
            fresh_start = False
        elif len(xt) > 1:
            xt = qr(xt)
            xt = xt[:40]  # 40 trial vectors at most

        axt = aop(xt)
        for k, xi in enumerate(xt):
            xs.append(xt[k])
            ax.append(axt[k])
        rnow = len(xt)
        head, space = space, space+rnow

        for i in range(rnow):
            for k in range(rnow):
                heff[head+k,head+i] = dot(xt[k].conj(), axt[i])
        for i in range(head):
            axi = ax[i]
            xi = xs[i]
            for k in range(rnow):
                heff[head+k,i] = dot(xt[k].conj(), axi)
                heff[i,head+k] = dot(xi.conj(), axt[k])

        w, v = scipy.linalg.eig(heff[:space,:space])
        idx = pick(w, v, nroots)

        if idx.size < nroots or e.size != nroots:
            de = w[idx].real
        else:
            de = w[idx].real - e
        e = w[idx].real
        v = v[:,idx].real

        x0 = []
        ax0 = []
        if lessio and not _incore:
            for k, ek in enumerate(e):
                x0.append(xs[space-1] * v[space-1,k])
            for i in reversed(range(space-1)):
                xsi = xs[i]
                for k, ek in enumerate(e):
                    x0[k] += v[i,k] * xsi
            ax0 = aop(x0)
        else:
            for k, ek in enumerate(e):
                x0 .append(xs[space-1] * v[space-1,k])
                ax0.append(ax[space-1] * v[space-1,k])
            for i in reversed(range(space-1)):
                xsi = xs[i]
                axi = ax[i]
                for k, ek in enumerate(e):
                    x0 [k] += v[i,k] * xsi
                    ax0[k] += v[i,k] * axi

        ide = numpy.argmax(abs(de))
        if abs(de[ide]) < tol:
            log.debug('converge %d %d  e= %s  max|de|= %4.3g',
                      icyc, space, e, de[ide])
            break

        dx_norm = []
        xt = []
        for k, ek in enumerate(e):
            dxtmp = ax0[k] - ek * x0[k]
            xt.append(dxtmp)
            dx_norm.append(numpy_helper.norm(dxtmp))

        if max(dx_norm) < toloose:
            log.debug('converge %d %d  |r|= %4.3g  e= %s  max|de|= %4.3g',
                      icyc, space, max(dx_norm), e, de[ide])
            break

        # remove subspace linear dependency
        for k, ek in enumerate(e):
            if dx_norm[k] > toloose:
                xt[k] = precond(xt[k], e[0], x0[k])
                xt[k] *= 1/numpy_helper.norm(xt[k])
            else:
                xt[k] = None
        xt = [xi for xi in xt if xi is not None]
        for i in range(space):
            xsi = xs[i]
            for xi in xt:
                xi -= xsi * numpy.dot(xi, xsi)
        norm_min = 1
        for i,xi in enumerate(xt):
            norm = numpy_helper.norm(xi)
            if norm > toloose:
                xt[i] *= 1/norm
                norm_min = min(norm_min, norm)
            else:
                xt[i] = None
        xt = [xi for xi in xt if xi is not None]
        if len(xt) == 0:
            log.debug('Linear dependency in trial subspace')
            break
        log.debug('davidson %d %d  |r|= %4.3g  e= %s  max|de|= %4.3g  lindep= %4.3g',
                  icyc, space, max(dx_norm), e, de[ide], norm_min)

        fresh_start = fresh_start or (space+len(xt) > max_space)

        if callable(callback):
            callback(locals())

    if nroots == 1:
        return e[0], x0[0]
    else:
        return e, x0

Example 11

Project: topical_word_embeddings Source File: lsimodel.py
    def merge(self, other, decay=1.0):
        """
        Merge this Projection with another.

        The content of `other` is destroyed in the process, so pass this function a
        copy of `other` if you need it further.
        """
        if other.u is None:
            # the other projection is empty => do nothing
            return
        if self.u is None:
            # we are empty => result of merge is the other projection, whatever it is
            self.u = other.u.copy()
            self.s = other.s.copy()
            return
        if self.m != other.m:
            raise ValueError("vector space mismatch: update is using %s features, expected %s" %
                             (other.m, self.m))
        logger.info("merging projections: %s + %s" % (str(self.u.shape), str(other.u.shape)))
        m, n1, n2 = self.u.shape[0], self.u.shape[1], other.u.shape[1]
        # TODO Maybe keep the bases as elementary reflectors, without
        # forming explicit matrices with ORGQR.
        # The only operation we ever need is basis^T*basis ond basis*component.
        # But how to do that in scipy? And is it fast(er)?

        # find component of u2 orthogonal to u1
        logger.debug("constructing orthogonal component")
        self.u = asfarray(self.u, 'self.u')
        c = numpy.dot(self.u.T, other.u)
        self.u = ascarray(self.u, 'self.u')
        other.u -= numpy.dot(self.u, c)

        other.u = [other.u] # do some reference magic and call qr_destroy, to save RAM
        q, r = matutils.qr_destroy(other.u) # q, r = QR(component)
        assert not other.u

        # find the rotation that diagonalizes r
        k = numpy.bmat([[numpy.diag(decay * self.s), numpy.multiply(c, other.s)],
                        [matutils.pad(numpy.array([]).reshape(0, 0), min(m, n2), n1), numpy.multiply(r, other.s)]])
        logger.debug("computing SVD of %s dense matrix" % str(k.shape))
        try:
            # in numpy < 1.1.0, running SVD sometimes results in "LinAlgError: SVD did not converge'.
            # for these early versions of numpy, catch the error and try to compute
            # SVD again, but over k*k^T.
            # see http://www.mail-archive.com/[email protected]/msg07224.html and
            # bug ticket http://projects.scipy.org/numpy/ticket/706
            # sdoering: replaced numpy's linalg.svd with scipy's linalg.svd:
            u_k, s_k, _ = scipy.linalg.svd(k, full_matrices=False) # TODO *ugly overkill*!! only need first self.k SVD factors... but there is no LAPACK wrapper for partial svd/eigendecomp in numpy :( //sdoering: maybe there is one in scipy?
        except scipy.linalg.LinAlgError:
            logger.error("SVD(A) failed; trying SVD(A * A^T)")
            u_k, s_k, _ = scipy.linalg.svd(numpy.dot(k, k.T), full_matrices=False) # if this fails too, give up with an exception
            s_k = numpy.sqrt(s_k) # go back from eigen values to singular values

        k = clip_spectrum(s_k**2, self.k)
        u1_k, u2_k, s_k = numpy.array(u_k[:n1, :k]), numpy.array(u_k[n1:, :k]), s_k[:k]

        # update & rotate current basis U = [U, U']*[U1_k, U2_k]
        logger.debug("updating orthonormal basis U")
        self.s = s_k
        self.u = ascarray(self.u, 'self.u')
        self.u = numpy.dot(self.u, u1_k)

        q = ascarray(q, 'q')
        q = numpy.dot(q, u2_k)
        self.u += q

        # make each column of U start with a non-negative number (to force canonical decomposition)
        if self.u.shape[0] > 0:
            for i in xrange(self.u.shape[1]):
                if self.u[0, i] < 0.0:
                    self.u[:, i] *= -1.0

Example 12

Project: STAMP Source File: PCA.py
def _symm_eig(a):
	"""Return the eigenvectors and eigenvalues of the symmetric matrix a'a. If
	a has more columns than rows, then that matrix will be rank-deficient,
	and the non-zero eigenvalues and eigenvectors can be more easily extracted
	from the matrix aa'.
	From the properties of the SVD:
		if a of shape (m,n) has SVD u*s*v', then:
			a'a = v*s's*v'
			aa' = u*ss'*u'
		let s_hat, an array of shape (m,n), be such that s * s_hat = I(m,m) 
		and s_hat * s = I(n,n). (Note that s_hat is just the elementwise 
		reciprocal of s, as s is zero except on the main diagonal.)
		
		Thus, we can solve for u or v in terms of the other:
			v = a'*u*s_hat'
			u = a*v*s_hat			
	"""
	m, n = a.shape
	if m >= n:
		# just return the eigenvalues and eigenvectors of a'a
		vecs, vals = _eigh(numpy.dot(a.transpose(), a))
		vecs = numpy.where(vecs < 0, 0, vecs)
		return vecs, vals
	else:
		# figure out the eigenvalues and vectors based on aa', which is smaller
		sst_diag, u = _eigh(numpy.dot(a, a.transpose()))
		# in case due to numerical instabilities we have sst_diag < 0 anywhere, 
		# peg them to zero
		sst_diag = numpy.where(sst_diag < 0, 0, sst_diag)
		# now get the inverse square root of the diagonal, which will form the
		# main diagonal of s_hat
		err = numpy.seterr(divide='ignore', invalid='ignore')
		s_hat_diag = 1/numpy.sqrt(sst_diag)
		numpy.seterr(**err)
		s_hat_diag[~numpy.isfinite(s_hat_diag)] = 0
		# s_hat_diag is a list of length m, a'u is (n,m), so we can just use
		# numpy's broadcasting instead of matrix multiplication, and only create
		# the upper mxm block of a'u, since that's all we'll use anyway...
		v = numpy.dot(a.transpose(), u[:,:m]) * s_hat_diag
		return sst_diag, v

Example 13

Project: krypy Source File: test_utils.py
def run_ritz(A, v, maxiter, ip_B, Aevals, An, with_V, hermitian, type):
    ortho = 'house' if ip_B is None else 'dmgs'
    V, H = krypy.utils.arnoldi(A, v, maxiter=maxiter, ortho=ortho, ip_B=ip_B)
    N = v.shape[0]
    n = H.shape[1]
    A = krypy.utils.get_linearoperator((N, N), A)

    Z = None
    if with_V:
        theta, U, resnorm, Z = krypy.utils.ritz(H, V=V, hermitian=hermitian,
                                                type=type)
    else:
        theta, U, resnorm = krypy.utils.ritz(H, hermitian=hermitian, type=type)
    # check Z
    if Z is not None:
        assert(numpy.linalg.norm(numpy.dot(V[:, :n], U) - Z, 2) <= 1e-14)
    else:
        Z = numpy.dot(V[:, :n], U)

    # check shapes
    assert(theta.shape == (n,))
    assert(U.shape == (n, n))
    assert(resnorm.shape == (n,))
    assert(Z.shape == (N, n))
    # check norm of Ritz vectors
    for i in range(n):
        assert(numpy.abs(numpy.linalg.norm(U[:, i], 2) - 1) <= 1e-14)
    # check residuals
    R = A*Z - numpy.dot(Z, numpy.diag(theta))
    for i in range(n):
        rnorm = krypy.utils.norm(R[:, [i]], ip_B=ip_B)
        assert(numpy.abs(rnorm - resnorm[i]) <= 1e-14*An)
    # check Ritz projection property
    if type == 'ritz':
        assert(numpy.linalg.norm(krypy.utils.inner(V[:, :n], R, ip_B=ip_B), 2)
               <= 1e-14*An)
    elif type == 'harmonic':
        AVortho = scipy.linalg.orth(A*V[:, :n])
        assert(numpy.linalg.norm(krypy.utils.inner(AVortho, R, ip_B=ip_B), 2)
               <= 1e-12*An)
    else:
        pass

    # compare Ritz values with eigenvalues if n==N
    if n == N:
        Aevals_sort = numpy.argsort(numpy.abs(Aevals))
        theta_sort = numpy.argsort(numpy.abs(theta))
        assert((numpy.abs(Aevals[Aevals_sort] - theta[theta_sort])
                <= 5e-14*An).all())

Example 14

Project: gensim Source File: lsimodel.py
    def merge(self, other, decay=1.0):
        """
        Merge this Projection with another.

        The content of `other` is destroyed in the process, so pass this function a
        copy of `other` if you need it further.
        """
        if other.u is None:
            # the other projection is empty => do nothing
            return
        if self.u is None:
            # we are empty => result of merge is the other projection, whatever it is
            self.u = other.u.copy()
            self.s = other.s.copy()
            return
        if self.m != other.m:
            raise ValueError("vector space mismatch: update is using %s features, expected %s" %
                             (other.m, self.m))
        logger.info("merging projections: %s + %s", str(self.u.shape), str(other.u.shape))
        m, n1, n2 = self.u.shape[0], self.u.shape[1], other.u.shape[1]
        # TODO Maybe keep the bases as elementary reflectors, without
        # forming explicit matrices with ORGQR.
        # The only operation we ever need is basis^T*basis ond basis*component.
        # But how to do that in scipy? And is it fast(er)?

        # find component of u2 orthogonal to u1
        logger.debug("constructing orthogonal component")
        self.u = asfarray(self.u, 'self.u')
        c = numpy.dot(self.u.T, other.u)
        self.u = ascarray(self.u, 'self.u')
        other.u -= numpy.dot(self.u, c)

        other.u = [other.u]  # do some reference magic and call qr_destroy, to save RAM
        q, r = matutils.qr_destroy(other.u)  # q, r = QR(component)
        assert not other.u

        # find the rotation that diagonalizes r
        k = numpy.bmat([[numpy.diag(decay * self.s), numpy.multiply(c, other.s)],
                        [matutils.pad(numpy.array([]).reshape(0, 0), min(m, n2), n1), numpy.multiply(r, other.s)]])
        logger.debug("computing SVD of %s dense matrix", k.shape)
        try:
            # in numpy < 1.1.0, running SVD sometimes results in "LinAlgError: SVD did not converge'.
            # for these early versions of numpy, catch the error and try to compute
            # SVD again, but over k*k^T.
            # see http://www.mail-archive.com/[email protected]/msg07224.html and
            # bug ticket http://projects.scipy.org/numpy/ticket/706
            # sdoering: replaced numpy's linalg.svd with scipy's linalg.svd:
            u_k, s_k, _ = scipy.linalg.svd(k, full_matrices=False)  # TODO *ugly overkill*!! only need first self.k SVD factors... but there is no LAPACK wrapper for partial svd/eigendecomp in numpy :( //sdoering: maybe there is one in scipy?
        except scipy.linalg.LinAlgError:
            logger.error("SVD(A) failed; trying SVD(A * A^T)")
            u_k, s_k, _ = scipy.linalg.svd(numpy.dot(k, k.T), full_matrices=False)  # if this fails too, give up with an exception
            s_k = numpy.sqrt(s_k)  # go back from eigen values to singular values

        k = clip_spectrum(s_k**2, self.k)
        u1_k, u2_k, s_k = numpy.array(u_k[:n1, :k]), numpy.array(u_k[n1:, :k]), s_k[:k]

        # update & rotate current basis U = [U, U']*[U1_k, U2_k]
        logger.debug("updating orthonormal basis U")
        self.s = s_k
        self.u = ascarray(self.u, 'self.u')
        self.u = numpy.dot(self.u, u1_k)

        q = ascarray(q, 'q')
        q = numpy.dot(q, u2_k)
        self.u += q

        # make each column of U start with a non-negative number (to force canonical decomposition)
        if self.u.shape[0] > 0:
            for i in xrange(self.u.shape[1]):
                if self.u[0, i] < 0.0:
                    self.u[:, i] *= -1.0

Example 15

Project: scikit-learn Source File: dict_learning.py
Function: sparse_encode
def sparse_encode(X, dictionary, gram=None, cov=None, algorithm='lasso_lars',
                  n_nonzero_coefs=None, alpha=None, copy_cov=True, init=None,
                  max_iter=1000, n_jobs=1, check_input=True, verbose=0):
    """Sparse coding

    Each row of the result is the solution to a sparse coding problem.
    The goal is to find a sparse array `code` such that::

        X ~= code * dictionary

    Read more in the :ref:`User Guide <SparseCoder>`.

    Parameters
    ----------
    X: array of shape (n_samples, n_features)
        Data matrix

    dictionary: array of shape (n_components, n_features)
        The dictionary matrix against which to solve the sparse coding of
        the data. Some of the algorithms assume normalized rows for meaningful
        output.

    gram: array, shape=(n_components, n_components)
        Precomputed Gram matrix, dictionary * dictionary'

    cov: array, shape=(n_components, n_samples)
        Precomputed covariance, dictionary' * X

    algorithm: {'lasso_lars', 'lasso_cd', 'lars', 'omp', 'threshold'}
        lars: uses the least angle regression method (linear_model.lars_path)
        lasso_lars: uses Lars to compute the Lasso solution
        lasso_cd: uses the coordinate descent method to compute the
        Lasso solution (linear_model.Lasso). lasso_lars will be faster if
        the estimated components are sparse.
        omp: uses orthogonal matching pursuit to estimate the sparse solution
        threshold: squashes to zero all coefficients less than alpha from
        the projection dictionary * X'

    n_nonzero_coefs: int, 0.1 * n_features by default
        Number of nonzero coefficients to target in each column of the
        solution. This is only used by `algorithm='lars'` and `algorithm='omp'`
        and is overridden by `alpha` in the `omp` case.

    alpha: float, 1. by default
        If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the
        penalty applied to the L1 norm.
        If `algorithm='threshold'`, `alpha` is the absolute value of the
        threshold below which coefficients will be squashed to zero.
        If `algorithm='omp'`, `alpha` is the tolerance parameter: the value of
        the reconstruction error targeted. In this case, it overrides
        `n_nonzero_coefs`.

    init: array of shape (n_samples, n_components)
        Initialization value of the sparse codes. Only used if
        `algorithm='lasso_cd'`.

    max_iter: int, 1000 by default
        Maximum number of iterations to perform if `algorithm='lasso_cd'`.

    copy_cov: boolean, optional
        Whether to copy the precomputed covariance matrix; if False, it may be
        overwritten.

    n_jobs: int, optional
        Number of parallel jobs to run.

    check_input: boolean, optional
        If False, the input arrays X and dictionary will not be checked.

    verbose : int, optional
        Controls the verbosity; the higher, the more messages. Defaults to 0.

    Returns
    -------
    code: array of shape (n_samples, n_components)
        The sparse codes

    See also
    --------
    sklearn.linear_model.lars_path
    sklearn.linear_model.orthogonal_mp
    sklearn.linear_model.Lasso
    SparseCoder
    """
    if check_input:
        if algorithm == 'lasso_cd':
            dictionary = check_array(dictionary, order='C', dtype='float64')
            X = check_array(X, order='C', dtype='float64')
        else:
            dictionary = check_array(dictionary)
            X = check_array(X)

    n_samples, n_features = X.shape
    n_components = dictionary.shape[0]

    if gram is None and algorithm != 'threshold':
        gram = np.dot(dictionary, dictionary.T)

    if cov is None and algorithm != 'lasso_cd':
        copy_cov = False
        cov = np.dot(dictionary, X.T)

    if algorithm in ('lars', 'omp'):
        regularization = n_nonzero_coefs
        if regularization is None:
            regularization = min(max(n_features / 10, 1), n_components)
    else:
        regularization = alpha
        if regularization is None:
            regularization = 1.

    if n_jobs == 1 or algorithm == 'threshold':
        code = _sparse_encode(X,
                              dictionary, gram, cov=cov,
                              algorithm=algorithm,
                              regularization=regularization, copy_cov=copy_cov,
                              init=init,
                              max_iter=max_iter,
                              check_input=False,
                              verbose=verbose)
        # This ensure that dimensionality of code is always 2,
        # consistant with the case n_jobs > 1
        if code.ndim == 1:
            code = code[np.newaxis, :]
        return code

    # Enter parallel code block
    code = np.empty((n_samples, n_components))
    slices = list(gen_even_slices(n_samples, _get_n_jobs(n_jobs)))

    code_views = Parallel(n_jobs=n_jobs, verbose=verbose)(
        delayed(_sparse_encode)(
            X[this_slice], dictionary, gram,
            cov[:, this_slice] if cov is not None else None,
            algorithm,
            regularization=regularization, copy_cov=copy_cov,
            init=init[this_slice] if init is not None else None,
            max_iter=max_iter,
            check_input=False)
        for this_slice in slices)
    for this_slice, this_view in zip(slices, code_views):
        code[this_slice] = this_view
    return code

Example 16

Project: RMG-Py Source File: mbo.py
Function: calculate
    def calculate(self, indices=None, fupdate=0.05):
        """Calculate Mayer's bond orders."""
    
        retval = super(MBO, self).calculate(fupdate)
        if not retval: #making density didn't work
            return False

        # Do we have the needed info in the ccData object?
        if not (hasattr(self.data, "aooverlaps")
                or hasattr(self.data, "fooverlaps")):
            self.logger.error("Missing overlap matrix")
            return False #let the caller of function know we didn't finish

        if not indices:

            # Build list of groups of orbitals in each atom for atomresults.
            if hasattr(self.data, "aonames"):
                names = self.data.aonames
                overlaps = self.data.aooverlaps
            elif hasattr(self.data, "fonames"):
                names = self.data.fonames
                overlaps = self.data.fooverlaps
            else:
                self.logger.error("Missing aonames or fonames")
                return False

            atoms = []
            indices = []

            name = names[0].split('_')[0]
            atoms.append(name)
            indices.append([0])

            for i in range(1, len(names)):
                name = names[i].split('_')[0]
                try:
                    index = atoms.index(name)
                except ValueError: #not found in atom list
                    atoms.append(name)
                    indices.append([i])
                else:
                    indices[index].append(i)

        self.logger.info("Creating attribute fragresults: array[3]")
        size = len(indices)

        # Determine number of steps, and whether process involves beta orbitals.
        PS = []
        PS.append(numpy.dot(self.density[0], overlaps))
        nstep = size**2 #approximately quadratic in size
        unrestricted = (len(self.data.mocoeffs) == 2)
        if unrestricted:
            self.fragresults = numpy.zeros([2, size, size], "d")
            PS.append(numpy.dot(self.density[1], overlaps))
        else:
            self.fragresults = numpy.zeros([1, size, size], "d")

        # Intialize progress if available.
        if self.progress:
            self.progress.initialize(nstep)

        step = 0
        for i in range(len(indices)):

            if self.progress and random.random() < fupdate:
                self.progress.update(step, "Mayer's Bond Order")

            for j in range(i+1, len(indices)):

                tempsumA = 0
                tempsumB = 0
                
                for a in indices[i]:

                    for b in indices[j]:

                        tempsumA += 2 * PS[0][a][b] * PS[0][b][a]
                        if unrestricted:
                            tempsumB += 2 * PS[1][a][b] * PS[1][b][a]

                self.fragresults[0][i, j] = tempsumA
                self.fragresults[0][j, i] = tempsumA

                if unrestricted:
                    self.fragresults[1][i, j] = tempsumB
                    self.fragresults[1][j, i] = tempsumB

        if self.progress:
            self.progress.update(nstep, "Done")

        return True

Example 17

Project: scikit-learn Source File: omp.py
Function: orthogonal_mp
def orthogonal_mp(X, y, n_nonzero_coefs=None, tol=None, precompute=False,
                  copy_X=True, return_path=False,
                  return_n_iter=False):
    """Orthogonal Matching Pursuit (OMP)

    Solves n_targets Orthogonal Matching Pursuit problems.
    An instance of the problem has the form:

    When parametrized by the number of non-zero coefficients using
    `n_nonzero_coefs`:
    argmin ||y - X\gamma||^2 subject to ||\gamma||_0 <= n_{nonzero coefs}

    When parametrized by error using the parameter `tol`:
    argmin ||\gamma||_0 subject to ||y - X\gamma||^2 <= tol

    Read more in the :ref:`User Guide <omp>`.

    Parameters
    ----------
    X : array, shape (n_samples, n_features)
        Input data. Columns are assumed to have unit norm.

    y : array, shape (n_samples,) or (n_samples, n_targets)
        Input targets

    n_nonzero_coefs : int
        Desired number of non-zero entries in the solution. If None (by
        default) this value is set to 10% of n_features.

    tol : float
        Maximum norm of the residual. If not None, overrides n_nonzero_coefs.

    precompute : {True, False, 'auto'},
        Whether to perform precomputations. Improves performance when n_targets
        or n_samples is very large.

    copy_X : bool, optional
        Whether the design matrix X must be copied by the algorithm. A false
        value is only helpful if X is already Fortran-ordered, otherwise a
        copy is made anyway.

    return_path : bool, optional. Default: False
        Whether to return every value of the nonzero coefficients along the
        forward path. Useful for cross-validation.

    return_n_iter : bool, optional default False
        Whether or not to return the number of iterations.

    Returns
    -------
    coef : array, shape (n_features,) or (n_features, n_targets)
        Coefficients of the OMP solution. If `return_path=True`, this contains
        the whole coefficient path. In this case its shape is
        (n_features, n_features) or (n_features, n_targets, n_features) and
        iterating over the last axis yields coefficients in increasing order
        of active features.

    n_iters : array-like or int
        Number of active features across every target. Returned only if
        `return_n_iter` is set to True.

    See also
    --------
    OrthogonalMatchingPursuit
    orthogonal_mp_gram
    lars_path
    decomposition.sparse_encode

    Notes
    -----
    Orthogonal matching pursuit was introduced in S. Mallat, Z. Zhang,
    Matching pursuits with time-frequency dictionaries, IEEE Transactions on
    Signal Processing, Vol. 41, No. 12. (December 1993), pp. 3397-3415.
    (http://blanche.polytechnique.fr/~mallat/papiers/MallatPursuit93.pdf)

    This implementation is based on Rubinstein, R., Zibulevsky, M. and Elad,
    M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal
    Matching Pursuit Technical Report - CS Technion, April 2008.
    http://www.cs.technion.ac.il/~ronrubin/Publications/KSVD-OMP-v2.pdf

    """
    X = check_array(X, order='F', copy=copy_X)
    copy_X = False
    if y.ndim == 1:
        y = y.reshape(-1, 1)
    y = check_array(y)
    if y.shape[1] > 1:  # subsequent targets will be affected
        copy_X = True
    if n_nonzero_coefs is None and tol is None:
        # default for n_nonzero_coefs is 0.1 * n_features
        # but at least one.
        n_nonzero_coefs = max(int(0.1 * X.shape[1]), 1)
    if tol is not None and tol < 0:
        raise ValueError("Epsilon cannot be negative")
    if tol is None and n_nonzero_coefs <= 0:
        raise ValueError("The number of atoms must be positive")
    if tol is None and n_nonzero_coefs > X.shape[1]:
        raise ValueError("The number of atoms cannot be more than the number "
                         "of features")
    if precompute == 'auto':
        precompute = X.shape[0] > X.shape[1]
    if precompute:
        G = np.dot(X.T, X)
        G = np.asfortranarray(G)
        Xy = np.dot(X.T, y)
        if tol is not None:
            norms_squared = np.sum((y ** 2), axis=0)
        else:
            norms_squared = None
        return orthogonal_mp_gram(G, Xy, n_nonzero_coefs, tol, norms_squared,
                                  copy_Gram=copy_X, copy_Xy=False,
                                  return_path=return_path)

    if return_path:
        coef = np.zeros((X.shape[1], y.shape[1], X.shape[1]))
    else:
        coef = np.zeros((X.shape[1], y.shape[1]))
    n_iters = []

    for k in range(y.shape[1]):
        out = _cholesky_omp(
            X, y[:, k], n_nonzero_coefs, tol,
            copy_X=copy_X, return_path=return_path)
        if return_path:
            _, idx, coefs, n_iter = out
            coef = coef[:, :, :len(idx)]
            for n_active, x in enumerate(coefs.T):
                coef[idx[:n_active + 1], k, n_active] = x[:n_active + 1]
        else:
            x, idx, n_iter = out
            coef[idx, k] = x
        n_iters.append(n_iter)

    if y.shape[1] == 1:
        n_iters = n_iters[0]

    if return_n_iter:
        return np.squeeze(coef), n_iters
    else:
        return np.squeeze(coef)

Example 18

Project: scipy Source File: matfuncs.py
def funm(A, func, disp=True):
    """
    Evaluate a matrix function specified by a callable.

    Returns the value of matrix-valued function ``f`` at `A`. The
    function ``f`` is an extension of the scalar-valued function `func`
    to matrices.

    Parameters
    ----------
    A : (N, N) array_like
        Matrix at which to evaluate the function
    func : callable
        Callable object that evaluates a scalar function f.
        Must be vectorized (eg. using vectorize).
    disp : bool, optional
        Print warning if error in the result is estimated large
        instead of returning estimated error. (Default: True)

    Returns
    -------
    funm : (N, N) ndarray
        Value of the matrix function specified by func evaluated at `A`
    errest : float
        (if disp == False)

        1-norm of the estimated error, ||err||_1 / ||A||_1

    Examples
    --------
    >>> from scipy.linalg import funm
    >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
    >>> funm(a, lambda x: x*x)
    array([[  4.,  15.],
           [  5.,  19.]])
    >>> a.dot(a)
    array([[  4.,  15.],
           [  5.,  19.]])

    Notes
    -----
    This function implements the general algorithm based on Schur decomposition
    (Algorithm 9.1.1. in [1]_).

    If the input matrix is known to be diagonalizable, then relying on the
    eigendecomposition is likely to be faster. For example, if your matrix is
    Hermitian, you can do

    >>> from scipy.linalg import eigh
    >>> def funm_herm(a, func, check_finite=False):
    ...     w, v = eigh(a, check_finite=check_finite)
    ...     ## if you further know that your matrix is positive semidefinite,
    ...     ## you can optionally guard against precision errors by doing
    ...     # w = np.maximum(w, 0)
    ...     w = func(w)
    ...     return (v * w).dot(v.conj().T)

    References
    ----------
    .. [1] Gene H. Golub, Charles F. van Loan, Matrix Computations 4th ed.

    """
    A = _asarray_square(A)
    # Perform Shur decomposition (lapack ?gees)
    T, Z = schur(A)
    T, Z = rsf2csf(T,Z)
    n,n = T.shape
    F = diag(func(diag(T)))  # apply function to diagonal elements
    F = F.astype(T.dtype.char)  # e.g. when F is real but T is complex

    minden = abs(T[0,0])

    # implement Algorithm 11.1.1 from Golub and Van Loan
    #                 "matrix Computations."
    for p in range(1,n):
        for i in range(1,n-p+1):
            j = i + p
            s = T[i-1,j-1] * (F[j-1,j-1] - F[i-1,i-1])
            ksl = slice(i,j-1)
            val = dot(T[i-1,ksl],F[ksl,j-1]) - dot(F[i-1,ksl],T[ksl,j-1])
            s = s + val
            den = T[j-1,j-1] - T[i-1,i-1]
            if den != 0.0:
                s = s / den
            F[i-1,j-1] = s
            minden = min(minden,abs(den))

    F = dot(dot(Z, F), transpose(conjugate(Z)))
    F = _maybe_real(A, F)

    tol = {0:feps, 1:eps}[_array_precision[F.dtype.char]]
    if minden == 0.0:
        minden = tol
    err = min(1, max(tol,(tol/minden)*norm(triu(T,1),1)))
    if product(ravel(logical_not(isfinite(F))),axis=0):
        err = Inf
    if disp:
        if err > 1000*tol:
            print("funm result may be inaccurate, approximate err =", err)
        return F
    else:
        return F, err

Example 19

Project: molecular-design-toolkit Source File: coords.py
@toplevel
#@singledispatch
def dihedral(a1, a2=None, a3=None, a4=None):
    """Twist angle of bonds a1-a2 and a4-a3 around around the central bond a2-a3

    Can be called as ``dihedral(a1, a2, a3, a4)``
              OR     ``dihedral(a2, a2)``
              OR     ``dihedral(bond)``

    Args:
        a1 (mdt.Bond): the central bond in the dihedral. OR
        a1,a2 (mdt.Atom): the atoms describing the dihedral
        a3,a4 (mdt.Atom): (optional) if not passed, ``a1`` and ``a2`` will be treated as the
            central atoms in this bond, and a3 and a4 will be inferred.

    Returns:
        (units.Scalar[angle]): angle -  [0, 2 pi) radians
    """
    if a3 is a4 is None:  # infer the first and last atoms
        a1, a2, a3, a4 = _infer_dihedral(a1, a2)

    r21 = (a1.position - a2.position).defunits_value()  # remove units immediately to improve speed
    r34 = (a4.position - a3.position).defunits_value()

    dim = len(r21.shape)
    center_bond = (a2.position - a3.position).defunits_value()
    plane_normal = normalized(center_bond)

    if dim == 1:
        r21_proj = r21 - plane_normal * r21.dot(plane_normal)
        r34_proj = r34 - plane_normal * r34.dot(plane_normal)
    else:
        assert dim == 2
        r21_proj = r21 - plane_normal * (r21*plane_normal).sum(axis=1)[:, None]
        r34_proj = r34 - plane_normal * (r34*plane_normal).sum(axis=1)[:, None]

    va = normalized(r21_proj)
    vb = normalized(r34_proj)

    if dim == 1:
        costheta = np.dot(va, vb)
        if np.allclose(costheta, 1.0):
            return 0.0 * u.radians
        elif np.allclose(costheta, -1.0):
            return u.pi * u.radians
    else:
        costheta = (va*vb).sum(axis=1)

    abstheta = safe_arccos(costheta)
    cross = np.cross(va, vb)

    if dim == 1:
        theta = abstheta * np.sign(np.dot(cross, plane_normal))
    else:
        theta = abstheta * np.sign((cross*plane_normal).sum(axis=1))
    return (theta * u.radians) % (2.0 * u.pi * u.radians)

Example 20

Project: python-control Source File: statesp.py
def _rss_generate(states, inputs, outputs, type):
    """Generate a random state space.

    This does the actual random state space generation expected from rss and
    drss.  type is 'c' for continuous systems and 'd' for discrete systems.

    """

    # Probability of repeating a previous root.
    pRepeat = 0.05
    # Probability of choosing a real root.  Note that when choosing a complex
    # root, the conjugate gets chosen as well.  So the expected proportion of
    # real roots is pReal / (pReal + 2 * (1 - pReal)).
    pReal = 0.6
    # Probability that an element in B or C will not be masked out.
    pBCmask = 0.8
    # Probability that an element in D will not be masked out.
    pDmask = 0.3
    # Probability that D = 0.
    pDzero = 0.5

    # Check for valid input arguments.
    if states < 1 or states % 1:
        raise ValueError("states must be a positive integer.  states = %g." %
            states)
    if inputs < 1 or inputs % 1:
        raise ValueError("inputs must be a positive integer.  inputs = %g." %
            inputs)
    if outputs < 1 or outputs % 1:
        raise ValueError("outputs must be a positive integer.  outputs = %g." %
            outputs)

    # Make some poles for A.  Preallocate a complex array.
    poles = zeros(states) + zeros(states) * 0.j
    i = 0

    while i < states:
        if rand() < pRepeat and i != 0 and i != states - 1:
            # Small chance of copying poles, if we're not at the first or last
            # element.
            if poles[i-1].imag == 0:
                # Copy previous real pole.
                poles[i] = poles[i-1]
                i += 1
            else:
                # Copy previous complex conjugate pair of poles.
                poles[i:i+2] = poles[i-2:i]
                i += 2
        elif rand() < pReal or i == states - 1:
            # No-oscillation pole.
            if type == 'c':
                poles[i] = -exp(randn()) + 0.j
            elif type == 'd':
                poles[i] = 2. * rand() - 1.
            i += 1
        else:
            # Complex conjugate pair of oscillating poles.
            if type == 'c':
                poles[i] = complex(-exp(randn()), 3. * exp(randn()))
            elif type == 'd':
                mag = rand()
                phase = 2. * pi * rand()
                poles[i] = complex(mag * cos(phase),
                    mag * sin(phase))
            poles[i+1] = complex(poles[i].real, -poles[i].imag)
            i += 2

    # Now put the poles in A as real blocks on the diagonal.
    A = zeros((states, states))
    i = 0
    while i < states:
        if poles[i].imag == 0:
            A[i, i] = poles[i].real
            i += 1
        else:
            A[i, i] = A[i+1, i+1] = poles[i].real
            A[i, i+1] = poles[i].imag
            A[i+1, i] = -poles[i].imag
            i += 2
    # Finally, apply a transformation so that A is not block-diagonal.
    while True:
        T = randn(states, states)
        try:
            A = dot(solve(T, A), T) # A = T \ A * T
            break
        except LinAlgError:
            # In the unlikely event that T is rank-deficient, iterate again.
            pass

    # Make the remaining matrices.
    B = randn(states, inputs)
    C = randn(outputs, states)
    D = randn(outputs, inputs)

    # Make masks to zero out some of the elements.
    while True:
        Bmask = rand(states, inputs) < pBCmask
        if any(Bmask): # Retry if we get all zeros.
            break
    while True:
        Cmask = rand(outputs, states) < pBCmask
        if any(Cmask): # Retry if we get all zeros.
            break
    if rand() < pDzero:
        Dmask = zeros((outputs, inputs))
    else:
        Dmask = rand(outputs, inputs) < pDmask

    # Apply masks.
    B = B * Bmask
    C = C * Cmask
    D = D * Dmask

    return StateSpace(A, B, C, D)

Example 21

Project: filterpy Source File: kalman_filter.py
    def test_matrix_dimensions(self, z=None, H=None, R=None, F=None, Q=None):
        """ Performs a series of asserts to check that the size of everything
        is what it should be. This can help you debug problems in your design.

        If you pass in H, R, F, Q those will be used instead of this object's
        value for those matrices.

        Testing `z` (the measurement) is problamatic. x is a vector, and can be
        implemented as either a 1D array or as a nx1 column vector. Thus Hx
        can be of different shapes. Then, if Hx is a single value, it can
        be either a 1D array or 2D vector. If either is true, z can reasonably
        be a scalar (either '3' or np.array('3') are scalars under this
        definition), a 1D, 1 element array, or a 2D, 1 element array. You are
        allowed to pass in any combination that works.
        """

        if H is None: H = self.H
        if R is None: R = self.R
        if F is None: F = self._F
        if Q is None: Q = self._Q
        x = self._x
        P = self._P

        assert x.ndim == 1 or x.ndim == 2, \
                "x must have one or two dimensions, but has {}".format(
                x.ndim)

        if x.ndim == 1:
            assert x.shape[0] == self.dim_x, \
                   "Shape of x must be ({},{}), but is {}".format(
                   self.dim_x, 1, x.shape)
        else:
            assert x.shape == (self.dim_x, 1), \
                   "Shape of x must be ({},{}), but is {}".format(
                   self.dim_x, 1, x.shape)

        assert P.shape == (self.dim_x, self.dim_x), \
               "Shape of P must be ({},{}), but is {}".format(
               self.dim_x, self.dim_x, P.shape)

        assert Q.shape == (self.dim_x, self.dim_x), \
               "Shape of P must be ({},{}), but is {}".format(
               self.dim_x, self.dim_x, P.shape)

        assert F.shape == (self.dim_x, self.dim_x), \
               "Shape of F must be ({},{}), but is {}".format(
               self.dim_x, self.dim_x, F.shape)


        assert np.ndim(H) == 2, \
               "Shape of H must be (dim_z, {}), but is {}".format(
               P.shape[0], shape(H))

        assert H.shape[1] == P.shape[0], \
               "Shape of H must be (dim_z, {}), but is {}".format(
               P.shape[0], H.shape)

        # shape of R must be the same as HPH'
        hph_shape = (H.shape[0], H.shape[0])
        r_shape = shape(R)

        if H.shape[0] == 1:
            # r can be scalar, 1D, or 2D in this case
            assert r_shape == () or r_shape == (1,) or r_shape == (1,1), \
            "R must be scalar or one element array, but is shaped {}".format(
            r_shape)
        else:
            assert r_shape == hph_shape, \
            "shape of R should be {} but it is {}".format(hph_shape, r_shape)


        if z is not None:
            z_shape = shape(z)
        else:
            z_shape = (self.dim_z, 1)

        # H@x must have shape of z
        Hx = dot(H, x)

        if z_shape == (): # scalar or np.array(scalar)
            assert Hx.ndim == 1 or shape(Hx) == (1,1), \
            "shape of z should be {}, not {} for the given H".format(
                   shape(Hx), z_shape)

        elif shape(Hx) == (1,):
            assert z_shape[0] == 1, 'Shape of z must be {} for the given H'.format(shape(Hx))

        else:
            assert (z_shape == shape(Hx) or
                    (len(z_shape) == 1 and shape(Hx) == (z_shape[0], 1))), \
                    "shape of z should be {}, not {} for the given H".format(
                    shape(Hx), z_shape)

        if np.ndim(Hx) > 1 and shape(Hx) != (1,1):
            assert shape(Hx) == z_shape, \
               'shape of z should be {} for the given H, but it is {}'.format(
               shape(Hx), z_shape)

Example 22

Project: filterpy Source File: kalman_filter.py
def update(x, P, z, R, H=None, return_all=False):
    """
    Add a new measurement (z) to the Kalman filter. If z is None, nothing
    is changed.

    This can handle either the multidimensional or unidimensional case. If
    all parameters are floats instead of arrays the filter will still work,
    and return floats for x, P as the result.

    update(1, 2, 1, 1, 1)  # univariate
    update(x, P, 1



    Parameters
    ----------

    x : numpy.array(dim_x, 1), or float
        State estimate vector

    P : numpy.array(dim_x, dim_x), or float
        Covariance matrix

    z : numpy.array(dim_z, 1), or float
        measurement for this update.

    R : numpy.array(dim_z, dim_z), or float
        Measurement noise matrix

    H : numpy.array(dim_x, dim_x), or float, optional
        Measurement function. If not provided, a value of 1 is assumed.

    return_all : bool, default False
        If true, y, K, S, and log_likelihood are returned, otherwise
        only x and P are returned.

    Returns
    -------

    x : numpy.array
        Posterior state estimate vector

    P : numpy.array
        Posterior covariance matrix

    y : numpy.array or scalar
        Residua. Difference between measurement and state in measurement space

    K : numpy.array
        Kalman gain

    S : numpy.array
        System uncertainty in measurement space

    log_likelihood : float
        log likelihood of the measurement
    """


    if z is None:
        if return_all:
            return x, P, None, None, None, None
        else:
            return x, P

    if H is None:
        H = np.array([1])

    if np.isscalar(H):
        H = np.array([H])

    if not np.isscalar(x):
        # handle special case: if z is in form [[z]] but x is not a column
        # vector dimensions will not match
        if x.ndim==1 and shape(z) == (1,1):
            z = z[0]

        if shape(z) == (): # is it scalar, e.g. z=3 or z=np.array(3)
            z = np.asarray([z])

    # error (residual) between measurement and prediction
    y = z - dot(H, x)

    # project system uncertainty into measurement space
    S = dot3(H, P, H.T) + R


    # map system uncertainty into kalman gain
    try:
        K = dot3(P, H.T, linalg.inv(S))
    except:
        K = dot3(P, H.T, 1/S)


    # predict new x with residual scaled by the kalman gain
    x = x + dot(K, y)

    # P = (I-KH)P(I-KH)' + KRK'
    KH = dot(K, H)


    try:
        I_KH = np.eye(KH.shape[0]) - KH
    except:
        I_KH = np.array(1 - KH)
    P = dot3(I_KH, P, I_KH.T) + dot3(K, R, K.T)

    # compute log likelihood
    mean = np.array(dot(H, x)).flatten()
    flatz = np.asarray(z).flatten()
    log_likelihood = multivariate_normal.logpdf(flatz, mean, cov=S, allow_singular=True)

    if return_all:
        return x, P, y, K, S, log_likelihood
    else:
        return x, P

Example 23

Project: QuantEcon.py Source File: lqnash.py
def nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,
          beta=1.0, tol=1e-8, max_iter=1000):
    r"""
    Compute the limit of a Nash linear quadratic dynamic game. In this
    problem, player i minimizes

    .. math::
        \sum_{t=0}^{\infty}
        \left\{
            x_t' r_i x_t + 2 x_t' w_i
            u_{it} +u_{it}' q_i u_{it} + u_{jt}' s_i u_{jt} + 2 u_{jt}'
            m_i u_{it}
        \right\}

    subject to the law of motion

    .. math::
        x_{t+1} = A x_t + b_1 u_{1t} + b_2 u_{2t}

    and a perceived control law :math:`u_j(t) = - f_j x_t` for the other
    player.

    The solution computed in this routine is the :math:`f_i` and
    :math:`p_i` of the associated double optimal linear regulator
    problem.

    Parameters
    ----------
    A : scalar(float) or array_like(float)
        Corresponds to the above equation, should be of size (n, n)
    B1 : scalar(float) or array_like(float)
        As above, size (n, k_1)
    B2 : scalar(float) or array_like(float)
        As above, size (n, k_2)
    R1 : scalar(float) or array_like(float)
        As above, size (n, n)
    R2 : scalar(float) or array_like(float)
        As above, size (n, n)
    Q1 : scalar(float) or array_like(float)
        As above, size (k_1, k_1)
    Q2 : scalar(float) or array_like(float)
        As above, size (k_2, k_2)
    S1 : scalar(float) or array_like(float)
        As above, size (k_1, k_1)
    S2 : scalar(float) or array_like(float)
        As above, size (k_2, k_2)
    W1 : scalar(float) or array_like(float)
        As above, size (n, k_1)
    W2 : scalar(float) or array_like(float)
        As above, size (n, k_2)
    M1 : scalar(float) or array_like(float)
        As above, size (k_2, k_1)
    M2 : scalar(float) or array_like(float)
        As above, size (k_1, k_2)
    beta : scalar(float), optional(default=1.0)
        Discount rate
    tol : scalar(float), optional(default=1e-8)
        This is the tolerance level for convergence
    max_iter : scalar(int), optional(default=1000)
        This is the maximum number of iteratiosn allowed

    Returns
    -------
    F1 : array_like, dtype=float, shape=(k_1, n)
        Feedback law for agent 1
    F2 : array_like, dtype=float, shape=(k_2, n)
        Feedback law for agent 2
    P1 : array_like, dtype=float, shape=(n, n)
        The steady-state solution to the associated discrete matrix
        Riccati equation for agent 1
    P2 : array_like, dtype=float, shape=(n, n)
        The steady-state solution to the associated discrete matrix
        Riccati equation for agent 2

    """
    # == Unload parameters and make sure everything is an array == #
    params = A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2
    params = map(np.asarray, params)
    A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2 = params

    # == Multiply A, B1, B2 by sqrt(beta) to enforce discounting == #
    A, B1, B2 = [np.sqrt(beta) * x for x in (A, B1, B2)]

    n = A.shape[0]

    if B1.ndim == 1:
        k_1 = 1
        B1 = np.reshape(B1, (n, 1))
    else:
        k_1 = B1.shape[1]

    if B2.ndim == 1:
        k_2 = 1
        B2 = np.reshape(B2, (n, 1))
    else:
        k_2 = B2.shape[1]

    v1 = eye(k_1)
    v2 = eye(k_2)
    P1 = np.zeros((n, n))
    P2 = np.zeros((n, n))
    F1 = np.random.randn(k_1, n)
    F2 = np.random.randn(k_2, n)

    for it in range(max_iter):
        # update
        F10 = F1
        F20 = F2

        G2 = solve(dot(B2.T, P2.dot(B2))+Q2, v2)
        G1 = solve(dot(B1.T, P1.dot(B1))+Q1, v1)
        H2 = dot(G2, B2.T.dot(P2))
        H1 = dot(G1, B1.T.dot(P1))

        # break up the computation of F1, F2
        F1_left = v1 - dot(H1.dot(B2)+G1.dot(M1.T),
                           H2.dot(B1)+G2.dot(M2.T))
        F1_right = H1.dot(A)+G1.dot(W1.T) - dot(H1.dot(B2)+G1.dot(M1.T),
                                                H2.dot(A)+G2.dot(W2.T))
        F1 = solve(F1_left, F1_right)
        F2 = H2.dot(A)+G2.dot(W2.T) - dot(H2.dot(B1)+G2.dot(M2.T), F1)

        Lambda1 = A - B2.dot(F2)
        Lambda2 = A - B1.dot(F1)
        Pi1 = R1 + dot(F2.T, S1.dot(F2))
        Pi2 = R2 + dot(F1.T, S2.dot(F1))

        P1 = dot(Lambda1.T, P1.dot(Lambda1)) + Pi1 - \
             dot(dot(Lambda1.T, P1.dot(B1)) + W1 - F2.T.dot(M1), F1)
        P2 = dot(Lambda2.T, P2.dot(Lambda2)) + Pi2 - \
             dot(dot(Lambda2.T, P2.dot(B2)) + W2 - F1.T.dot(M2), F2)

        dd = np.max(np.abs(F10 - F1)) + np.max(np.abs(F20 - F2))

        if dd < tol:  # success!
            break

    else:
        msg = 'No convergence: Iteration limit of {0} reached in nnash'
        raise ValueError(msg.format(max_iter))

    return F1, F2, P1, P2

Example 24

Project: filterpy Source File: UKF.py
    def rts_smoother(self, Xs, Ps, Qs=None, dt=None):
        """ Runs the Rauch-Tung-Striebal Kalman smoother on a set of
        means and covariances computed by the UKF. The usual input
        would come from the output of `batch_filter()`.

        Parameters
        ----------

        Xs : numpy.array
           array of the means (state variable x) of the output of a Kalman
           filter.

        Ps : numpy.array
            array of the covariances of the output of a kalman filter.

        Qs: list-like collection of numpy.array, optional
            Process noise of the Kalman filter at each time step. Optional,
            if not provided the filter's self.Q will be used

        dt : optional, float or array-like of float
            If provided, specifies the time step of each step of the filter.
            If float, then the same time step is used for all steps. If
            an array, then each element k contains the time  at step k.
            Units are seconds.

        Returns
        -------

        x : numpy.ndarray
           smoothed means

        P : numpy.ndarray
           smoothed state covariances

        K : numpy.ndarray
            smoother gain at each step

        Examples
        --------

        .. code-block:: Python

            zs = [t + random.randn()*4 for t in range (40)]

            (mu, cov, _, _) = kalman.batch_filter(zs)
            (x, P, K) = rts_smoother(mu, cov, fk.F, fk.Q)
        """

        assert len(Xs) == len(Ps)
        n, dim_x = Xs.shape

        if dt is None:
            dt = [self._dt] * n
        elif isscalar(dt):
            dt = [dt] * n

        if Qs is None:
            Qs = [self.Q] * n

        # smoother gain
        Ks = zeros((n,dim_x,dim_x))

        num_sigmas = self._num_sigmas

        xs, ps = Xs.copy(), Ps.copy()
        sigmas_f = zeros((num_sigmas, dim_x))

        for k in range(n-2,-1,-1):
            # create sigma points from state estimate, pass through state func
            sigmas = self.points_fn.sigma_points(xs[k], ps[k])
            for i in range(num_sigmas):
                sigmas_f[i] = self.fx(sigmas[i], dt[k])

            # compute backwards prior state and covariance
            xb = dot(self.Wm, sigmas_f)
            Pb = 0
            x = Xs[k]
            for i in range(num_sigmas):
                y = self.residual_x(sigmas_f[i], x)
                Pb += self.Wc[i] * outer(y, y)
            Pb += Qs[k]

            # compute cross variance
            Pxb = 0
            for i in range(num_sigmas):
                z = self.residual_x(sigmas[i], Xs[k])
                y = self.residual_x(sigmas_f[i], xb)
                Pxb += self.Wc[i] * outer(z, y)

            # compute gain
            K = dot(Pxb, inv(Pb))

            # update the smoothed estimates
            xs[k] += dot (K, self.residual_x(xs[k+1], xb))
            ps[k] += dot3(K, ps[k+1] - Pb, K.T)
            Ks[k] = K

        return (xs, ps, Ks)

Example 25

Project: filterpy Source File: unscented_transform.py
Function: unscented_transform
def unscented_transform(sigmas, Wm, Wc, noise_cov=None,
                        mean_fn=None, residual_fn=None):
    """ Computes unscented transform of a set of sigma points and weights.
    returns the mean and covariance in a tuple.

    Parameters
    ----------

    sigmas: ndarray [#sigmas per dimension, dimension]
        2D array of sigma points.

    Wm : ndarray [# sigmas per dimension]
        Weights for the mean. Must sum to 1.


    Wc : ndarray [# sigmas per dimension]
        Weights for the covariance. Must sum to 1.

    noise_cov : ndarray, optional
        noise matrix added to the final computed covariance matrix.

    mean_fn : callable (sigma_points, weights), optional
        Function that computes the mean of the provided sigma points
        and weights. Use this if your state variable contains nonlinear
        values such as angles which cannot be summed.

        .. code-block:: Python

            def state_mean(sigmas, Wm):
                x = np.zeros(3)
                sum_sin, sum_cos = 0., 0.

                for i in range(len(sigmas)):
                    s = sigmas[i]
                    x[0] += s[0] * Wm[i]
                    x[1] += s[1] * Wm[i]
                    sum_sin += sin(s[2])*Wm[i]
                    sum_cos += cos(s[2])*Wm[i]
                x[2] = atan2(sum_sin, sum_cos)
                return x

    residual_fn : callable (x, y), optional

        Function that computes the residual (difference) between x and y.
        You will have to supply this if your state variable cannot support
        subtraction, such as angles (359-1 degreees is 2, not 358). x and y
        are state vectors, not scalars.

        .. code-block:: Python

            def residual(a, b):
                y = a[0] - b[0]
                y = y % (2 * np.pi)
                if y > np.pi:
                    y -= 2*np.pi
                return y


    Returns
    -------

    x : ndarray [dimension]
        Mean of the sigma points after passing through the transform.

    P : ndarray
        covariance of the sigma points after passing throgh the transform.
    """

    kmax, n = sigmas.shape


    if mean_fn is None:
        # new mean is just the sum of the sigmas * weight
        x = np.dot(Wm, sigmas)    # dot = \Sigma^n_1 (W[k]*Xi[k])
    else:
        x = mean_fn(sigmas, Wm)


    # new covariance is the sum of the outer product of the residuals
    # times the weights

    # this is the fast way to do this - see 'else' for the slow way
    if residual_fn is None:
        y = sigmas - x[np.newaxis,:]
        P = y.T.dot(np.diag(Wc)).dot(y)
    else:
        P = np.zeros((n, n))
        for k in range(kmax):
            y = residual_fn(sigmas[k], x)
            P += Wc[k] * np.outer(y, y)

    if noise_cov is not None:
        P += noise_cov

    return (x, P)

Example 26

Project: CorpusTools Source File: representations.py
def to_mfcc(filename, freq_lims,num_coeffs,win_len,time_step,num_filters = 26, use_power = False,debug=False):
    """Generate MFCCs in the style of HTK from a full path to a .wav file.

    Parameters
    ----------
    filename : str
        Full path to .wav file to process.
    freq_lims : tuple
        Minimum and maximum frequencies in Hertz to use.
    num_coeffs : int
        Number of coefficients of the cepstrum to return.
    win_len : float
        Window length in seconds to use for FFT.
    time_step : float
        Time step in seconds for windowing.
    num_filters : int
        Number of mel filters to use in the filter bank, defaults to 26.
    use_power : bool
        If true, use the first coefficient of the cepstrum, which is power
        based.  Defaults to false.

    Returns
    -------
    2D array
        MFCCs for each frame.  The first dimension is the time in frames,
        the second dimension is the MFCC values.

    """
    sr, proc = preproc(filename,alpha=0.97)

    minHz = freq_lims[0]
    maxHz = freq_lims[1]

    L = 22
    n = arange(num_filters)
    lift = 1+ (L/2)*sin(pi*n/L)
    lift = diag(lift)

    pspec = to_powerspec(proc,sr,win_len,time_step)

    filterbank = filter_bank((pspec.shape[1]-1) * 2,num_filters,minHz,maxHz,sr)


    num_frames = pspec.shape[0]

    mfccs = zeros((num_frames,num_coeffs))
    aspec =zeros((num_frames,num_filters))
    for k in range(num_frames):
        filteredSpectrum = dot(sqrt(pspec[k,:]), filterbank)**2
        aspec[k,:] = filteredSpectrum
        dctSpectrum = dct_spectrum(filteredSpectrum)
        dctSpectrum = dot(dctSpectrum , lift)
        if not use_power:
            dctSpectrum = dctSpectrum[1:]
        mfccs[k,:] = dctSpectrum[:num_coeffs]
    if debug:
        return mfccs,pspec,aspec
    return mfccs

Example 27

Project: pybrain Source File: cmaes.py
Function: learnstep
    def _learnStep(self):
        # Generate and evaluate lambda offspring
        arz = randn(self.numParameters, self.batchSize)
        arx = tile(self.center.reshape(self.numParameters, 1), (1, self.batchSize))\
                        + self.stepSize * dot(dot(self.B, self.D), arz)
        arfitness = zeros(self.batchSize)
        for k in range(self.batchSize):
            arfitness[k] = self._oneEvaluation(arx[:, k])

        # Sort by fitness and compute weighted mean into center
        arfitness, arindex = sorti(arfitness)  # minimization
        arz = arz[:, arindex]
        arx = arx[:, arindex]
        arzsel = arz[:, range(self.mu)]
        arxsel = arx[:, range(self.mu)]
        arxmut = arxsel - tile(self.center.reshape(self.numParameters, 1), (1, self.mu))

        zmean = dot(arzsel, self.weights)
        self.center = dot(arxsel, self.weights)

        if self.storeAllCenters:
            self._allCenters.append(self.center)

        # Cumulation: Update evolution paths
        self.stepPath = (1 - self.cuemStep) * self.stepPath \
                + sqrt(self.cuemStep * (2 - self.cuemStep) * self.muEff) * dot(self.B, zmean)         # Eq. (4)
        hsig = norm(self.stepPath) / sqrt(1 - (1 - self.cuemStep) ** (2 * self.numEvaluations / float(self.batchSize))) / self.chiN \
                    < 1.4 + 2. / (self.numParameters + 1)
        self.covPath = (1 - self.cuemCov) * self.covPath + hsig * \
                sqrt(self.cuemCov * (2 - self.cuemCov) * self.muEff) * dot(dot(self.B, self.D), zmean) # Eq. (2)

        # Adapt covariance matrix C
        self.C = ((1 - self.covLearningRate) * self.C                    # regard old matrix   % Eq. (3)
             + self.covLearningRate * (1 / self.muCov) * (outer(self.covPath, self.covPath) # plus rank one update
                                   + (1 - hsig) * self.cuemCov * (2 - self.cuemCov) * self.C)
             + self.covLearningRate * (1 - 1 / self.muCov)                 # plus rank mu update
             * dot(dot(arxmut, diag(self.weights)), arxmut.T)
            )

        # Adapt step size self.stepSize
        self.stepSize *= exp((self.cuemStep / self.dampings) * (norm(self.stepPath) / self.chiN - 1)) # Eq. (5)

        # Update B and D from C
        # This is O(n^3). When strategy internal CPU-time is critical, the
        # next three lines should be executed only every (alpha/covLearningRate/N)-th
        # iteration, where alpha is e.g. between 0.1 and 10
        self.C = (self.C + self.C.T) / 2 # enforce symmetry
        Ev, self.B = eig(self.C)          # eigen decomposition, B==normalized eigenvectors
        Ev = real(Ev)       # enforce real value
        self.D = diag(sqrt(Ev))      #diag(ravel(sqrt(Ev))) # D contains standard deviations now
        self.B = real(self.B)

        # convergence is reached
        if arfitness[0] == arfitness[-1] or (abs(arfitness[0] - arfitness[-1]) /
                                             (abs(arfitness[0]) + abs(arfitness[-1]))) <= self.stopPrecision:
            if self.verbose:
                print("Converged.")
            self.maxLearningSteps = self.numLearningSteps

        # or diverged, unfortunately
        if min(Ev) > 1e5:
            if self.verbose:
                print("Diverged.")
            self.maxLearningSteps = self.numLearningSteps

Example 28

Project: QuantEcon.py Source File: matrix_eqn.py
def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500):
    """
    Solves the discrete-time algebraic Riccati equation

        X = A'XA - (N + B'XA)'(B'XB + R)^{-1}(N + B'XA) + Q

    via a modified structured doubling algorithm.  An explanation of the
    algorithm can be found in the reference below.

    Note that SciPy also has a discrete riccati equation solver.  However it
    cannot handle the case where R is not invertible, or when N is nonzero.
    Both of these cases can be handled in the algorithm implemented below.

    Parameters
    ----------
    A : array_like(float, ndim=2)
        k x k array.
    B : array_like(float, ndim=2)
        k x n array
    Q : array_like(float, ndim=2)
        k x k, should be symmetric and non-negative definite
    R : array_like(float, ndim=2)
        n x n, should be symmetric and positive definite
    N : array_like(float, ndim=2)
        n x k array
    tolerance : scalar(float), optional(default=1e-10)
        The tolerance level for convergence
    max_iter : scalar(int), optional(default=500)
        The maximum number of iterations allowed

    Returns
    -------
    X : array_like(float, ndim=2)
        The fixed point of the Riccati equation; a  k x k array
        representing the approximate solution

    References
    ----------
    Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. "STRUCTURED DOUBLING
    ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR
    CONTROL WEIGHTING MATRICES." Taiwanese Journal of Mathematics 14, no. 3A
    (2010): pp-935.

    """
    # == Set up == #
    error = tolerance + 1
    fail_msg = "Convergence failed after {} iterations."

    # == Make sure that all array_likes are np arrays, two-dimensional == #
    A, B, Q, R = np.atleast_2d(A, B, Q, R)
    n, k = R.shape[0], Q.shape[0]
    I = np.identity(k)
    if N is None:
        N = np.zeros((n, k))
    else:
        N = np.atleast_2d(N)

    # == Choose optimal value of gamma in R_hat = R + gamma B'B == #
    current_min = np.inf
    candidates = (0.0, 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 10.0, 100.0, 10e5)
    BB = dot(B.T, B)
    BTA = dot(B.T, A)
    for gamma in candidates:
        Z = R + gamma * BB
        cn = np.linalg.cond(Z)
        if np.isfinite(cn):
            Q_tilde = - Q + dot(N.T, solve(Z, N + gamma * BTA)) + gamma * I
            G0 = dot(B, solve(Z, B.T))
            A0 = dot(I - gamma * G0, A) - dot(B, solve(Z, N))
            H0 = gamma * dot(A.T, A0) - Q_tilde
            f1 = np.linalg.cond(Z, np.inf)
            f2 = gamma * f1
            f3 = np.linalg.cond(I + dot(G0, H0))
            f_gamma = max(f1, f2, f3)
            if f_gamma < current_min:
                best_gamma = gamma
                current_min = f_gamma

    # == If no candidate successful then fail == #
    if current_min == np.inf:
        msg = "Unable to initialize routine due to ill conditioned arguments"
        raise ValueError(msg)

    gamma = best_gamma
    R_hat = R + gamma * BB

    # == Initial conditions == #
    Q_tilde = - Q + dot(N.T, solve(R_hat, N + gamma * BTA)) + gamma * I
    G0 = dot(B, solve(R_hat, B.T))
    A0 = dot(I - gamma * G0, A) - dot(B, solve(R_hat, N))
    H0 = gamma * dot(A.T, A0) - Q_tilde
    i = 1

    # == Main loop == #
    while error > tolerance:

        if i > max_iter:
            raise ValueError(fail_msg.format(i))

        else:
            A1 = dot(A0, solve(I + dot(G0, H0), A0))
            G1 = G0 + dot(dot(A0, G0), solve(I + dot(H0, G0), A0.T))
            H1 = H0 + dot(A0.T, solve(I + dot(H0, G0), dot(H0, A0)))

            error = np.max(np.abs(H1 - H0))
            A0 = A1
            G0 = G1
            H0 = H1
            i += 1

    return H1 + gamma * I  # Return X

Example 29

Project: rlpy Source File: Greedy_GQ.py
Function: learn
    def learn(self, s, p_actions, a, r, ns, np_actions, na, terminal):
        self.representation.pre_discover(s, False, a, ns, terminal)
        discount_factor = self.discount_factor
        weight_vec = self.representation.weight_vec
        phi_s = self.representation.phi(s, False)
        phi = self.representation.phi_sa(s, False, a, phi_s)
        phi_prime_s = self.representation.phi(ns, terminal)
        na = self.representation.bestAction(
            ns,
            terminal,
            np_actions,
            phi_prime_s)  # Switch na to the best possible action
        phi_prime = self.representation.phi_sa(
            ns,
            terminal,
            na,
            phi_prime_s)
        nnz = count_nonzero(phi_s)    # Number of non-zero elements

        expanded = (- len(self.GQWeight) + len(phi)) / self.representation.actions_num
        if expanded:
            self._expand_vectors(expanded)
        # Set eligibility traces:
        if self.lambda_:
            self.eligibility_trace *= discount_factor * self.lambda_
            self.eligibility_trace += phi

            self.eligibility_trace_s *= discount_factor * self.lambda_
            self.eligibility_trace_s += phi_s

            # Set max to 1
            self.eligibility_trace[self.eligibility_trace > 1] = 1
            self.eligibility_trace_s[self.eligibility_trace_s > 1] = 1
        else:
            self.eligibility_trace = phi
            self.eligibility_trace_s = phi_s

        td_error                     = r + \
            np.dot(discount_factor * phi_prime - phi, weight_vec)
        self.updateLearnRate(
            phi_s,
            phi_prime_s,
            self.eligibility_trace_s,
            discount_factor,
            nnz,
            terminal)

        if nnz > 0:  # Phi has some nonzero elements, proceed with update
            td_error_estimate_now = np.dot(phi, self.GQWeight)
            Delta_weight_vec                 = td_error * self.eligibility_trace - \
                discount_factor * td_error_estimate_now * phi_prime
            weight_vec += self.learn_rate * Delta_weight_vec
            Delta_GQWeight = (
                td_error - td_error_estimate_now) * phi
            self.GQWeight               += self.learn_rate * \
                self.secondLearningRateCoef * Delta_GQWeight

        expanded = self.representation.post_discover(
            s,
            False,
            a,
            td_error,
            phi_s)
        if expanded:
            self._expand_vectors(expanded)
        if terminal:
            self.episodeTerminated()

Example 30

Project: scikit-learn Source File: test_pls.py
def test_pls():
    d = load_linnerud()
    X = d.data
    Y = d.target
    # 1) Canonical (symmetric) PLS (PLS 2 blocks canonical mode A)
    # ===========================================================
    # Compare 2 algo.: nipals vs. svd
    # ------------------------------
    pls_bynipals = pls_.PLSCanonical(n_components=X.shape[1])
    pls_bynipals.fit(X, Y)
    pls_bysvd = pls_.PLSCanonical(algorithm="svd", n_components=X.shape[1])
    pls_bysvd.fit(X, Y)
    # check equalities of loading (up to the sign of the second column)
    assert_array_almost_equal(
        pls_bynipals.x_loadings_,
        pls_bysvd.x_loadings_, decimal=5,
        err_msg="nipals and svd implementations lead to different x loadings")

    assert_array_almost_equal(
        pls_bynipals.y_loadings_,
        pls_bysvd.y_loadings_, decimal=5,
        err_msg="nipals and svd implementations lead to different y loadings")

    # Check PLS properties (with n_components=X.shape[1])
    # ---------------------------------------------------
    plsca = pls_.PLSCanonical(n_components=X.shape[1])
    plsca.fit(X, Y)
    T = plsca.x_scores_
    P = plsca.x_loadings_
    Wx = plsca.x_weights_
    U = plsca.y_scores_
    Q = plsca.y_loadings_
    Wy = plsca.y_weights_

    def check_ortho(M, err_msg):
        K = np.dot(M.T, M)
        assert_array_almost_equal(K, np.diag(np.diag(K)), err_msg=err_msg)

    # Orthogonality of weights
    # ~~~~~~~~~~~~~~~~~~~~~~~~
    check_ortho(Wx, "x weights are not orthogonal")
    check_ortho(Wy, "y weights are not orthogonal")

    # Orthogonality of latent scores
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    check_ortho(T, "x scores are not orthogonal")
    check_ortho(U, "y scores are not orthogonal")

    # Check X = TP' and Y = UQ' (with (p == q) components)
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # center scale X, Y
    Xc, Yc, x_mean, y_mean, x_std, y_std =\
        pls_._center_scale_xy(X.copy(), Y.copy(), scale=True)
    assert_array_almost_equal(Xc, np.dot(T, P.T), err_msg="X != TP'")
    assert_array_almost_equal(Yc, np.dot(U, Q.T), err_msg="Y != UQ'")

    # Check that rotations on training data lead to scores
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    Xr = plsca.transform(X)
    assert_array_almost_equal(Xr, plsca.x_scores_,
                              err_msg="rotation on X failed")
    Xr, Yr = plsca.transform(X, Y)
    assert_array_almost_equal(Xr, plsca.x_scores_,
                              err_msg="rotation on X failed")
    assert_array_almost_equal(Yr, plsca.y_scores_,
                              err_msg="rotation on Y failed")

    # "Non regression test" on canonical PLS
    # --------------------------------------
    # The results were checked against the R-package plspm
    pls_ca = pls_.PLSCanonical(n_components=X.shape[1])
    pls_ca.fit(X, Y)

    x_weights = np.array(
        [[-0.61330704,  0.25616119, -0.74715187],
         [-0.74697144,  0.11930791,  0.65406368],
         [-0.25668686, -0.95924297, -0.11817271]])
    # x_weights_sign_flip holds columns of 1 or -1, depending on sign flip
    # between R and python
    x_weights_sign_flip = pls_ca.x_weights_ / x_weights

    x_rotations = np.array(
        [[-0.61330704,  0.41591889, -0.62297525],
         [-0.74697144,  0.31388326,  0.77368233],
         [-0.25668686, -0.89237972, -0.24121788]])
    x_rotations_sign_flip = pls_ca.x_rotations_ / x_rotations

    y_weights = np.array(
        [[+0.58989127,  0.7890047,   0.1717553],
         [+0.77134053, -0.61351791,  0.16920272],
         [-0.23887670, -0.03267062,  0.97050016]])
    y_weights_sign_flip = pls_ca.y_weights_ / y_weights

    y_rotations = np.array(
        [[+0.58989127,  0.7168115,  0.30665872],
         [+0.77134053, -0.70791757,  0.19786539],
         [-0.23887670, -0.00343595,  0.94162826]])
    y_rotations_sign_flip = pls_ca.y_rotations_ / y_rotations

    # x_weights = X.dot(x_rotation)
    # Hence R/python sign flip should be the same in x_weight and x_rotation
    assert_array_almost_equal(x_rotations_sign_flip, x_weights_sign_flip)
    # This test that R / python give the same result up to column
    # sign indeterminacy
    assert_array_almost_equal(np.abs(x_rotations_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(x_weights_sign_flip), 1, 4)


    assert_array_almost_equal(y_rotations_sign_flip, y_weights_sign_flip)
    assert_array_almost_equal(np.abs(y_rotations_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(y_weights_sign_flip), 1, 4)

    # 2) Regression PLS (PLS2): "Non regression test"
    # ===============================================
    # The results were checked against the R-packages plspm, misOmics and pls
    pls_2 = pls_.PLSRegression(n_components=X.shape[1])
    pls_2.fit(X, Y)

    x_weights = np.array(
        [[-0.61330704, -0.00443647,  0.78983213],
         [-0.74697144, -0.32172099, -0.58183269],
         [-0.25668686,  0.94682413, -0.19399983]])
    x_weights_sign_flip = pls_2.x_weights_ / x_weights

    x_loadings = np.array(
        [[-0.61470416, -0.24574278,  0.78983213],
         [-0.65625755, -0.14396183, -0.58183269],
         [-0.51733059,  1.00609417, -0.19399983]])
    x_loadings_sign_flip = pls_2.x_loadings_ / x_loadings

    y_weights = np.array(
        [[+0.32456184,  0.29892183,  0.20316322],
         [+0.42439636,  0.61970543,  0.19320542],
         [-0.13143144, -0.26348971, -0.17092916]])
    y_weights_sign_flip = pls_2.y_weights_ / y_weights

    y_loadings = np.array(
        [[+0.32456184,  0.29892183,  0.20316322],
         [+0.42439636,  0.61970543,  0.19320542],
         [-0.13143144, -0.26348971, -0.17092916]])
    y_loadings_sign_flip = pls_2.y_loadings_ / y_loadings

    # x_loadings[:, i] = Xi.dot(x_weights[:, i]) \forall i
    assert_array_almost_equal(x_loadings_sign_flip, x_weights_sign_flip, 4)
    assert_array_almost_equal(np.abs(x_loadings_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(x_weights_sign_flip), 1, 4)

    assert_array_almost_equal(y_loadings_sign_flip, y_weights_sign_flip, 4)
    assert_array_almost_equal(np.abs(y_loadings_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(y_weights_sign_flip), 1, 4)

    # 3) Another non-regression test of Canonical PLS on random dataset
    # =================================================================
    # The results were checked against the R-package plspm
    n = 500
    p_noise = 10
    q_noise = 5
    # 2 latents vars:
    np.random.seed(11)
    l1 = np.random.normal(size=n)
    l2 = np.random.normal(size=n)
    latents = np.array([l1, l1, l2, l2]).T
    X = latents + np.random.normal(size=4 * n).reshape((n, 4))
    Y = latents + np.random.normal(size=4 * n).reshape((n, 4))
    X = np.concatenate(
        (X, np.random.normal(size=p_noise * n).reshape(n, p_noise)), axis=1)
    Y = np.concatenate(
        (Y, np.random.normal(size=q_noise * n).reshape(n, q_noise)), axis=1)
    np.random.seed(None)
    pls_ca = pls_.PLSCanonical(n_components=3)
    pls_ca.fit(X, Y)

    x_weights = np.array(
        [[0.65803719,  0.19197924,  0.21769083],
         [0.7009113,  0.13303969, -0.15376699],
         [0.13528197, -0.68636408,  0.13856546],
         [0.16854574, -0.66788088, -0.12485304],
         [-0.03232333, -0.04189855,  0.40690153],
         [0.1148816, -0.09643158,  0.1613305],
         [0.04792138, -0.02384992,  0.17175319],
         [-0.06781, -0.01666137, -0.18556747],
         [-0.00266945, -0.00160224,  0.11893098],
         [-0.00849528, -0.07706095,  0.1570547],
         [-0.00949471, -0.02964127,  0.34657036],
         [-0.03572177,  0.0945091,  0.3414855],
         [0.05584937, -0.02028961, -0.57682568],
         [0.05744254, -0.01482333, -0.17431274]])
    x_weights_sign_flip = pls_ca.x_weights_ / x_weights


    x_loadings = np.array(
        [[0.65649254,  0.1847647,  0.15270699],
         [0.67554234,  0.15237508, -0.09182247],
         [0.19219925, -0.67750975,  0.08673128],
         [0.2133631, -0.67034809, -0.08835483],
         [-0.03178912, -0.06668336,  0.43395268],
         [0.15684588, -0.13350241,  0.20578984],
         [0.03337736, -0.03807306,  0.09871553],
         [-0.06199844,  0.01559854, -0.1881785],
         [0.00406146, -0.00587025,  0.16413253],
         [-0.00374239, -0.05848466,  0.19140336],
         [0.00139214, -0.01033161,  0.32239136],
         [-0.05292828,  0.0953533,  0.31916881],
         [0.04031924, -0.01961045, -0.65174036],
         [0.06172484, -0.06597366, -0.1244497]])
    x_loadings_sign_flip = pls_ca.x_loadings_ / x_loadings

    y_weights = np.array(
        [[0.66101097,  0.18672553,  0.22826092],
         [0.69347861,  0.18463471, -0.23995597],
         [0.14462724, -0.66504085,  0.17082434],
         [0.22247955, -0.6932605, -0.09832993],
         [0.07035859,  0.00714283,  0.67810124],
         [0.07765351, -0.0105204, -0.44108074],
         [-0.00917056,  0.04322147,  0.10062478],
         [-0.01909512,  0.06182718,  0.28830475],
         [0.01756709,  0.04797666,  0.32225745]])
    y_weights_sign_flip = pls_ca.y_weights_ / y_weights

    y_loadings = np.array(
        [[0.68568625,  0.1674376,  0.0969508],
         [0.68782064,  0.20375837, -0.1164448],
         [0.11712173, -0.68046903,  0.12001505],
         [0.17860457, -0.6798319, -0.05089681],
         [0.06265739, -0.0277703,  0.74729584],
         [0.0914178,  0.00403751, -0.5135078],
         [-0.02196918, -0.01377169,  0.09564505],
         [-0.03288952,  0.09039729,  0.31858973],
         [0.04287624,  0.05254676,  0.27836841]])
    y_loadings_sign_flip = pls_ca.y_loadings_ / y_loadings

    assert_array_almost_equal(x_loadings_sign_flip, x_weights_sign_flip, 4)
    assert_array_almost_equal(np.abs(x_weights_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(x_loadings_sign_flip), 1, 4)

    assert_array_almost_equal(y_loadings_sign_flip, y_weights_sign_flip, 4)
    assert_array_almost_equal(np.abs(y_weights_sign_flip), 1, 4)
    assert_array_almost_equal(np.abs(y_loadings_sign_flip), 1, 4)

    # Orthogonality of weights
    # ~~~~~~~~~~~~~~~~~~~~~~~~
    check_ortho(pls_ca.x_weights_, "x weights are not orthogonal")
    check_ortho(pls_ca.y_weights_, "y weights are not orthogonal")

    # Orthogonality of latent scores
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    check_ortho(pls_ca.x_scores_, "x scores are not orthogonal")
    check_ortho(pls_ca.y_scores_, "y scores are not orthogonal")

Example 31

Project: orbital Source File: utilities.py
def elements_from_state_vector(r, v, mu):
    h = angular_momentum(r, v)
    n = node_vector(h)

    ev = eccentricity_vector(r, v, mu)

    E = specific_orbital_energy(r, v, mu)

    a = -mu / (2 * E)
    e = norm(ev)

    SMALL_NUMBER = 1e-15

    # Inclination is the angle between the angular
    # momentum vector and its z component.
    i = acos(h.z / norm(h))

    if abs(i - 0) < SMALL_NUMBER:
        # For non-inclined orbits, raan is undefined;
        # set to zero by convention
        raan = 0
        if abs(e - 0) < SMALL_NUMBER:
            # For circular orbits, place periapsis
            # at ascending node by convention
            arg_pe = 0
        else:
            # Argument of periapsis is the angle between
            # eccentricity vector and its x component.
            arg_pe = acos(ev.x / norm(ev))
    else:
        # Right ascension of ascending node is the angle
        # between the node vector and its x component.
        raan = acos(n.x / norm(n))
        if n.y < 0:
            raan = 2 * pi - raan

        # Argument of periapsis is angle between
        # node and eccentricity vectors.
        arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev)))

    if abs(e - 0) < SMALL_NUMBER:
        if abs(i - 0) < SMALL_NUMBER:
            # True anomaly is angle between position
            # vector and its x component.
            f = acos(r.x / norm(r))
            if v.x > 0:
                f = 2 * pi - f
        else:
            # True anomaly is angle between node
            # vector and position vector.
            f = acos(dot(n, r) / (norm(n) * norm(r)))
            if dot(n, v) > 0:
                f = 2 * pi - f
    else:
        if ev.z < 0:
            arg_pe = 2 * pi - arg_pe

        # True anomaly is angle between eccentricity
        # vector and position vector.
        f = acos(dot(ev, r) / (norm(ev) * norm(r)))

        if dot(r, v) < 0:
            f = 2 * pi - f

    return OrbitalElements(a=a, e=e, i=i, raan=raan, arg_pe=arg_pe, f=f)

Example 32

Project: scikit-learn Source File: nmf.py
def _nls_subproblem(V, W, H, tol, max_iter, alpha=0., l1_ratio=0.,
                    sigma=0.01, beta=0.1):
    """Non-negative least square solver

    Solves a non-negative least squares subproblem using the projected
    gradient descent algorithm.

    Parameters
    ----------
    V : array-like, shape (n_samples, n_features)
        Constant matrix.

    W : array-like, shape (n_samples, n_components)
        Constant matrix.

    H : array-like, shape (n_components, n_features)
        Initial guess for the solution.

    tol : float
        Tolerance of the stopping condition.

    max_iter : int
        Maximum number of iterations before timing out.

    alpha : double, default: 0.
        Constant that multiplies the regularization terms. Set it to zero to
        have no regularization.

    l1_ratio : double, default: 0.
        The regularization mixing parameter, with 0 <= l1_ratio <= 1.
        For l1_ratio = 0 the penalty is an L2 penalty.
        For l1_ratio = 1 it is an L1 penalty.
        For 0 < l1_ratio < 1, the penalty is a combination of L1 and L2.

    sigma : float
        Constant used in the sufficient decrease condition checked by the line
        search.  Smaller values lead to a looser sufficient decrease condition,
        thus reducing the time taken by the line search, but potentially
        increasing the number of iterations of the projected gradient
        procedure. 0.01 is a commonly used value in the optimization
        literature.

    beta : float
        Factor by which the step size is decreased (resp. increased) until
        (resp. as long as) the sufficient decrease condition is satisfied.
        Larger values allow to find a better step size but lead to longer line
        search. 0.1 is a commonly used value in the optimization literature.

    Returns
    -------
    H : array-like, shape (n_components, n_features)
        Solution to the non-negative least squares problem.

    grad : array-like, shape (n_components, n_features)
        The gradient.

    n_iter : int
        The number of iterations done by the algorithm.

    References
    ----------
    C.-J. Lin. Projected gradient methods for non-negative matrix
    factorization. Neural Computation, 19(2007), 2756-2779.
    http://www.csie.ntu.edu.tw/~cjlin/nmf/
    """
    WtV = safe_sparse_dot(W.T, V)
    WtW = fast_dot(W.T, W)

    # values justified in the paper (alpha is renamed gamma)
    gamma = 1
    for n_iter in range(1, max_iter + 1):
        grad = np.dot(WtW, H) - WtV
        if alpha > 0 and l1_ratio == 1.:
            grad += alpha
        elif alpha > 0:
            grad += alpha * (l1_ratio + (1 - l1_ratio) * H)

        # The following multiplication with a boolean array is more than twice
        # as fast as indexing into grad.
        if norm(grad * np.logical_or(grad < 0, H > 0)) < tol:
            break

        Hp = H

        for inner_iter in range(20):
            # Gradient step.
            Hn = H - gamma * grad
            # Projection step.
            Hn *= Hn > 0
            d = Hn - H
            gradd = np.dot(grad.ravel(), d.ravel())
            dQd = np.dot(np.dot(WtW, d).ravel(), d.ravel())
            suff_decr = (1 - sigma) * gradd + 0.5 * dQd < 0
            if inner_iter == 0:
                decr_gamma = not suff_decr

            if decr_gamma:
                if suff_decr:
                    H = Hn
                    break
                else:
                    gamma *= beta
            elif not suff_decr or (Hp == Hn).all():
                H = Hp
                break
            else:
                gamma /= beta
                Hp = Hn

    if n_iter == max_iter:
        warnings.warn("Iteration limit reached in nls subproblem.")

    return H, grad, n_iter

Example 33

Project: RLScore Source File: test_rls.py
Function: test_linear
    def test_linear(self):
        #Test that learning with linear kernel works correctly both
        #with low and high-dimensional data
        for X in [self.Xtrain1, self.Xtrain2]:
            for Y in [self.Ytrain1, self.Ytrain2]:
                #Basic case
                primal_rls = RLS(X, Y, regparam=1.0, bias=0.)
                W = primal_rls.predictor.W
                d = X.shape[1]
                W2 = np.linalg.solve(np.dot(X.T, X) + np.eye(d), np.dot(X.T, Y))
                assert_allclose(W, W2)
                #Fast regularization algorithm
                primal_rls.solve(10.)
                W = primal_rls.predictor.W
                W2 = np.linalg.solve(np.dot(X.T, X) + 10.*np.eye(d), np.dot(X.T, Y))
                assert_allclose(W, W2)
                #Bias term included
                primal_rls = RLS(X, Y, regparam=1.0, bias=2.)
                O = np.sqrt(2.) * np.ones((X.shape[0],1))
                X_new = np.hstack((X, O))
                W = primal_rls.predictor.W
                W2 = np.linalg.solve(np.dot(X_new.T, X_new) + np.eye(d+1), np.dot(X_new.T, Y))
                b = primal_rls.predictor.b
                b2 = W2[-1]
                W2 = W2[:-1]
                assert_allclose(W, W2)
                assert_allclose(b, np.sqrt(2) * b2)
                #reduced set approximation
                primal_rls = RLS(X, Y, basis_vectors = X[self.bvectors], regparam=5.0, bias=2.)
                W = primal_rls.predictor.W
                b = primal_rls.predictor.b
                K = np.dot(X_new, X_new.T)
                Kr = K[:, self.bvectors]
                Krr = K[np.ix_(self.bvectors, self.bvectors)]
                A = np.linalg.solve(np.dot(Kr.T, Kr)+ 5.0 * Krr, np.dot(Kr.T, Y))
                W2 = np.dot(X_new[self.bvectors].T, A)
                b2 = W2[-1]
                W2 = W2[:-1]
                assert_allclose(W, W2)
                assert_allclose(b, np.sqrt(2) * b2)
                #Using pre-computed linear kernel matrix
                kernel = LinearKernel(X, bias = 2.)
                K = kernel.getKM(X)
                dual_rls = RLS(K, Y, kernel = "PrecomputedKernel", regparam=0.01)
                W = np.dot(X_new.T, dual_rls.predictor.W)
                b = W[-1]
                W = W[:-1]
                W2 = np.linalg.solve(np.dot(X_new.T, X_new) + 0.01 * np.eye(d+1), np.dot(X_new.T, Y))
                b2 = W2[-1]
                W2 = W2[:-1]
                assert_allclose(W, W2)
                assert_allclose(b, b2)
                #Pre-computed linear kernel, reduced set approximation
                kernel = LinearKernel(X[self.bvectors], bias = 2.)
                dual_rls = RLS(kernel.getKM(X), Y, kernel="PrecomputedKernel", basis_vectors = kernel.getKM(X[self.bvectors]), regparam=5.0)
                W = np.dot(X_new[self.bvectors].T, dual_rls.predictor.W)
                b = W[-1]
                W = W[:-1]
                K = np.dot(X_new, X_new.T)
                Kr = K[:, self.bvectors]
                Krr = K[np.ix_(self.bvectors, self.bvectors)]
                A = np.linalg.solve(np.dot(Kr.T, Kr)+ 5.0 * Krr, np.dot(Kr.T, Y))
                W2 = np.dot(X_new[self.bvectors].T, A)
                b2 = W2[-1]
                W2 = W2[:-1]
                assert_allclose(W, W2)
                assert_allclose(b, b2)

Example 34

Project: scikit-learn Source File: extmath.py
Function: randomized_svd
def randomized_svd(M, n_components, n_oversamples=10, n_iter='auto',
                   power_iteration_normalizer='auto', transpose='auto',
                   flip_sign=True, random_state=0):
    """Computes a truncated randomized SVD

    Parameters
    ----------
    M: ndarray or sparse matrix
        Matrix to decompose

    n_components: int
        Number of singular values and vectors to extract.

    n_oversamples: int (default is 10)
        Additional number of random vectors to sample the range of M so as
        to ensure proper conditioning. The total number of random vectors
        used to find the range of M is n_components + n_oversamples. Smaller
        number can improve speed but can negatively impact the quality of
        approximation of singular vectors and singular values.

    n_iter: int or 'auto' (default is 'auto')
        Number of power iterations. It can be used to deal with very noisy
        problems. When 'auto', it is set to 4, unless `n_components` is small
        (< .1 * min(X.shape)) `n_iter` in which case is set to 7.
        This improves precision with few components.

        .. versionchanged:: 0.18

    power_iteration_normalizer: 'auto' (default), 'QR', 'LU', 'none'
        Whether the power iterations are normalized with step-by-step
        QR factorization (the slowest but most accurate), 'none'
        (the fastest but numerically unstable when `n_iter` is large, e.g.
        typically 5 or larger), or 'LU' factorization (numerically stable
        but can lose slightly in accuracy). The 'auto' mode applies no
        normalization if `n_iter`<=2 and switches to LU otherwise.

        .. versionadded:: 0.18

    transpose: True, False or 'auto' (default)
        Whether the algorithm should be applied to M.T instead of M. The
        result should approximately be the same. The 'auto' mode will
        trigger the transposition if M.shape[1] > M.shape[0] since this
        implementation of randomized SVD tend to be a little faster in that
        case.

        .. versionchanged:: 0.18

    flip_sign: boolean, (True by default)
        The output of a singular value decomposition is only unique up to a
        permutation of the signs of the singular vectors. If `flip_sign` is
        set to `True`, the sign ambiguity is resolved by making the largest
        loadings for each component in the left singular vectors positive.

    random_state: RandomState or an int seed (0 by default)
        A random number generator instance to make behavior

    Notes
    -----
    This algorithm finds a (usually very good) approximate truncated
    singular value decomposition using randomization to speed up the
    computations. It is particularly fast on large matrices on which
    you wish to extract only a small number of components. In order to
    obtain further speed up, `n_iter` can be set <=2 (at the cost of
    loss of precision).

    References
    ----------
    * Finding structure with randomness: Stochastic algorithms for constructing
      approximate matrix decompositions
      Halko, et al., 2009 http://arxiv.org/abs/arXiv:0909.4061

    * A randomized algorithm for the decomposition of matrices
      Per-Gunnar Martinsson, Vladimir Rokhlin and Mark Tygert

    * An implementation of a randomized algorithm for principal component
      analysis
      A. Szlam et al. 2014
    """
    random_state = check_random_state(random_state)
    n_random = n_components + n_oversamples
    n_samples, n_features = M.shape

    if n_iter == 'auto':
        # Checks if the number of iterations is explicitely specified
        # Adjust n_iter. 7 was found a good compromise for PCA. See #5299
        n_iter = 7 if n_components < .1 * min(M.shape) else 4

    if transpose == 'auto':
        transpose = n_samples < n_features
    if transpose:
        # this implementation is a bit faster with smaller shape[1]
        M = M.T

    Q = randomized_range_finder(M, n_random, n_iter,
                                power_iteration_normalizer, random_state)

    # project M to the (k + p) dimensional space using the basis vectors
    B = safe_sparse_dot(Q.T, M)

    # compute the SVD on the thin matrix: (k + p) wide
    Uhat, s, V = linalg.svd(B, full_matrices=False)
    del B
    U = np.dot(Q, Uhat)

    if flip_sign:
        if not transpose:
            U, V = svd_flip(U, V)
        else:
            # In case of transpose u_based_decision=false
            # to actually flip based on u and not v.
            U, V = svd_flip(U, V, u_based_decision=False)

    if transpose:
        # transpose back the results according to the input convention
        return V[:n_components, :].T, s[:n_components], U[:, :n_components].T
    else:
        return U[:, :n_components], s[:n_components], V[:n_components, :]

Example 35

Project: RMG-Py Source File: network.py
    def solveFullME(self, tlist, x0):
        """
        Directly solve the full master equation using a stiff ODE solver. Pass the
        reaction `network` to solve, the temperature `T` in K and pressure `P` in
        Pa to solve at, the energies `Elist` in J/mol to use, the output time
        points `tlist` in s, the initial total populations `x0`, the full master
        equation matrix `M`, the accounting matrix `indices` relating isomer and
        energy grain indices to indices of the master equation matrix, and the
        densities of states `densStates` in mol/J of each isomer.
        Returns the times in s, population distributions for each isomer, and total
        population profiles for each configuration.
        """
        import scipy.integrate
    
        Elist = self.Elist
        Jlist = self.Jlist
        densStates = self.densStates
        
        Nisom = self.Nisom
        Nreac = self.Nreac
        Nprod = self.Nprod
        Ngrains = len(Elist)
        NJ = len(Jlist)
        Ntime = len(tlist)
        
        def residual(t, y, K):
            return numpy.dot(K, y)
        
        def jacobian(t, y, K):
            return K
    
        ymB = self.P / constants.R / self.T
        M, indices = self.generateFullMEMatrix()
        Nrows = M.shape[0]
        M[:,Nrows-Nreac-Nprod:] *= ymB
        
        if self.ymB is not None:
            if isinstance(self.ymB, float):
                assert Nreac <= 1
                M[:,Nrows-Nreac-Nprod:] *= self.ymB
            else:
                for n in range(Nreac+Nprod):
                    M[:,Nrows-Nreac-Nprod+n] *= self.ymB[n]
        
        # Get equilibrium distributions
        eqDist = numpy.zeros_like(densStates)
        for i in range(Nisom):
            for s in range(NJ):
                eqDist[i,:,s] = densStates[i,:,s] * (2*Jlist[s]+1) * numpy.exp(-Elist / constants.R / self.T)
            eqDist[i,:,:] /= sum(eqDist[i,:,:])
    
        # Set initial conditions
        p0 = numpy.zeros([M.shape[0]], float)
        for i in range(Nisom):
            for r in range(Ngrains):
                for s in range(NJ):
                    index = indices[i,r,s]
                    if indices[i,r,s] > 0:
                        p0[index] = x0[i] * eqDist[i,r,s]
        for i in range(Nreac+Nprod):
            p0[-Nreac-Nprod + i] = x0[i+Nisom]
    
        # Set up ODEs
        ode = scipy.integrate.ode(residual, jacobian).set_integrator('vode', method='bdf', with_jacobian=True, atol=1e-16, rtol=1e-8)
        ode.set_initial_value(p0, 0.0).set_f_params(M).set_jac_params(M)
    
        # Generate solution
        t = numpy.zeros([Ntime], float)
        p = numpy.zeros([Ntime, Nisom, Ngrains, NJ], float)
        x = numpy.zeros([Ntime, Nisom+Nreac+Nprod], float)
        for m in range(Ntime):
            ode.integrate(tlist[m])
            t[m] = ode.t
            for r in range(Ngrains):
                for s in range(NJ):
                    for i in range(0, Nisom):
                        index = indices[i,r,s]
                        if index > 0:
                            p[m,i,r,s] += ode.y[index]
                            x[m,i] += ode.y[index]
            for n in range(Nisom, Nisom+Nreac+Nprod):
                x[m,n] = ode.y[-(Nisom+Nreac+Nprod)+n]
    
        return t, p, x

Example 36

Project: director Source File: polarisplatformplanner.py
    def fitRunningBoardAtFeet(self):

        # get stance frame
        startPose = self.getPlanningStartPose()
        stanceFrame = self.robotSystem.footstepsDriver.getFeetMidPoint(self.robotSystem.robotStateModel, useWorldZ=False)
        stanceFrameAxes = transformUtils.getAxesFromTransform(stanceFrame)

        # get pointcloud and extract search region covering the running board
        polyData = segmentation.getCurrentRevolutionData()
        polyData = segmentation.applyVoxelGrid(polyData, leafSize=0.01)
        _, polyData = segmentation.removeGround(polyData)
        polyData = segmentation.cropToBox(polyData, stanceFrame, [1.0, 1.0, 0.1])

        if not polyData.GetNumberOfPoints():
            print 'empty search region point cloud'
            return

        vis.updatePolyData(polyData, 'running board search points', parent=segmentation.getDebugFolder(), color=[0,1,0], visible=False)

        # extract maximal points along the stance x axis
        perpAxis = stanceFrameAxes[0]
        edgeAxis = stanceFrameAxes[1]
        edgePoints = segmentation.computeEdge(polyData, edgeAxis, perpAxis)
        edgePoints = vnp.getVtkPolyDataFromNumpyPoints(edgePoints)
        vis.updatePolyData(edgePoints, 'edge points', parent=segmentation.getDebugFolder(), visible=True)

        # ransac fit a line to the edge points
        linePoint, lineDirection, fitPoints = segmentation.applyLineFit(edgePoints)
        if np.dot(lineDirection, stanceFrameAxes[1]) < 0:
            lineDirection = -lineDirection

        linePoints = segmentation.thresholdPoints(fitPoints, 'ransac_labels', [1.0, 1.0])
        dists = np.dot(vnp.getNumpyFromVtk(linePoints, 'Points')-linePoint, lineDirection)
        p1 = linePoint + lineDirection*np.min(dists)
        p2 = linePoint + lineDirection*np.max(dists)
        vis.updatePolyData(fitPoints, 'line fit points', parent=segmentation.getDebugFolder(), colorByName='ransac_labels', visible=False)


        # compute a new frame that is in plane with the stance frame
        # and matches the orientation and position of the detected edge
        origin = np.array(stanceFrame.GetPosition())
        normal = np.array(stanceFrameAxes[2])

        # project stance origin to edge, then back to foot frame
        originProjectedToEdge = linePoint + lineDirection*np.dot(origin - linePoint, lineDirection)
        originProjectedToPlane = segmentation.projectPointToPlane(originProjectedToEdge, origin, normal)
        zaxis = np.array(stanceFrameAxes[2])
        yaxis = np.array(lineDirection)
        xaxis = np.cross(yaxis, zaxis)
        xaxis /= np.linalg.norm(xaxis)
        yaxis = np.cross(zaxis, xaxis)
        yaxis /= np.linalg.norm(yaxis)

        d = DebugData()
        d.addSphere(p1, radius=0.005)
        d.addSphere(p2, radius=0.005)
        d.addLine(p1, p2)
        d.addSphere(originProjectedToEdge, radius=0.001, color=[1,0,0])
        d.addSphere(originProjectedToPlane, radius=0.001, color=[0,1,0])
        d.addLine(originProjectedToPlane, origin, color=[0,1,0])
        d.addLine(originProjectedToEdge, origin, color=[1,0,0])
        vis.updatePolyData(d.getPolyData(), 'running board edge', parent=segmentation.getDebugFolder(), colorByName='RGB255', visible=False)

        # update the running board box affordance position and orientation to
        # fit the detected edge
        box = self.spawnRunningBoardAffordance()
        boxDimensions = box.getProperty('Dimensions')
        t = transformUtils.getTransformFromAxesAndOrigin(xaxis, yaxis, zaxis, originProjectedToPlane)
        t.PreMultiply()
        t.Translate(-boxDimensions[0]/2.0, 0.0, -boxDimensions[2]/2.0)
        box.getChildFrame().copyFrame(t)

        self.initialize()

Example 37

Project: qutip Source File: random_objects.py
def rand_super_bcsz(N=2, enforce_tp=True, rank=None, dims=None):
    """
    Returns a random superoperator drawn from the Bruzda
    et al ensemble for CPTP maps [BCSZ08]_. Note that due to
    finite numerical precision, for ranks less than full-rank,
    zero eigenvalues may become slightly negative, such that the
    returned operator is not actually completely positive.


    Parameters
    ----------
    N : int
        Square root of the dimension of the superoperator to be returned.
    enforce_tp : bool
        If True, the trace-preserving condition of [BCSZ08]_ is enforced;
        otherwise only complete positivity is enforced.
    rank : int or None
        Rank of the sampled superoperator. If None, a full-rank
        superoperator is generated.
    dims : list
        Dimensions of quantum object.  Used for specifying
        tensor structure. Default is dims=[[[N],[N]], [[N],[N]]].

    Returns
    -------
    rho : Qobj
        A superoperator acting on vectorized dim × dim density operators,
        sampled from the BCSZ distribution.
    """
    if dims is not None:
        # TODO: check!
        pass
    else:
        dims = [[[N],[N]], [[N],[N]]]

    if rank is None:
        rank = N**2
    if rank > N**2:
        raise ValueError("Rank cannot exceed superoperator dimension.")

    # We use mainly dense matrices here for speed in low
    # dimensions. In the future, it would likely be better to switch off
    # between sparse and dense matrices as the dimension grows.

    # We start with a Ginibre uniform matrix X of the appropriate rank,
    # and use it to construct a positive semidefinite matrix X X⁺.
    X = randnz((N**2, rank), norm='ginibre')
    
    # Precompute X X⁺, as we'll need it in two different places.
    XXdag = np.dot(X, X.T.conj())
    
    if enforce_tp:
        # We do the partial trace over the first index by using dense reshape
        # operations, so that we can avoid bouncing to a sparse representation
        # and back.
        Y = np.einsum('ijik->jk', XXdag.reshape((N, N, N, N)))

        # Now we have the matrix 𝟙 ⊗ Y^{-1/2}, which we can find by doing
        # the square root and the inverse separately. As a possible improvement,
        # iterative methods exist to find inverse square root matrices directly,
        # as this is important in statistics.
        Z = np.kron(
            np.eye(N),
            sqrtm(la.inv(Y))
        )

        # Finally, we dot everything together and pack it into a Qobj,
        # marking the dimensions as that of a type=super (that is,
        # with left and right compound indices, each representing
        # left and right indices on the underlying Hilbert space).
        D = Qobj(np.dot(Z, np.dot(XXdag, Z)))
    else:
        D = N * Qobj(XXdag / np.trace(XXdag))

    D.dims = [
        # Left dims
        [[N], [N]],
        # Right dims
        [[N], [N]]
    ]

    # Since [BCSZ08] gives a row-stacking Choi matrix, but QuTiP
    # expects a column-stacking Choi matrix, we must permute the indices.
    D = D.permute([[1], [0]])

    D.dims = dims

    # Mark that we've made a Choi matrix.
    D.superrep = 'choi'

    return sr.to_super(D)

Example 38

Project: pylon Source File: pips.py
def pips(f_fcn, x0, A=None, l=None, u=None, xmin=None, xmax=None,
         gh_fcn=None, hess_fcn=None, opt=None):
    """Primal-dual interior point method for NLP (non-linear programming).
    Minimize a function F(X) beginning from a starting point M{x0}, subject to
    optional linear and non-linear constraints and variable bounds::

            min f(x)
             x

    subject to::

            g(x) = 0            (non-linear equalities)
            h(x) <= 0           (non-linear inequalities)
            l <= A*x <= u       (linear constraints)
            xmin <= x <= xmax   (variable bounds)

    Note: The calling syntax is almost identical to that of FMINCON from
    MathWorks' Optimization Toolbox. The main difference is that the linear
    constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
    C{Aeq}, C{Beq}. The functions for evaluating the objective function,
    constraints and Hessian are identical.

    Example from U{http://en.wikipedia.org/wiki/Nonlinear_programming}:
        >>> from numpy import array, r_, float64, dot
        >>> from scipy.sparse import csr_matrix
        >>> def f2(x):
        ...     f = -x[0] * x[1] - x[1] * x[2]
        ...     df = -r_[x[1], x[0] + x[2], x[1]]
        ...     # actually not used since 'hess_fcn' is provided
        ...     d2f = -array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], float64)
        ...     return f, df, d2f
        >>> def gh2(x):
        ...     h = dot(array([[1, -1, 1],
        ...                    [1,  1, 1]]), x**2) + array([-2.0, -10.0])
        ...     dh = 2 * csr_matrix(array([[ x[0], x[0]],
        ...                                [-x[1], x[1]],
        ...                                [ x[2], x[2]]]))
        ...     g = array([])
        ...     dg = None
        ...     return h, g, dh, dg
        >>> def hess2(x, lam):
        ...     mu = lam["ineqnonlin"]
        ...     a = r_[dot(2 * array([1, 1]), mu), -1, 0]
        ...     b = r_[-1, dot(2 * array([-1, 1]),mu),-1]
        ...     c = r_[0, -1, dot(2 * array([1, 1]),mu)]
        ...     Lxx = csr_matrix(array([a, b, c]))
        ...     return Lxx
        >>> x0 = array([1, 1, 0], float64)
        >>> solution = pips(f2, x0, gh_fcn=gh2, hess_fcn=hess2)
        >>> round(solution["f"], 11) == -7.07106725919
        True
        >>> solution["output"]["iterations"]
        8

    Ported by Richard Lincoln from the MATLAB Interior Point Solver (MIPS)
    (v1.9) by Ray Zimmerman.  MIPS is distributed as part of the MATPOWER
    project, developed at the Power System Engineering Research Center (PSERC),
    Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
    MIPS was ported by Ray Zimmerman from C code written by H. Wang for his
    PhD dissertation:
      - "On the Computation and Application of Multi-period
        Security-Constrained Optimal Power Flow for Real-time
        Electricity Market Operations", Cornell University, May 2007.

    See also:
      - H. Wang, C. E. Murillo-Sanchez, R. D. Zimmerman, R. J. Thomas,
        "On Computational Issues of Market-Based Optimal Power Flow",
        IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
        pp. 1185-1193.

    All parameters are optional except C{f_fcn} and C{x0}.
    @param f_fcn: Function that evaluates the objective function, its gradients
                  and Hessian for a given value of M{x}. If there are
                  non-linear constraints, the Hessian information is provided
                  by the 'hess_fcn' argument and is not required here.
    @type f_fcn: callable
    @param x0: Starting value of optimization vector M{x}.
    @type x0: array
    @param A: Optional linear constraints.
    @type A: csr_matrix
    @param l: Optional linear constraints. Default values are M{-Inf}.
    @type l: array
    @param u: Optional linear constraints. Default values are M{Inf}.
    @type u: array
    @param xmin: Optional lower bounds on the M{x} variables, defaults are
                 M{-Inf}.
    @type xmin: array
    @param xmax: Optional upper bounds on the M{x} variables, defaults are
                 M{Inf}.
    @type xmax: array
    @param gh_fcn: Function that evaluates the optional non-linear constraints
                   and their gradients for a given value of M{x}.
    @type gh_fcn: callable
    @param hess_fcn: Handle to function that computes the Hessian of the
                     Lagrangian for given values of M{x}, M{lambda} and M{mu},
                     where M{lambda} and M{mu} are the multipliers on the
                     equality and inequality constraints, M{g} and M{h},
                     respectively.
    @type hess_fcn: callable
    @param opt: optional options dictionary with the following keys, all of
                which are also optional (default values shown in parentheses)
                  - C{verbose} (False) - Controls level of progress output
                    displayed
                  - C{feastol} (1e-6) - termination tolerance for feasibility
                    condition
                  - C{gradtol} (1e-6) - termination tolerance for gradient
                    condition
                  - C{comptol} (1e-6) - termination tolerance for
                    complementarity condition
                  - C{costtol} (1e-6) - termination tolerance for cost
                    condition
                  - C{max_it} (150) - maximum number of iterations
                  - C{step_control} (False) - set to True to enable step-size
                    control
                  - C{max_red} (20) - maximum number of step-size reductions if
                    step-control is on
                  - C{cost_mult} (1.0) - cost multiplier used to scale the
                    objective function for improved conditioning. Note: The
                    same value must also be passed to the Hessian evaluation
                    function so that it can appropriately scale the objective
                    function term in the Hessian of the Lagrangian.
    @type opt: dict

    @rtype: dict
    @return: The solution dictionary has the following keys:
               - C{x} - solution vector
               - C{f} - final objective function value
               - C{converged} - exit status
                   - True = first order optimality conditions satisfied
                   - False = maximum number of iterations reached
                   - None = numerically failed
               - C{output} - output dictionary with keys:
                   - C{iterations} - number of iterations performed
                   - C{hist} - dictionary of arrays with trajectories of the
                     following: feascond, gradcond, compcond, costcond, gamma,
                     stepsize, obj, alphap, alphad
                   - C{message} - exit message
               - C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
                 multipliers on the constraints, with keys:
                   - C{eqnonlin} - non-linear equality constraints
                   - C{ineqnonlin} - non-linear inequality constraints
                   - C{mu_l} - lower (left-hand) limit on linear constraints
                   - C{mu_u} - upper (right-hand) limit on linear constraints
                   - C{lower} - lower bound on optimization variables
                   - C{upper} - upper bound on optimization variables

    @license: Apache License version 2.0
    """
    nx = x0.shape[0]                        # number of variables
    nA = A.shape[0] if A is not None else 0 # number of original linear constr

    # default argument values
#    l = array([]) if A is None else l
#    u = array([]) if A is None else u
    l = -Inf * ones(nA) if l is None else l
    u =  Inf * ones(nA) if u is None else u
    xmin = -Inf * ones(x0.shape[0]) if xmin is None else xmin
    xmax =  Inf * ones(x0.shape[0]) if xmax is None else xmax
    if gh_fcn is None:
        nonlinear = False
        gn = array([])
        hn = array([])
    else:
        nonlinear = True

    opt = {} if opt is None else opt
    # options
    if not opt.has_key("feastol"):
        opt["feastol"] = 1e-06
    if not opt.has_key("gradtol"):
        opt["gradtol"] = 1e-06
    if not opt.has_key("comptol"):
        opt["comptol"] = 1e-06
    if not opt.has_key("costtol"):
        opt["costtol"] = 1e-06
    if not opt.has_key("max_it"):
        opt["max_it"] = 150
    if not opt.has_key("max_red"):
        opt["max_red"] = 20
    if not opt.has_key("step_control"):
        opt["step_control"] = False
    if not opt.has_key("cost_mult"):
        opt["cost_mult"] = 1
    if not opt.has_key("verbose"):
        opt["verbose"] = False

    # initialize history
    hist = {}

    # constants
    xi = 0.99995
    sigma = 0.1
    z0 = 1
    alpha_min = 1e-8
#    rho_min = 0.95
#    rho_max = 1.05
    mu_threshold = 1e-5

    # initialize
    i = 0                       # iteration counter
    converged = False           # flag
    eflag = False               # exit flag

    # add var limits to linear constraints
    eyex = eye(nx, nx, format="csr")
    AA = eyex if A is None else vstack([eyex, A], "csr")
    ll = r_[xmin, l]
    uu = r_[xmax, u]

    # split up linear constraints
    ieq = flatnonzero( absolute(uu - ll) <= EPS )
    igt = flatnonzero( (uu >=  1e10) & (ll > -1e10) )
    ilt = flatnonzero( (ll <= -1e10) & (uu <  1e10) )
    ibx = flatnonzero( (absolute(uu - ll) > EPS) & (uu < 1e10) & (ll > -1e10) )
    # zero-sized sparse matrices unsupported
    Ae = AA[ieq, :] if len(ieq) else None
    if len(ilt) or len(igt) or len(ibx):
        idxs = [(1, ilt), (-1, igt), (1, ibx), (-1, ibx)]
        Ai = vstack([sig * AA[idx, :] for sig, idx in idxs if len(idx)])
    else:
        Ai = None
    be = uu[ieq, :]
    bi = r_[uu[ilt], -ll[igt], uu[ibx], -ll[ibx]]

    # evaluate cost f(x0) and constraints g(x0), h(x0)
    x = x0
    f, df, _ = f_fcn(x)                 # cost
    f = f * opt["cost_mult"]
    df = df * opt["cost_mult"]
    if nonlinear:
        hn, gn, dhn, dgn = gh_fcn(x)        # non-linear constraints
        h = hn if Ai is None else r_[hn, Ai * x - bi] # inequality constraints
        g = gn if Ae is None else r_[gn, Ae * x - be] # equality constraints

        if (dhn is None) and (Ai is None):
            dh = None
        elif dhn is None:
            dh = Ai.T
        elif Ae is None:
            dh = dhn
        else:
            dh = hstack([dhn, Ai.T])

        if (dgn is None) and (Ae is None):
            dg = None
        elif dgn is None:
            dg = Ae.T
        elif Ae is None:
            dg = dgn
        else:
            dg = hstack([dgn, Ae.T])
    else:
        h = -bi if Ai is None else Ai * x - bi        # inequality constraints
        g = -be if Ae is None else Ae * x - be        # equality constraints
        dh = None if Ai is None else Ai.T     # 1st derivative of inequalities
        dg = None if Ae is None else Ae.T     # 1st derivative of equalities

    # some dimensions
    neq = g.shape[0]           # number of equality constraints
    niq = h.shape[0]           # number of inequality constraints
    neqnln = gn.shape[0]       # number of non-linear equality constraints
    niqnln = hn.shape[0]       # number of non-linear inequality constraints
    nlt = len(ilt)             # number of upper bounded linear inequalities
    ngt = len(igt)             # number of lower bounded linear inequalities
    nbx = len(ibx)             # number of doubly bounded linear inequalities

    # initialize gamma, lam, mu, z, e
    gamma = 1                  # barrier coefficient
    lam = zeros(neq)
    z = z0 * ones(niq)
    mu = z0 * ones(niq)
    k = flatnonzero(h < -z0)
    z[k] = -h[k]
    k = flatnonzero((gamma / z) > z0)
    mu[k] = gamma / z[k]
    e = ones(niq)

    # check tolerance
    f0 = f
#    if opt["step_control"]:
#        L = f + lam.T * g + mu.T * (h + z) - gamma * sum(log(z))

    Lx = df
    Lx = Lx + dg * lam if dg is not None else Lx
    Lx = Lx + dh * mu  if dh is not None else Lx

    gnorm = norm(g, Inf) if len(g) else 0.0
    lam_norm = norm(lam, Inf) if len(lam) else 0.0
    mu_norm = norm(mu, Inf) if len(mu) else 0.0
    feascond = \
        max([gnorm, max(h)]) / (1 + max([norm(x, Inf), norm(z, Inf)]))
    gradcond = \
        norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
    compcond = dot(z, mu) / (1 + norm(x, Inf))
    costcond = absolute(f - f0) / (1 + absolute(f0))

    # save history
    hist[i] = {'feascond': feascond, 'gradcond': gradcond,
        'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
        'stepsize': 0, 'obj': f / opt["cost_mult"], 'alphap': 0, 'alphad': 0}

    if opt["verbose"]:
#        s = '-sc' if opt["step_control"] else ''
#        version, date = '1.0b2', '24-Mar-2010'
#        print 'Python Interior Point Solver - PIPS%s, Version %s, %s' % \
#                    (s, version, date)
        print " it    objective   step size   feascond     gradcond     " \
              "compcond     costcond  "
        print "----  ------------ --------- ------------ ------------ " \
              "------------ ------------"
        print "%3d  %12.8g %10s %12g %12g %12g %12g" % \
            (i, (f / opt["cost_mult"]), "",
             feascond, gradcond, compcond, costcond)

    if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
        compcond < opt["comptol"] and costcond < opt["costtol"]:
        converged = True
        if opt["verbose"]:
            print "Converged!"

    # do Newton iterations
    while (not converged and i < opt["max_it"]):
        # update iteration counter
        i += 1

        # compute update step
        lmbda = {"eqnonlin": lam[range(neqnln)],
                 "ineqnonlin": mu[range(niqnln)]}
        if nonlinear:
            if hess_fcn is None:
                print "pips: Hessian evaluation via finite differences " \
                      "not yet implemented.\nPlease provide " \
                      "your own hessian evaluation function."
            Lxx = hess_fcn(x, lmbda)
        else:
            _, _, d2f = f_fcn(x)      # cost
            Lxx = d2f * opt["cost_mult"]
        rz = range(len(z))
        zinvdiag = csr_matrix((1.0 / z, (rz, rz))) if len(z) else None
        rmu = range(len(mu))
        mudiag = csr_matrix((mu, (rmu, rmu))) if len(mu) else None
        dh_zinv = None if dh is None else dh * zinvdiag
        M = Lxx if dh is None else Lxx + dh_zinv * mudiag * dh.T
        N = Lx if dh is None else Lx + dh_zinv * (mudiag * h + gamma * e)

        Ab = M if dg is None else vstack([
            hstack([M, dg]),
            hstack([dg.T, csr_matrix((neq, neq))])
        ])
        bb = r_[-N, -g]

        dxdlam = spsolve(Ab.tocsr(), bb)

        dx = dxdlam[:nx]
        dlam = dxdlam[nx:nx + neq]
        dz = -h - z if dh is None else -h - z - dh.T * dx
        dmu = -mu if dh is None else -mu + zinvdiag * (gamma * e - mudiag * dz)

        # optional step-size control
#        sc = False
        if opt["step_control"]:
            raise NotImplementedError
#            x1 = x + dx
#
#            # evaluate cost, constraints, derivatives at x1
#            f1, df1 = ipm_f(x1)          # cost
#            f1 = f1 * opt["cost_mult"]
#            df1 = df1 * opt["cost_mult"]
#            gn1, hn1, dgn1, dhn1 = ipm_gh(x1) # non-linear constraints
#            g1 = gn1 if Ai is None else r_[gn1, Ai * x1 - bi] # ieq constraints
#            h1 = hn1 if Ae is None else r_[hn1, Ae * x1 - be] # eq constraints
#            dg1 = dgn1 if Ai is None else r_[dgn1, Ai.T]      # 1st der of ieq
#            dh1 = dhn1 if Ae is None else r_[dhn1, Ae.T]      # 1st der of eqs
#
#            # check tolerance
#            Lx1 = df1 + dh1 * lam + dg1 * mu
#            feascond1 = max([ norm(h1, Inf), max(g1) ]) / \
#                (1 + max([ norm(x1, Inf), norm(z, Inf) ]))
#            gradcond1 = norm(Lx1, Inf) / \
#                (1 + max([ norm(lam, Inf), norm(mu, Inf) ]))
#
#            if feascond1 > feascond and gradcond1 > gradcond:
#                sc = True
#        if sc:
#            alpha = 1.0
#            for j in range(opt["max_red"]):
#                dx1 = alpha * dx
#                x1 = x + dx1
#                f1 = ipm_f(x1)             # cost
#                f1 = f1 * opt["cost_mult"]
#                gn1, hn1 = ipm_gh(x1)              # non-linear constraints
#                g1 = r_[gn1, Ai * x1 - bi]         # inequality constraints
#                h1 = r_[hn1, Ae * x1 - be]         # equality constraints
#                L1 = f1 + lam.H * h1 + mu.H * (g1 + z) - gamma * sum(log(z))
#                if opt["verbose"]:
#                    logger.info("\n   %3d            %10.f" % (-j, norm(dx1)))
#                rho = (L1 - L) / (Lx.H * dx1 + 0.5 * dx1.H * Lxx * dx1)
#                if rho > rho_min and rho < rho_max:
#                    break
#                else:
#                    alpha = alpha / 2.0
#            dx = alpha * dx
#            dz = alpha * dz
#            dlam = alpha * dlam
#            dmu = alpha * dmu

        # do the update
        k = flatnonzero(dz < 0.0)
        alphap = min([xi * min(z[k] / -dz[k]), 1]) if len(k) else 1.0
        k = flatnonzero(dmu < 0.0)
        alphad = min([xi * min(mu[k] / -dmu[k]), 1]) if len(k) else 1.0
        x = x + alphap * dx
        z = z + alphap * dz
        lam = lam + alphad * dlam
        mu = mu + alphad * dmu
        if niq > 0:
            gamma = sigma * dot(z, mu) / niq

        # evaluate cost, constraints, derivatives
        f, df, _ = f_fcn(x)             # cost
        f = f * opt["cost_mult"]
        df = df * opt["cost_mult"]
        if nonlinear:
            hn, gn, dhn, dgn = gh_fcn(x)                   # nln constraints
#            g = gn if Ai is None else r_[gn, Ai * x - bi] # ieq constraints
#            h = hn if Ae is None else r_[hn, Ae * x - be] # eq constraints
            h = hn if Ai is None else r_[hn, Ai * x - bi] # ieq constr
            g = gn if Ae is None else r_[gn, Ae * x - be]  # eq constr

            if (dhn is None) and (Ai is None):
                dh = None
            elif dhn is None:
                dh = Ai.T
            elif Ae is None:
                dh = dhn
            else:
                dh = hstack([dhn, Ai.T])

            if (dgn is None) and (Ae is None):
                dg = None
            elif dgn is None:
                dg = Ae.T
            elif Ae is None:
                dg = dgn
            else:
                dg = hstack([dgn, Ae.T])
        else:
            h = -bi if Ai is None else Ai * x - bi    # inequality constraints
            g = -be if Ae is None else Ae * x - be    # equality constraints
            # 1st derivatives are constant, still dh = Ai.T, dg = Ae.T

        Lx = df
        Lx = Lx + dg * lam if dg is not None else Lx
        Lx = Lx + dh * mu  if dh is not None else Lx

        gnorm = norm(g, Inf) if len(g) else 0.0
        lam_norm = norm(lam, Inf) if len(lam) else 0.0
        mu_norm = norm(mu, Inf) if len(mu) else 0.0
        feascond = \
            max([gnorm, max(h)]) / (1+max([norm(x, Inf), norm(z, Inf)]))
        gradcond = \
            norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
        compcond = dot(z, mu) / (1 + norm(x, Inf))
        costcond = float(absolute(f - f0) / (1 + absolute(f0)))

        hist[i] = {'feascond': feascond, 'gradcond': gradcond,
            'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
            'stepsize': norm(dx), 'obj': f / opt["cost_mult"],
            'alphap': alphap, 'alphad': alphad}

        if opt["verbose"]:
            print "%3d  %12.8g %10.5g %12g %12g %12g %12g" % \
                (i, (f / opt["cost_mult"]), norm(dx), feascond, gradcond,
                 compcond, costcond)

        if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
            compcond < opt["comptol"] and costcond < opt["costtol"]:
            converged = True
            if opt["verbose"]:
                print "Converged!"
        else:
            if any(isnan(x)) or (alphap < alpha_min) or \
                (alphad < alpha_min) or (gamma < EPS) or (gamma > 1.0 / EPS):
                if opt["verbose"]:
                    print "Numerically failed."
                eflag = -1
                break
            f0 = f

#            if opt["step_control"]:
#                L = f + dot(lam, g) + dot(mu * (h + z)) - gamma * sum(log(z))

    if opt["verbose"]:
        if not converged:
            print "Did not converge in %d iterations." % i

    # package results
    if eflag != -1:
        eflag = converged

    if eflag == 0:
        message = 'Did not converge'
    elif eflag == 1:
        message = 'Converged'
    elif eflag == -1:
        message = 'Numerically failed'
    else:
        raise

    output = {"iterations": i, "history": hist, "message": message}

    # zero out multipliers on non-binding constraints
    mu[flatnonzero( (h < -opt["feastol"]) & (mu < mu_threshold) )] = 0.0

    # un-scale cost and prices
    f = f / opt["cost_mult"]
    lam = lam / opt["cost_mult"]
    mu = mu / opt["cost_mult"]

    # re-package multipliers into struct
    lam_lin = lam[neqnln:neq]           # lambda for linear constraints
    mu_lin = mu[niqnln:niq]             # mu for linear constraints
    kl = flatnonzero(lam_lin < 0.0)     # lower bound binding
    ku = flatnonzero(lam_lin > 0.0)     # upper bound binding

    mu_l = zeros(nx + nA)
    mu_l[ieq[kl]] = -lam_lin[kl]
    mu_l[igt] = mu_lin[nlt:nlt + ngt]
    mu_l[ibx] = mu_lin[nlt + ngt + nbx:nlt + ngt + nbx + nbx]

    mu_u = zeros(nx + nA)
    mu_u[ieq[ku]] = lam_lin[ku]
    mu_u[ilt] = mu_lin[:nlt]
    mu_u[ibx] = mu_lin[nlt + ngt:nlt + ngt + nbx]

    lmbda = {'mu_l': mu_l[nx:], 'mu_u': mu_u[nx:],
             'lower': mu_l[:nx], 'upper': mu_u[:nx]}

    if niqnln > 0:
        lmbda['ineqnonlin'] = mu[:niqnln]
    if neqnln > 0:
        lmbda['eqnonlin'] = lam[:neqnln]

#    lmbda = {"eqnonlin": lam[:neqnln], 'ineqnonlin': mu[:niqnln],
#             "mu_l": mu_l[nx:], "mu_u": mu_u[nx:],
#             "lower": mu_l[:nx], "upper": mu_u[:nx]}

    solution =  {"x": x, "f": f, "converged": converged,
                 "lmbda": lmbda, "output": output}

    return solution

Example 39

Project: pylon Source File: pips.py
Function: qps_pips
def qps_pips(H, c, A, l, u, xmin=None, xmax=None, x0=None, opt=None):
    """Uses the Python Interior Point Solver (PIPS) to solve the following
    QP (quadratic programming) problem::

            min 1/2 x'*H*x + C'*x
             x

    subject to::

            l <= A*x <= u       (linear constraints)
            xmin <= x <= xmax   (variable bounds)

    Note the calling syntax is almost identical to that of QUADPROG from
    MathWorks' Optimization Toolbox. The main difference is that the linear
    constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
    C{Aeq}, C{Beq}.

    See also L{pips}.

    Example from U{http://www.uc.edu/sashtml/iml/chap8/sect12.htm}:

        >>> from numpy import array, zeros, Inf
        >>> from scipy.sparse import csr_matrix
        >>> H = csr_matrix(array([[1003.1,  4.3,     6.3,     5.9],
        ...                       [4.3,     2.2,     2.1,     3.9],
        ...                       [6.3,     2.1,     3.5,     4.8],
        ...                       [5.9,     3.9,     4.8,     10 ]]))
        >>> c = zeros(4)
        >>> A = csr_matrix(array([[1,       1,       1,       1   ],
        ...                       [0.17,    0.11,    0.10,    0.18]]))
        >>> l = array([1, 0.10])
        >>> u = array([1, Inf])
        >>> xmin = zeros(4)
        >>> xmax = None
        >>> x0 = array([1, 0, 0, 1])
        >>> solution = qps_pips(H, c, A, l, u, xmin, xmax, x0)
        >>> round(solution["f"], 11) == 1.09666678128
        True
        >>> solution["converged"]
        True
        >>> solution["output"]["iterations"]
        10

    All parameters are optional except C{H}, C{C}, C{A} and C{L}.
    @param H: Quadratic cost coefficients.
    @type H: csr_matrix
    @param c: vector of linear cost coefficients
    @type c: array
    @param A: Optional linear constraints.
    @type A: csr_matrix
    @param l: Optional linear constraints. Default values are M{-Inf}.
    @type l: array
    @param u: Optional linear constraints. Default values are M{Inf}.
    @type u: array
    @param xmin: Optional lower bounds on the M{x} variables, defaults are
                 M{-Inf}.
    @type xmin: array
    @param xmax: Optional upper bounds on the M{x} variables, defaults are
                 M{Inf}.
    @type xmax: array
    @param x0: Starting value of optimization vector M{x}.
    @type x0: array
    @param opt: optional options dictionary with the following keys, all of
                which are also optional (default values shown in parentheses)
                  - C{verbose} (False) - Controls level of progress output
                    displayed
                  - C{feastol} (1e-6) - termination tolerance for feasibility
                    condition
                  - C{gradtol} (1e-6) - termination tolerance for gradient
                    condition
                  - C{comptol} (1e-6) - termination tolerance for
                    complementarity condition
                  - C{costtol} (1e-6) - termination tolerance for cost
                    condition
                  - C{max_it} (150) - maximum number of iterations
                  - C{step_control} (False) - set to True to enable step-size
                    control
                  - C{max_red} (20) - maximum number of step-size reductions if
                    step-control is on
                  - C{cost_mult} (1.0) - cost multiplier used to scale the
                    objective function for improved conditioning. Note: The
                    same value must also be passed to the Hessian evaluation
                    function so that it can appropriately scale the objective
                    function term in the Hessian of the Lagrangian.
    @type opt: dict

    @rtype: dict
    @return: The solution dictionary has the following keys:
               - C{x} - solution vector
               - C{f} - final objective function value
               - C{converged} - exit status
                   - True = first order optimality conditions satisfied
                   - False = maximum number of iterations reached
                   - None = numerically failed
               - C{output} - output dictionary with keys:
                   - C{iterations} - number of iterations performed
                   - C{hist} - dictionary of arrays with trajectories of the
                     following: feascond, gradcond, compcond, costcond, gamma,
                     stepsize, obj, alphap, alphad
                   - C{message} - exit message
               - C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
                 multipliers on the constraints, with keys:
                   - C{eqnonlin} - non-linear equality constraints
                   - C{ineqnonlin} - non-linear inequality constraints
                   - C{mu_l} - lower (left-hand) limit on linear constraints
                   - C{mu_u} - upper (right-hand) limit on linear constraints
                   - C{lower} - lower bound on optimization variables
                   - C{upper} - upper bound on optimization variables

    @license: Apache License version 2.0
    """
    if H is None or H.nnz == 0:
        if A is None or A.nnz == 0 and \
           xmin is None or len(xmin) == 0 and \
           xmax is None or len(xmax) == 0:
            print 'qps_pips: LP problem must include constraints or variable bounds'
            return
        else:
            if A is not None and A.nnz >= 0:
                nx = A.shape[1]
            elif xmin is not None and len(xmin) > 0:
                nx = xmin.shape[0]
            elif xmax is not None and len(xmax) > 0:
                nx = xmax.shape[0]
        H = csr_matrix((nx, nx))
    else:
        nx = H.shape[0]

    xmin = -Inf * ones(nx) if xmin is None else xmin
    xmax =  Inf * ones(nx) if xmax is None else xmax

    c = zeros(nx) if c is None else c

#    if x0 is None:
#        x0 = zeros(nx)
#        k = flatnonzero( (VUB < 1e10) & (VLB > -1e10) )
#        x0[k] = ((VUB[k] + VLB[k]) / 2)
#        k = flatnonzero( (VUB < 1e10) & (VLB <= -1e10) )
#        x0[k] = VUB[k] - 1
#        k = flatnonzero( (VUB >= 1e10) & (VLB > -1e10) )
#        x0[k] = VLB[k] + 1

    x0 = zeros(nx) if x0 is None else x0

    opt = {} if opt is None else opt
    if not opt.has_key("cost_mult"):
        opt["cost_mult"] = 1

    def qp_f(x):
        f = 0.5 * dot(x.T * H, x) + dot(c.T, x)
        df = H * x + c
        d2f = H
        return f, df, d2f

#    def qp_gh(x):
#        g = array([])
#        h = array([])
#        dg = None
#        dh = None
#        return g, h, dg, dh
#
#    def qp_hessian(x, lmbda):
#        Lxx = H * opt["cost_mult"]
#        return Lxx

#    l = -Inf * ones(b.shape[0])
#    l[:N] = b[:N]

    return pips(qp_f, x0, A, l, u, xmin, xmax, opt=opt)

Example 40

Project: pyspace Source File: SOR.py
Function: iteration_loop
    def iteration_loop(self, M, reduced_indices=[]):
        """ The algorithm is calling the :func:`reduced_descent<pySPACE.missions.nodes.classifiers.ada_SVM.SORSVMNode.reduced_descent>` method in loops over alpha

        In the first step it uses a complete loop over all components of alpha
        and in the second inner loop only the non zero alpha are observed till
        come convergence criterion is reached.

        *reduced_indices* will be skipped in observation.
        """
        ## Definition of tracking variables ##
        self.iterations = 0
        self.difference = numpy.inf
        ## outer iteration loop ##
        while (self.difference > self.tolerance and
               self.iterations <= self.max_iterations):
            # inner iteration loop only on active vectors/alpha (non zero) ##
            self.sub_iterations = 0
            # sorting or randomizing non zero indices
            # arrays are mapped to lists for later iteration
            sort_dual = self.dual_solution

            num_non_zeros = len(map(list,sort_dual.nonzero())[0])
            max_values = len(map(list,
                                 numpy.where(sort_dual == sort_dual.max()))[0])
            # sort the entries of the current dual
            # and get the corresponding indices
            sorted_indices = map(list,[numpy.argsort(sort_dual)])[0]
            if num_non_zeros == 0 or num_non_zeros==max_values:
                # skip sub iteration if everything is zero or maximal
                non_zero_indices = []
            else:
                non_zero_indices = sorted_indices[-num_non_zeros:-max_values]
            for index in reduced_indices:
                try:
                    non_zero_indices.remove(index)
                except ValueError:
                    pass
            if self.random:
                random.shuffle(non_zero_indices)
            self.max_sub_iterations = self.max_iterations_factor * \
                len(non_zero_indices) * 0.5
            while (self.difference > self.tolerance and
                   self.sub_iterations < self.max_sub_iterations
                   and self.iterations < self.max_iterations):
                ## iteration step ##
                self.reduced_descent(self.dual_solution, M, non_zero_indices)
                ## outer loop ##
            if not (self.iterations < self.max_iterations):
                break
            # For the first run, the previous reduced descent is skipped
            # but for retraining it is important
            # to have first the small loop, since normally, this is sufficient.
            # Furthermore having it at the end simplifies the stop criterion
            self.max_sub_iterations = numpy.inf
            self.total_descent(self.dual_solution, M, reduced_indices)
            ## Final solution ##
        # in the case without kernels, we have to calculate the result
        # by hand new for each incoming sample
        if self.version == "matrix":
            self.b = self.offset_factor * dot(self.dual_solution, self.bi)
            # self.w = self.samples[0]*self.dual_solution[0]*self.bi[0]
            # for i in range(self.num_samples-1):
            #     self.w = self.w + self.bi[i+1] * self.samples[i+1] *
            #         self.dual_solution[i+1]
            if self.kernel_type == "LINEAR":
                self.w = numpy.array([dot(dot(self.A.T, self.D),
                                          self.dual_solution)]).T
        elif self.version == "samples" and self.kernel_type == "LINEAR":
            # w and b are pre-computed in the loop
            # transferring of 1-d array to 2d array
            # self.w = numpy.array([self.w]).T
            pass

Example 41

Project: pyoptools Source File: calc.py
def pupil_location(opsys,ccds,opaxis):
    '''
    Function to find the optical system pupils position
    
    .. note:
        For this function to operate, the system should have a rotational
        symmetry around the optical axis. 
       
    **Parameters:**
        
        opsys   Optical system to use.
        opaxis  Ray representing the optical axis
        ccds    Surface that represents a detector in the aperture plane
    
    **Return Value**
        
        (enpl,expl)
        
        enpl tuple (xen,yen,zen) containing the entrance pupil coordinates
        expl tuple (xex,yex,zex) containing the exit pupil coordinates
        
    
    '''
    
    #log.info("Propagate Optical axis ray")
    opsys.clear_ray_list()
    opsys.reset()
    #opsys.ray_add(cray)
    opsys.ray_add(opaxis)

    opsys.propagate()
    
        
    if (len(ccds.hit_list)==0):
        raise Exception("The optical axis did not intersect the aperture") 
    if(len(ccds.hit_list)>1):
        raise Exception("The optical axis intersected the aperture more than once") 
        
    aip=ccds.hit_list[0][0]
    air=ccds.hit_list[0][1]
    
    #log.info("Optical Axis Intersection point= "+str(aip))
    #log.info("Intersection Ray= "+str(air))
    
    #Getting Intersection point in global coordinates
    
    if(len(air.childs)!=1):
        raise Exception("The intersected ray can only have one child")
    
    ip=air.childs[0].pos
    d=air.childs[0].dir
    #log.info("Intersection point in world coordinates= "+str(ip))
    #log.info("Direction of the optical axis at the intersection point"+str(d))
    
    #Todo: Check if the optical axis and the aperture are perpendicular
    
    # Calculate vectors perpendicular to the optical axis and to the XYZ axes
    pv1= cross(d,(0,0,1))
    pv2= cross(d,(0,1,0))
    pv3= cross(d,(1,0,0))
    
    pv=[pv1,pv2,pv3]
    
    # Search for the longest pv
    pvn=array((dot(pv1,pv1),dot(pv2,pv2),dot(pv3,pv3)))
    
    pvm=pv[pvn.argmax()]
    
    #log.info("Displacement vector found: "+str(pvm))

    # Create ray to calculate the exit pupil
    expuray=air.childs[0].copy()
    expuray.dir=expuray.dir+pvm*.0001
    
    # Create the ray to calculate the entrance pupil
    enpuray=expuray.reverse()
    
    opsys.clear_ray_list()
    opsys.reset()
    opsys.ray_add(enpuray)
    opsys.ray_add(expuray)

    opsys.propagate()
    
    enp=enpuray.get_final_rays(inc_zeros = False)
    exp=expuray.get_final_rays(inc_zeros = False)
    oax=opaxis.get_final_rays(inc_zeros = False)
    #log.info("enp="+str(enp))
    #log.info("exp="+str(exp))
    #log.info("oax="+str(oax))
    if len(enp)!=1 or len(exp)!=1 or len(oax)!=1: 
        raise Exception("The principal ray or the optical axis ray have more"
        " than one final ray")
    #log.info("Calculating entrance pupil location")
    
    # Find the nearest points between the rays. 
    # Some times because of numerical errors, or some aberrations in the optical
    # system, the rays do not trully intersect.
    # Use instead the nearest points and issue a warning when the rays do not trully
    # intersect.
    
    enpl=intersection(opaxis,enp[0])[0]
    if (isnan(enpl)).all():
        p1, p2, d, rv =nearest_points(opaxis,enp[0])
        print("Warning: The optical axis does not intersect the principal ray at the entrance")
        print("pupil. The minimum distance is:", d)
        enpl=(p1+p2)/2
        
    #log.info("Calculating exit pupil location")
    expl=intersection(oax[0],exp[0])[0]
    if (isnan(expl)).all():
        p1, p2, d, rv =nearest_points(oax[0],exp[0])
        print("Warning: The optical axis does not intersect the principal ray at the exit")
        print("pupil. The minimum distance is:", d)
        expl=(p1+p2)/2
    
    return enpl,expl

Example 42

Project: PYPOWER Source File: dcopf_solver.py
def dcopf_solver(om, ppopt, out_opt=None):
    """Solves a DC optimal power flow.

    Inputs are an OPF model object, a PYPOWER options dict and
    a dict containing fields (can be empty) for each of the desired
    optional output fields.

    Outputs are a C{results} dict, C{success} flag and C{raw} output dict.

    C{results} is a PYPOWER case dict (ppc) with the usual baseMVA, bus
    branch, gen, gencost fields, along with the following additional
    fields:
        - C{order}      see 'help ext2int' for details of this field
        - C{x}          final value of optimization variables (internal order)
        - C{f}          final objective function value
        - C{mu}         shadow prices on ...
            - C{var}
                - C{l}  lower bounds on variables
                - C{u}  upper bounds on variables
            - C{lin}
                - C{l}  lower bounds on linear constraints
                - C{u}  upper bounds on linear constraints
        - C{g}          (optional) constraint values
        - C{dg}         (optional) constraint 1st derivatives
        - C{df}         (optional) obj fun 1st derivatives (not yet implemented)
        - C{d2f}        (optional) obj fun 2nd derivatives (not yet implemented)

    C{success} is C{True} if solver converged successfully, C{False} otherwise.

    C{raw} is a raw output dict in form returned by MINOS
        - C{xr}     final value of optimization variables
        - C{pimul}  constraint multipliers
        - C{info}   solver specific termination code
        - C{output} solver specific output information

    @see: L{opf}, L{qps_pypower}

    @author: Ray Zimmerman (PSERC Cornell)
    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
    Autonoma de Manizales)
    """
    if out_opt is None:
        out_opt = {}

    ## options
    verbose = ppopt['VERBOSE']
    alg     = ppopt['OPF_ALG_DC']

    if alg == 0:
        if have_fcn('cplex'):        ## use CPLEX by default, if available
            alg = 500
        elif have_fcn('mosek'):      ## if not, then MOSEK, if available
            alg = 600
        elif have_fcn('gurobi'):     ## if not, then Gurobi, if available
            alg = 700
        else:                        ## otherwise PIPS
            alg = 200

    ## unpack data
    ppc = om.get_ppc()
    baseMVA, bus, gen, branch, gencost = \
        ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"], ppc["gencost"]
    cp = om.get_cost_params()
    N, H, Cw = cp["N"], cp["H"], cp["Cw"]
    fparm = array(c_[cp["dd"], cp["rh"], cp["kk"], cp["mm"]])
    Bf = om.userdata('Bf')
    Pfinj = om.userdata('Pfinj')
    vv, ll, _, _ = om.get_idx()

    ## problem dimensions
    ipol = find(gencost[:, MODEL] == POLYNOMIAL) ## polynomial costs
    ipwl = find(gencost[:, MODEL] == PW_LINEAR)  ## piece-wise linear costs
    nb = bus.shape[0]              ## number of buses
    nl = branch.shape[0]           ## number of branches
    nw = N.shape[0]                ## number of general cost vars, w
    ny = om.getN('var', 'y')       ## number of piece-wise linear costs
    nxyz = om.getN('var')          ## total number of control vars of all types

    ## linear constraints & variable bounds
    A, l, u = om.linear_constraints()
    x0, xmin, xmax = om.getv()

    ## set up objective function of the form: f = 1/2 * X'*HH*X + CC'*X
    ## where X = [x;y;z]. First set up as quadratic function of w,
    ## f = 1/2 * w'*HHw*w + CCw'*w, where w = diag(M) * (N*X - Rhat). We
    ## will be building on the (optionally present) user supplied parameters.

    ## piece-wise linear costs
    any_pwl = int(ny > 0)
    if any_pwl:
        # Sum of y vars.
        Npwl = sparse((ones(ny), (zeros(ny), arange(vv["i1"]["y"], vv["iN"]["y"]))), (1, nxyz))
        Hpwl = sparse((1, 1))
        Cpwl = array([1])
        fparm_pwl = array([[1, 0, 0, 1]])
    else:
        Npwl = None#zeros((0, nxyz))
        Hpwl = None#array([])
        Cpwl = array([])
        fparm_pwl = zeros((0, 4))

    ## quadratic costs
    npol = len(ipol)
    if any(find(gencost[ipol, NCOST] > 3)):
        stderr.write('DC opf cannot handle polynomial costs with higher '
                     'than quadratic order.\n')
    iqdr = find(gencost[ipol, NCOST] == 3)
    ilin = find(gencost[ipol, NCOST] == 2)
    polycf = zeros((npol, 3))         ## quadratic coeffs for Pg
    if len(iqdr) > 0:
        polycf[iqdr, :] = gencost[ipol[iqdr], COST:COST + 3]
    if npol:
        polycf[ilin, 1:3] = gencost[ipol[ilin], COST:COST + 2]
    polycf = dot(polycf, diag([ baseMVA**2, baseMVA, 1]))     ## convert to p.u.
    if npol:
        Npol = sparse((ones(npol), (arange(npol), vv["i1"]["Pg"] + ipol)),
                      (npol, nxyz))  # Pg vars
        Hpol = sparse((2 * polycf[:, 0], (arange(npol), arange(npol))),
                      (npol, npol))
    else:
        Npol = None
        Hpol = None
    Cpol = polycf[:, 1]
    fparm_pol = ones((npol, 1)) * array([[1, 0, 0, 1]])

    ## combine with user costs
    NN = vstack([n for n in [Npwl, Npol, N] if n is not None and n.shape[0] > 0], "csr")
    # FIXME: Zero dimension sparse matrices.
    if (Hpwl is not None) and any_pwl and (npol + nw):
        Hpwl = hstack([Hpwl, sparse((any_pwl, npol + nw))])
    if Hpol is not None:
        if any_pwl and npol:
            Hpol = hstack([sparse((npol, any_pwl)), Hpol])
        if npol and nw:
            Hpol = hstack([Hpol, sparse((npol, nw))])
    if (H is not None) and nw and (any_pwl + npol):
        H = hstack([sparse((nw, any_pwl + npol)), H])
    HHw = vstack([h for h in [Hpwl, Hpol, H] if h is not None and h.shape[0] > 0], "csr")
    CCw = r_[Cpwl, Cpol, Cw]
    ffparm = r_[fparm_pwl, fparm_pol, fparm]

    ## transform quadratic coefficients for w into coefficients for X
    nnw = any_pwl + npol + nw
    M = sparse((ffparm[:, 3], (range(nnw), range(nnw))))
    MR = M * ffparm[:, 1]
    HMR = HHw * MR
    MN = M * NN
    HH = MN.T * HHw * MN
    CC = MN.T * (CCw - HMR)
    C0 = 0.5 * dot(MR, HMR) + sum(polycf[:, 2])  # Constant term of cost.

    ## set up input for QP solver
    opt = {'alg': alg, 'verbose': verbose}
    if (alg == 200) or (alg == 250):
        ## try to select an interior initial point
        Varefs = bus[bus[:, BUS_TYPE] == REF, VA] * (pi / 180.0)

        lb, ub = xmin.copy(), xmax.copy()
        lb[xmin == -Inf] = -1e10   ## replace Inf with numerical proxies
        ub[xmax ==  Inf] =  1e10
        x0 = (lb + ub) / 2;
        # angles set to first reference angle
        x0[vv["i1"]["Va"]:vv["iN"]["Va"]] = Varefs[0]
        if ny > 0:
            ipwl = find(gencost[:, MODEL] == PW_LINEAR)
            # largest y-value in CCV data
            c = gencost.flatten('F')[sub2ind(gencost.shape, ipwl,
                                NCOST + 2 * gencost[ipwl, NCOST])]
            x0[vv["i1"]["y"]:vv["iN"]["y"]] = max(c) + 0.1 * abs(max(c))

        ## set up options
        feastol = ppopt['PDIPM_FEASTOL']
        gradtol = ppopt['PDIPM_GRADTOL']
        comptol = ppopt['PDIPM_COMPTOL']
        costtol = ppopt['PDIPM_COSTTOL']
        max_it  = ppopt['PDIPM_MAX_IT']
        max_red = ppopt['SCPDIPM_RED_IT']
        if feastol == 0:
            feastol = ppopt['OPF_VIOLATION']    ## = OPF_VIOLATION by default
        opt["pips_opt"] = {  'feastol': feastol,
                             'gradtol': gradtol,
                             'comptol': comptol,
                             'costtol': costtol,
                             'max_it':  max_it,
                             'max_red': max_red,
                             'cost_mult': 1  }
    elif alg == 400:
        opt['ipopt_opt'] = ipopt_options([], ppopt)
    elif alg == 500:
        opt['cplex_opt'] = cplex_options([], ppopt)
    elif alg == 600:
        opt['mosek_opt'] = mosek_options([], ppopt)
    elif alg == 700:
        opt['grb_opt'] = gurobi_options([], ppopt)
    else:
        raise ValueError("Unrecognised solver [%d]." % alg)

    ##-----  run opf  -----
    x, f, info, output, lmbda = \
            qps_pypower(HH, CC, A, l, u, xmin, xmax, x0, opt)
    success = (info == 1)

    ##-----  calculate return values  -----
    if not any(isnan(x)):
        ## update solution data
        Va = x[vv["i1"]["Va"]:vv["iN"]["Va"]]
        Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]
        f = f + C0

        ## update voltages & generator outputs
        bus[:, VA] = Va * 180 / pi
        gen[:, PG] = Pg * baseMVA

        ## compute branch flows
        branch[:, [QF, QT]] = zeros((nl, 2))
        branch[:, PF] = (Bf * Va + Pfinj) * baseMVA
        branch[:, PT] = -branch[:, PF]

    ## package up results
    mu_l = lmbda["mu_l"]
    mu_u = lmbda["mu_u"]
    muLB = lmbda["lower"]
    muUB = lmbda["upper"]

    ## update Lagrange multipliers
    il = find((branch[:, RATE_A] != 0) & (branch[:, RATE_A] < 1e10))
    bus[:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]] = zeros((nb, 4))
    gen[:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]] = zeros((gen.shape[0], 4))
    branch[:, [MU_SF, MU_ST]] = zeros((nl, 2))
    bus[:, LAM_P]       = (mu_u[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]] -
                           mu_l[ll["i1"]["Pmis"]:ll["iN"]["Pmis"]]) / baseMVA
    branch[il, MU_SF]   = mu_u[ll["i1"]["Pf"]:ll["iN"]["Pf"]] / baseMVA
    branch[il, MU_ST]   = mu_u[ll["i1"]["Pt"]:ll["iN"]["Pt"]] / baseMVA
    gen[:, MU_PMIN]     = muLB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA
    gen[:, MU_PMAX]     = muUB[vv["i1"]["Pg"]:vv["iN"]["Pg"]] / baseMVA

    pimul = r_[
      mu_l - mu_u,
     -ones((ny > 0)), ## dummy entry corresponding to linear cost row in A
      muLB - muUB
    ]

    mu = { 'var': {'l': muLB, 'u': muUB},
           'lin': {'l': mu_l, 'u': mu_u} }

    results = deepcopy(ppc)
    results["bus"], results["branch"], results["gen"], \
        results["om"], results["x"], results["mu"], results["f"] = \
            bus, branch, gen, om, x, mu, f

    raw = {'xr': x, 'pimul': pimul, 'info': info, 'output': output}

    return results, success, raw

Example 43

Project: PYPOWER Source File: opf_costfcn.py
def opf_costfcn(x, om, return_hessian=False):
    """Evaluates objective function, gradient and Hessian for OPF.

    Objective function evaluation routine for AC optimal power flow,
    suitable for use with L{pips}. Computes objective function value,
    gradient and Hessian.

    @param x: optimization vector
    @param om: OPF model object

    @return: C{F} - value of objective function. C{df} - (optional) gradient
    of objective function (column vector). C{d2f} - (optional) Hessian of
    objective function (sparse matrix).

    @see: L{opf_consfcn}, L{opf_hessfcn}

    @author: Carlos E. Murillo-Sanchez (PSERC Cornell & Universidad
    Autonoma de Manizales)
    @author: Ray Zimmerman (PSERC Cornell)
    """
    ##----- initialize -----
    ## unpack data
    ppc = om.get_ppc()
    baseMVA, gen, gencost = ppc["baseMVA"], ppc["gen"], ppc["gencost"]
    cp = om.get_cost_params()
    N, Cw, H, dd, rh, kk, mm = \
        cp["N"], cp["Cw"], cp["H"], cp["dd"], cp["rh"], cp["kk"], cp["mm"]
    vv, _, _, _ = om.get_idx()

    ## problem dimensions
    ng = gen.shape[0]          ## number of dispatchable injections
    ny = om.getN('var', 'y')   ## number of piece-wise linear costs
    nxyz = len(x)              ## total number of control vars of all types

    ## grab Pg & Qg
    Pg = x[vv["i1"]["Pg"]:vv["iN"]["Pg"]]  ## active generation in p.u.
    Qg = x[vv["i1"]["Qg"]:vv["iN"]["Qg"]]  ## reactive generation in p.u.

    ##----- evaluate objective function -----
    ## polynomial cost of P and Q
    # use totcost only on polynomial cost in the minimization problem
    # formulation, pwl cost is the sum of the y variables.
    ipol = find(gencost[:, MODEL] == POLYNOMIAL)   ## poly MW and MVAr costs
    xx = r_[ Pg, Qg ] * baseMVA
    if any(ipol):
        f = sum( totcost(gencost[ipol, :], xx[ipol]) )  ## cost of poly P or Q
    else:
        f = 0

    ## piecewise linear cost of P and Q
    if ny > 0:
        ccost = sparse((ones(ny),
                        (zeros(ny), arange(vv["i1"]["y"], vv["iN"]["y"]))),
                       (1, nxyz)).toarray().flatten()
        f = f + dot(ccost, x)
    else:
        ccost = zeros(nxyz)

    ## generalized cost term
    if issparse(N) and N.nnz > 0:
        nw = N.shape[0]
        r = N * x - rh                   ## Nx - rhat
        iLT = find(r < -kk)              ## below dead zone
        iEQ = find((r == 0) & (kk == 0)) ## dead zone doesn't exist
        iGT = find(r > kk)               ## above dead zone
        iND = r_[iLT, iEQ, iGT]          ## rows that are Not in the Dead region
        iL = find(dd == 1)           ## rows using linear function
        iQ = find(dd == 2)           ## rows using quadratic function
        LL = sparse((ones(len(iL)), (iL, iL)), (nw, nw))
        QQ = sparse((ones(len(iQ)), (iQ, iQ)), (nw, nw))
        kbar = sparse((r_[ones(len(iLT)), zeros(len(iEQ)), -ones(len(iGT))],
                       (iND, iND)), (nw, nw)) * kk
        rr = r + kbar                  ## apply non-dead zone shift
        M = sparse((mm[iND], (iND, iND)), (nw, nw))  ## dead zone or scale
        diagrr = sparse((rr, (arange(nw), arange(nw))), (nw, nw))

        ## linear rows multiplied by rr(i), quadratic rows by rr(i)^2
        w = M * (LL + QQ * diagrr) * rr

        f = f + dot(w * H, w) / 2 + dot(Cw, w)

    ##----- evaluate cost gradient -----
    ## index ranges
    iPg = range(vv["i1"]["Pg"], vv["iN"]["Pg"])
    iQg = range(vv["i1"]["Qg"], vv["iN"]["Qg"])

    ## polynomial cost of P and Q
    df_dPgQg = zeros(2 * ng)        ## w.r.t p.u. Pg and Qg
    df_dPgQg[ipol] = baseMVA * polycost(gencost[ipol, :], xx[ipol], 1)
    df = zeros(nxyz)
    df[iPg] = df_dPgQg[:ng]
    df[iQg] = df_dPgQg[ng:ng + ng]

    ## piecewise linear cost of P and Q
    df = df + ccost  # The linear cost row is additive wrt any nonlinear cost.

    ## generalized cost term
    if issparse(N) and N.nnz > 0:
        HwC = H * w + Cw
        AA = N.T * M * (LL + 2 * QQ * diagrr)
        df = df + AA * HwC

        ## numerical check
        if 0:    ## 1 to check, 0 to skip check
            ddff = zeros(df.shape)
            step = 1e-7
            tol  = 1e-3
            for k in range(len(x)):
                xx = x
                xx[k] = xx[k] + step
                ddff[k] = (opf_costfcn(xx, om) - f) / step
            if max(abs(ddff - df)) > tol:
                idx = find(abs(ddff - df) == max(abs(ddff - df)))
                print('Mismatch in gradient')
                print('idx             df(num)         df              diff')
                print('%4d%16g%16g%16g' %
                      (range(len(df)), ddff.T, df.T, abs(ddff - df).T))
                print('MAX')
                print('%4d%16g%16g%16g' %
                      (idx.T, ddff[idx].T, df[idx].T,
                       abs(ddff[idx] - df[idx]).T))

    if not return_hessian:
        return f, df

    ## ---- evaluate cost Hessian -----
    pcost = gencost[range(ng), :]
    if gencost.shape[0] > ng:
        qcost = gencost[ng + 1:2 * ng, :]
    else:
        qcost = array([])

    ## polynomial generator costs
    d2f_dPg2 = zeros(ng)               ## w.r.t. p.u. Pg
    d2f_dQg2 = zeros(ng)               ## w.r.t. p.u. Qg
    ipolp = find(pcost[:, MODEL] == POLYNOMIAL)
    d2f_dPg2[ipolp] = \
            baseMVA**2 * polycost(pcost[ipolp, :], Pg[ipolp]*baseMVA, 2)
    if any(qcost):          ## Qg is not free
        ipolq = find(qcost[:, MODEL] == POLYNOMIAL)
        d2f_dQg2[ipolq] = \
                baseMVA**2 * polycost(qcost[ipolq, :], Qg[ipolq] * baseMVA, 2)
    i = r_[iPg, iQg].T
    d2f = sparse((r_[d2f_dPg2, d2f_dQg2], (i, i)), (nxyz, nxyz))

    ## generalized cost
    if N is not None and issparse(N):
        d2f = d2f + AA * H * AA.T + 2 * N.T * M * QQ * \
                sparse((HwC, (range(nw), range(nw))), (nw, nw)) * N

    return f, df, d2f

Example 44

Project: pyoptools Source File: calc.py
def find_ppp(opsys, opaxis):
    """Function to find the primary principal plane location of a lens or an 
    optical component
    
    Arguments:

    
    opsys
        Optical system or optical component whose principal planes are to be 
        found
    opaxis
        Ray defining the optical axis of the system
    
    For this function to operate, the system should have a rotational symetry
    around the optical axis. 
    
    Note: 
        This function is returns the intersection point of the optical axis and
        the principal plane.
    """
    
    # Create a system with the component
    if isinstance(opsys,(Component)):
        c=opsys
        opsys=System(complist=[(c,(0,0,0),(0,0,0)),
                    ],n=1)
                    
    # To create a ray parallel to the optical axis, find a displacement vector
    # perpendicular to the optical axis, and to the XYZ axes
    
    d=opaxis.dir
    pv1= cross(d,(0,0,1))
    pv2= cross(d,(0,1,0))
    pv3= cross(d,(1,0,0))
    
    pv=[pv1,pv2,pv3]
    
    # Search for the longest pv
    pvn=array((dot(pv1,pv1),dot(pv2,pv2),dot(pv3,pv3)))
    
    pvm=pv[pvn.argmax()]

    # Create parallel ray
    
    par_ray=opaxis.copy()
    par_ray.pos=par_ray.pos+pvm*.0001
    
    opsys.clear_ray_list()
    opsys.ray_add([opaxis, par_ray])
    opsys.propagate()
    par_ray_end=par_ray.get_final_rays(inc_zeros = False)
    
    if len(par_ray_end)!=1: 
        raise Exception("The paraxial ray has more than one final ray")

    pppl=intersection(par_ray,par_ray_end[0])
    
    #Move the intersection point toward the optical axis
     
    ppp=pppl[0]-pvm*.0001
    return ppp #, pppl[1])

Example 45

Project: PYPOWER Source File: pips.py
def pips(f_fcn, x0=None, A=None, l=None, u=None, xmin=None, xmax=None,
         gh_fcn=None, hess_fcn=None, opt=None):
    """Primal-dual interior point method for NLP (nonlinear programming).
    Minimize a function F(X) beginning from a starting point M{x0}, subject to
    optional linear and nonlinear constraints and variable bounds::

            min f(x)
             x

    subject to::

            g(x) = 0            (nonlinear equalities)
            h(x) <= 0           (nonlinear inequalities)
            l <= A*x <= u       (linear constraints)
            xmin <= x <= xmax   (variable bounds)

    Note: The calling syntax is almost identical to that of FMINCON from
    MathWorks' Optimization Toolbox. The main difference is that the linear
    constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
    C{Aeq}, C{Beq}. The functions for evaluating the objective function,
    constraints and Hessian are identical.

    Example from U{http://en.wikipedia.org/wiki/Nonlinear_programming}:
        >>> from numpy import array, r_, float64, dot
        >>> from scipy.sparse import csr_matrix
        >>> def f2(x):
        ...     f = -x[0] * x[1] - x[1] * x[2]
        ...     df = -r_[x[1], x[0] + x[2], x[1]]
        ...     # actually not used since 'hess_fcn' is provided
        ...     d2f = -array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], float64)
        ...     return f, df, d2f
        >>> def gh2(x):
        ...     h = dot(array([[1, -1, 1],
        ...                    [1,  1, 1]]), x**2) + array([-2.0, -10.0])
        ...     dh = 2 * csr_matrix(array([[ x[0], x[0]],
        ...                                [-x[1], x[1]],
        ...                                [ x[2], x[2]]]))
        ...     g = array([])
        ...     dg = None
        ...     return h, g, dh, dg
        >>> def hess2(x, lam, cost_mult=1):
        ...     mu = lam["ineqnonlin"]
        ...     a = r_[dot(2 * array([1, 1]), mu), -1, 0]
        ...     b = r_[-1, dot(2 * array([-1, 1]), mu),-1]
        ...     c = r_[0, -1, dot(2 * array([1, 1]), mu)]
        ...     Lxx = csr_matrix(array([a, b, c]))
        ...     return Lxx
        >>> x0 = array([1, 1, 0], float64)
        >>> solution = pips(f2, x0, gh_fcn=gh2, hess_fcn=hess2)
        >>> round(solution["f"], 11) == -7.07106725919
        True
        >>> solution["output"]["iterations"]
        8

    Ported by Richard Lincoln from the MATLAB Interior Point Solver (MIPS)
    (v1.9) by Ray Zimmerman.  MIPS is distributed as part of the MATPOWER
    project, developed at the Power System Engineering Research Center (PSERC) (PSERC),
    Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
    MIPS was ported by Ray Zimmerman from C code written by H. Wang for his
    PhD dissertation:
      - "On the Computation and Application of Multi-period
        Security-Constrained Optimal Power Flow for Real-time
        Electricity Market Operations", Cornell University, May 2007.

    See also:
      - H. Wang, C. E. Murillo-Sanchez, R. D. Zimmerman, R. J. Thomas,
        "On Computational Issues of Market-Based Optimal Power Flow",
        IEEE Transactions on Power Systems, Vol. 22, No. 3, Aug. 2007,
        pp. 1185-1193.

    All parameters are optional except C{f_fcn} and C{x0}.
    @param f_fcn: Function that evaluates the objective function, its gradients
                  and Hessian for a given value of M{x}. If there are
                  nonlinear constraints, the Hessian information is provided
                  by the 'hess_fcn' argument and is not required here.
    @type f_fcn: callable
    @param x0: Starting value of optimization vector M{x}.
    @type x0: array
    @param A: Optional linear constraints.
    @type A: csr_matrix
    @param l: Optional linear constraints. Default values are M{-Inf}.
    @type l: array
    @param u: Optional linear constraints. Default values are M{Inf}.
    @type u: array
    @param xmin: Optional lower bounds on the M{x} variables, defaults are
                 M{-Inf}.
    @type xmin: array
    @param xmax: Optional upper bounds on the M{x} variables, defaults are
                 M{Inf}.
    @type xmax: array
    @param gh_fcn: Function that evaluates the optional nonlinear constraints
                   and their gradients for a given value of M{x}.
    @type gh_fcn: callable
    @param hess_fcn: Handle to function that computes the Hessian of the
                     Lagrangian for given values of M{x}, M{lambda} and M{mu},
                     where M{lambda} and M{mu} are the multipliers on the
                     equality and inequality constraints, M{g} and M{h},
                     respectively.
    @type hess_fcn: callable
    @param opt: optional options dictionary with the following keys, all of
                which are also optional (default values shown in parentheses)
                  - C{verbose} (False) - Controls level of progress output
                    displayed
                  - C{feastol} (1e-6) - termination tolerance for feasibility
                    condition
                  - C{gradtol} (1e-6) - termination tolerance for gradient
                    condition
                  - C{comptol} (1e-6) - termination tolerance for
                    complementarity condition
                  - C{costtol} (1e-6) - termination tolerance for cost
                    condition
                  - C{max_it} (150) - maximum number of iterations
                  - C{step_control} (False) - set to True to enable step-size
                    control
                  - C{max_red} (20) - maximum number of step-size reductions if
                    step-control is on
                  - C{cost_mult} (1.0) - cost multiplier used to scale the
                    objective function for improved conditioning. Note: This
                    value is also passed as the 3rd argument to the Hessian
                    evaluation function so that it can appropriately scale the
                    objective function term in the Hessian of the Lagrangian.
    @type opt: dict

    @rtype: dict
    @return: The solution dictionary has the following keys:
               - C{x} - solution vector
               - C{f} - final objective function value
               - C{converged} - exit status
                   - True = first order optimality conditions satisfied
                   - False = maximum number of iterations reached
                   - None = numerically failed
               - C{output} - output dictionary with keys:
                   - C{iterations} - number of iterations performed
                   - C{hist} - list of arrays with trajectories of the
                     following: feascond, gradcond, compcond, costcond, gamma,
                     stepsize, obj, alphap, alphad
                   - C{message} - exit message
               - C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
                 multipliers on the constraints, with keys:
                   - C{eqnonlin} - nonlinear equality constraints
                   - C{ineqnonlin} - nonlinear inequality constraints
                   - C{mu_l} - lower (left-hand) limit on linear constraints
                   - C{mu_u} - upper (right-hand) limit on linear constraints
                   - C{lower} - lower bound on optimization variables
                   - C{upper} - upper bound on optimization variables

    @see: U{http://www.pserc.cornell.edu/matpower/}

    @author: Ray Zimmerman (PSERC Cornell)
    """
    if isinstance(f_fcn, dict):  ## problem dict
        p = f_fcn
        f_fcn = p['f_fcn']
        x0 = p['x0']
        if 'opt' in p: opt = p['opt']
        if 'hess_fcn' in p: hess_fcn = p['hess_fcn']
        if 'gh_fcn' in p: gh_fcn = p['gh_fcn']
        if 'xmax' in p: xmax = p['xmax']
        if 'xmin' in p: xmin = p['xmin']
        if 'u' in p: u = p['u']
        if 'l' in p: l = p['l']
        if 'A' in p: A = p['A']

    nx = x0.shape[0]                        # number of variables
    nA = A.shape[0] if A is not None else 0 # number of original linear constr

    # default argument values
    if l is None or len(l) == 0: l = -Inf * ones(nA)
    if u is None or len(u) == 0: u =  Inf * ones(nA)
    if xmin is None or len(xmin) == 0: xmin = -Inf * ones(x0.shape[0])
    if xmax is None or len(xmax) == 0: xmax =  Inf * ones(x0.shape[0])
    if gh_fcn is None:
        nonlinear = False
        gn = array([])
        hn = array([])
    else:
        nonlinear = True

    if opt is None: opt = {}
    # options
    if "feastol" not in opt:
        opt["feastol"] = 1e-06
    if "gradtol" not in opt:
        opt["gradtol"] = 1e-06
    if "comptol" not in opt:
        opt["comptol"] = 1e-06
    if "costtol" not in opt:
        opt["costtol"] = 1e-06
    if "max_it" not in opt:
        opt["max_it"] = 150
    if "max_red" not in opt:
        opt["max_red"] = 20
    if "step_control" not in opt:
        opt["step_control"] = False
    if "cost_mult" not in opt:
        opt["cost_mult"] = 1
    if "verbose" not in opt:
        opt["verbose"] = 0

    # initialize history
    hist = []

    # constants
    xi = 0.99995
    sigma = 0.1
    z0 = 1
    alpha_min = 1e-8
    rho_min = 0.95
    rho_max = 1.05
    mu_threshold = 1e-5

    # initialize
    i = 0                       # iteration counter
    converged = False           # flag
    eflag = False               # exit flag

    # add var limits to linear constraints
    eyex = eye(nx, nx, format="csr")
    AA = eyex if A is None else vstack([eyex, A], "csr")
    ll = r_[xmin, l]
    uu = r_[xmax, u]

    # split up linear constraints
    ieq = find( absolute(uu - ll) <= EPS )
    igt = find( (uu >=  1e10) & (ll > -1e10) )
    ilt = find( (ll <= -1e10) & (uu <  1e10) )
    ibx = find( (absolute(uu - ll) > EPS) & (uu < 1e10) & (ll > -1e10) )
    # zero-sized sparse matrices unsupported
    Ae = AA[ieq, :] if len(ieq) else None
    if len(ilt) or len(igt) or len(ibx):
        idxs = [(1, ilt), (-1, igt), (1, ibx), (-1, ibx)]
        Ai = vstack([sig * AA[idx, :] for sig, idx in idxs if len(idx)], 'csr')
    else:
        Ai = None
    be = uu[ieq]
    bi = r_[uu[ilt], -ll[igt], uu[ibx], -ll[ibx]]

    # evaluate cost f(x0) and constraints g(x0), h(x0)
    x = x0
    f, df = f_fcn(x)                 # cost
    f = f * opt["cost_mult"]
    df = df * opt["cost_mult"]
    if nonlinear:
        hn, gn, dhn, dgn = gh_fcn(x)        # nonlinear constraints
        h = hn if Ai is None else r_[hn, Ai * x - bi] # inequality constraints
        g = gn if Ae is None else r_[gn, Ae * x - be] # equality constraints

        if (dhn is None) and (Ai is None):
            dh = None
        elif dhn is None:
            dh = Ai.T
        elif Ai is None:
            dh = dhn
        else:
            dh = hstack([dhn, Ai.T])

        if (dgn is None) and (Ae is None):
            dg = None
        elif dgn is None:
            dg = Ae.T
        elif Ae is None:
            dg = dgn
        else:
            dg = hstack([dgn, Ae.T])
    else:
        h = -bi if Ai is None else Ai * x - bi        # inequality constraints
        g = -be if Ae is None else Ae * x - be        # equality constraints
        dh = None if Ai is None else Ai.T     # 1st derivative of inequalities
        dg = None if Ae is None else Ae.T     # 1st derivative of equalities

    # some dimensions
    neq = g.shape[0]           # number of equality constraints
    niq = h.shape[0]           # number of inequality constraints
    neqnln = gn.shape[0]       # number of nonlinear equality constraints
    niqnln = hn.shape[0]       # number of nonlinear inequality constraints
    nlt = len(ilt)             # number of upper bounded linear inequalities
    ngt = len(igt)             # number of lower bounded linear inequalities
    nbx = len(ibx)             # number of doubly bounded linear inequalities

    # initialize gamma, lam, mu, z, e
    gamma = 1                  # barrier coefficient
    lam = zeros(neq)
    z = z0 * ones(niq)
    mu = z0 * ones(niq)
    k = find(h < -z0)
    z[k] = -h[k]
    k = find((gamma / z) > z0)
    mu[k] = gamma / z[k]
    e = ones(niq)

    # check tolerance
    f0 = f
    if opt["step_control"]:
        L = f + dot(lam, g) + dot(mu, h + z) - gamma * sum(log(z))

    Lx = df.copy()
    Lx = Lx + dg * lam if dg is not None else Lx
    Lx = Lx + dh * mu  if dh is not None else Lx

    maxh = zeros(1) if len(h) == 0 else max(h)

    gnorm = norm(g, Inf) if len(g) else 0.0
    lam_norm = norm(lam, Inf) if len(lam) else 0.0
    mu_norm = norm(mu, Inf) if len(mu) else 0.0
    znorm = norm(z, Inf) if len(z) else 0.0
    feascond = \
        max([gnorm, maxh]) / (1 + max([norm(x, Inf), znorm]))
    gradcond = \
        norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
    compcond = dot(z, mu) / (1 + norm(x, Inf))
    costcond = absolute(f - f0) / (1 + absolute(f0))

    # save history
    hist.append({'feascond': feascond, 'gradcond': gradcond,
        'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
        'stepsize': 0, 'obj': f / opt["cost_mult"], 'alphap': 0, 'alphad': 0})

    if opt["verbose"]:
        s = '-sc' if opt["step_control"] else ''
        v = pipsver('all')
        print('Python Interior Point Solver - PIPS%s, Version %s, %s' %
                    (s, v['Version'], v['Date']))
        if opt['verbose'] > 1:
            print(" it    objective   step size   feascond     gradcond     "
                  "compcond     costcond  ")
            print("----  ------------ --------- ------------ ------------ "
                  "------------ ------------")
            print("%3d  %12.8g %10s %12g %12g %12g %12g" %
                (i, (f / opt["cost_mult"]), "",
                 feascond, gradcond, compcond, costcond))

    if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
        compcond < opt["comptol"] and costcond < opt["costtol"]:
        converged = True
        if opt["verbose"]:
            print("Converged!")

    # do Newton iterations
    while (not converged) and (i < opt["max_it"]):
        # update iteration counter
        i += 1

        # compute update step
        lmbda = {"eqnonlin": lam[range(neqnln)],
                 "ineqnonlin": mu[range(niqnln)]}
        if nonlinear:
            if hess_fcn is None:
                print("pips: Hessian evaluation via finite differences "
                      "not yet implemented.\nPlease provide "
                      "your own hessian evaluation function.")
            Lxx = hess_fcn(x, lmbda, opt["cost_mult"])
        else:
            _, _, d2f = f_fcn(x, True)      # cost
            Lxx = d2f * opt["cost_mult"]
        rz = range(len(z))
        zinvdiag = sparse((1.0 / z, (rz, rz))) if len(z) else None
        rmu = range(len(mu))
        mudiag = sparse((mu, (rmu, rmu))) if len(mu) else None
        dh_zinv = None if dh is None else dh * zinvdiag
        M = Lxx if dh is None else Lxx + dh_zinv * mudiag * dh.T
        N = Lx if dh is None else Lx + dh_zinv * (mudiag * h + gamma * e)

        Ab = sparse(M) if dg is None else vstack([
            hstack([M, dg]),
            hstack([dg.T, sparse((neq, neq))])
        ])
        bb = r_[-N, -g]

        dxdlam = spsolve(Ab.tocsr(), bb)

        if any(isnan(dxdlam)):
            if opt["verbose"]:
                print('\nNumerically Failed\n')
            eflag = -1
            break

        dx = dxdlam[:nx]
        dlam = dxdlam[nx:nx + neq]
        dz = -h - z if dh is None else -h - z - dh.T * dx
        dmu = -mu if dh is None else -mu + zinvdiag * (gamma * e - mudiag * dz)

        # optional step-size control
        sc = False
        if opt["step_control"]:
            x1 = x + dx

            # evaluate cost, constraints, derivatives at x1
            f1, df1 = f_fcn(x1)          # cost
            f1 = f1 * opt["cost_mult"]
            df1 = df1 * opt["cost_mult"]
            if nonlinear:
                hn1, gn1, dhn1, dgn1 = gh_fcn(x1) # nonlinear constraints

                h1 = hn1 if Ai is None else r_[hn1, Ai * x1 - bi] # ieq constraints
                g1 = gn1 if Ae is None else r_[gn1, Ae * x1 - be] # eq constraints

                # 1st der of ieq
                if (dhn1 is None) and (Ai is None):
                    dh1 = None
                elif dhn1 is None:
                    dh1 = Ai.T
                elif Ai is None:
                    dh1 = dhn1
                else:
                    dh1 = hstack([dhn1, Ai.T])

                # 1st der of eqs
                if (dgn1 is None) and (Ae is None):
                    dg1 = None
                elif dgn is None:
                    dg1 = Ae.T
                elif Ae is None:
                    dg1 = dgn1
                else:
                    dg1 = hstack([dgn1, Ae.T])
            else:
                h1 = -bi if Ai is None else Ai * x1 - bi    # inequality constraints
                g1 = -be if Ae is None else Ae * x1 - be    # equality constraints

                dh1 = dh                       ## 1st derivative of inequalities
                dg1 = dg                       ## 1st derivative of equalities

            # check tolerance
            Lx1 = df1
            Lx1 = Lx1 + dg1 * lam if dg1 is not None else Lx1
            Lx1 = Lx1 + dh1 * mu  if dh1 is not None else Lx1

            maxh1 = zeros(1) if len(h1) == 0 else max(h1)

            g1norm = norm(g1, Inf) if len(g1) else 0.0
            lam1_norm = norm(lam, Inf) if len(lam) else 0.0
            mu1_norm = norm(mu, Inf) if len(mu) else 0.0
            z1norm = norm(z, Inf) if len(z) else 0.0

            feascond1 = max([ g1norm, maxh1 ]) / \
                (1 + max([ norm(x1, Inf), z1norm ]))
            gradcond1 = norm(Lx1, Inf) / (1 + max([ lam1_norm, mu1_norm ]))

            if (feascond1 > feascond) and (gradcond1 > gradcond):
                sc = True
        if sc:
            alpha = 1.0
            for j in range(opt["max_red"]):
                dx1 = alpha * dx
                x1 = x + dx1
                f1, _ = f_fcn(x1)             # cost
                f1 = f1 * opt["cost_mult"]
                if nonlinear:
                    hn1, gn1, _, _ = gh_fcn(x1)              # nonlinear constraints
                    h1 = hn1 if Ai is None else r_[hn1, Ai * x1 - bi]         # inequality constraints
                    g1 = gn1 if Ae is None else r_[gn1, Ae * x1 - be]         # equality constraints
                else:
                    h1 = -bi if Ai is None else Ai * x1 - bi    # inequality constraints
                    g1 = -be if Ae is None else Ae * x1 - be    # equality constraints

                L1 = f1 + dot(lam, g1) + dot(mu, h1 + z) - gamma * sum(log(z))

                if opt["verbose"] > 2:
                    print("   %3d            %10.5f" % (-j, norm(dx1)))

                rho = (L1 - L) / (dot(Lx, dx1) + 0.5 * dot(dx1, Lxx * dx1))

                if (rho > rho_min) and (rho < rho_max):
                    break
                else:
                    alpha = alpha / 2.0
            dx = alpha * dx
            dz = alpha * dz
            dlam = alpha * dlam
            dmu = alpha * dmu

        # do the update
        k = find(dz < 0.0)
        alphap = min([xi * min(z[k] / -dz[k]), 1]) if len(k) else 1.0
        k = find(dmu < 0.0)
        alphad = min([xi * min(mu[k] / -dmu[k]), 1]) if len(k) else 1.0
        x = x + alphap * dx
        z = z + alphap * dz
        lam = lam + alphad * dlam
        mu = mu + alphad * dmu
        if niq > 0:
            gamma = sigma * dot(z, mu) / niq

        # evaluate cost, constraints, derivatives
        f, df = f_fcn(x)             # cost
        f = f * opt["cost_mult"]
        df = df * opt["cost_mult"]
        if nonlinear:
            hn, gn, dhn, dgn = gh_fcn(x)                   # nln constraints
#            g = gn if Ai is None else r_[gn, Ai * x - bi] # ieq constraints
#            h = hn if Ae is None else r_[hn, Ae * x - be] # eq constraints
            h = hn if Ai is None else r_[hn, Ai * x - bi] # ieq constr
            g = gn if Ae is None else r_[gn, Ae * x - be]  # eq constr

            if (dhn is None) and (Ai is None):
                dh = None
            elif dhn is None:
                dh = Ai.T
            elif Ai is None:
                dh = dhn
            else:
                dh = hstack([dhn, Ai.T])

            if (dgn is None) and (Ae is None):
                dg = None
            elif dgn is None:
                dg = Ae.T
            elif Ae is None:
                dg = dgn
            else:
                dg = hstack([dgn, Ae.T])
        else:
            h = -bi if Ai is None else Ai * x - bi    # inequality constraints
            g = -be if Ae is None else Ae * x - be    # equality constraints
            # 1st derivatives are constant, still dh = Ai.T, dg = Ae.T

        Lx = df
        Lx = Lx + dg * lam if dg is not None else Lx
        Lx = Lx + dh * mu  if dh is not None else Lx

        if len(h) == 0:
            maxh = zeros(1)
        else:
            maxh = max(h)

        gnorm = norm(g, Inf) if len(g) else 0.0
        lam_norm = norm(lam, Inf) if len(lam) else 0.0
        mu_norm = norm(mu, Inf) if len(mu) else 0.0
        znorm = norm(z, Inf) if len(z) else 0.0
        feascond = \
            max([gnorm, maxh]) / (1 + max([norm(x, Inf), znorm]))
        gradcond = \
            norm(Lx, Inf) / (1 + max([lam_norm, mu_norm]))
        compcond = dot(z, mu) / (1 + norm(x, Inf))
        costcond = float(absolute(f - f0) / (1 + absolute(f0)))

        hist.append({'feascond': feascond, 'gradcond': gradcond,
            'compcond': compcond, 'costcond': costcond, 'gamma': gamma,
            'stepsize': norm(dx), 'obj': f / opt["cost_mult"],
            'alphap': alphap, 'alphad': alphad})

        if opt["verbose"] > 1:
            print("%3d  %12.8g %10.5g %12g %12g %12g %12g" %
                (i, (f / opt["cost_mult"]), norm(dx), feascond, gradcond,
                 compcond, costcond))

        if feascond < opt["feastol"] and gradcond < opt["gradtol"] and \
            compcond < opt["comptol"] and costcond < opt["costtol"]:
            converged = True
            if opt["verbose"]:
                print("Converged!")
        else:
            if any(isnan(x)) or (alphap < alpha_min) or \
                (alphad < alpha_min) or (gamma < EPS) or (gamma > 1.0 / EPS):
                if opt["verbose"]:
                    print("Numerically failed.")
                eflag = -1
                break
            f0 = f

            if opt["step_control"]:
                L = f + dot(lam, g) + dot(mu, (h + z)) - gamma * sum(log(z))

    if opt["verbose"]:
        if not converged:
            print("Did not converge in %d iterations." % i)

    # package results
    if eflag != -1:
        eflag = converged

    if eflag == 0:
        message = 'Did not converge'
    elif eflag == 1:
        message = 'Converged'
    elif eflag == -1:
        message = 'Numerically failed'
    else:
        raise

    output = {"iterations": i, "hist": hist, "message": message}

    # zero out multipliers on non-binding constraints
    mu[find( (h < -opt["feastol"]) & (mu < mu_threshold) )] = 0.0

    # un-scale cost and prices
    f = f / opt["cost_mult"]
    lam = lam / opt["cost_mult"]
    mu = mu / opt["cost_mult"]

    # re-package multipliers into struct
    lam_lin = lam[neqnln:neq]           # lambda for linear constraints
    mu_lin = mu[niqnln:niq]             # mu for linear constraints
    kl = find(lam_lin < 0.0)     # lower bound binding
    ku = find(lam_lin > 0.0)     # upper bound binding

    mu_l = zeros(nx + nA)
    mu_l[ieq[kl]] = -lam_lin[kl]
    mu_l[igt] = mu_lin[nlt:nlt + ngt]
    mu_l[ibx] = mu_lin[nlt + ngt + nbx:nlt + ngt + nbx + nbx]

    mu_u = zeros(nx + nA)
    mu_u[ieq[ku]] = lam_lin[ku]
    mu_u[ilt] = mu_lin[:nlt]
    mu_u[ibx] = mu_lin[nlt + ngt:nlt + ngt + nbx]

    lmbda = {'mu_l': mu_l[nx:], 'mu_u': mu_u[nx:],
             'lower': mu_l[:nx], 'upper': mu_u[:nx]}

    if niqnln > 0:
        lmbda['ineqnonlin'] = mu[:niqnln]
    if neqnln > 0:
        lmbda['eqnonlin'] = lam[:neqnln]

#    lmbda = {"eqnonlin": lam[:neqnln], 'ineqnonlin': mu[:niqnln],
#             "mu_l": mu_l[nx:], "mu_u": mu_u[nx:],
#             "lower": mu_l[:nx], "upper": mu_u[:nx]}

    solution =  {"x": x, "f": f, "eflag": converged,
                 "output": output, "lmbda": lmbda}

    return solution

Example 46

Project: pele Source File: _lbfgs_py.py
    def adjustStepSize(self, X, E, G, stp):
        """
        We now have a proposed step.  This function will make sure it is 
        a good step and then take it.  This is known as a Backtracking linesearch
        
        http://en.wikipedia.org/wiki/Backtracking_line_search
        
        1) if the step is not anti-aligned with the gradient (i.e. downhill), 
           then reverse the step
        
        2) if the step is larger than maxstep, then rescale the step
        
        3) calculate the energy and gradient of the new position
        
        4) if the step increases the energy by more than maxErise, 
            then reduce the step size and go to 3)
        
        5) if the step is reduced more than 10 times and the energy is still 
           not acceptable, then increment nfail, reset the lbfgs optimizer and 
           continue
            
        6) if nfail is greater than 5 abort the quench
                
        """
        f = 1.
        X0 = X.copy()
        G0 = G.copy()
        E0 = E

        if np.dot(G, stp) > 0:
            if self.debug:
                overlap = np.dot(G, stp) / np.linalg.norm(G) / np.linalg.norm(stp)
                self.logger.warn("LBFGS returned uphill step, reversing step direction: overlap {}".format(overlap))
            stp = -stp

        stepsize = np.linalg.norm(stp)

        if f * stepsize > self.maxstep:
            f = self.maxstep / stepsize

        nincrease = 0
        while True:
            X = X0 + f * stp
            E, G = self.pot.getEnergyGradient(X)
            self.funcalls += 1

            # if the increase is greater than maxErise reduce the step size
            if self._accept_step(E, E0, G, G0, f * stp):
                break
            else:
                if self.debug:
                    self.logger.warn("energy increased, trying a smaller step %s %s %s %s", E, E0, f * stepsize,
                                     nincrease)
                f /= 10.
                nincrease += 1
                if nincrease > 10:
                    break

        if nincrease > 10:
            self.nfailed += 1
            if self.nfailed > 10:
                raise (LineSearchError("lbfgs: too many failures in adjustStepSize, exiting"))

            # abort the linesearch, reset the memory and reset the coordinates            
            self.logger.warning("lbfgs: having trouble finding a good step size. %s %s, resetting the minimizer",
                                f * stepsize, stepsize)
            self.reset()
            E = E0
            G = G0
            X = X0
            f = 0.

        self.stepsize = f * stepsize
        return X, E, G

Example 47

Project: pyconsensus Source File: __init__.py
    def consensus(self):
        # Handle missing values
        reports_filled = self.interpolate(self.reports)

        # Consensus - Row Players
        player_info = self.lie_detector(reports_filled)

        # Column Players (The Event Creators)
        outcomes_raw = np.dot(player_info['smooth_rep'], reports_filled)
        if outcomes_raw.shape != (1,):
            outcomes_raw = outcomes_raw.squeeze()

        # Discriminate Based on Contract Type
        if self.event_bounds is not None:
            for i in range(reports_filled.shape[1]):

                # Our Current best-guess for this Scaled Event (weighted median)
                if self.event_bounds[i]["scaled"]:
                    outcomes_raw[i] = weighted_median(
                        reports_filled[:,i],
                        weights=player_info["smooth_rep"].ravel(),
                    )

        # The Outcome (Discriminate Based on Contract Type)
        outcomes_adj = []
        for i, raw in enumerate(outcomes_raw):
            if self.event_bounds is not None and self.event_bounds[i]["scaled"]:
                outcomes_adj.append(raw)
            else:
                outcomes_adj.append(self.catch(raw))

        outcomes_final = []
        for i, raw in enumerate(outcomes_raw):
            outcomes_final.append(outcomes_adj[i])
            if self.event_bounds is not None and self.event_bounds[i]["scaled"]:
                outcomes_final[i] *= self.event_bounds[i]["max"] - self.event_bounds[i]["min"]
                outcomes_final[i] += self.event_bounds[i]["min"]

        certainty = []
        for i, adj in enumerate(outcomes_adj):
            certainty.append(sum(player_info["smooth_rep"][reports_filled[:,i] == adj]))

        certainty = np.array(certainty)
        consensus_reward = self.normalize(certainty)
        avg_certainty = np.mean(certainty)

        # Participation: information about missing values
        na_mat = self.reports * 0
        na_mat[np.isnan(self.reports)] = 1  # indicator matrix for missing
        na_mat[self.reports == NA] = 1
        if self.verbose:
            print "NA Mat:"
            print na_mat
            print

        # Participation Within Events (Columns)
        # % of reputation that answered each Event
        participation_columns = 1 - np.dot(player_info['smooth_rep'], na_mat)

        # Participation Within Agents (Rows)
        # Democracy Option - all Events treated equally.
        if self.verbose:
            print "Sum:"
            print na_mat.sum(axis=1)
            print
        participation_rows = 1 - na_mat.sum(axis=1) / na_mat.shape[1]

        # General Participation
        percent_na = 1 - np.mean(participation_columns)
        if self.verbose:
            print percent_na

        # Combine Information
        # Row
        na_bonus_reporters = self.normalize(participation_rows)
        reporter_bonus = na_bonus_reporters * percent_na + player_info['smooth_rep'] * (1 - percent_na)

        # Column
        na_bonus_events = self.normalize(participation_columns)
        author_bonus = na_bonus_events * percent_na + consensus_reward * (1 - percent_na)

        return {
            'original': self.reports.data,
            'filled': reports_filled.data,
            'agents': {
                'old_rep': player_info['old_rep'],
                'this_rep': player_info['this_rep'],
                'smooth_rep': player_info['smooth_rep'],
                'na_row': na_mat.sum(axis=1).data.tolist(),
                'participation_rows': participation_rows.data.tolist(),
                'relative_part': na_bonus_reporters.data.tolist(),
                'reporter_bonus': reporter_bonus.data.tolist(),
                'scores': player_info['scores'],
            },
            'events': {
                'adj_first_loadings': player_info['first_loading'].data.tolist(),
                'outcomes_raw': outcomes_raw.tolist(),
                'consensus_reward': consensus_reward,
                'certainty': certainty,
                'NAs Filled': na_mat.sum(axis=0).data.tolist(),
                'participation_columns': participation_columns.data.tolist(),
                'author_bonus': author_bonus.data.tolist(),
                'outcomes_adjusted': outcomes_adj,
                'outcomes_final': outcomes_final,
            },
            'participation': 1 - percent_na,
            'avg_certainty': avg_certainty,
            'convergence': self.convergence,
            'components': self.num_components,
        }

Example 48

Project: pyspace Source File: SOR.py
    def reduced_descent(self, current_dual, M, relevant_indices):
        """ Basic iteration step over a set of indices, possibly subset of all

        The main principle is to make a descent step with just one index,
        while fixing the other dual_solutions.

        The main formula comes from *M&M*:

        .. math::

            d        = \\alpha_i - \\frac{\\omega}{M[i][i]}(M[i]\\alpha-1)

            \\text{with } M[i][j]  = y_i y_j(<x_i,x_j>+1)

            \\text{and final projection: }\\alpha_i = \\max(0,\\min(d,c_i)).

        Here we use c for the weights for each sample in the loss term,
        which is normally complexity times corresponding class weight.
        y is used for the labels, which have to be 1 or -1.

        In the *sample* version only the diagonal of M is used.
        The sum with the alpha is tracked by using the classification vector w
        and the offset b.

        .. math::

            o        = \\alpha_i

            d        = \\alpha_i - \\frac{\\omega}{M[i][i]}(y_i(<w,x_i>+b)-1)

            \\text{with projection: }\\alpha_i = \\max(0,\\min(d,c_i)),

            b=b+(\\alpha_i-o)y_i

            w=w+(\\alpha_i-o)y_i x_i
        """
        self.irrelevant_indices = []
        self.difference = 0
        for i in relevant_indices:
            old_dual = current_dual[i]
            ### Main Function ###
            ### elemental update step of SOR algorithm ###

            if self.version == "matrix":
                # this step is kernel independent
                x = old_dual - self.omega / (
                    M[i][i] + self.squ_factor/(2 * self.ci[i])) * \
                    (dot(M[i], current_dual) - 1)
            elif self.version == "samples":
                xi = self.samples[i]
                bi = self.bi[i]
                x = old_dual - self.omega * (M[i]) * \
                    (bi * (dot(xi.T, self.w) + self.b) - 1 +
                     self.squ_factor * old_dual / (2 * self.ci[i]))
            # map dual solution to the interval [0,C]
            if x <= 0:
                self.irrelevant_indices.append(i)
                current_dual[i] = 0
            elif not self.squ_factor:
                current_dual[i] = min(x, self.ci[i])
            else:
                current_dual[i] = x
            if self.version == "matrix":
                delta = (current_dual[i] - old_dual)
                # update w and b in samples case
            if self.version == "samples":
                delta = (current_dual[i] - old_dual) * bi
                # update classification function parameter w and b
                # self.update_classification_function(delta=delta, index=i)
                self.b = self.b + self.offset_factor * delta
                self.w = self.w + delta * xi
            current_difference = numpy.abs(delta)
            if current_difference > self.difference:
                self.difference = current_difference
            self.sub_iterations += 1
            self.iterations += 1
            
            if not (self.sub_iterations < self.max_sub_iterations
                    and self.iterations < self.max_iterations):
                break
        if self.reduce_non_zeros:
            for index in self.irrelevant_indices:
                try:
                    relevant_indices.remove(index)
                except:
                    # special mapping for RMM case
                    if index < self.num_samples:
                        relevant_indices.remove(index+self.num_samples)
                    else:
                        relevant_indices.remove(index-self.num_samples)
        if self.random:
            random.shuffle(relevant_indices)

Example 49

Project: PYPOWER Source File: qps_pips.py
def qps_pips(H, c, A, l, u, xmin=None, xmax=None, x0=None, opt=None):
    """Uses the Python Interior Point Solver (PIPS) to solve the following
    QP (quadratic programming) problem::

            min 1/2 x'*H*x + C'*x
             x

    subject to::

            l <= A*x <= u       (linear constraints)
            xmin <= x <= xmax   (variable bounds)

    Note the calling syntax is almost identical to that of QUADPROG from
    MathWorks' Optimization Toolbox. The main difference is that the linear
    constraints are specified with C{A}, C{L}, C{U} instead of C{A}, C{B},
    C{Aeq}, C{Beq}.

    Example from U{http://www.uc.edu/sashtml/iml/chap8/sect12.htm}:

        >>> from numpy import array, zeros, Inf
        >>> from scipy.sparse import csr_matrix
        >>> H = csr_matrix(array([[1003.1,  4.3,     6.3,     5.9],
        ...                       [4.3,     2.2,     2.1,     3.9],
        ...                       [6.3,     2.1,     3.5,     4.8],
        ...                       [5.9,     3.9,     4.8,     10 ]]))
        >>> c = zeros(4)
        >>> A = csr_matrix(array([[1,       1,       1,       1   ],
        ...                       [0.17,    0.11,    0.10,    0.18]]))
        >>> l = array([1, 0.10])
        >>> u = array([1, Inf])
        >>> xmin = zeros(4)
        >>> xmax = None
        >>> x0 = array([1, 0, 0, 1])
        >>> solution = qps_pips(H, c, A, l, u, xmin, xmax, x0)
        >>> round(solution["f"], 11) == 1.09666678128
        True
        >>> solution["converged"]
        True
        >>> solution["output"]["iterations"]
        10

    All parameters are optional except C{H}, C{c}, C{A} and C{l} or C{u}.
    @param H: Quadratic cost coefficients.
    @type H: csr_matrix
    @param c: vector of linear cost coefficients
    @type c: array
    @param A: Optional linear constraints.
    @type A: csr_matrix
    @param l: Optional linear constraints. Default values are M{-Inf}.
    @type l: array
    @param u: Optional linear constraints. Default values are M{Inf}.
    @type u: array
    @param xmin: Optional lower bounds on the M{x} variables, defaults are
                 M{-Inf}.
    @type xmin: array
    @param xmax: Optional upper bounds on the M{x} variables, defaults are
                 M{Inf}.
    @type xmax: array
    @param x0: Starting value of optimization vector M{x}.
    @type x0: array
    @param opt: optional options dictionary with the following keys, all of
                which are also optional (default values shown in parentheses)
                  - C{verbose} (False) - Controls level of progress output
                    displayed
                  - C{feastol} (1e-6) - termination tolerance for feasibility
                    condition
                  - C{gradtol} (1e-6) - termination tolerance for gradient
                    condition
                  - C{comptol} (1e-6) - termination tolerance for
                    complementarity condition
                  - C{costtol} (1e-6) - termination tolerance for cost
                    condition
                  - C{max_it} (150) - maximum number of iterations
                  - C{step_control} (False) - set to True to enable step-size
                    control
                  - C{max_red} (20) - maximum number of step-size reductions if
                    step-control is on
                  - C{cost_mult} (1.0) - cost multiplier used to scale the
                    objective function for improved conditioning. Note: The
                    same value must also be passed to the Hessian evaluation
                    function so that it can appropriately scale the objective
                    function term in the Hessian of the Lagrangian.
    @type opt: dict

    @rtype: dict
    @return: The solution dictionary has the following keys:
               - C{x} - solution vector
               - C{f} - final objective function value
               - C{converged} - exit status
                   - True = first order optimality conditions satisfied
                   - False = maximum number of iterations reached
                   - None = numerically failed
               - C{output} - output dictionary with keys:
                   - C{iterations} - number of iterations performed
                   - C{hist} - dictionary of arrays with trajectories of the
                     following: feascond, gradcond, coppcond, costcond, gamma,
                     stepsize, obj, alphap, alphad
                   - C{message} - exit message
               - C{lmbda} - dictionary containing the Langrange and Kuhn-Tucker
                 multipliers on the constraints, with keys:
                   - C{eqnonlin} - nonlinear equality constraints
                   - C{ineqnonlin} - nonlinear inequality constraints
                   - C{mu_l} - lower (left-hand) limit on linear constraints
                   - C{mu_u} - upper (right-hand) limit on linear constraints
                   - C{lower} - lower bound on optimization variables
                   - C{upper} - upper bound on optimization variables

    @see: L{pips}

    @author: Ray Zimmerman (PSERC Cornell)
    """
    if isinstance(H, dict):
        p = H
    else:
        p = {'H': H, 'c': c, 'A': A, 'l': l, 'u': u}
        if xmin is not None: p['xmin'] = xmin
        if xmax is not None: p['xmax'] = xmax
        if x0 is not None: p['x0'] = x0
        if opt is not None: p['opt'] = opt

    if 'H' not in p or p['H'] == None:#p['H'].nnz == 0:
        if p['A'] is None or p['A'].nnz == 0 and \
           'xmin' not in p and \
           'xmax' not in p:
#           'xmin' not in p or len(p['xmin']) == 0 and \
#           'xmax' not in p or len(p['xmax']) == 0:
            print('qps_pips: LP problem must include constraints or variable bounds')
            return
        else:
            if p['A'] is not None and p['A'].nnz >= 0:
                nx = p['A'].shape[1]
            elif 'xmin' in p and len(p['xmin']) > 0:
                nx = p['xmin'].shape[0]
            elif 'xmax' in p and len(p['xmax']) > 0:
                nx = p['xmax'].shape[0]
        p['H'] = sparse((nx, nx))
    else:
        nx = p['H'].shape[0]

    p['xmin'] = -Inf * ones(nx) if 'xmin' not in p else p['xmin']
    p['xmax'] =  Inf * ones(nx) if 'xmax' not in p else p['xmax']

    p['c'] = zeros(nx) if p['c'] is None else p['c']

    p['x0'] = zeros(nx) if 'x0' not in p else p['x0']

    def qp_f(x, return_hessian=False):
        f = 0.5 * dot(x * p['H'], x) + dot(p['c'], x)
        df = p['H'] * x + p['c']
        if not return_hessian:
            return f, df
        d2f = p['H']
        return f, df, d2f

    p['f_fcn'] = qp_f

    sol = pips(p)

    return sol["x"], sol["f"], sol["eflag"], sol["output"], sol["lmbda"]

Example 50

Project: scikit-image Source File: _geometric.py
def _umeyama(src, dst, estimate_scale):
    """Estimate N-D similarity transformation with or without scaling.

    Parameters
    ----------
    src : (M, N) array
        Source coordinates.
    dst : (M, N) array
        Destination coordinates.
    estimate_scale : bool
        Whether to estimate scaling factor.

    Returns
    -------
    T : (N + 1, N + 1)
        The humogeneous similarity transformation matrix. The matrix contains
        NaN values only if the problem is not well-conditioned.

    References
    ----------
    .. [1] "Least-squares estimation of transformation parameters between two
            point patterns", Shinji Umeyama, PAMI 1991, DOI: 10.1109/34.88573

    """

    num = src.shape[0]
    dim = src.shape[1]

    # Compute mean of src and dst.
    src_mean = src.mean(axis=0)
    dst_mean = dst.mean(axis=0)

    # Subtract mean from src and dst.
    src_demean = src - src_mean
    dst_demean = dst - dst_mean

    # Eq. (38).
    A = np.dot(dst_demean.T, src_demean) / num

    # Eq. (39).
    d = np.ones((dim,), dtype=np.double)
    if np.linalg.det(A) < 0:
        d[dim - 1] = -1

    T = np.eye(dim + 1, dtype=np.double)

    U, S, V = np.linalg.svd(A)

    # Eq. (40) and (43).
    rank = np.linalg.matrix_rank(A)
    if rank == 0:
        return np.nan * T
    elif rank == dim - 1:
        if np.linalg.det(U) * np.linalg.det(V) > 0:
            T[:dim, :dim] = np.dot(U, V)
        else:
            s = d[dim - 1]
            d[dim - 1] = -1
            T[:dim, :dim] = np.dot(U, np.dot(np.diag(d), V))
            d[dim - 1] = s
    else:
        T[:dim, :dim] = np.dot(U, np.dot(np.diag(d), V.T))

    if estimate_scale:
        # Eq. (41) and (42).
        scale = 1.0 / src_demean.var(axis=0).sum() * np.dot(S, d)
    else:
        scale = 1.0

    T[:dim, dim] = dst_mean - scale * np.dot(T[:dim, :dim], src_mean.T)
    T[:dim, :dim] *= scale

    return T
See More Examples - Go to Next Page
Page 1 Selected Page 2 Page 3 Page 4