@@ -1221,3 +1221,226 @@ bsxfun(f, a, b) = f(a, b)
1221
1221
bsxfun (f, a:: AbstractArray , b) = f (a, b)
1222
1222
bsxfun (f, a, b:: AbstractArray ) = f (a, b)
1223
1223
bsxfun (f, a, b, c... ) = bsxfun (f, bsxfun (f, a, b), c... )
1224
+
1225
+ # Basic AbstractArray functions
1226
+
1227
+ max {T} (A:: AbstractArray{T} , b:: () , region) = reducedim (max,A,region,typemin (T))
1228
+ min {T} (A:: AbstractArray{T} , b:: () , region) = reducedim (min,A,region,typemax (T))
1229
+ sum {T} (A:: AbstractArray{T} , region) = reducedim (+ ,A,region,zero (T))
1230
+ prod {T} (A:: AbstractArray{T} , region) = reducedim (* ,A,region,one (T))
1231
+
1232
+ all (A:: AbstractArray{Bool} , region) = reducedim (all,A,region,true )
1233
+ any (A:: AbstractArray{Bool} , region) = reducedim (any,A,region,false )
1234
+ sum (A:: AbstractArray{Bool} , region) = reducedim (+ ,A,region,0 ,similar (A,Int,reduced_dims (A,region)))
1235
+ sum (A:: AbstractArray{Bool} ) = sum (A, [1 : ndims (A)])[1 ]
1236
+ prod (A:: AbstractArray{Bool} ) =
1237
+ error (" use all() instead of prod() for boolean arrays" )
1238
+ prod (A:: AbstractArray{Bool} , region) =
1239
+ error (" use all() instead of prod() for boolean arrays" )
1240
+
1241
+ function sum {T} (A:: AbstractArray{T} )
1242
+ if isempty (A)
1243
+ return zero (T)
1244
+ end
1245
+ v = A[1 ]
1246
+ for i= 2 : length (A)
1247
+ v += A[i]
1248
+ end
1249
+ v
1250
+ end
1251
+
1252
+ function sum_kbn {T<:FloatingPoint} (A:: AbstractArray{T} )
1253
+ n = length (A)
1254
+ if (n == 0 )
1255
+ return zero (T)
1256
+ end
1257
+ s = A[1 ]
1258
+ c = zero (T)
1259
+ for i in 2 : n
1260
+ Ai = A[i]
1261
+ t = s + Ai
1262
+ if abs (s) >= abs (Ai)
1263
+ c += ((s- t) + Ai)
1264
+ else
1265
+ c += ((Ai- t) + s)
1266
+ end
1267
+ s = t
1268
+ end
1269
+
1270
+ s + c
1271
+ end
1272
+
1273
+ # Uses K-B-N summation
1274
+ function cumsum_kbn {T<:FloatingPoint} (v:: AbstractVector{T} )
1275
+ n = length (v)
1276
+ r = similar (v, n)
1277
+ if n == 0 ; return r; end
1278
+
1279
+ s = r[1 ] = v[1 ]
1280
+ c = zero (T)
1281
+ for i= 2 : n
1282
+ vi = v[i]
1283
+ t = s + vi
1284
+ if abs (s) >= abs (vi)
1285
+ c += ((s- t) + vi)
1286
+ else
1287
+ c += ((vi- t) + s)
1288
+ end
1289
+ s = t
1290
+ r[i] = s+ c
1291
+ end
1292
+ return r
1293
+ end
1294
+
1295
+ # Uses K-B-N summation
1296
+ function cumsum_kbn {T<:FloatingPoint} (A:: AbstractArray{T} , axis:: Integer )
1297
+ dimsA = size (A)
1298
+ ndimsA = ndims (A)
1299
+ axis_size = dimsA[axis]
1300
+ axis_stride = 1
1301
+ for i = 1 : (axis- 1 )
1302
+ axis_stride *= size (A,i)
1303
+ end
1304
+
1305
+ if axis_size <= 1
1306
+ return A
1307
+ end
1308
+
1309
+ B = similar (A)
1310
+ C = similar (A)
1311
+
1312
+ for i = 1 : length (A)
1313
+ if div (i- 1 , axis_stride) % axis_size == 0
1314
+ B[i] = A[i]
1315
+ C[i] = zero (T)
1316
+ else
1317
+ s = B[i- axis_stride]
1318
+ Ai = A[i]
1319
+ B[i] = t = s + Ai
1320
+ if abs (s) >= abs (Ai)
1321
+ C[i] = C[i- axis_stride] + ((s- t) + Ai)
1322
+ else
1323
+ C[i] = C[i- axis_stride] + ((Ai- t) + s)
1324
+ end
1325
+ end
1326
+ end
1327
+
1328
+ return B + C
1329
+ end
1330
+
1331
+ function prod {T} (A:: AbstractArray{T} )
1332
+ if isempty (A)
1333
+ return one (T)
1334
+ end
1335
+ v = A[1 ]
1336
+ for i= 2 : length (A)
1337
+ v *= A[i]
1338
+ end
1339
+ v
1340
+ end
1341
+
1342
+ function min {T<:Integer} (A:: AbstractArray{T} )
1343
+ v = typemax (T)
1344
+ for i= 1 : length (A)
1345
+ x = A[i]
1346
+ if x < v
1347
+ v = x
1348
+ end
1349
+ end
1350
+ v
1351
+ end
1352
+
1353
+ function max {T<:Integer} (A:: AbstractArray{T} )
1354
+ v = typemin (T)
1355
+ for i= 1 : length (A)
1356
+ x = A[i]
1357
+ if x > v
1358
+ v = x
1359
+ end
1360
+ end
1361
+ v
1362
+ end
1363
+
1364
+ # # map over arrays ##
1365
+
1366
+ # # along an axis
1367
+ function amap (f:: Function , A:: AbstractArray , axis:: Integer )
1368
+ dimsA = size (A)
1369
+ ndimsA = ndims (A)
1370
+ axis_size = dimsA[axis]
1371
+
1372
+ if axis_size == 0
1373
+ return f (A)
1374
+ end
1375
+
1376
+ idx = ntuple (ndimsA, j -> j == axis ? 1 : 1 : dimsA[j])
1377
+ r = f (sub (A, idx))
1378
+ R = Array (typeof (r), axis_size)
1379
+ R[1 ] = r
1380
+
1381
+ for i = 2 : axis_size
1382
+ idx = ntuple (ndimsA, j -> j == axis ? i : 1 : dimsA[j])
1383
+ R[i] = f (sub (A, idx))
1384
+ end
1385
+
1386
+ return R
1387
+ end
1388
+
1389
+
1390
+ # # 1 argument
1391
+ function map_to2 (f, first, dest:: AbstractArray , A:: AbstractArray )
1392
+ dest[1 ] = first
1393
+ for i= 2 : length (A)
1394
+ dest[i] = f (A[i])
1395
+ end
1396
+ return dest
1397
+ end
1398
+
1399
+ function map (f, A:: AbstractArray )
1400
+ if isempty (A); return A; end
1401
+ first = f (A[1 ])
1402
+ dest = similar (A, typeof (first))
1403
+ return map_to2 (f, first, dest, A)
1404
+ end
1405
+
1406
+ # # 2 argument
1407
+ function map_to2 (f, first, dest:: AbstractArray , A:: AbstractArray , B:: AbstractArray )
1408
+ dest[1 ] = first
1409
+ for i= 2 : length (A)
1410
+ dest[i] = f (A[i], B[i])
1411
+ end
1412
+ return dest
1413
+ end
1414
+
1415
+ function map (f, A:: AbstractArray , B:: AbstractArray )
1416
+ shp = promote_shape (size (A),size (B))
1417
+ if isempty (A)
1418
+ return similar (A, eltype (A), shp)
1419
+ end
1420
+ first = f (A[1 ], B[1 ])
1421
+ dest = similar (A, typeof (first), shp)
1422
+ return map_to2 (f, first, dest, A, B)
1423
+ end
1424
+
1425
+ # # N argument
1426
+ function map_to2 (f, first, dest:: AbstractArray , As:: AbstractArray... )
1427
+ n = length (As[1 ])
1428
+ i = 1
1429
+ ith = a-> a[i]
1430
+ dest[1 ] = first
1431
+ for i= 2 : n
1432
+ dest[i] = f (map (ith, As)... )
1433
+ end
1434
+ return dest
1435
+ end
1436
+
1437
+ function map (f, As:: AbstractArray... )
1438
+ shape = mapreduce (size, promote_shape, As)
1439
+ if prod (shape) == 0
1440
+ return similar (As[1 ], eltype (As[1 ]), shape)
1441
+ end
1442
+ first = f (map (a-> a[1 ], As)... )
1443
+ dest = similar (As[1 ], typeof (first), shape)
1444
+ return map_to2 (f, first, dest, As... )
1445
+ end
1446
+
0 commit comments