Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

# -*- coding: utf-8 -*- 

r""" 

Basic Linear Algebra Subroutines on dictionaries 

  

This module provides functions for "level 1 basic linear algebra 

subroutines" on sparse vectors represented as dictionaries. The API is 

inspired from the standard BLAS API 

:wikipedia:`Basic_Linear_Algebra_Subprograms`, but does not try to 

follow it exactly. 

  

For ``a, b, ...`` hashable objects, the dictionary ``{a: 2, b:3, ...}`` 

represents the formal linear combination `2 \cdot a + 3 \cdot b + \cdots`. 

The values of the dictionary should all lie in the same parent `K`. For 

simplicity, we call this formal linear combination a vector, as the 

typical use case is `K` being a field. However the routines here can 

be used with a ring `K` to represent free module elements, or a 

semiring like `\NN` to represent elements of a commutative monoid, or 

even just an additive semigroup. Of course not all operations are 

meaningful in those cases. We are also assuming that ``-1 * x = -x`` 

and ``bool(x) == bool(-x)`` for all ``x`` in `K`. 

  

Unless stated overwise, all values `v` in the dictionaries should be 

non zero (as tested with `bool(v)`). 

  

This is mostly used by :class:`CombinatorialFreeModule`. 

""" 

#***************************************************************************** 

# Copyright (C) 2010 Christian Stump <christian.stump@univie.ac.at> 

# 2016 Travis Scrimshaw <tscrimsh@umn.edu> 

# 2016 Nicolas M. Thiéry <nthiery at users.sf.net> 

# 

# This program is free software: you can redistribute it and/or modify 

# it under the terms of the GNU General Public License as published by 

# the Free Software Foundation, either version 2 of the License, or 

# (at your option) any later version. 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

  

cpdef int iaxpy(a, dict X, dict Y, bint remove_zeros=True, bint factor_on_left=True) except -1: 

r""" 

Mutate `Y` to represent `a X + Y`. 

  

INPUT: 

  

- ``a`` -- element of a parent `K` or `±1` 

- ``X,Y`` -- dictionaries representing a vector `X` over `K` 

- ``remove_zeros`` -- boolean (default: ``True``); whether to 

remove the keys whose values are zero after the addition has 

been performed 

- ``factor_on_left`` -- boolean (default: ``False``); 

whether to compute `a X + Y` or `X a + Y` 

  

OUTPUT: 

  

``0`` when successful and ``-1`` on error. This is done because 

returning an ``int`` is faster than a (Python) ``None``. 

  

``Y`` has been mutated to represent the vector `a \cdot X + Y`. 

If ``remove_zeros=True`` (and the input have no zero values 

themselves), then the output is guaranteed to have no zero 

values. Otherwise some zero values may be left in the output. 

Which ones is voluntarily left undefined. Use this in cases 

where you want to postpone clearing until the end of a long 

series of operations. 

  

The parent `K` should support addition unless `a = -1`, negation 

if `a = -1`, and multiplication if `a \neq ±1`. We are also 

assuming that ``-1 * x = -x`` and ``bool(x) == bool(-x)`` in `K`. 

  

See :mod:`sage.data_structures.blas_dict` for an overview. 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: X = {0: -1, 1: 1} 

sage: Y = {0: 1, 1: 1} 

  

Computing `X+Y`:: 

  

sage: blas.iaxpy(1, X, Y) 

0 

sage: Y 

{1: 2} 

  

Reseting `Y` and computing `-X+Y`:: 

  

sage: Y = {0: 1, 1: 1} 

sage: blas.iaxpy(-1, X, Y) 

0 

sage: Y 

{0: 2} 

  

Reseting `Y` and computing `2X+Y`:: 

  

sage: Y = {0: 1, 1: 1} 

sage: blas.iaxpy(2, X, Y) 

0 

sage: Y 

{0: -1, 1: 3} 

  

Reseting `Y` and computing `-X+Y` without removing zeros:: 

  

sage: Y = {0: 1, 1: 1} 

sage: blas.iaxpy(-1, X, Y, remove_zeros=False) 

0 

sage: Y 

{0: 2, 1: 0} 

""" 

cdef int flag = 0 

if a == 1: 

flag = 1 

elif a == -1: 

flag = -1 

elif not a: 

return 0 

for (key, value) in X.iteritems(): 

if flag == -1: 

if key in Y: 

Y[key] -= value 

else: 

Y[key] = -value 

continue # no need to check for zero 

else: 

if flag != 1: 

# If we had the guarantee that `a` is in the base 

# ring, we could use a._mul_(value), as suggested in 

# the documentation of sage.structure.element. However 

# we want to support elements that can act on values 

# in the base ring as well 

if factor_on_left: 

value = a*value 

else: 

value = value*a 

if not value: 

continue # a is a zero divisor 

# assert value 

if key in Y: 

Y[key] += value 

else: 

Y[key] = value 

continue # no need to check for zero 

if remove_zeros and not Y[key]: 

del Y[key] 

return 0 

  

cpdef dict axpy(a, dict X, dict Y, bint factor_on_left=True): 

""" 

Return `a X + Y`. 

  

All input dictionaries are supposed to have values in the same 

base ring `K` and have no nonzero values. 

  

INPUT: 

  

- ``a`` -- an element of `K` or `±1` 

- ``X``, ``Y`` -- dictionaries representing respectively vectors `X` and `Y` 

- ``factor_on_left`` -- boolean (default: ``False``); 

whether to compute `a X + Y` or `X a + Y` 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: Y = {0: 1, 1: 1} 

sage: X = {0: -1, 1: 1} 

  

Computing `X + Y`:: 

  

sage: blas.axpy(1, X, Y) 

{1: 2} 

  

`X` and `Y` are not mutated: 

  

sage: X, Y 

({0: -1, 1: 1}, {0: 1, 1: 1}) 

  

Computing `-X + Y`:: 

  

sage: blas.axpy(-1, X, Y) 

{0: 2} 

  

Computing `2X + Y`:: 

  

sage: blas.axpy(2, X, Y) 

{0: -1, 1: 3} 

  

Testing the special cases of empty dictionaries:: 

  

sage: Y = {} 

sage: Z = blas.axpy(2, X, Y) 

sage: Z is Y 

False 

sage: X, Y 

({0: -1, 1: 1}, {}) 

  

sage: Z = blas.axpy(2, Y, X) 

sage: Z is X 

True 

sage: X, Y 

({0: -1, 1: 1}, {}) 

""" 

if X: 

Y = Y.copy() 

iaxpy(a, X, Y, True, factor_on_left) 

return Y 

  

cpdef dict negate(dict D): 

r""" 

Return a dictionary representing the vector `-X`. 

  

INPUT: 

  

- ``X`` -- a dictionary representing a vector `X` 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: D1 = {0: 1, 1: 1} 

sage: blas.negate(D1) 

{0: -1, 1: -1} 

""" 

return { key: -value for key, value in D.iteritems() } 

  

cpdef dict scal(a, dict D, bint factor_on_left=True): 

r""" 

Return a dictionary representing the vector `a*X`. 

  

INPUT: 

  

- ``a`` -- an element of the base ring `K` 

- ``X`` -- a dictionary representing a vector `X` 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: R = IntegerModRing(12) # a ring with zero divisors 

sage: D = {0: R(1), 1: R(2), 2: R(4)} 

sage: blas.scal(R(3), D) 

{0: 3, 1: 6} 

""" 

# We could use a dict comprehension like for negate, but care 

# needs to be taken to remove products that cancel out. 

# So for now we just delegate to axpy. 

return axpy(a, D, {}, factor_on_left=factor_on_left) 

  

cpdef dict add(dict D, dict D2): 

r""" 

Return the pointwise addition of dictionaries ``D`` and ``D2``. 

  

INPUT: 

  

- ``D``, ``D2`` -- dictionaries whose values are in a common ring 

and all values are non-zero 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: D = {0: 1, 1: 1} 

sage: D2 = {0: -1, 1: 1} 

sage: blas.add(D, D2) 

{1: 2} 

  

`D` and `D2` are left unchanged:: 

  

sage: D, D2 

({0: 1, 1: 1}, {0: -1, 1: 1}) 

""" 

# Optimization: swap the two dicts to ensure D is the largest 

if len(D) < len(D2): 

D, D2 = D2, D 

return axpy(1, D2, D) 

  

cpdef dict sum(dict_iter): 

r""" 

Return the pointwise addition of dictionaries with coefficients. 

  

INPUT: 

  

- ``dict_iter`` -- iterator of dictionaries whose values are in 

a common ring and all values are non-zero 

  

OUTPUT: 

  

- a dictionary containing all keys of the dictionaries in ``dict_list`` 

with values being the sum of the values in the different dictionaries 

(keys with zero value are omitted) 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: D = {0: 1, 1: 1}; D 

{0: 1, 1: 1} 

sage: blas.sum(D for x in range(5)) 

{0: 5, 1: 5} 

  

sage: D1 = {0: 1, 1: 1}; D2 = {0: -1, 1: 1} 

sage: blas.sum([D1, D2]) 

{1: 2} 

  

sage: blas.sum([{}, {}, D1, D2, {}]) 

{1: 2} 

""" 

cdef dict result = {} 

cdef D 

cdef list for_removal 

  

for D in dict_iter: 

if result: 

iaxpy(1, D, result, remove_zeros=False) 

elif D: 

result = D.copy() 

  

for_removal = [key for key in result if not result[key]] 

for key in for_removal: 

del result[key] 

return result 

  

cpdef dict linear_combination(dict_factor_iter, bint factor_on_left=True): 

r""" 

Return the pointwise addition of dictionaries with coefficients. 

  

INPUT: 

  

- ``dict_factor_iter`` -- iterator of pairs ``D``, ``coeff``, where 

  

* the ``D``'s are dictionaries with values in a common ring 

* the ``coeff``'s are coefficients in this ring 

  

- ``factor_on_left`` -- boolean (default: ``True``); if ``True``, 

the coefficients are multiplied on the left, otherwise they are 

multiplied on the right 

  

OUTPUT: 

  

- a dictionary containing all keys of dictionaries in ``dict_list`` 

with values being the sum of the values in the different 

dictionaries and each one first multiplied by the given factor 

(keys with zero value are omitted) 

  

EXAMPLES:: 

  

sage: import sage.data_structures.blas_dict as blas 

sage: D = { 0:1, 1:1 }; D 

{0: 1, 1: 1} 

sage: blas.linear_combination( (D,i) for i in range(5) ) 

{0: 10, 1: 10} 

sage: blas.linear_combination( [(D,1), (D,-1)] ) 

{} 

""" 

cdef dict result = {} 

cdef dict D 

cdef list for_removal 

  

for D, a in dict_factor_iter: 

if not a: # We multiply by 0, so nothing to do 

continue 

if not result and a == 1: 

result = D.copy() 

else: 

iaxpy(a, D, result, remove_zeros=False) 

  

for_removal = [key for key in result if not result[key]] 

for key in for_removal: 

del result[key] 

  

return result