Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Problems.GaussianProcessTuning/ILNumerics.2.14.4735.573/Functions/builtin/distL1.cs @ 11194

Last change on this file since 11194 was 9102, checked in by gkronber, 12 years ago

#1967: ILNumerics source for experimentation
File size: 68.2 KB
Line 
1///
2///    This file is part of ILNumerics Community Edition.
3///
4///    ILNumerics Community Edition - high performance computing for applications.
5///    Copyright (C) 2006 - 2012 Haymo Kutschbach, http://ilnumerics.net
6///
7///    ILNumerics Community Edition is free software: you can redistribute it and/or modify
8///    it under the terms of the GNU General Public License version 3 as published by
9///    the Free Software Foundation.
10///
11///    ILNumerics Community Edition is distributed in the hope that it will be useful,
12///    but WITHOUT ANY WARRANTY; without even the implied warranty of
13///    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14///    GNU General Public License for more details.
15///
16///    You should have received a copy of the GNU General Public License
17///    along with ILNumerics Community Edition. See the file License.txt in the root
18///    of your distribution package. If not, see <http://www.gnu.org/licenses/>.
19///
20///    In addition this software uses the following components and/or licenses:
21///
22///    =================================================================================
23///    The Open Toolkit Library License
24///   
25///    Copyright (c) 2006 - 2009 the Open Toolkit library.
26///   
27///    Permission is hereby granted, free of charge, to any person obtaining a copy
28///    of this software and associated documentation files (the "Software"), to deal
29///    in the Software without restriction, including without limitation the rights to
30///    use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
31///    the Software, and to permit persons to whom the Software is furnished to do
32///    so, subject to the following conditions:
33///
34///    The above copyright notice and this permission notice shall be included in all
35///    copies or substantial portions of the Software.
36///
37///    =================================================================================
38///   
39
40using System;
41using ILNumerics;
42using ILNumerics.Misc;
43using ILNumerics.Storage;
44using ILNumerics.Native;
45using ILNumerics.Exceptions;
46
47
48
49namespace ILNumerics {
50
51    public partial class ILMath {
52
53
54
55        /// <summary>
56        /// pairwise L1 distance
57        /// </summary>
58        /// <param name="A">input points (matrix)</param>
59        /// <param name="B">input point (vector)</param>
60        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
61        /// <remarks>
62        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
63        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
64        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
65        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
66        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
67        /// </remarks>
68        public static unsafe ILRetArray<double> distL1(ILInArray<double> A, ILInArray<double> B) {
69            using (ILScope.Enter(A, B)) {
70
71                #region parameter checking
72                if (isnull(A) || isnull(B))
73                    return empty<double>(ILSize.Empty00);
74                if (A.IsEmpty) {
75                    return empty<double>(B.S);
76                } else if (B.IsEmpty) {
77                    return empty<double>(A.S);
78                }
79                // early exit: make the function cummutative
80                if (A.IsVector && !B.IsVector) {
81                    return distL1(B, A);
82                }
83                int dim = -1;
84                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
85                    if (A.S[l] != B.S[l]) {
86                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
87                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
88                        }
89                        dim = l;
90                    }
91                }
92               
93                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
94                #endregion
95
96                #region parameter preparation
97                ILSize outDims;
98                if (dim == 0) {
99                    outDims = size(1,A.S[1]);   
100                } else if (dim == 1) {
101                    outDims = size(A.S[0], 1);
102                } else if (A.S.IsSameSize(B.S)) {
103                    outDims = A.S;
104                    dim = 0;
105                } else
106                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
107
108               
109                double[] retArr = ILMemoryPool.Pool.New<double>(outDims.NumberOfElements);
110
111               
112                double[] arrA = A.GetArrayForRead();
113
114               
115                double[] arrB = B.GetArrayForRead();
116                int itLen = A.S[0];
117                #endregion
118
119                #region worker loops definition
120                ILDenseStorage<double> retStorage = new ILDenseStorage<double>(retArr, outDims);
121                int workerCount = 1;
122                Action<object> worker = data => {
123                    // expects: iStart, iLen, ap, bp, cp
124                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
125                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
126
127                   
128                    double* ap;
129
130                   
131                    double* bp;
132
133                   
134                    double* cp;
135                    if (dim == 0) {
136                        ap = (double*)range.Item2;
137                        cp = (double*)range.Item4;
138                        for (int s = 0; s < range.Item1; s++) {
139                            bp = (double*)range.Item3;
140                           
141                            double sum = 0;
142                            int l = itLen;
143                            while (l > 20) {
144                                sum +=  Math.Abs(ap[0] - bp[0])
145                                    +   Math.Abs(ap[1] - bp[1])
146                                    +   Math.Abs(ap[2] - bp[2])
147                                    +   Math.Abs(ap[3] - bp[3])
148                                    +   Math.Abs(ap[4] - bp[4])
149                                    +   Math.Abs(ap[5] - bp[5])
150                                    +   Math.Abs(ap[6] - bp[6])
151                                    +   Math.Abs(ap[7] - bp[7])
152                                    +   Math.Abs(ap[8] - bp[8])
153                                    +   Math.Abs(ap[9] - bp[9])
154                                    +   Math.Abs(ap[10] - bp[10])
155                                    +   Math.Abs(ap[11] - bp[11])
156                                    +   Math.Abs(ap[12] - bp[12])
157                                    +   Math.Abs(ap[13] - bp[13])
158                                    +   Math.Abs(ap[14] - bp[14])
159                                    +   Math.Abs(ap[15] - bp[15])
160                                    +   Math.Abs(ap[16] - bp[16])
161                                    +   Math.Abs(ap[17] - bp[17])
162                                    +   Math.Abs(ap[18] - bp[18])
163                                    +   Math.Abs(ap[19] - bp[19])
164                                    +   Math.Abs(ap[20] - bp[20]);
165                                ap += 21;
166                                bp += 21;
167                                l -= 21;
168                            }
169                            while (l-- > 0) {
170                                sum +=  Math.Abs(*ap - *bp);
171                                ap++;
172                                bp++;
173                            }
174                            *cp++ = sum;
175                        }
176                    } else {
177                        // dim = 1
178                        ap = (double*)range.Item2;
179                        bp = (double*)range.Item3;
180                        cp = (double*)range.Item4;
181                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
182                           
183                            double val = *bp++;
184                            int l = itLen;
185                            while (l > 10) {
186                                cp[0] +=   Math.Abs(ap[0] - val);
187                                cp[1] +=   Math.Abs(ap[1] - val);
188                                cp[2] +=   Math.Abs(ap[2] - val);
189                                cp[3] +=   Math.Abs(ap[3] - val);
190                                cp[4] +=   Math.Abs(ap[4] - val);
191                                cp[5] +=   Math.Abs(ap[5] - val);
192                                cp[6] +=   Math.Abs(ap[6] - val);
193                                cp[7] +=   Math.Abs(ap[7] - val);
194                                cp[8] +=   Math.Abs(ap[8] - val);
195                                cp[9] +=   Math.Abs(ap[9] - val);
196                                cp[10] +=  Math.Abs(ap[10] - val);
197                                cp[11] +=  Math.Abs(ap[11] - val);
198                                cp[12] +=  Math.Abs(ap[12] - val);
199                                cp[13] +=  Math.Abs(ap[13] - val);
200                                cp[14] +=  Math.Abs(ap[14] - val);
201                                cp[15] +=  Math.Abs(ap[15] - val);
202                                cp[16] +=  Math.Abs(ap[16] - val);
203                                cp[17] +=  Math.Abs(ap[17] - val);
204                                cp[18] +=  Math.Abs(ap[18] - val);
205                                cp[19] +=  Math.Abs(ap[19] - val);
206                                cp[20] +=  Math.Abs(ap[20] - val);
207                                ap += 21;
208                                cp += 21;
209                                l -= 21;
210                            }
211                            while (l-- > 0) {
212                                *cp +=  Math.Abs(*ap - val);
213                                ap++;
214                                cp++;
215                            }
216                        }
217                    }
218                    System.Threading.Interlocked.Decrement(ref workerCount);
219                };
220                #endregion
221
222                #region work distribution
223                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
224                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
225                    && outDims[1] > 1) {
226                        if (outDims[1] > workItemCount) {
227                        workItemLength = outDims[1] / workItemCount;
228                    } else {
229                        workItemLength = outDims[1] / 2;
230                        workItemCount = 2;
231                    }
232                } else {
233                    workItemLength = outDims[1];
234                    workItemCount = 1;
235                }
236
237                //int[] partResults = new int[workItemCount];
238                fixed ( double* arrAP = arrA)
239                fixed ( double* arrBP = arrB)
240                fixed ( double* retArrP = retArr) {
241
242                    for (; i < workItemCount - 1; i++) {
243                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
244                               (workItemLength
245                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
246                               , (IntPtr)(arrBP + i * dim * workItemLength)
247                               , (IntPtr)(retArrP + i * workItemLength));
248                        System.Threading.Interlocked.Increment(ref workerCount);
249                        ILThreadPool.QueueUserWorkItem(i, worker, range);
250                    }
251                    // the last (or may the only) chunk is done right here
252                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
253                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
254                                (outDims[1] - i * workItemLength
255                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
256                                , (IntPtr)(arrBP + i * dim * workItemLength)
257                                , (IntPtr)(retArrP + i * workItemLength)));
258
259                    ILThreadPool.Wait4Workers(ref workerCount);
260                }
261                #endregion
262
263                return new ILRetArray<double>(retStorage);
264            }
265        }
266
267#region HYCALPER AUTO GENERATED CODE
268
269        /// <summary>
270        /// pairwise L1 distance
271        /// </summary>
272        /// <param name="A">input points (matrix)</param>
273        /// <param name="B">input point (vector)</param>
274        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
275        /// <remarks>
276        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
277        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
278        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
279        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
280        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
281        /// </remarks>
282        public static unsafe ILRetArray<Int64> distL1(ILInArray<Int64> A, ILInArray<Int64> B) {
283            using (ILScope.Enter(A, B)) {
284
285                #region parameter checking
286                if (isnull(A) || isnull(B))
287                    return empty<Int64>(ILSize.Empty00);
288                if (A.IsEmpty) {
289                    return empty<Int64>(B.S);
290                } else if (B.IsEmpty) {
291                    return empty<Int64>(A.S);
292                }
293                // early exit: make the function cummutative
294                if (A.IsVector && !B.IsVector) {
295                    return distL1(B, A);
296                }
297                int dim = -1;
298                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
299                    if (A.S[l] != B.S[l]) {
300                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
301                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
302                        }
303                        dim = l;
304                    }
305                }
306               
307                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
308                #endregion
309
310                #region parameter preparation
311                ILSize outDims;
312                if (dim == 0) {
313                    outDims = size(1,A.S[1]);   
314                } else if (dim == 1) {
315                    outDims = size(A.S[0], 1);
316                } else if (A.S.IsSameSize(B.S)) {
317                    outDims = A.S;
318                    dim = 0;
319                } else
320                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
321
322               
323                Int64[] retArr = ILMemoryPool.Pool.New<Int64>(outDims.NumberOfElements);
324
325               
326                Int64[] arrA = A.GetArrayForRead();
327
328               
329                Int64[] arrB = B.GetArrayForRead();
330                int itLen = A.S[0];
331                #endregion
332
333                #region worker loops definition
334                ILDenseStorage<Int64> retStorage = new ILDenseStorage<Int64>(retArr, outDims);
335                int workerCount = 1;
336                Action<object> worker = data => {
337                    // expects: iStart, iLen, ap, bp, cp
338                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
339                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
340
341                   
342                    Int64* ap;
343
344                   
345                    Int64* bp;
346
347                   
348                    Int64* cp;
349                    if (dim == 0) {
350                        ap = (Int64*)range.Item2;
351                        cp = (Int64*)range.Item4;
352                        for (int s = 0; s < range.Item1; s++) {
353                            bp = (Int64*)range.Item3;
354                           
355                            Int64 sum = 0;
356                            int l = itLen;
357                            while (l > 20) {
358                                sum +=  Math.Abs(ap[0] - bp[0])
359                                    +  Math.Abs(ap[1] - bp[1])
360                                    +  Math.Abs(ap[2] - bp[2])
361                                    +  Math.Abs(ap[3] - bp[3])
362                                    +  Math.Abs(ap[4] - bp[4])
363                                    +  Math.Abs(ap[5] - bp[5])
364                                    +  Math.Abs(ap[6] - bp[6])
365                                    +  Math.Abs(ap[7] - bp[7])
366                                    +  Math.Abs(ap[8] - bp[8])
367                                    +  Math.Abs(ap[9] - bp[9])
368                                    +  Math.Abs(ap[10] - bp[10])
369                                    +  Math.Abs(ap[11] - bp[11])
370                                    +  Math.Abs(ap[12] - bp[12])
371                                    +  Math.Abs(ap[13] - bp[13])
372                                    +  Math.Abs(ap[14] - bp[14])
373                                    +  Math.Abs(ap[15] - bp[15])
374                                    +  Math.Abs(ap[16] - bp[16])
375                                    +  Math.Abs(ap[17] - bp[17])
376                                    +  Math.Abs(ap[18] - bp[18])
377                                    +  Math.Abs(ap[19] - bp[19])
378                                    +  Math.Abs(ap[20] - bp[20]);
379                                ap += 21;
380                                bp += 21;
381                                l -= 21;
382                            }
383                            while (l-- > 0) {
384                                sum +=  Math.Abs(*ap - *bp);
385                                ap++;
386                                bp++;
387                            }
388                            *cp++ = sum;
389                        }
390                    } else {
391                        // dim = 1
392                        ap = (Int64*)range.Item2;
393                        bp = (Int64*)range.Item3;
394                        cp = (Int64*)range.Item4;
395                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
396                           
397                            Int64 val = *bp++;
398                            int l = itLen;
399                            while (l > 10) {
400                                cp[0] +=  Math.Abs(ap[0] - val);
401                                cp[1] +=  Math.Abs(ap[1] - val);
402                                cp[2] +=  Math.Abs(ap[2] - val);
403                                cp[3] +=  Math.Abs(ap[3] - val);
404                                cp[4] +=  Math.Abs(ap[4] - val);
405                                cp[5] +=  Math.Abs(ap[5] - val);
406                                cp[6] +=  Math.Abs(ap[6] - val);
407                                cp[7] +=  Math.Abs(ap[7] - val);
408                                cp[8] +=  Math.Abs(ap[8] - val);
409                                cp[9] +=  Math.Abs(ap[9] - val);
410                                cp[10] +=  Math.Abs(ap[10] - val);
411                                cp[11] +=  Math.Abs(ap[11] - val);
412                                cp[12] +=  Math.Abs(ap[12] - val);
413                                cp[13] +=  Math.Abs(ap[13] - val);
414                                cp[14] +=  Math.Abs(ap[14] - val);
415                                cp[15] +=  Math.Abs(ap[15] - val);
416                                cp[16] +=  Math.Abs(ap[16] - val);
417                                cp[17] +=  Math.Abs(ap[17] - val);
418                                cp[18] +=  Math.Abs(ap[18] - val);
419                                cp[19] +=  Math.Abs(ap[19] - val);
420                                cp[20] +=  Math.Abs(ap[20] - val);
421                                ap += 21;
422                                cp += 21;
423                                l -= 21;
424                            }
425                            while (l-- > 0) {
426                                *cp +=  Math.Abs(*ap - val);
427                                ap++;
428                                cp++;
429                            }
430                        }
431                    }
432                    System.Threading.Interlocked.Decrement(ref workerCount);
433                };
434                #endregion
435
436                #region work distribution
437                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
438                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
439                    && outDims[1] > 1) {
440                        if (outDims[1] > workItemCount) {
441                        workItemLength = outDims[1] / workItemCount;
442                    } else {
443                        workItemLength = outDims[1] / 2;
444                        workItemCount = 2;
445                    }
446                } else {
447                    workItemLength = outDims[1];
448                    workItemCount = 1;
449                }
450
451                //int[] partResults = new int[workItemCount];
452                fixed ( Int64* arrAP = arrA)
453                fixed ( Int64* arrBP = arrB)
454                fixed ( Int64* retArrP = retArr) {
455
456                    for (; i < workItemCount - 1; i++) {
457                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
458                               (workItemLength
459                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
460                               , (IntPtr)(arrBP + i * dim * workItemLength)
461                               , (IntPtr)(retArrP + i * workItemLength));
462                        System.Threading.Interlocked.Increment(ref workerCount);
463                        ILThreadPool.QueueUserWorkItem(i, worker, range);
464                    }
465                    // the last (or may the only) chunk is done right here
466                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
467                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
468                                (outDims[1] - i * workItemLength
469                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
470                                , (IntPtr)(arrBP + i * dim * workItemLength)
471                                , (IntPtr)(retArrP + i * workItemLength)));
472
473                    ILThreadPool.Wait4Workers(ref workerCount);
474                }
475                #endregion
476
477                return new ILRetArray<Int64>(retStorage);
478            }
479        }
480        /// <summary>
481        /// pairwise L1 distance
482        /// </summary>
483        /// <param name="A">input points (matrix)</param>
484        /// <param name="B">input point (vector)</param>
485        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
486        /// <remarks>
487        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
488        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
489        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
490        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
491        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
492        /// </remarks>
493        public static unsafe ILRetArray<Int32> distL1(ILInArray<Int32> A, ILInArray<Int32> B) {
494            using (ILScope.Enter(A, B)) {
495
496                #region parameter checking
497                if (isnull(A) || isnull(B))
498                    return empty<Int32>(ILSize.Empty00);
499                if (A.IsEmpty) {
500                    return empty<Int32>(B.S);
501                } else if (B.IsEmpty) {
502                    return empty<Int32>(A.S);
503                }
504                // early exit: make the function cummutative
505                if (A.IsVector && !B.IsVector) {
506                    return distL1(B, A);
507                }
508                int dim = -1;
509                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
510                    if (A.S[l] != B.S[l]) {
511                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
512                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
513                        }
514                        dim = l;
515                    }
516                }
517               
518                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
519                #endregion
520
521                #region parameter preparation
522                ILSize outDims;
523                if (dim == 0) {
524                    outDims = size(1,A.S[1]);   
525                } else if (dim == 1) {
526                    outDims = size(A.S[0], 1);
527                } else if (A.S.IsSameSize(B.S)) {
528                    outDims = A.S;
529                    dim = 0;
530                } else
531                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
532
533               
534                Int32[] retArr = ILMemoryPool.Pool.New<Int32>(outDims.NumberOfElements);
535
536               
537                Int32[] arrA = A.GetArrayForRead();
538
539               
540                Int32[] arrB = B.GetArrayForRead();
541                int itLen = A.S[0];
542                #endregion
543
544                #region worker loops definition
545                ILDenseStorage<Int32> retStorage = new ILDenseStorage<Int32>(retArr, outDims);
546                int workerCount = 1;
547                Action<object> worker = data => {
548                    // expects: iStart, iLen, ap, bp, cp
549                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
550                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
551
552                   
553                    Int32* ap;
554
555                   
556                    Int32* bp;
557
558                   
559                    Int32* cp;
560                    if (dim == 0) {
561                        ap = (Int32*)range.Item2;
562                        cp = (Int32*)range.Item4;
563                        for (int s = 0; s < range.Item1; s++) {
564                            bp = (Int32*)range.Item3;
565                           
566                            Int32 sum = 0;
567                            int l = itLen;
568                            while (l > 20) {
569                                sum +=  Math.Abs(ap[0] - bp[0])
570                                    +  Math.Abs(ap[1] - bp[1])
571                                    +  Math.Abs(ap[2] - bp[2])
572                                    +  Math.Abs(ap[3] - bp[3])
573                                    +  Math.Abs(ap[4] - bp[4])
574                                    +  Math.Abs(ap[5] - bp[5])
575                                    +  Math.Abs(ap[6] - bp[6])
576                                    +  Math.Abs(ap[7] - bp[7])
577                                    +  Math.Abs(ap[8] - bp[8])
578                                    +  Math.Abs(ap[9] - bp[9])
579                                    +  Math.Abs(ap[10] - bp[10])
580                                    +  Math.Abs(ap[11] - bp[11])
581                                    +  Math.Abs(ap[12] - bp[12])
582                                    +  Math.Abs(ap[13] - bp[13])
583                                    +  Math.Abs(ap[14] - bp[14])
584                                    +  Math.Abs(ap[15] - bp[15])
585                                    +  Math.Abs(ap[16] - bp[16])
586                                    +  Math.Abs(ap[17] - bp[17])
587                                    +  Math.Abs(ap[18] - bp[18])
588                                    +  Math.Abs(ap[19] - bp[19])
589                                    +  Math.Abs(ap[20] - bp[20]);
590                                ap += 21;
591                                bp += 21;
592                                l -= 21;
593                            }
594                            while (l-- > 0) {
595                                sum +=  Math.Abs(*ap - *bp);
596                                ap++;
597                                bp++;
598                            }
599                            *cp++ = sum;
600                        }
601                    } else {
602                        // dim = 1
603                        ap = (Int32*)range.Item2;
604                        bp = (Int32*)range.Item3;
605                        cp = (Int32*)range.Item4;
606                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
607                           
608                            Int32 val = *bp++;
609                            int l = itLen;
610                            while (l > 10) {
611                                cp[0] +=  Math.Abs(ap[0] - val);
612                                cp[1] +=  Math.Abs(ap[1] - val);
613                                cp[2] +=  Math.Abs(ap[2] - val);
614                                cp[3] +=  Math.Abs(ap[3] - val);
615                                cp[4] +=  Math.Abs(ap[4] - val);
616                                cp[5] +=  Math.Abs(ap[5] - val);
617                                cp[6] +=  Math.Abs(ap[6] - val);
618                                cp[7] +=  Math.Abs(ap[7] - val);
619                                cp[8] +=  Math.Abs(ap[8] - val);
620                                cp[9] +=  Math.Abs(ap[9] - val);
621                                cp[10] +=  Math.Abs(ap[10] - val);
622                                cp[11] +=  Math.Abs(ap[11] - val);
623                                cp[12] +=  Math.Abs(ap[12] - val);
624                                cp[13] +=  Math.Abs(ap[13] - val);
625                                cp[14] +=  Math.Abs(ap[14] - val);
626                                cp[15] +=  Math.Abs(ap[15] - val);
627                                cp[16] +=  Math.Abs(ap[16] - val);
628                                cp[17] +=  Math.Abs(ap[17] - val);
629                                cp[18] +=  Math.Abs(ap[18] - val);
630                                cp[19] +=  Math.Abs(ap[19] - val);
631                                cp[20] +=  Math.Abs(ap[20] - val);
632                                ap += 21;
633                                cp += 21;
634                                l -= 21;
635                            }
636                            while (l-- > 0) {
637                                *cp +=  Math.Abs(*ap - val);
638                                ap++;
639                                cp++;
640                            }
641                        }
642                    }
643                    System.Threading.Interlocked.Decrement(ref workerCount);
644                };
645                #endregion
646
647                #region work distribution
648                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
649                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
650                    && outDims[1] > 1) {
651                        if (outDims[1] > workItemCount) {
652                        workItemLength = outDims[1] / workItemCount;
653                    } else {
654                        workItemLength = outDims[1] / 2;
655                        workItemCount = 2;
656                    }
657                } else {
658                    workItemLength = outDims[1];
659                    workItemCount = 1;
660                }
661
662                //int[] partResults = new int[workItemCount];
663                fixed ( Int32* arrAP = arrA)
664                fixed ( Int32* arrBP = arrB)
665                fixed ( Int32* retArrP = retArr) {
666
667                    for (; i < workItemCount - 1; i++) {
668                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
669                               (workItemLength
670                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
671                               , (IntPtr)(arrBP + i * dim * workItemLength)
672                               , (IntPtr)(retArrP + i * workItemLength));
673                        System.Threading.Interlocked.Increment(ref workerCount);
674                        ILThreadPool.QueueUserWorkItem(i, worker, range);
675                    }
676                    // the last (or may the only) chunk is done right here
677                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
678                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
679                                (outDims[1] - i * workItemLength
680                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
681                                , (IntPtr)(arrBP + i * dim * workItemLength)
682                                , (IntPtr)(retArrP + i * workItemLength)));
683
684                    ILThreadPool.Wait4Workers(ref workerCount);
685                }
686                #endregion
687
688                return new ILRetArray<Int32>(retStorage);
689            }
690        }
691        /// <summary>
692        /// pairwise L1 distance
693        /// </summary>
694        /// <param name="A">input points (matrix)</param>
695        /// <param name="B">input point (vector)</param>
696        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
697        /// <remarks>
698        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
699        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
700        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
701        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
702        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
703        /// </remarks>
704        public static unsafe ILRetArray<float> distL1(ILInArray<float> A, ILInArray<float> B) {
705            using (ILScope.Enter(A, B)) {
706
707                #region parameter checking
708                if (isnull(A) || isnull(B))
709                    return empty<float>(ILSize.Empty00);
710                if (A.IsEmpty) {
711                    return empty<float>(B.S);
712                } else if (B.IsEmpty) {
713                    return empty<float>(A.S);
714                }
715                // early exit: make the function cummutative
716                if (A.IsVector && !B.IsVector) {
717                    return distL1(B, A);
718                }
719                int dim = -1;
720                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
721                    if (A.S[l] != B.S[l]) {
722                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
723                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
724                        }
725                        dim = l;
726                    }
727                }
728               
729                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
730                #endregion
731
732                #region parameter preparation
733                ILSize outDims;
734                if (dim == 0) {
735                    outDims = size(1,A.S[1]);   
736                } else if (dim == 1) {
737                    outDims = size(A.S[0], 1);
738                } else if (A.S.IsSameSize(B.S)) {
739                    outDims = A.S;
740                    dim = 0;
741                } else
742                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
743
744               
745                float[] retArr = ILMemoryPool.Pool.New<float>(outDims.NumberOfElements);
746
747               
748                float[] arrA = A.GetArrayForRead();
749
750               
751                float[] arrB = B.GetArrayForRead();
752                int itLen = A.S[0];
753                #endregion
754
755                #region worker loops definition
756                ILDenseStorage<float> retStorage = new ILDenseStorage<float>(retArr, outDims);
757                int workerCount = 1;
758                Action<object> worker = data => {
759                    // expects: iStart, iLen, ap, bp, cp
760                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
761                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
762
763                   
764                    float* ap;
765
766                   
767                    float* bp;
768
769                   
770                    float* cp;
771                    if (dim == 0) {
772                        ap = (float*)range.Item2;
773                        cp = (float*)range.Item4;
774                        for (int s = 0; s < range.Item1; s++) {
775                            bp = (float*)range.Item3;
776                           
777                            float sum = 0;
778                            int l = itLen;
779                            while (l > 20) {
780                                sum +=  Math.Abs(ap[0] - bp[0])
781                                    +  Math.Abs(ap[1] - bp[1])
782                                    +  Math.Abs(ap[2] - bp[2])
783                                    +  Math.Abs(ap[3] - bp[3])
784                                    +  Math.Abs(ap[4] - bp[4])
785                                    +  Math.Abs(ap[5] - bp[5])
786                                    +  Math.Abs(ap[6] - bp[6])
787                                    +  Math.Abs(ap[7] - bp[7])
788                                    +  Math.Abs(ap[8] - bp[8])
789                                    +  Math.Abs(ap[9] - bp[9])
790                                    +  Math.Abs(ap[10] - bp[10])
791                                    +  Math.Abs(ap[11] - bp[11])
792                                    +  Math.Abs(ap[12] - bp[12])
793                                    +  Math.Abs(ap[13] - bp[13])
794                                    +  Math.Abs(ap[14] - bp[14])
795                                    +  Math.Abs(ap[15] - bp[15])
796                                    +  Math.Abs(ap[16] - bp[16])
797                                    +  Math.Abs(ap[17] - bp[17])
798                                    +  Math.Abs(ap[18] - bp[18])
799                                    +  Math.Abs(ap[19] - bp[19])
800                                    +  Math.Abs(ap[20] - bp[20]);
801                                ap += 21;
802                                bp += 21;
803                                l -= 21;
804                            }
805                            while (l-- > 0) {
806                                sum +=  Math.Abs(*ap - *bp);
807                                ap++;
808                                bp++;
809                            }
810                            *cp++ = sum;
811                        }
812                    } else {
813                        // dim = 1
814                        ap = (float*)range.Item2;
815                        bp = (float*)range.Item3;
816                        cp = (float*)range.Item4;
817                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
818                           
819                            float val = *bp++;
820                            int l = itLen;
821                            while (l > 10) {
822                                cp[0] +=  Math.Abs(ap[0] - val);
823                                cp[1] +=  Math.Abs(ap[1] - val);
824                                cp[2] +=  Math.Abs(ap[2] - val);
825                                cp[3] +=  Math.Abs(ap[3] - val);
826                                cp[4] +=  Math.Abs(ap[4] - val);
827                                cp[5] +=  Math.Abs(ap[5] - val);
828                                cp[6] +=  Math.Abs(ap[6] - val);
829                                cp[7] +=  Math.Abs(ap[7] - val);
830                                cp[8] +=  Math.Abs(ap[8] - val);
831                                cp[9] +=  Math.Abs(ap[9] - val);
832                                cp[10] +=  Math.Abs(ap[10] - val);
833                                cp[11] +=  Math.Abs(ap[11] - val);
834                                cp[12] +=  Math.Abs(ap[12] - val);
835                                cp[13] +=  Math.Abs(ap[13] - val);
836                                cp[14] +=  Math.Abs(ap[14] - val);
837                                cp[15] +=  Math.Abs(ap[15] - val);
838                                cp[16] +=  Math.Abs(ap[16] - val);
839                                cp[17] +=  Math.Abs(ap[17] - val);
840                                cp[18] +=  Math.Abs(ap[18] - val);
841                                cp[19] +=  Math.Abs(ap[19] - val);
842                                cp[20] +=  Math.Abs(ap[20] - val);
843                                ap += 21;
844                                cp += 21;
845                                l -= 21;
846                            }
847                            while (l-- > 0) {
848                                *cp +=  Math.Abs(*ap - val);
849                                ap++;
850                                cp++;
851                            }
852                        }
853                    }
854                    System.Threading.Interlocked.Decrement(ref workerCount);
855                };
856                #endregion
857
858                #region work distribution
859                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
860                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
861                    && outDims[1] > 1) {
862                        if (outDims[1] > workItemCount) {
863                        workItemLength = outDims[1] / workItemCount;
864                    } else {
865                        workItemLength = outDims[1] / 2;
866                        workItemCount = 2;
867                    }
868                } else {
869                    workItemLength = outDims[1];
870                    workItemCount = 1;
871                }
872
873                //int[] partResults = new int[workItemCount];
874                fixed ( float* arrAP = arrA)
875                fixed ( float* arrBP = arrB)
876                fixed ( float* retArrP = retArr) {
877
878                    for (; i < workItemCount - 1; i++) {
879                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
880                               (workItemLength
881                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
882                               , (IntPtr)(arrBP + i * dim * workItemLength)
883                               , (IntPtr)(retArrP + i * workItemLength));
884                        System.Threading.Interlocked.Increment(ref workerCount);
885                        ILThreadPool.QueueUserWorkItem(i, worker, range);
886                    }
887                    // the last (or may the only) chunk is done right here
888                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
889                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
890                                (outDims[1] - i * workItemLength
891                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
892                                , (IntPtr)(arrBP + i * dim * workItemLength)
893                                , (IntPtr)(retArrP + i * workItemLength)));
894
895                    ILThreadPool.Wait4Workers(ref workerCount);
896                }
897                #endregion
898
899                return new ILRetArray<float>(retStorage);
900            }
901        }
902        /// <summary>
903        /// pairwise L1 distance
904        /// </summary>
905        /// <param name="A">input points (matrix)</param>
906        /// <param name="B">input point (vector)</param>
907        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
908        /// <remarks>
909        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
910        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
911        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
912        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
913        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
914        /// </remarks>
915        public static unsafe ILRetArray<fcomplex> distL1(ILInArray<fcomplex> A, ILInArray<fcomplex> B) {
916            using (ILScope.Enter(A, B)) {
917
918                #region parameter checking
919                if (isnull(A) || isnull(B))
920                    return empty<fcomplex>(ILSize.Empty00);
921                if (A.IsEmpty) {
922                    return empty<fcomplex>(B.S);
923                } else if (B.IsEmpty) {
924                    return empty<fcomplex>(A.S);
925                }
926                // early exit: make the function cummutative
927                if (A.IsVector && !B.IsVector) {
928                    return distL1(B, A);
929                }
930                int dim = -1;
931                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
932                    if (A.S[l] != B.S[l]) {
933                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
934                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
935                        }
936                        dim = l;
937                    }
938                }
939               
940                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
941                #endregion
942
943                #region parameter preparation
944                ILSize outDims;
945                if (dim == 0) {
946                    outDims = size(1,A.S[1]);   
947                } else if (dim == 1) {
948                    outDims = size(A.S[0], 1);
949                } else if (A.S.IsSameSize(B.S)) {
950                    outDims = A.S;
951                    dim = 0;
952                } else
953                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
954
955               
956                fcomplex[] retArr = ILMemoryPool.Pool.New<fcomplex>(outDims.NumberOfElements);
957
958               
959                fcomplex[] arrA = A.GetArrayForRead();
960
961               
962                fcomplex[] arrB = B.GetArrayForRead();
963                int itLen = A.S[0];
964                #endregion
965
966                #region worker loops definition
967                ILDenseStorage<fcomplex> retStorage = new ILDenseStorage<fcomplex>(retArr, outDims);
968                int workerCount = 1;
969                Action<object> worker = data => {
970                    // expects: iStart, iLen, ap, bp, cp
971                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
972                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
973
974                   
975                    fcomplex* ap;
976
977                   
978                    fcomplex* bp;
979
980                   
981                    fcomplex* cp;
982                    if (dim == 0) {
983                        ap = (fcomplex*)range.Item2;
984                        cp = (fcomplex*)range.Item4;
985                        for (int s = 0; s < range.Item1; s++) {
986                            bp = (fcomplex*)range.Item3;
987                           
988                            fcomplex sum = 0;
989                            int l = itLen;
990                            while (l > 20) {
991                                sum +=  fcomplex.Abs(ap[0] - bp[0])
992                                    +  fcomplex.Abs(ap[1] - bp[1])
993                                    +  fcomplex.Abs(ap[2] - bp[2])
994                                    +  fcomplex.Abs(ap[3] - bp[3])
995                                    +  fcomplex.Abs(ap[4] - bp[4])
996                                    +  fcomplex.Abs(ap[5] - bp[5])
997                                    +  fcomplex.Abs(ap[6] - bp[6])
998                                    +  fcomplex.Abs(ap[7] - bp[7])
999                                    +  fcomplex.Abs(ap[8] - bp[8])
1000                                    +  fcomplex.Abs(ap[9] - bp[9])
1001                                    +  fcomplex.Abs(ap[10] - bp[10])
1002                                    +  fcomplex.Abs(ap[11] - bp[11])
1003                                    +  fcomplex.Abs(ap[12] - bp[12])
1004                                    +  fcomplex.Abs(ap[13] - bp[13])
1005                                    +  fcomplex.Abs(ap[14] - bp[14])
1006                                    +  fcomplex.Abs(ap[15] - bp[15])
1007                                    +  fcomplex.Abs(ap[16] - bp[16])
1008                                    +  fcomplex.Abs(ap[17] - bp[17])
1009                                    +  fcomplex.Abs(ap[18] - bp[18])
1010                                    +  fcomplex.Abs(ap[19] - bp[19])
1011                                    +  fcomplex.Abs(ap[20] - bp[20]);
1012                                ap += 21;
1013                                bp += 21;
1014                                l -= 21;
1015                            }
1016                            while (l-- > 0) {
1017                                sum +=  fcomplex.Abs(*ap - *bp);
1018                                ap++;
1019                                bp++;
1020                            }
1021                            *cp++ = sum;
1022                        }
1023                    } else {
1024                        // dim = 1
1025                        ap = (fcomplex*)range.Item2;
1026                        bp = (fcomplex*)range.Item3;
1027                        cp = (fcomplex*)range.Item4;
1028                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
1029                           
1030                            fcomplex val = *bp++;
1031                            int l = itLen;
1032                            while (l > 10) {
1033                                cp[0] +=  fcomplex.Abs(ap[0] - val);
1034                                cp[1] +=  fcomplex.Abs(ap[1] - val);
1035                                cp[2] +=  fcomplex.Abs(ap[2] - val);
1036                                cp[3] +=  fcomplex.Abs(ap[3] - val);
1037                                cp[4] +=  fcomplex.Abs(ap[4] - val);
1038                                cp[5] +=  fcomplex.Abs(ap[5] - val);
1039                                cp[6] +=  fcomplex.Abs(ap[6] - val);
1040                                cp[7] +=  fcomplex.Abs(ap[7] - val);
1041                                cp[8] +=  fcomplex.Abs(ap[8] - val);
1042                                cp[9] +=  fcomplex.Abs(ap[9] - val);
1043                                cp[10] +=  fcomplex.Abs(ap[10] - val);
1044                                cp[11] +=  fcomplex.Abs(ap[11] - val);
1045                                cp[12] +=  fcomplex.Abs(ap[12] - val);
1046                                cp[13] +=  fcomplex.Abs(ap[13] - val);
1047                                cp[14] +=  fcomplex.Abs(ap[14] - val);
1048                                cp[15] +=  fcomplex.Abs(ap[15] - val);
1049                                cp[16] +=  fcomplex.Abs(ap[16] - val);
1050                                cp[17] +=  fcomplex.Abs(ap[17] - val);
1051                                cp[18] +=  fcomplex.Abs(ap[18] - val);
1052                                cp[19] +=  fcomplex.Abs(ap[19] - val);
1053                                cp[20] +=  fcomplex.Abs(ap[20] - val);
1054                                ap += 21;
1055                                cp += 21;
1056                                l -= 21;
1057                            }
1058                            while (l-- > 0) {
1059                                *cp +=  fcomplex.Abs(*ap - val);
1060                                ap++;
1061                                cp++;
1062                            }
1063                        }
1064                    }
1065                    System.Threading.Interlocked.Decrement(ref workerCount);
1066                };
1067                #endregion
1068
1069                #region work distribution
1070                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
1071                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
1072                    && outDims[1] > 1) {
1073                        if (outDims[1] > workItemCount) {
1074                        workItemLength = outDims[1] / workItemCount;
1075                    } else {
1076                        workItemLength = outDims[1] / 2;
1077                        workItemCount = 2;
1078                    }
1079                } else {
1080                    workItemLength = outDims[1];
1081                    workItemCount = 1;
1082                }
1083
1084                //int[] partResults = new int[workItemCount];
1085                fixed ( fcomplex* arrAP = arrA)
1086                fixed ( fcomplex* arrBP = arrB)
1087                fixed ( fcomplex* retArrP = retArr) {
1088
1089                    for (; i < workItemCount - 1; i++) {
1090                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
1091                               (workItemLength
1092                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
1093                               , (IntPtr)(arrBP + i * dim * workItemLength)
1094                               , (IntPtr)(retArrP + i * workItemLength));
1095                        System.Threading.Interlocked.Increment(ref workerCount);
1096                        ILThreadPool.QueueUserWorkItem(i, worker, range);
1097                    }
1098                    // the last (or may the only) chunk is done right here
1099                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
1100                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
1101                                (outDims[1] - i * workItemLength
1102                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
1103                                , (IntPtr)(arrBP + i * dim * workItemLength)
1104                                , (IntPtr)(retArrP + i * workItemLength)));
1105
1106                    ILThreadPool.Wait4Workers(ref workerCount);
1107                }
1108                #endregion
1109
1110                return new ILRetArray<fcomplex>(retStorage);
1111            }
1112        }
1113        /// <summary>
1114        /// pairwise L1 distance
1115        /// </summary>
1116        /// <param name="A">input points (matrix)</param>
1117        /// <param name="B">input point (vector)</param>
1118        /// <returns>pairwise L1 distances between the data point provided in the input vector <paramref name="B"/> and the data points stored in the matrix <paramref name="A"/>.</returns>
1119        /// <remarks>
1120        /// <para>If <paramref name="B"/> is a colum vector, the distances between <paramref name="B"/> and the columns of <paramref name="A"/> are calculated. The number of rows of <paramref name="A"/>
1121        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of columns of <paramref name="A"/>: <code>A.S[1]</code>.</para>
1122        /// <para>If <paramref name="B"/> is a row vector, the distances between <paramref name="B"/> and the rows of <paramref name="A"/> are calculated. The number of columns of <paramref name="A"/>
1123        /// must match the length of vector <paramref name="B"/> than. Therefore, the length of the returned vector of distances matches the number of rows of <paramref name="A"/>: <code>A.S[0]</code>.</para>
1124        /// <para>This function is cummulative, the single data point may be provided in <paramref name="A"/> and the data point matrix may be provided in <paramref name="B"/> as well.</para>
1125        /// </remarks>
1126        public static unsafe ILRetArray<complex> distL1(ILInArray<complex> A, ILInArray<complex> B) {
1127            using (ILScope.Enter(A, B)) {
1128
1129                #region parameter checking
1130                if (isnull(A) || isnull(B))
1131                    return empty<complex>(ILSize.Empty00);
1132                if (A.IsEmpty) {
1133                    return empty<complex>(B.S);
1134                } else if (B.IsEmpty) {
1135                    return empty<complex>(A.S);
1136                }
1137                // early exit: make the function cummutative
1138                if (A.IsVector && !B.IsVector) {
1139                    return distL1(B, A);
1140                }
1141                int dim = -1;
1142                for (int l = 0; l < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); l++) {
1143                    if (A.S[l] != B.S[l]) {
1144                        if (dim >= 0 || (A.S[l] != 1 && B.S[l] != 1)) {
1145                            throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
1146                        }
1147                        dim = l;
1148                    }
1149                }
1150               
1151                dim = -(dim - 1);  // 0 -> 1, 1 -> 0
1152                #endregion
1153
1154                #region parameter preparation
1155                ILSize outDims;
1156                if (dim == 0) {
1157                    outDims = size(1,A.S[1]);   
1158                } else if (dim == 1) {
1159                    outDims = size(A.S[0], 1);
1160                } else if (A.S.IsSameSize(B.S)) {
1161                    outDims = A.S;
1162                    dim = 0;
1163                } else
1164                    throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
1165
1166               
1167                complex[] retArr = ILMemoryPool.Pool.New<complex>(outDims.NumberOfElements);
1168
1169               
1170                complex[] arrA = A.GetArrayForRead();
1171
1172               
1173                complex[] arrB = B.GetArrayForRead();
1174                int itLen = A.S[0];
1175                #endregion
1176
1177                #region worker loops definition
1178                ILDenseStorage<complex> retStorage = new ILDenseStorage<complex>(retArr, outDims);
1179                int workerCount = 1;
1180                Action<object> worker = data => {
1181                    // expects: iStart, iLen, ap, bp, cp
1182                    Tuple<int, IntPtr, IntPtr, IntPtr> range =
1183                        (Tuple<int, IntPtr, IntPtr, IntPtr>)data;
1184
1185                   
1186                    complex* ap;
1187
1188                   
1189                    complex* bp;
1190
1191                   
1192                    complex* cp;
1193                    if (dim == 0) {
1194                        ap = (complex*)range.Item2;
1195                        cp = (complex*)range.Item4;
1196                        for (int s = 0; s < range.Item1; s++) {
1197                            bp = (complex*)range.Item3;
1198                           
1199                            complex sum = 0;
1200                            int l = itLen;
1201                            while (l > 20) {
1202                                sum +=  complex.Abs(ap[0] - bp[0])
1203                                    +  complex.Abs(ap[1] - bp[1])
1204                                    +  complex.Abs(ap[2] - bp[2])
1205                                    +  complex.Abs(ap[3] - bp[3])
1206                                    +  complex.Abs(ap[4] - bp[4])
1207                                    +  complex.Abs(ap[5] - bp[5])
1208                                    +  complex.Abs(ap[6] - bp[6])
1209                                    +  complex.Abs(ap[7] - bp[7])
1210                                    +  complex.Abs(ap[8] - bp[8])
1211                                    +  complex.Abs(ap[9] - bp[9])
1212                                    +  complex.Abs(ap[10] - bp[10])
1213                                    +  complex.Abs(ap[11] - bp[11])
1214                                    +  complex.Abs(ap[12] - bp[12])
1215                                    +  complex.Abs(ap[13] - bp[13])
1216                                    +  complex.Abs(ap[14] - bp[14])
1217                                    +  complex.Abs(ap[15] - bp[15])
1218                                    +  complex.Abs(ap[16] - bp[16])
1219                                    +  complex.Abs(ap[17] - bp[17])
1220                                    +  complex.Abs(ap[18] - bp[18])
1221                                    +  complex.Abs(ap[19] - bp[19])
1222                                    +  complex.Abs(ap[20] - bp[20]);
1223                                ap += 21;
1224                                bp += 21;
1225                                l -= 21;
1226                            }
1227                            while (l-- > 0) {
1228                                sum +=  complex.Abs(*ap - *bp);
1229                                ap++;
1230                                bp++;
1231                            }
1232                            *cp++ = sum;
1233                        }
1234                    } else {
1235                        // dim = 1
1236                        ap = (complex*)range.Item2;
1237                        bp = (complex*)range.Item3;
1238                        cp = (complex*)range.Item4;
1239                        for (int s = 0; s < range.Item1; s++) { // TODO: further potential speedup: unroll in s direction!
1240                           
1241                            complex val = *bp++;
1242                            int l = itLen;
1243                            while (l > 10) {
1244                                cp[0] +=  complex.Abs(ap[0] - val);
1245                                cp[1] +=  complex.Abs(ap[1] - val);
1246                                cp[2] +=  complex.Abs(ap[2] - val);
1247                                cp[3] +=  complex.Abs(ap[3] - val);
1248                                cp[4] +=  complex.Abs(ap[4] - val);
1249                                cp[5] +=  complex.Abs(ap[5] - val);
1250                                cp[6] +=  complex.Abs(ap[6] - val);
1251                                cp[7] +=  complex.Abs(ap[7] - val);
1252                                cp[8] +=  complex.Abs(ap[8] - val);
1253                                cp[9] +=  complex.Abs(ap[9] - val);
1254                                cp[10] +=  complex.Abs(ap[10] - val);
1255                                cp[11] +=  complex.Abs(ap[11] - val);
1256                                cp[12] +=  complex.Abs(ap[12] - val);
1257                                cp[13] +=  complex.Abs(ap[13] - val);
1258                                cp[14] +=  complex.Abs(ap[14] - val);
1259                                cp[15] +=  complex.Abs(ap[15] - val);
1260                                cp[16] +=  complex.Abs(ap[16] - val);
1261                                cp[17] +=  complex.Abs(ap[17] - val);
1262                                cp[18] +=  complex.Abs(ap[18] - val);
1263                                cp[19] +=  complex.Abs(ap[19] - val);
1264                                cp[20] +=  complex.Abs(ap[20] - val);
1265                                ap += 21;
1266                                cp += 21;
1267                                l -= 21;
1268                            }
1269                            while (l-- > 0) {
1270                                *cp +=  complex.Abs(*ap - val);
1271                                ap++;
1272                                cp++;
1273                            }
1274                        }
1275                    }
1276                    System.Threading.Interlocked.Decrement(ref workerCount);
1277                };
1278                #endregion
1279
1280                #region work distribution
1281                int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
1282                if (Settings.s_maxNumberThreads > 1 && outDims.NumberOfElements >= Settings.s_minParallelElement1Count
1283                    && outDims[1] > 1) {
1284                        if (outDims[1] > workItemCount) {
1285                        workItemLength = outDims[1] / workItemCount;
1286                    } else {
1287                        workItemLength = outDims[1] / 2;
1288                        workItemCount = 2;
1289                    }
1290                } else {
1291                    workItemLength = outDims[1];
1292                    workItemCount = 1;
1293                }
1294
1295                //int[] partResults = new int[workItemCount];
1296                fixed ( complex* arrAP = arrA)
1297                fixed ( complex* arrBP = arrB)
1298                fixed ( complex* retArrP = retArr) {
1299
1300                    for (; i < workItemCount - 1; i++) {
1301                        Tuple<int, IntPtr, IntPtr, IntPtr> range = new Tuple<int, IntPtr, IntPtr, IntPtr>
1302                               (workItemLength
1303                               , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
1304                               , (IntPtr)(arrBP + i * dim * workItemLength)
1305                               , (IntPtr)(retArrP + i * workItemLength));
1306                        System.Threading.Interlocked.Increment(ref workerCount);
1307                        ILThreadPool.QueueUserWorkItem(i, worker, range);
1308                    }
1309                    // the last (or may the only) chunk is done right here
1310                    //System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
1311                    worker(new Tuple<int, IntPtr, IntPtr, IntPtr>
1312                                (outDims[1] - i * workItemLength
1313                                , (IntPtr)(arrAP + i * outDims[0] * workItemLength)
1314                                , (IntPtr)(arrBP + i * dim * workItemLength)
1315                                , (IntPtr)(retArrP + i * workItemLength)));
1316
1317                    ILThreadPool.Wait4Workers(ref workerCount);
1318                }
1319                #endregion
1320
1321                return new ILRetArray<complex>(retStorage);
1322            }
1323        }
1324
1325#endregion HYCALPER AUTO GENERATED CODE
1326
1327    }
1328}
Note: See TracBrowser for help on using the repository browser.