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
5
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)
5
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)
5
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
2
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
2
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()
2
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
2
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()
2
Example 8
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
2
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
2
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
2
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
2
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
2
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())
2
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
0
Example 15
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
0
Example 16
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
0
Example 17
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)
0
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
0
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)
0
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)
0
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)
0
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
0
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
0
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)
0
Example 25
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)
0
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
0
Example 27
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
0
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
0
Example 29
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()
0
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")
0
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)
0
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
0
Example 33
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)
0
Example 34
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, :]
0
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
0
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()
0
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)
0
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
0
Example 39
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)
0
Example 40
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
0
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
0
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
0
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
0
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])
0
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
0
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
0
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,
}
0
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)
0
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"]
0
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