/*
* Copyright (c) 2002-2006 Samit Basu
*
* 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.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "Array.hpp"
#include "Exception.hpp"
#include "Data.hpp"
#include "Malloc.hpp"
#include "IEEEFP.hpp"
#include "Interpreter.hpp"
#include "Sparse.hpp"
#include <math.h>
#include <stdio.h>
#include <set>
#include "FunctionDef.hpp"
#include "NumericArray.hpp"
#include "LAPACK.hpp"
static int objectBalance;
#define MSGBUFLEN 2048
typedef std::set<uint32, std::less<uint32> > intSet;
intSet addresses;
bool isColonOperator(Array& a) {
return ((a.dataClass() == FM_STRING) &&
(a.getLength() == 1) &&
(((const char*) a.getDataPointer())[0] == ':'));
}
ArrayVector singleArrayVector(Array a) {
ArrayVector retval;
retval.push_back(a);
return retval;
}
ArrayVector operator+(Array a, Array b) {
ArrayVector retval;
retval.push_back(a);
retval.push_back(b);
return retval;
}
ArrayVector operator+(ArrayVector a, Array b) {
a.push_back(b);
return a;
}
ArrayVector operator+(Array a, ArrayVector b) {
b.push_front(a);
return b;
}
void outputDoublePrecisionFloat(char *buf, double num) {
char temp_buf[100];
char *tbuf;
int len;
tbuf = temp_buf;
if (num>=0) {
sprintf(tbuf," ");
tbuf++;
}
if (IsNaN(num))
sprintf(tbuf," nan");
else if (fabs(num)>=0.1f && fabs(num)<1.0f || num == 0.0f)
sprintf(tbuf," %0.15f",num);
else if (fabs(num)>=0.01f && fabs(num)<0.1f)
sprintf(tbuf," %0.16f",num);
else if (fabs(num)>=0.001f && fabs(num)<0.01f)
sprintf(tbuf," %0.17f",num);
else if (fabs(num)>=1.0f && fabs(num)<10.0f)
sprintf(tbuf," %1.15f",num);
else if (fabs(num)>=10.0f && fabs(num)<100.0f)
sprintf(tbuf," %2.13f",num);
else if (fabs(num)>=100.0f && fabs(num)<1000.0f)
sprintf(tbuf,"%3.12f",num);
else
sprintf(tbuf," %1.14e",num);
len = strlen(temp_buf);
memcpy(buf,temp_buf,len);
memset(buf+len,' ',24-len);
buf[24] = 0;
}
void outputSinglePrecisionFloat(char *buf, float num) {
char temp_buf[100];
char *tbuf;
int len;
tbuf = temp_buf;
if (num>=0) {
sprintf(tbuf," ");
tbuf++;
}
if (IsNaN(num))
sprintf(tbuf," nan");
else if (fabs(num)>=0.1f && fabs(num)<1.0f || num == 0.0f)
sprintf(tbuf," %0.8f",num);
else if (fabs(num)>=0.01f && fabs(num)<0.1f)
sprintf(tbuf," %0.9f",num);
else if (fabs(num)>=0.001f && fabs(num)<0.01f)
sprintf(tbuf," %0.10f",num);
else if (fabs(num)>=1.0f && fabs(num)<10.0f)
sprintf(tbuf," %1.7f",num);
else if (fabs(num)>=10.0f && fabs(num)<100.0f)
sprintf(tbuf," %2.6f",num);
else if (fabs(num)>=100.0f && fabs(num)<1000.0f)
sprintf(tbuf,"%3.5f",num);
else
sprintf(tbuf," %1.7e",num);
len = strlen(temp_buf);
memcpy(buf,temp_buf,len);
memset(buf+len,' ',17-len);
buf[17] = 0;
}
void* Array::allocateArray(Class type, uint32 length, rvstring names) {
if (length == 0) return NULL;
switch(type) {
case FM_FUNCPTR_ARRAY: {
FuncPtr *dp = new FuncPtr[length];
for (int i=0;i<length;i++)
dp[i] = FuncPtr((FunctionDef*) NULL);
return dp;
}
case FM_CELL_ARRAY: {
Array *dp = new Array[length];
for (int i=0;i<length;i++)
dp[i] = Array(FM_DOUBLE);
return dp;
}
case FM_STRUCT_ARRAY: {
Array *dp = new Array[length*names.size()];
for (int i=0;i<length*names.size();i++)
dp[i] = Array(FM_DOUBLE);
return dp;
}
case FM_LOGICAL:
return Malloc(sizeof(logical)*length);
case FM_UINT8:
return Malloc(sizeof(uint8)*length);
case FM_INT8:
return Malloc(sizeof(int8)*length);
case FM_UINT16:
return Malloc(sizeof(uint16)*length);
case FM_INT16:
return Malloc(sizeof(int16)*length);
case FM_UINT32:
return Malloc(sizeof(uint32)*length);
case FM_INT32:
return Malloc(sizeof(int32)*length);
case FM_UINT64:
return Malloc(sizeof(uint64)*length);
case FM_INT64:
return Malloc(sizeof(int64)*length);
case FM_FLOAT:
return Malloc(sizeof(float)*length);
case FM_DOUBLE:
return Malloc(sizeof(double)*length);
case FM_COMPLEX:
return Malloc(2*sizeof(float)*length);
case FM_DCOMPLEX:
return Malloc(2*sizeof(double)*length);
case FM_STRING:
return Malloc(sizeof(char)*length);
}
return 0;
}
bool* Array::getBinaryMap(uint32 maxD) {
bool* map = (bool*) Malloc(maxD*sizeof(bool));
int N = getLength();
constIndexPtr rp = (constIndexPtr) data();
for (int i=0;i<N;i++) {
indexType n = (rp[i]-1);
if (n < 0 || n >= maxD) {
Free(map);
throw Exception("Array index exceeds bounds");
}
map[n] = true;
}
return map;
}
uint32 Array::getMaxAsIndex() {
indexType maxval;
int k;
constIndexPtr rp = (constIndexPtr) data();
int K = getLength();
maxval = rp[0];
for (k=1;k<K;k++)
if (rp[k] > maxval) maxval = rp[k];
if (maxval <= 0) throw Exception("Illegal zero or negative index");
return maxval;
}
void Array::toOrdinalType(Interpreter *m_eval) {
// Special case : sparse matrices with logical type can be efficiently
// converted to ordinal types
if (sparse() && dataClass() == FM_LOGICAL) {
int outcount;
uint32 *sp = SparseLogicalToOrdinal(dimensions().get(0),dimensions().get(1),data(),outcount);
setData(FM_UINT32,Dimensions(outcount,1),sp);
return;
}
if (sparse())
makeDense();
switch(dataClass()) {
case FM_LOGICAL:
{
// We make a first pass through the array, and count the number of
// non-zero entries.
const logical *rp = (const logical *) data();
int indexCount = 0;
int len = getLength();
int i;
for (i=0;i<len;i++)
if (rp[i] != 0) indexCount++;
// Allocate space to hold the new type.
indexType *lp = (indexType *) Malloc(indexCount*sizeof(indexType));
indexType *qp = lp;
for (i=0;i<len;i++)
if (rp[i] != 0) *qp++ = (indexType) (i+1);
// Reset our data pointer to the new vector.
Dimensions dimensions;
dimensions.set(1,1);
dimensions.set(0,indexCount);
// Change the class to an FM_UINT32.
setData(FM_UINT32,dimensions,lp);
}
break;
case FM_STRING:
{
throw Exception("Cannot convert string data types to indices.");
}
break;
case FM_DCOMPLEX:
{
m_eval->warningMessage("Imaginary part of complex index ignored.\n");
// We convert complex values into real values
const double *rp = (const double *) data();
int len = getLength();
indexType ndx;
bool fractionalWarning = false;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = (indexType) rp[i<<1];
if (!fractionalWarning && (double) ndx != rp[i<<1]) {
m_eval->warningMessage("Fractional part of index ignored.\n");
fractionalWarning = true;
}
if (ndx <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_COMPLEX:
{
m_eval->warningMessage("Imaginary part of complex index ignored.\n");
// We convert complex values into real values
const float *rp = (const float *) data();
int len = getLength();
indexType ndx;
bool fractionalWarning = false;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = (indexType) rp[i<<1];
if (!fractionalWarning && (double) ndx != rp[i<<1]) {
m_eval->warningMessage("Fractional part of index ignored.\n");
fractionalWarning = true;
}
if (ndx <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_DOUBLE:
{
const double *rp = (const double *) data();
int len = getLength();
indexType ndx;
bool fractionalWarning = false;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = (indexType) rp[i];
if (!fractionalWarning && (double) ndx != rp[i]) {
m_eval->warningMessage("Fractional part of index ignored.\n");
fractionalWarning = true;
}
if (ndx <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_FLOAT:
{
const float *rp = (const float *) data();
int len = getLength();
indexType ndx;
bool fractionalWarning = false;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = (indexType) rp[i];
if (!fractionalWarning && (double) ndx != rp[i]) {
m_eval->warningMessage("Fractional part of index ignored.\n");
fractionalWarning = true;
}
if (ndx <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_INT64:
{
const int64 *rp = (const int64 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT64,dimensions(),lp);
}
break;
case FM_UINT64:
{
const uint64 *rp = (const uint64 *) data();
int len = getLength();
// Allocate space to hold the new type
for (int i=0;i<len;i++)
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
}
break;
case FM_INT32:
{
const int32 *rp = (const int32 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_UINT32:
{
const uint32 *rp = (const uint32 *) data();
int len = getLength();
// Allocate space to hold the new type
for (int i=0;i<len;i++)
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
}
break;
case FM_INT16:
{
const int16 *rp = (const int16 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_UINT16:
{
const uint16 *rp = (const uint16 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_INT8:
{
const int8 *rp = (const int8 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_UINT8:
{
const uint8 *rp = (const uint8 *) data();
int len = getLength();
indexType ndx;
// Allocate space to hold the new type
indexType *lp = (indexType *) Malloc(len*sizeof(indexType));
for (int i=0;i<len;i++) {
ndx = rp[i];
if (rp[i] <= 0)
throw Exception("Zero or negative index encountered.");
lp[i] = ndx;
}
setData(FM_UINT32,dimensions(),lp);
}
break;
case FM_CELL_ARRAY:
{
throw Exception("Cannot convert cell arrays to indices.");
}
break;
case FM_STRUCT_ARRAY:
{
throw Exception("Cannot convert structure arrays to indices.");
}
break;
case FM_FUNCPTR_ARRAY:
{
throw Exception("Cannot convert function pointer arrays to indices.");
}
break;
}
}
/**
* Compute the ordinal index for a given field name from a
* structure when the list of field names is given.
*/
int32 getFieldIndexFromList(std::string fName,rvstring fieldNames) {
bool foundName = false;
uint32 i;
i = 0;
while (i<fieldNames.size() && !foundName) {
foundName = (fieldNames.at(i) == fName);
if (!foundName) i++;
}
if (foundName)
return i;
else
return -1;
}
/**
* Compute the ordinal index for a given field name from a
* structure using the current set of field names.
*/
int32 Array::getFieldIndex(std::string fName) {
return getFieldIndexFromList(fName,fieldNames());
}
//!
//@Module RESIZE Resizing an Array
//@@Section Array
//@@Usage
//Arrays in FreeMat will resize themselves automatically
//as required in order to accomodate assignments. The rules
//for resizing are as follows. If an assignment is made to
//an n-dimensional array (where n >= 2) that is outside the
//current dimension bounds of the array, then the array is
//zero padded until the it is large enough for the assignment
//to work. If the array is a scalar, and an assignment is
//made to the non-unity element, such as:
//@[
//a = 1;
//a(3) = 4;
//@]
//then the result will be a row vector (in this case, of size 3).
//Row and column vectors will be resized so as to preserve their
//orientation. And if an n-dimensional array is forced to resize
//using the vector notation, then the result is a row vector.
//@@Tests
//@{ test_resize1.m
//function test_val = test_resize1
//% Check a normal resize in 2D
//a = [];
//a = [1,2;3,4];
//a(3,3) = 1;
//test_val = all(a == [1,2,0;3,4,0;0,0,1]);
//@}
//@{ test_resize2.m
//function test_val = test_resize2
//% Check a scalar -> vector resize
//a = 1;
//a(3) = 1;
//test_val = all(a == [1,0,1]);
//@}
//@{ test_resize3.m
//function test_val = test_resize3
//% Check a row-vector -> vector resize
//a = [1 2];
//a(3) = 1;
//test_val = all(a == [1,2,1]);
//@}
//@{ test_resize4.m
//function test_val = test_resize4
//% Check a column-vector -> vector resize
//a = [1;2];
//a(3) = 1;
//test_val = all(a == [1;2;1]);
//@}
//@{ test_resize5.m
//function test_val = test_resize5
//a = 0;
//a(2,1,2) = 0;
//a(5) = 1;
//test_val = all(a == [0,0,0,0,1]);
//@}
//@{ test_delete1.m
//% Check the delete functionality in a vector setting
//function test_val = test_delete1
//a = [1,2,3,4,5];
//a([2,3,4]) = [];
//test_val = test(a == [1,5]);
//@}
//@{ test_delete2.m
//% Check the delete functionality in an ndim setting
//function test_val = test_delete2
//a = [1,2,3;4,5,6];
//a(:,2) = [];
//test_val = test(a == [1,3;4,6]);
//@}
//@{ test_delete3.m
//% Check the delete functionality in an ndim setting
//function test_val = test_delete3
//a = [1,2;3,4];
//a(2,2,2) = 7;
//a(:,:,2) = [];
//test_val = test(a == [1,2;3,4]);
//@}
//@{ test_delete4.m
//% Check the delete-all functionality in an ndim setting
//function test_val = test_delete4
//a = [1,2,3;4,5,6];
//a(:,:) = [];
//test_val = test(isempty(a));
//@}
//@{ test_delete5.m
//% Check the delete-all functionality in a vector setting
//function test_val = test_delete5
//a = [1,2,3;4,5,6];
//a(:) = [];
//test_val = test(isempty(a));
//@}
//@{ test_delete6.m
//% Check the multi-dim delete function with an invalid argument
//function test_val = test_delete6
//a = [1,2,3;4,5,6;7,8,9];
//a(3,3,2) = 10;
//test_val = 0;
//try
// a(1,:,2) = [];
//catch
// test_val = 1;
//end
//@}
//!
void Array::resize(Dimensions& a) {
Dimensions newSize;
// Make a copy of the current dimension vector, and
// compute the new dimension size.
newSize = dimensions();
newSize.expandToCover(a);
// Check to see if the dimensions are unchanged.
if (newSize.equals(dimensions())) return;
// Check to see if the total number of elements is unchanged.
if (newSize.getElementCount() == getLength()) {
reshape(newSize);
return;
}
if (!dp) {
void *dst_data = allocateArray(dataClass(),newSize.getElementCount(),
fieldNames());
setData(dataClass(),newSize,dst_data,sparse(),fieldNames(),className());
return;
}
if (sparse()) {
setData(dataClass(),newSize,
CopyResizeSparseMatrix(dataClass(),
dimensions().get(0),
dimensions().get(1),
data(),
newSize.get(0),
newSize.get(1)),true);
return;
}
// Allocate space for our new size.
void *dst_data = allocateArray(dataClass(),newSize.getElementCount(),
fieldNames());
if (!isEmpty()) {
// Initialize a pointer to zero.
Dimensions curPos(dimensions().getLength());
// Because we copy & convert data a column at a time, we retrieve
// the number of rows in each column.
int rowCount = dimensions().get(0);
// Track our offset into the original data.
int srcIndex = 0;
// Loop until we have exhausted the original data.
int dstIndex;
while (curPos.inside(dimensions())) {
// Get the destination index for the current source position.
dstIndex = newSize.mapPoint(curPos);
// Copy the data from our original data structure to the
// new data structure, starting from the source index
// srcIndex, and moving to dstIndex.
copyElements(srcIndex,dst_data,dstIndex,rowCount);
// Update the column number (as we have just copied an
// entire column).
curPos.incrementModulo(dimensions(),1);
// Advance the source data pointer so that it points to the
// start of the next column.
srcIndex += rowCount;
}
}
setData(dataClass(),newSize,dst_data,
sparse(),fieldNames(),className());
}
void Array::permute(const int32 *permutation) {
// Check for an identity permutation
int Adims = dimensions().getLength();
bool id_perm = true;
for (int i=0;i<Adims;i++)
if ((permutation[i]-1) != i) id_perm = false;
// It is an identity permutation - do nothing
if (id_perm) return;
// In 2D, the only non-identity permutation is a transpose
if ((dimensions().getLength() == 2) && (is2D())) {
transpose();
return;
}
if (isEmpty()) return;
// Allocate a space to store the permuted data
void *dst_data = allocateArray(dataClass(),getLength(),fieldNames());
// Initialize a pointer to zero.
Dimensions curPos(dimensions().getLength());
// Create a permuted dimensions vector
Dimensions newdims(dimensions().permute(permutation));
newdims.updateCacheVariables();
int dstIndex = 0;
// Loop over all points
for (int srcIndex=0;srcIndex < getLength();srcIndex++) {
// Permute the current point into the new array
// and then map it to a linear index
dstIndex = newdims.mapPoint(curPos.permute(permutation));
// Copy this one element
copyElements(srcIndex,dst_data,dstIndex,1);
// Move to the next point in the source array
curPos.incrementModulo(dimensions(),0);
}
// Replace our data with the new array.
setData(dataClass(),newdims,dst_data,sparse(),fieldNames(),className());
}
void Array::vectorResize(int max_index) {
if (max_index > getLength()) {
Dimensions newDim;
if (isEmpty() || dimensions().isScalar()) {
newDim.reset();
newDim = Dimensions(1,max_index);
} else if (dimensions().isVector()) {
newDim = dimensions();
if (dimensions().get(0) != 1)
newDim.set(0,max_index);
else
newDim.set(1,max_index);
} else {
// First reshape it
Dimensions tDim(1,getLength());
reshape(tDim);
newDim = Dimensions(1,max_index);
}
resize(newDim);
}
}
/**
* Reshape an array. This is only legal if the number of
* elements remains the same after reshaping.
*/
void Array::reshape(Dimensions& a) {
if (a.getElementCount() != getLength())
throw Exception("Reshape operation cannot change the number of elements in array.");
if (sparse()) {
a.simplify();
if (a.getLength() > 2)
throw Exception("Cannot reshape sparse matrix to an N-dimensional array - FreeMat does not support N-dimensional sparse arrays");
setData(dataClass(),a,
ReshapeSparseMatrix(dataClass(),
dimensions().get(0),
dimensions().get(1),
data(),
a.get(0),
a.get(1)),true);
} else {
dp->setDimensions(a);
}
}
/**
* Hermitian transpose our array. By default, this is just a transpose
* operation. The distinction is made in complex classed by overriding
* this method.
*/
void Array::hermitian() {
if (!is2D())
throw Exception("Cannot apply Hermitian transpose operation to multi-dimensional array.");
if (isEmpty()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
dp->setDimensions(Dimensions(cols,rows));
return;
}
if (!isComplex())
transpose();
else {
if (sparse()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
void *qp = SparseArrayHermitian(dataClass(), rows, cols, data());
setData(dataClass(),Dimensions(cols,rows),qp,true);
return;
}
if (dataClass() == FM_COMPLEX) {
// Allocate space for our transposed array
void *dstPtr = allocateArray(dataClass(),getLength());
float *qp, *sp;
int i, j;
int rowCount;
int colCount;
rowCount = dimensions().get(0);
colCount = dimensions().get(1);
int ptr;
qp = (float*) dstPtr;
sp = (float*) data();
ptr = 0;
for (i=0;i<rowCount;i++)
for (j=0;j<colCount;j++) {
qp[2*ptr] = sp[2*(i + j*rowCount)];
qp[2*ptr+1] = -sp[2*(i + j*rowCount) + 1];
ptr++;
}
setData(FM_COMPLEX,Dimensions(colCount,rowCount),dstPtr);
} else {
// Allocate space for our transposed array
void *dstPtr = allocateArray(dataClass(),getLength());
double *qp, *sp;
int i, j;
int rowCount;
int colCount;
rowCount = dimensions().get(0);
colCount = dimensions().get(1);
int ptr;
qp = (double*) dstPtr;
sp = (double*) data();
ptr = 0;
for (i=0;i<rowCount;i++)
for (j=0;j<colCount;j++) {
qp[2*ptr] = sp[2*(i + j*rowCount)];
qp[2*ptr+1] = -sp[2*(i + j*rowCount) + 1];
ptr++;
}
setData(FM_DCOMPLEX,Dimensions(colCount,rowCount),dstPtr);
}
}
}
/**
* Transpose our array.
*/
void Array::transpose() {
if (!is2D())
throw Exception("Cannot apply transpose operation to multi-dimensional array.");
if (isEmpty()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
dp->setDimensions(Dimensions(cols,rows));
return;
}
if (sparse()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
void *qp = SparseArrayTranspose(dataClass(), rows, cols, data());
setData(dataClass(),Dimensions(cols,rows),qp,true);
return;
}
// Allocate space for our transposed array
void *dstPtr = allocateArray(dataClass(),getLength(),fieldNames());
int i, j;
int rowCount;
int colCount;
rowCount = dimensions().get(0);
colCount = dimensions().get(1);
int ptr;
ptr = 0;
for (i=0;i<rowCount;i++)
for (j=0;j<colCount;j++) {
copyElements(i+j*rowCount,dstPtr,ptr,1);
ptr++;
}
setData(dataClass(),Dimensions(colCount,rowCount),dstPtr,
sparse(),fieldNames(),className());
}
/**
* Calculate the size of each element in this array.
*/
int Array::getElementSize() const {
switch(dataClass()) {
case FM_FUNCPTR_ARRAY:
return sizeof(FuncPtr);
case FM_CELL_ARRAY:
return sizeof(Array);
case FM_STRUCT_ARRAY:
return (sizeof(Array)*fieldNames().size());
case FM_LOGICAL:
return sizeof(logical);
case FM_UINT8:
return sizeof(uint8);
case FM_INT8:
return sizeof(int8);
case FM_UINT16:
return sizeof(uint16);
case FM_INT16:
return sizeof(int16);
case FM_UINT32:
return sizeof(uint32);
case FM_INT32:
return sizeof(int32);
case FM_UINT64:
return sizeof(uint64);
case FM_INT64:
return sizeof(int64);
case FM_FLOAT:
return sizeof(float);
case FM_DOUBLE:
return sizeof(double);
case FM_COMPLEX:
return sizeof(float)*2;
case FM_DCOMPLEX:
return sizeof(double)*2;
case FM_STRING:
return sizeof(char);
}
return 0;
}
#define caseReal(caseLabel,dpType) \
case caseLabel:\
{\
const dpType* qp = (const dpType*) data();\
bool allPositive = true;\
int N = dimensions().get(0); \
int i, j;\
for (i=0;i<N;i++)\
for (j=i+1;j<N;j++)\
if (qp[i+j*N] != qp[j+i*N]) return false;\
return true;\
break;\
}
#define caseComplex(caseLabel,dpType) \
case caseLabel:\
{\
const dpType* qp = (const dpType*) data();\
bool allPositive = true;\
int N = dimensions().get(0); \
int i, j;\
for (i=0;i<N;i++)\
for (j=i+1;j<N;j++) {\
if (qp[2*(i+j*N)] != qp[2*(j+i*N)]) return false;\
if (qp[2*(i+j*N)+1] != -qp[2*(j+i*N)+1]) return false;\
}\
return true;\
break;\
}
const bool Array::isSymmetric() const {
if (!is2D()) return false;
if (isReferenceType()) return false;
if (sparse())
throw Exception("Cannot determine symmetry of sparse arrays");
switch(dataClass()) {
caseReal(FM_INT8,int8);
caseReal(FM_INT16,int16);
caseReal(FM_INT32,int32);
caseReal(FM_INT64,int32);
caseReal(FM_UINT8,uint8);
caseReal(FM_UINT16,uint16);
caseReal(FM_UINT32,uint32);
caseReal(FM_UINT64,uint32);
caseReal(FM_FLOAT,float);
caseReal(FM_DOUBLE,double);
caseComplex(FM_COMPLEX,float);
caseComplex(FM_DCOMPLEX,double);
}
}
#undef caseReal
#undef caseComplex
/**
* Returns true if we are positive.
*/
#define caseMacro(caseLabel,dpType) \
case caseLabel:\
{\
const dpType* qp = (const dpType*) data();\
bool allPositive = true;\
int len = getLength();\
int i = 0;\
while (allPositive && (i<len)) {\
allPositive = allPositive && (qp[i] >= 0);\
i++;\
}\
return allPositive;\
}
const bool Array::isPositive() const {
if (dataClass() == FM_UINT8 ||
dataClass() == FM_UINT16 ||
dataClass() == FM_UINT32 ||
dataClass() == FM_UINT64)
return true;
if (dataClass() == FM_COMPLEX ||
dataClass() == FM_DCOMPLEX)
return false;
if (sparse())
return SparseIsPositive(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
getSparseDataPointer());
switch (dataClass()) {
caseMacro(FM_FLOAT,float);
caseMacro(FM_DOUBLE,double);
caseMacro(FM_INT8,int8);
caseMacro(FM_INT16,int16);
caseMacro(FM_INT32,int32);
caseMacro(FM_INT64,int64);
}
return false;
}
#undef caseMacro
const bool Array::isRealAllZeros() const {
bool allZeros;
int len = getLength();
int i;
#define caseMacro(caseLabel,dpType,testcode) \
case caseLabel:\
{\
const dpType* qp = (const dpType*) data();\
while (allZeros && (i<len)) {\
allZeros = allZeros && (testcode);\
i++;\
}\
return allZeros;\
}
allZeros = true;
i = 0;
if (sparse())
throw Exception("isPositive not supported for sparse arrays.");
switch (dataClass()) {
caseMacro(FM_LOGICAL,logical,qp[i]==0);
caseMacro(FM_UINT8,uint8,qp[i]==0);
caseMacro(FM_INT8,int8,qp[i]==0);
caseMacro(FM_UINT16,uint16,qp[i]==0);
caseMacro(FM_INT16,int16,qp[i]==0);
caseMacro(FM_UINT32,uint32,qp[i]==0);
caseMacro(FM_INT32,int32,qp[i]==0);
caseMacro(FM_UINT64,uint64,qp[i]==0);
caseMacro(FM_INT64,int64,qp[i]==0);
caseMacro(FM_FLOAT,float,qp[i]==0.0f);
caseMacro(FM_DOUBLE,double,qp[i]==0.0);
caseMacro(FM_COMPLEX,float,qp[i<<1]==0.0f);
caseMacro(FM_DCOMPLEX,double,qp[i<<1]==0.0);
default:
throw Exception("Unable to convert variable type to test for if/while statement");
}
#undef caseMacro
}
#define caseMacroReal(caseLabel,type) \
case caseLabel:\
retval = (*((const type*) x_dp) == *((const type*) y_dp)); \
break;
#define caseMacroComplex(caseLabel,type) \
case caseLabel:\
retval = (((const type*) x_dp)[0] == ((const type*) y_dp)[0]) && \
(((const type*) x_dp)[1] == ((const type*) y_dp)[1]); \
break;
const bool Array::testCaseMatchScalar(Array x) const {
if (sparse())
throw Exception("isPositive not supported for sparse arrays.");
// Now we have to compare ourselves to the argument. Check for the
// case that we are a string type
if (isString()) {
// If x is not a string, we cannot match
if (!x.isString())
return false;
// if x is a string do a string, string compare.
string s1 = getContentsAsString();
string s2 = x.getContentsAsString();
bool retval = (s1 == s2);
return retval;
}
if (!x.isScalar())
return false;
// OK - we are not a string, so we have a numerical comparison. To do this,
// we have to make both objects the same type.
Array y = *this;
if (x.dataClass() > y.dataClass())
y.promoteType(x.dataClass());
else
x.promoteType(y.dataClass());
// Finally, we can do a compare....
const void *x_dp = x.data();
const void *y_dp = y.data();
bool retval;
switch(x.dataClass()) {
caseMacroReal(FM_LOGICAL,logical);
caseMacroReal(FM_UINT8,uint8);
caseMacroReal(FM_INT8,int8);
caseMacroReal(FM_UINT16,uint16);
caseMacroReal(FM_INT16,int16);
caseMacroReal(FM_UINT32,uint32);
caseMacroReal(FM_INT32,int32);
caseMacroReal(FM_UINT64,uint64);
caseMacroReal(FM_INT64,int64);
caseMacroReal(FM_FLOAT,float);
caseMacroReal(FM_DOUBLE,double);
caseMacroComplex(FM_COMPLEX,float);
caseMacroComplex(FM_DCOMPLEX,double);
}
return retval;
}
#undef caseMacroReal
#undef caseMacroComplex
const bool Array::testForCaseMatch(Array x) const {
if (sparse())
throw Exception("isPositive not supported for sparse arrays.");
// We had better be a scalar
if (!(isScalar() || isString()))
throw Exception("Switch argument must be a scalar or a string");
// And we had better not be a reference type
if (isReferenceType())
throw Exception("Switch argument cannot be a reference type (struct or cell array)");
// If x is a scalar, we just need to call the scalar version
if (((x.dataClass() != FM_CELL_ARRAY) && x.isScalar()) || x.isString())
return testCaseMatchScalar(x);
if (x.dataClass() != FM_CELL_ARRAY)
throw Exception("Case arguments must either be a scalar or a cell array");
const Array* qp = (const Array*) x.data();
int len;
len = x.getLength();
bool foundMatch = false;
int i = 0;
while (i<len && !foundMatch) {
foundMatch = testCaseMatchScalar(qp[i]);
i++;
}
return foundMatch;
}
void Array::copyElements(int srcIndex, void* dstPtr, int dstIndex,
int count) const {
int elSize(getElementSize());
if (sparse())
throw Exception("copyElements not supported for sparse arrays.");
switch(dataClass()) {
case FM_CELL_ARRAY:
{
const Array* sp = (const Array*) data();
Array* qp = (Array*) dstPtr;
for (int i=0;i<count;i++)
qp[dstIndex+i] = sp[srcIndex+i];
}
break;
case FM_STRUCT_ARRAY:
{
const Array* sp = (const Array*) data();
Array* qp = (Array*) dstPtr;
int fieldCount(fieldNames().size());
for (int i=0;i<count;i++)
for (int j=0;j<fieldCount;j++)
qp[(dstIndex+i)*fieldCount+j] = sp[(srcIndex+i)*fieldCount+j];
}
break;
default:
{
const char* sp = (const char*) data();
char* qp = (char *) dstPtr;
memcpy(qp + dstIndex*elSize, sp + srcIndex*elSize, count*elSize);
}
break;
}
}
/**
* Promote our data to a new type.
*
* Copy data from our data array to the specified
* array, converting the data as we go. We can only
* convert data to or from base types. So if the source
* or destination types are reference types, we cannot
* perform the conversion.
*
* For the remaining types, we have a matrix of
* possibilities. Here we list the conversion rules.
*
* Source type
* - string
* - logical dest = (source == 0) ? 0 : 1
* - real dest = (double) source
* - complex dest = (double) source
* - logical
* - string dest = (char) source
* - real dest = (double) source
* - complex dest = (double) source
* - real
* - string dest = (char) source
* - logical dest = (source == 0) ? 0 : 1
* - complex dest = (double) source
* - complex
* - string dest = (char) real(source)
* - logical dest = (real(source) == 0 && imag(source) == 0) ? 0:1
* - real dest = real(source)
*/
void Array::promoteType(Class dstClass, rvstring fNames) {
int elCount;
int elSize;
void *dstPtr;
if (!dp) return;
if (isEmpty()) {
setData(dstClass,dimensions(),NULL,false,fNames);
return;
}
// Handle the reference types.
// Cell arrays can be promoted with no effort to cell arrays.
if (dataClass() == FM_FUNCPTR_ARRAY)
if (dstClass == FM_FUNCPTR_ARRAY)
return;
else
throw Exception("Cannot convert function pointer arrays to any other type.");
if (dataClass() == FM_CELL_ARRAY)
if (dstClass == FM_CELL_ARRAY)
return;
else
throw Exception("Cannot convert cell-arrays to any other type.");
// Structure arrays can be promoted to structure arrays with different
// field structures, but have to be rearranged.
if (dataClass() == FM_STRUCT_ARRAY)
if (dstClass == FM_STRUCT_ARRAY) {
// TODO: Generalize this code to allow for one more field in destination
// than in source...
if (fieldNames().size() > fNames.size())
throw Exception("Cannot combine structures with different fields if the combination requires fields to be deleted from one of the structures.");
// We are promoting a struct array to a struct array.
// To do so, we have to make sure that the field names work out.
// The only thing we must check for is that every field name
// in fieldnames is present in fnames.
int extraCount = 0;
int matchCount = 0;
int ndx;
int i;
for (i=0;i<fNames.size();i++) {
ndx = getFieldIndex(fNames.at(i));
if (ndx == -1)
extraCount++;
else
matchCount++;
}
// Now, matchCount should be equal to the size of fieldNames
if (matchCount != fieldNames().size())
throw Exception("Cannot combine structures with different fields if the combination requires fields to be deleted from one of the structures.");
void *dstPtr = allocateArray(dataClass(),getLength(),fNames);
const Array *src_rp = (const Array*) data();
Array * dst_rp = (Array*) dstPtr;
int elCount(getLength());
int fieldCount(fieldNames().size());
int newFieldCount(fNames.size());;
// Now we have to copy our existing fields into the new order...
for (i=0;i<fieldCount;i++) {
int newNdx = getFieldIndexFromList(fieldNames().at(i),fNames);
for (int j=0;j<elCount;j++)
dst_rp[j*newFieldCount + newNdx] = src_rp[j*fieldCount + i];
}
setData(dataClass(),dimensions(),dstPtr,false,fNames);
return;
}
else
throw Exception("Cannot convert struct-arrays to any other type.");
// Catch attempts to convert data types to reference types.
if ((dstClass == FM_CELL_ARRAY) || (dstClass == FM_STRUCT_ARRAY)
|| (dstClass == FM_FUNCPTR_ARRAY))
throw Exception("Cannot convert base types to reference types.");
// Do nothing for promoting to same class (no-op).
if (dstClass == dataClass()) return;
if (sparse()) {
setData(dstClass,dimensions(),
TypeConvertSparse(dataClass(),
dimensions().get(0),
dimensions().get(1),
data(),
dstClass),
true);
return;
}
elCount = getLength();
// We have to promote...
dstPtr = allocateArray(dstClass,elCount);
int count = elCount;
switch (dataClass()) {
#define caseMacro(caseLabel,dpType,convCode) \
case caseLabel: \
{ dpType* qp = (dpType*) dstPtr; \
for (int i=0;i<count;i++) convCode; \
} \
break;
case FM_STRING:
{
const char* sp = (const char *) data();
switch (dstClass) {
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_LOGICAL:
{
const logical* sp = (const logical *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_UINT8:
{
const uint8* sp = (const uint8 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_INT8:
{
const int8* sp = (const int8 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_UINT16:
{
const uint16* sp = (const uint16 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_INT16:
{
const int16* sp = (const int16 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_UINT32:
{
const uint32* sp = (const uint32 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_INT32:
{
const int32* sp = (const int32 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_UINT64:
{
const uint64* sp = (const uint64 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_INT64:
{
const int64* sp = (const int64 *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_FLOAT:
{
const float* sp = (const float *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_DOUBLE:
{
const double* sp = (const double *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i]);
caseMacro(FM_LOGICAL,logical,qp[i] = (sp[i]==0) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i]);
caseMacro(FM_COMPLEX,float,qp[i<<1] = (float) sp[i]);
caseMacro(FM_DCOMPLEX,double,qp[i<<1] = (double) sp[i]);
}
}
break;
case FM_COMPLEX:
{
const float* sp = (const float *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i<<1]);
caseMacro(FM_LOGICAL,logical,qp[i] = ((sp[i<<1]==0.0) && (sp[(i<<1) + 1] == 0.0)) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i<<1]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i<<1]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i<<1]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i<<1]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i<<1]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i<<1]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i<<1]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i<<1]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i<<1]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i<<1]);
caseMacro(FM_DCOMPLEX,double,{qp[i<<1]=(double)sp[i<<1];qp[(i<<1)+1]=(double)sp[(i<<1)+1];});
}
}
break;
case FM_DCOMPLEX:
{
const double* sp = (const double *) data();
switch (dstClass) {
caseMacro(FM_STRING,char,qp[i] = (char) sp[i<<1]);
caseMacro(FM_LOGICAL,logical,qp[i] = ((sp[i<<1]==0.0) && (sp[(i<<1) + 1] == 0.0)) ? 0 : 1);
caseMacro(FM_UINT8,uint8,qp[i] = (uint8) sp[i<<1]);
caseMacro(FM_INT8,int8,qp[i] = (int8) sp[i<<1]);
caseMacro(FM_UINT16,uint16,qp[i] = (uint16) sp[i<<1]);
caseMacro(FM_INT16,int16,qp[i] = (int16) sp[i<<1]);
caseMacro(FM_UINT32,uint32,qp[i] = (uint32) sp[i<<1]);
caseMacro(FM_INT32,int32,qp[i] = (int32) sp[i<<1]);
caseMacro(FM_UINT64,uint64,qp[i] = (uint64) sp[i<<1]);
caseMacro(FM_INT64,int64,qp[i] = (int64) sp[i<<1]);
caseMacro(FM_FLOAT,float,qp[i] = (float) sp[i<<1]);
caseMacro(FM_DOUBLE,double,qp[i] = (double) sp[i<<1]);
caseMacro(FM_COMPLEX,float,{qp[i<<1]=(float)sp[i<<1];qp[(i<<1)+1]=(float)sp[(i<<1)+1];});
}
}
break;
}
setData(dstClass,dimensions(),dstPtr);
}
#undef caseMacro
void Array::promoteType(Class dstClass) {
promoteType(dstClass,rvstring());
}
/********************************************************************************
* Constructors *
********************************************************************************/
Array Array::diagonalConstructor(Array src, int diagonalOrder) {
Array retval;
if (!src.isVector())
throw Exception("Argument to diagonal constructor must by a vector!\n");
int length;
length = src.getLength();
int M;
// Calculate the size of the output matrix (square of size outLen + abs(diagonalOrder)).
M = length + abs(diagonalOrder);
Dimensions dims(M,M);
// Allocate space for the output
void *rp = allocateArray(src.dataClass(),dims.getElementCount(),src.fieldNames());
int i;
int dstIndex;
if (diagonalOrder < 0) {
for (i=0;i<length;i++) {
dstIndex = -diagonalOrder + i * (M+1);
src.copyElements(i,rp,dstIndex,1);
}
} else {
for (i=0;i<length;i++) {
dstIndex = diagonalOrder*M + i * (M+1);
src.copyElements(i,rp,dstIndex,1);
}
}
return Array(src.dataClass(),dims,rp,false,src.fieldNames(),
src.className());
}
Array Array::logicalConstructor(bool aval) {
logical *data = (logical *) allocateArray(FM_LOGICAL,1);
*data = (logical) aval;
return Array(FM_LOGICAL,Dimensions(1,1),data);
}
Array Array::uint8Constructor(uint8 aval) {
uint8 *data = (uint8 *) allocateArray(FM_UINT8,1);
*data = aval;
return Array(FM_UINT8,Dimensions(1,1),data);
}
Array Array::int8Constructor(int8 aval) {
int8 *data = (int8 *) allocateArray(FM_INT8,1);
*data = aval;
return Array(FM_INT8,Dimensions(1,1),data);
}
Array Array::uint16Constructor(uint16 aval) {
uint16 *data = (uint16 *) allocateArray(FM_UINT16,1);
*data = aval;
return Array(FM_UINT16,Dimensions(1,1),data);
}
Array Array::int16Constructor(int16 aval) {
int16 *data = (int16 *) allocateArray(FM_INT16,1);
*data = aval;
return Array(FM_INT16,Dimensions(1,1),data);
}
Array Array::uint32Constructor(uint32 aval) {
uint32 *data = (uint32 *) allocateArray(FM_UINT32,1);
*data = aval;
return Array(FM_UINT32,Dimensions(1,1),data);
}
Array Array::int32Constructor(int32 aval) {
int32 *data = (int32 *) allocateArray(FM_INT32,1);
*data = aval;
return Array(FM_INT32,Dimensions(1,1),data);
}
Array Array::uint64Constructor(uint64 aval) {
uint64 *data = (uint64 *) allocateArray(FM_UINT64,1);
*data = aval;
return Array(FM_UINT64,Dimensions(1,1),data);
}
Array Array::int64Constructor(int64 aval) {
int64 *data = (int64 *) allocateArray(FM_INT64,1);
*data = aval;
return Array(FM_INT64,Dimensions(1,1),data);
}
Array Array::floatConstructor(float aval) {
float *data = (float *) allocateArray(FM_FLOAT,1);
*data = aval;
return Array(FM_FLOAT,Dimensions(1,1),data);
}
Array Array::doubleConstructor(double aval) {
double *data = (double *) allocateArray(FM_DOUBLE,1);
*data = aval;
return Array(FM_DOUBLE,Dimensions(1,1),data);
}
Array Array::floatVectorConstructor(int len) {
float *data = (float*) allocateArray(FM_FLOAT,len);
return Array(FM_FLOAT,Dimensions(1,len),data);
}
Array Array::uint32VectorConstructor(int len) {
uint32 *data = (uint32*) allocateArray(FM_UINT32,len);
return Array(FM_UINT32,Dimensions(1,len),data);
}
Array Array::doubleVectorConstructor(int len) {
double *data = (double*) allocateArray(FM_DOUBLE,len);
return Array(FM_DOUBLE,Dimensions(1,len),data);
}
Array Array::complexConstructor(float aval, float bval) {
float *data = (float *) allocateArray(FM_COMPLEX,1);
data[0] = aval;
data[1] = bval;
return Array(FM_COMPLEX,Dimensions(1,1),data);
}
Array Array::dcomplexConstructor(double aval, double bval) {
double *data = (double *) allocateArray(FM_DCOMPLEX,1);
data[0] = aval;
data[1] = bval;
return Array(FM_DCOMPLEX,Dimensions(1,1),data);
}
Array Array::emptyConstructor() {
return Array(FM_DOUBLE,Dimensions(0,0),NULL);
}
Array Array::funcPtrConstructor(FuncPtr fptr) {
FuncPtr *data = (FuncPtr*) allocateArray(FM_FUNCPTR_ARRAY,1);
data[0] = fptr;
return Array(FM_FUNCPTR_ARRAY, Dimensions(1,1), data);
}
//!
//@Module STRING String Arrays
//@@Section VARIABLES
//@@Usage
//FreeMat supports a @|string| array type that operates very
//much as you would expect. Strings are stored internally as
//8-bit values, and are otherwise similar to numerical arrays in
//all respects. In some respects, this makes strings arrays
//less useful than one might imagine. For example, numerical
//arrays in 2-D are rectangular. Thus, each row in the array
//must have the same number of columns. This requirement is
//natural for numerical arrays and matrices, but consider a
//string array. If one wants to store multiple strings in
//a single data structure, they must all be the
//same length (unlikely). The alternative is to use a cell
//array of strings, in which case, each string can be of
//arbitrary length. Most of the functions that support
//strings in a set-theoretic way, like @|unique| and @|sort|
//operate on cell-arrays of strings instead of string arrays.
//Just to make the example concrete, here is the old way of
//storing several strings in an array:
//@<1
//% This is an error
//A = ['hello';'bye']
//% This is OK, but awkward
//A = ['hello';'bye ']
//% This is the right way to do it
//A = {'hello','bye'}
//@>
//
//One important (tricky) point in FreeMat is the treatment
//of escape sequences. Recall that in @|C| programming, an
//escape sequence is a special character that causes the output
//to do something unusual. FreeMat supports the following
//escape sequences:
//\begin{itemize}
//\item @|\t| - causes a tab to be output
//\item @|\r| - causes a carriage return (return to the beginning
// of the line of output, and overwrite the text)
//\item @|\n| - causes a linefeed (advance to next line)
//\end{itemize}
//FreeMat follows the @|Unix/Linux| convention, that a @|\n|
//causes both a carriage return and a linefeed.
//To put a single quote into a string use the MATLAB convention
//of two single quotes, not the @|\'| sequence. Here is an
//example of a string containing some escape sequences:
//@<
//a = 'I can''t use contractions\n\tOr can I?\n'
//@>
//Now, note that the string itself still contains the @|\n|
//characters. With the exception of the @|\'|, the escape
//sequences do not affect the output unless the strings are
//put through @|printf| or @|fprintf|. For example, if we
//@|printf| the variable @|a|, we see the @|\n| and @|\t|
//take effect:
//@<
//printf(a);
//@>
//The final complicating factor is on @|MSWin| systems. There,
//filenames regularly contain @|\| characters. Thus, if you
//try to print a string containing the filename
//@|C:\redball\timewarp\newton.txt|, the output will be mangled
//because FreeMat thinks the @|\r|, @|\t| and @|\n| are escape
//sequences. You have two options. You can use @|disp| to show
//the filename (@|disp| does not do escape translation to be
//compatible with MATLAB). The second option is to escape the
//backslashes in the string, so that the string you send to
//@|printf| contains @|C:\\redball\\timewarp\\newton.txt|.
//@<
//% disp displays it ok
//a = 'C:\redball\timewarp\newton.txt'
//% printf makes a mess
//printf(a)
//% If we double up the slashes it works fine
//a = 'C:\\redball\\timewarp\\newton.txt'
//printf(a)
//@>
//!
Array Array::stringConstructor(std::string astr) {
int length;
char *cp;
length = astr.length();
cp = (char *) allocateArray(FM_STRING,length);
memcpy(cp,astr.c_str(),length);
return Array(FM_STRING,Dimensions(1,length),cp);
}
Array Array::int32RangeConstructor(int32 minval, int32 stepsize,
int32 maxval, bool vert) {
Dimensions dim;
int32 *rp = NULL;
if (stepsize == 0) throw Exception("step size must be nonzero in colon expression");
int scount = 0;
int accum = minval;
if (stepsize > 0) {
while (accum <= maxval) {
accum += stepsize;
scount++;
}
} else {
while (accum >= maxval) {
accum += stepsize;
scount++;
}
}
if (vert) {
dim = Dimensions(scount,1);
} else {
dim = Dimensions(1,scount);
}
rp = (int32 *) allocateArray(FM_INT32,scount);
for (int i=0;i<scount;i++)
rp[i] = minval + i*stepsize;
return Array(FM_INT32,dim,rp);
}
Array Array::int64RangeConstructor(int64 minval, int64 stepsize,
int64 maxval, bool vert) {
Dimensions dim;
int64 *rp = NULL;
if (stepsize == 0) throw Exception("step size must be nonzero in colon expression");
int scount = 0;
int accum = minval;
if (stepsize > 0) {
while (accum <= maxval) {
accum += stepsize;
scount++;
}
} else {
while (accum >= maxval) {
accum += stepsize;
scount++;
}
}
if (vert) {
dim = Dimensions(scount,1);
} else {
dim = Dimensions(1,scount);
}
rp = (int64 *) allocateArray(FM_INT64,scount);
for (int i=0;i<scount;i++)
rp[i] = minval + i*stepsize;
return Array(FM_INT64,dim,rp);
}
void do_single_sided_algo_float(float a, float b,float *pvec, int adder, int count) {
double d = a;
for (int i=0;i<count;i++) {
pvec[i*adder] = (float) d;
d += b;
}
}
void do_double_sided_algo_float(float a, float b, float c, float *pvec, int adder, int count) {
if (count%2) {
do_single_sided_algo_float(a,b,pvec,adder,count/2);
do_single_sided_algo_float(c,-b,pvec+(count-1)*adder,-adder,count/2+1);
} else {
do_single_sided_algo_float(a,b,pvec,adder,count/2);
do_single_sided_algo_float(c,-b,pvec+(count-1)*adder,-adder,count/2);
}
}
Array Array::floatRangeConstructor(float minval, float stepsize,
float maxval, bool vert) {
Dimensions dim;
float *rp = NULL;
if (stepsize == 0) throw Exception("step size must be nonzero in colon expression");
//ideally, n = (c-a)/b
// But this really defines an interval... we let
// n_min = min(c-a)/max(b)
// n_max = max(c-a)/min(b)
// where min(x) = {y \in fp | |y| is max, |y| < |x|, sign(y) = sign(x)}
// max(x) = {y \in fp | |y| is min, |y| > |x|, sign(y) = sign(x)}
float ntest_min = nextafterf(maxval-minval,0)/nextafterf(stepsize,stepsize+stepsize);
float ntest_max = nextafterf(maxval-minval,maxval-minval+stepsize)/nextafterf(stepsize,0);
int npts = (int) floor(ntest_max);
bool use_double_sided = (ntest_min <= npts) && (npts <= ntest_max);
npts++;
if (npts < 0) npts = 0;
if (vert) {
dim = Dimensions(npts,1);
} else {
dim = Dimensions(1,npts);
}
rp = (float *) allocateArray(FM_FLOAT,npts);
if (use_double_sided)
do_double_sided_algo_float(minval,stepsize,maxval,rp,1,npts);
else
do_single_sided_algo_float(minval,stepsize,rp,1,npts);
return Array(FM_FLOAT,dim,rp);
}
void do_single_sided_algo_double(double a, double b,double *pvec, int adder, int count) {
double d = a;
for (int i=0;i<count;i++) {
pvec[i*adder] = (double) d;
d += b;
}
}
void do_double_sided_algo_double(double a, double b, double c, double *pvec, int adder, int count) {
if (count%2) {
do_single_sided_algo_double(a,b,pvec,adder,count/2);
do_single_sided_algo_double(c,-b,pvec+(count-1)*adder,-adder,count/2+1);
} else {
do_single_sided_algo_double(a,b,pvec,adder,count/2);
do_single_sided_algo_double(c,-b,pvec+(count-1)*adder,-adder,count/2);
}
}
Array Array::doubleRangeConstructor(double minval, double stepsize,
double maxval, bool vert) {
Dimensions dim;
double *rp = NULL;
if (stepsize == 0) throw Exception("step size must be nonzero in colon expression");
//ideally, n = (c-a)/b
// But this really defines an interval... we let
// n_min = min(c-a)/max(b)
// n_max = max(c-a)/min(b)
// where min(x) = {y \in fp | |y| is max, |y| < |x|, sign(y) = sign(x)}
// max(x) = {y \in fp | |y| is min, |y| > |x|, sign(y) = sign(x)}
double ntest_min = nextafter(maxval-minval,0)/nextafter(stepsize,stepsize+stepsize);
double ntest_max = nextafter(maxval-minval,maxval-minval+stepsize)/nextafter(stepsize,0);
int npts = (int) floor(ntest_max);
bool use_double_sided = (ntest_min <= npts) && (npts <= ntest_max);
npts++;
if (npts < 0) npts = 0;
if (vert) {
dim = Dimensions(npts,1);
} else {
dim = Dimensions(1,npts);
}
rp = (double *) allocateArray(FM_DOUBLE,npts);
if (use_double_sided)
do_double_sided_algo_double(minval,stepsize,maxval,rp,1,npts);
else
do_single_sided_algo_double(minval,stepsize,rp,1,npts);
return Array(FM_DOUBLE,dim,rp);
}
Array Array::matrixConstructor(ArrayMatrix& m) {
Dimensions mat_dims;
Dimensions row_dims;
Class maxType, minType;
Class retType;
rvstring retNames;
Dimensions retDims;
void *dstPtr = NULL;
bool sparseArg = false;
try {
maxType = FM_FUNCPTR_ARRAY;
minType = FM_STRING;
ArrayMatrix::iterator i = m.begin();
bool firstNonzeroColumn = true;
bool firstNonzeroRow = true;
while (i != m.end()) {
ArrayVector ptr = (ArrayVector) *i;
for (int j=0;j<ptr.size();j++) {
const Array& d = ptr[j];
if (d.sparse())
sparseArg = true;
if (!d.isEmpty()) {
if (maxType < d.dataClass()) maxType = d.dataClass();
if (minType > d.dataClass()) minType = d.dataClass();
if (firstNonzeroColumn) {
/**
* For the first element in the row, we copy
* the dimensions into the row buffer.
*/
row_dims = d.dimensions();
/**
* Because we are concatenating in columns,
* we must have at least 1 column in the data.
* The dimension array can be thought of as being
* infinitely extendable with 1's for all higher
* dimensions. We will take advantage of this to
* force the elements to have at least one column.
*/
firstNonzeroColumn = false;
} else {
/**
* Cycle through the elements of the row.
* Check these elements against the dimension
* buffer, and update the second dimension (column
* count).
*/
if ((d.dimensions().getLength() != row_dims.getLength()) &&
(d.dimensions().getLength()>1))
throw Exception("Number of dimensions must match for each element in a row definition");
if (d.dimensions().get(0) != row_dims.get(0))
throw Exception("Mismatch in first dimension for elements in row definition");
for (int k=2;k<row_dims.getLength();k++)
if (d.dimensions().get(k) != row_dims.get(k))
throw Exception("Mismatch in dimension for elements in row definition");
row_dims.set(1,row_dims.get(1)+d.dimensions().get(1));
}
}
}
/**
* For the first row, we copy the row dimensions into
* the column buffer.
*/
if (!firstNonzeroColumn) {
if (firstNonzeroRow) {
mat_dims = row_dims;
firstNonzeroRow = false;
} else {
if (mat_dims.getLength() != row_dims.getLength())
throw Exception("Number of dimensions must match for each row in a matrix definition");
for (int k=1;k<row_dims.getLength();k++)
if (row_dims.get(k) != mat_dims.get(k))
throw Exception("Mismatch in dimension for rows in matrix definition");
mat_dims.set(0,mat_dims.get(0)+row_dims.get(0));
}
firstNonzeroColumn = true;
}
i++;
}
retType = maxType;
retDims = mat_dims;
/**
* Check for the special case of a struct array - if the min_type
* is not also a struct array, then we have an error condition.
* if min_type and max_type are the same, we have a special case.
*/
if ((maxType == FM_STRUCT_ARRAY) && (minType != maxType))
throw Exception("Cannot convert cell arrays to struct arrays");
/**
* If the final type is a struct array, we can get the field names from
* any of the elements in the ArrayMatrix pointer argument.
*/
if (maxType == FM_STRUCT_ARRAY) {
ArrayMatrix::iterator i = m.begin();
ArrayVector ptr = *i;
const Array& d = ptr.front();
retNames = d.fieldNames();
}
/**
* Check for the sparse case - if any of the elements are sparse,
* the output is sparse.
*/
if (sparseArg) {
if (retType < FM_INT32) retType = FM_INT32;
return Array(retType,retDims,
SparseMatrixConstructor(retType,
retDims.get(1),
m),
true);
}
/**
* At this point, we now have the dimensions of the output
* dataset, as well as the type of the output data.
*
* Next, we allocate space for the output array.
*/
dstPtr = allocateArray(retType,retDims.getElementCount(),retNames);
/**
* Now, we have to copy the data from the source
* vector to the destination vector. Consider first
* the matrix case. We have a matrix
* A = [B,C;D,E] - we start at (0,0), and copy in
* B, one column at a time.
*
* Start with a pointer to the origin.
*/
/**
* We only need to track the row and column displacement for each
* block.
*/
int row_corner = 0;
int column_corner = 0;
int dim_count(mat_dims.getLength());
Dimensions aptr(dim_count);
Dimensions bptr(dim_count);
int el_size;
int row_count;
i = m.begin();
while (i != m.end()) {
ArrayVector ptr = *i;
row_count = 0;
for (int j=0;j<ptr.size();j++) {
Array d(ptr[j]);
if (!d.isEmpty()) {
/**
* Initialize the aptr to [0,...,0]
*/
aptr.zeroOut();
/**
* Get the source pointer and a destination pointer
*/
int dstIndex;
int srcIndex;
bptr = d.dimensions();
row_count = d.dimensions().get(0);
// Promote d to our ultimate type
d.promoteType(retType,retNames);
srcIndex = 0;
while(aptr.inside(bptr)) {
/**
* Convert aptr to an index into the destination array.
* We do so by using the mapping vector to map aptr
* to an index, and then adjust the index using the row
* and column corners.
*/
dstIndex = mat_dims.mapPoint(aptr);
dstIndex += row_corner + column_corner*mat_dims.get(0);
/**
* Copy the column from src to dst
*/
d.copyElements(srcIndex,dstPtr,dstIndex,row_count);
/**
* Advance the column index
*/
aptr.incrementModulo(bptr,1);
/**
* Advance the source pointer
*/
srcIndex += row_count;
}
/**
* Update the column corner pointer by the width of d.
*/
column_corner += bptr.get(1);
}
}
/**
* Reset the column corner
*/
column_corner = 0;
/**
* Advance the row corner
*/
row_corner += row_count;
i++;
}
return Array(retType,retDims,dstPtr,false,retNames);
} catch (Exception &e) {
// Memory leak?
// freeArray(dstPtr);
throw e;
}
}
Array Array::cellConstructor(Array m) {
return cellConstructor(ArrayVector() << m);
}
Array Array::cellConstructor(ArrayVector& m) {
return cellConstructor(ArrayMatrix() << m);
}
Array Array::cellConstructor(ArrayMatrix& m) {
int columnCount, rowCount;
Array* qp = NULL;
try {
ArrayMatrix::iterator i = m.begin();
while (i != m.end()) {
ArrayVector ptr = *i;
/**
* If this is the first row in the matrix def, then we
* record its size in columnCount.
*/
if (i == m.begin())
columnCount = ptr.size();
else {
/**
* Otherwise, make sure the column counts are all the same...
*/
if (ptr.size() != columnCount)
throw Exception("Cell definition must have same number of elements in each row");
}
i++;
}
/**
* At this point, we know how many columns our cell array has,
* and the number of rows is also known (size of m). So, set
* up our dimensions, and allocate the output.
*/
rowCount = m.size();
Dimensions retDims(rowCount,columnCount);
/**
* Allocate storage space for the contents.
*/
qp = (Array *) allocateArray(FM_CELL_ARRAY,retDims.getElementCount());
Array *cp;
Array *sp;
/**
* Loop through the rows.
*/
sp = qp;
i = m.begin();
while (i != m.end()) {
ArrayVector ptr = *i;
cp = sp;
for (int j=0;j<ptr.size();j++) {
*cp = ptr[j];
cp += rowCount;
}
i++;
sp++;
}
return Array(FM_CELL_ARRAY,retDims,qp);
} catch (Exception &e) {
Free(qp);
throw e;
}
}
Array Array::structConstructor(rvstring fNames, ArrayVector& values) {
const Array* rptr;
Dimensions dims;
bool nonSingularFound;
int i, j;
Array *qp = NULL;
try {
if (fNames.size() != values.size())
throw Exception("Number of field names must match number of values in structure constructor.");
/**
* First, we have to make sure that each entry of "values" have
* 1. cell arrays of the same size,
* 2. single element cell arrays,
* 3. single values.
*/
nonSingularFound = false;
for (i=0;i<values.size();i++) {
/**
* Check the type of the entry. If its a non-cell array, then
* then ignore this entry.
*/
if (values[i].dataClass() == FM_CELL_ARRAY) {
/**
* This is a cell-array, so look for non-scalar cell-arrays.
*/
if (!values[i].isScalar()) {
if (!nonSingularFound) {
nonSingularFound = true;
dims = values[i].dimensions();
} else
if (!dims.equals(values[i].dimensions()))
throw Exception("Array dimensions of non-scalar entries must agree in structure construction.");
}
}
}
/**
* At this point we can construct the dimensions of the output.
*/
if (!nonSingularFound) {
dims = Dimensions(1,1);
}
/**
* The dimensions of the object have been identified. Set the
* dimensions of the object and the field names. Then allocate
* the space.
*/
qp = (Array*) allocateArray(FM_STRUCT_ARRAY,dims.getElementCount(),fNames);
/**
* Work through the values, and copy the values back one at a time.
*/
int length = dims.getElementCount();
int offset;
offset = 0;
for (j=0;j<length;j++)
for (i=0;i<fNames.size();i++) {
Array rval = values[i];
if (!rval.isEmpty() && rval.sparse())
throw Exception("sparse arrays not supported for struct constructor.");
rptr = (const Array*) rval.data();
if (rval.isEmpty())
qp[offset] = rval;
else {
if (rval.dataClass() == FM_CELL_ARRAY) {
if (rval.isScalar())
qp[offset] = rptr[0];
else
qp[offset] = rptr[j];
} else
qp[offset] = rval;
}
offset++;
}
return Array(FM_STRUCT_ARRAY,dims,qp,false,fNames);
} catch (Exception &e) {
Free(qp);
throw;
}
}
/********************************************************************************
* Get functions *
********************************************************************************/
/**
* Take the current variable, and return a new array consisting of
* the elements in source indexed by the index argument. Indexing
* is done using vector ordinals.
*/
Array Array::getVectorSubset(Array& index, Interpreter* m_eval) {
void *qp = NULL;
try {
if (index.isEmpty()) {
return Array(dataClass(),index.dimensions(),
NULL,false,fieldNames(),className());
}
if (isColonOperator(index) && isEmpty())
return Array::emptyConstructor();
if (isEmpty())
throw Exception("Cannot index into empty variable.");
// Check for a(:), which is mapped to a reshape operation
if (isColonOperator(index)) {
Array ret(*this);
Dimensions retDims(getLength(),1);
ret.reshape(retDims);
return ret;
}
index.toOrdinalType(m_eval);
Dimensions retdims;
if (isColumnVector() && index.isRowVector())
retdims = Dimensions(index.getLength(),1);
else
retdims = index.dimensions();
retdims.simplify();
if (sparse()) {
if (index.getLength() == 1) {
int indx;
indx = index.getContentsAsIntegerScalar()-1;
int row = indx % getDimensionLength(0);
int col = indx / getDimensionLength(0);
return Array(dataClass(),retdims,
GetSparseScalarElement(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
data(),
row+1, col+1),
false);
} else
return Array(dataClass(),retdims,
GetSparseVectorSubsets(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
data(),
(const indexType*)
index.data(),
index.getDimensionLength(0),
index.getDimensionLength(1)),
true);
}
//
// The output is the same size as the _index_, not the
// source variable (neat, huh?). But it inherits the
// type of the source variable.
//
// Bug 1221845 - there is an anomaly in the treatment of vectors
// If the source is a column vector, and the index is a row vector
// then the output is a column vector. This behaviour is an
// inconsistency in M that we will reproduce here.
//
int length = index.getLength();
qp = allocateArray(dataClass(),index.getLength(),fieldNames());
// Get a pointer to the index data set
const indexType *index_p = (const indexType *) index.data();
int bound = getLength();
int ndx;
for (int i=0;i<length;i++) {
ndx = (int) index_p[i] - 1;
if (ndx < 0 || ndx >= bound)
throw Exception("Index exceeds variable dimensions");
copyElements(ndx,qp,i,1);
}
return Array(dataClass(),retdims,qp,sparse(),fieldNames(),
className());
} catch (Exception &e) {
throw;
}
}
/**
* Given a vector of indexing arrays, convert them into
* index pointers. If a colon is encountered, it is
* preserved (the first one -- the remaining colon expressions
* are expanded out into vectors).
*/
constIndexPtr* ProcessNDimIndexes(bool preserveColons,
Dimensions dims,
ArrayVector& index,
bool& anyEmpty,
int& colonIndex,
Dimensions& outDims,
bool argCheck,
Interpreter* m_eval) {
int L = index.size();
constIndexPtr* outndx = (constIndexPtr *) Malloc(sizeof(constIndexPtr*)*L);
bool colonFound = false;
anyEmpty = false;
colonIndex = -1;
for (int i=0;i<index.size();i++) {
bool isColon = isColonOperator(index[i]);
if (!colonFound && isColon && preserveColons) {
colonFound = true;
colonIndex = i;
outndx[i] = NULL;
outDims.set(i,dims.get(i));
} else if (isColon) {
indexType* buildcolon = (indexType*) Malloc(sizeof(indexType)*dims.get(i));
for (int j=1;j<=dims.get(i);j++)
buildcolon[j-1] = (indexType) j;
outndx[i] = buildcolon;
outDims.set(i,dims.get(i));
} else if (index[i].isEmpty()) {
anyEmpty = true;
outndx[i] = NULL;
outDims.set(i,0);
} else {
index[i].toOrdinalType(m_eval);
if (argCheck && (index[i].getMaxAsIndex() > dims.get(i)))
throw Exception("index exceeds array bounds");
outndx[i] = (constIndexPtr) index[i].getDataPointer();
outDims.set(i,index[i].getLength());
}
}
return outndx;
}
bool allScalars(const ArrayVector& index) {
for (int i=0;i<index.size();i++)
if (!(index[i].isScalar() && (!index[i].isReferenceType()) &&
(!index[i].isString())))
return false;
return true;
}
static indexType scalarIndex(const Array& a, Interpreter* m_eval) {
switch(a.dataClass()) {
case FM_LOGICAL:
{
const logical *rp = (const logical *) a.data();
return (indexType) rp[0];
}
break;
case FM_STRING:
{
throw Exception("Cannot convert string data types to indices.");
}
break;
case FM_DCOMPLEX:
{
m_eval->warningMessage("Imaginary part of complex index ignored.\n");
// We convert complex values into real values
const double *rp = (const double *) a.data();
return (indexType) rp[0];
}
break;
case FM_COMPLEX:
{
m_eval->warningMessage("Imaginary part of complex index ignored.\n");
// We convert complex values into real values
const float *rp = (const float *) a.data();
return (indexType) rp[0];
}
break;
case FM_DOUBLE:
{
const double *rp = (const double *) a.data();
if (rp[0] != (indexType) rp[0])
m_eval->warningMessage("Fractional part of index ignored.\n");
return (indexType) rp[0];
}
break;
case FM_FLOAT:
{
const float *rp = (const float *) a.data();
if (rp[0] != (indexType) rp[0])
m_eval->warningMessage("Fractional part of index ignored.\n");
return (indexType) rp[0];
}
break;
case FM_INT64:
{
const int64 *rp = (const int64 *) a.data();
return (indexType) rp[0];
}
break;
case FM_UINT64:
{
const uint64 *rp = (const uint64 *) a.data();
return (indexType) rp[0];
}
break;
case FM_INT32:
{
const int32 *rp = (const int32 *) a.data();
return (indexType) rp[0];
}
break;
case FM_UINT32:
{
const uint32 *rp = (const uint32 *) a.data();
return (indexType) rp[0];
}
break;
case FM_INT16:
{
const int16 *rp = (const int16 *) a.data();
return (indexType) rp[0];
}
break;
case FM_UINT16:
{
const uint16 *rp = (const uint16 *) a.data();
return (indexType) rp[0];
}
break;
case FM_INT8:
{
const int8 *rp = (const int8 *) a.data();
return (indexType) rp[0];
}
break;
case FM_UINT8:
{
const uint8 *rp = (const uint8 *) a.data();
return (indexType) rp[0];
}
break;
case FM_CELL_ARRAY:
{
throw Exception("Cannot convert cell arrays to indices.");
}
break;
case FM_STRUCT_ARRAY:
{
throw Exception("Cannot convert structure arrays to indices.");
}
break;
case FM_FUNCPTR_ARRAY:
{
throw Exception("Cannot convert function pointer arrays to indices.");
}
break;
}
}
Array Array::getNDimSubsetScalars(ArrayVector& index, Interpreter* m_eval) {
int ndx = 0;
int pagesize = 1;
for (int i=0;i<index.size();i++) {
indexType q(scalarIndex(index[i],m_eval));
if ((q < 1) || (q > getDimensionLength(i)))
throw Exception("index exceeds array bounds");
ndx += pagesize*(q-1);
pagesize *= getDimensionLength(i);
}
void *qp = allocateArray(dataClass(),1,fieldNames());
copyElements(ndx,qp,0,1);
return Array(dataClass(),Dimensions(1,1),qp,sparse(),fieldNames(),className());
}
void Array::setNDimSubsetScalars(ArrayVector& index, Array &rdata, Interpreter* m_eval) {
if (!rdata.isScalar()) throw Exception("rhs must be scalar for scalar set!");
if (sparse() || rdata.sparse()) throw Exception("sparse case not allowed for scalar set!");
if (dataClass() == FM_STRUCT_ARRAY) throw Exception("structure arrays not allowed for scalar set!");
if (dataClass() < rdata.dataClass()) throw Exception("type mismatch not allowed for scalar set!");
rdata.promoteType(dataClass());
int ndx = 0;
int pagesize = 1;
for (int i=0;i<index.size();i++) {
indexType q(scalarIndex(index[i],m_eval));
if ((q < 1) || (q > getDimensionLength(i)))
throw Exception("index exceeds array bounds");
ndx += pagesize*(q-1);
pagesize *= getDimensionLength(i);
}
rdata.copyElements(0,getReadWriteDataPointer(),ndx,1);
}
/**
* Take the current variable, and return a new array consisting of
* the elements in source indexed by the index argument. Indexing
* is done using ndimensional indices.
*/
Array Array::getNDimSubset(ArrayVector& index, Interpreter* m_eval) {
constIndexPtr* indx = NULL;
void *qp = NULL;
int i;
bool anyEmpty;
int colonIndex;
Dimensions myDims(dimensions());
if (!sparse() && allScalars(index))
return getNDimSubsetScalars(index, m_eval);
// if (isEmpty())
// throw Exception("Cannot index into empty variable.");
try {
int L = index.size();
Dimensions outDims(L);
indx = ProcessNDimIndexes(true,myDims,
index, anyEmpty,
colonIndex, outDims,
true, m_eval);
if (anyEmpty) {
Free(indx);
return Array::emptyConstructor();
}
if (sparse()) {
if (L > 2)
throw Exception("multidimensional indexing (more than 2 dimensions) not legal for sparse arrays");
if ((outDims.get(0) == 1) && (outDims.get(1) == 1))
return Array(dataClass(),outDims,
GetSparseScalarElement(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
data(),
*((const indexType*) indx[0]),
*((const indexType*) indx[1])),
false);
else
return Array(dataClass(),outDims,
GetSparseNDimSubsets(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
data(),
(const indexType*) indx[0],
outDims.get(0),
(const indexType*) indx[1],
outDims.get(1)),
true);
}
qp = allocateArray(dataClass(),outDims.getElementCount(),fieldNames());
int outDimsInt[maxDims];
int srcDimsInt[maxDims];
for (int i=0;i<L;i++) {
outDimsInt[i] = outDims.get(i);
srcDimsInt[i] = myDims.get(i);
}
outDims.simplify();
switch (dataClass()) {
case FM_COMPLEX:
getNDimSubsetNumericDispatchBurst<float>(colonIndex,(const float*) getDataPointer(),
(float*) qp,outDimsInt,srcDimsInt,
indx, L, 2);
break;
case FM_DCOMPLEX:
getNDimSubsetNumericDispatchBurst<double>(colonIndex,(const double*) getDataPointer(),
(double*) qp,outDimsInt,srcDimsInt,
indx, L, 2);
break;
case FM_LOGICAL:
getNDimSubsetNumericDispatchReal<logical>(colonIndex,(const logical*) getDataPointer(),
(logical*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_FLOAT:
getNDimSubsetNumericDispatchReal<float>(colonIndex,(const float*) getDataPointer(),
(float*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_DOUBLE:
getNDimSubsetNumericDispatchReal<double>(colonIndex,(const double*) getDataPointer(),
(double*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_INT8:
getNDimSubsetNumericDispatchReal<int8>(colonIndex,(const int8*) getDataPointer(),
(int8*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_UINT8:
getNDimSubsetNumericDispatchReal<uint8>(colonIndex,(const uint8*) getDataPointer(),
(uint8*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_INT16:
getNDimSubsetNumericDispatchReal<int16>(colonIndex,(const int16*) getDataPointer(),
(int16*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_UINT16:
getNDimSubsetNumericDispatchReal<uint16>(colonIndex,(const uint16*) getDataPointer(),
(uint16*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_INT32:
getNDimSubsetNumericDispatchReal<int32>(colonIndex,(const int32*) getDataPointer(),
(int32*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_UINT32:
getNDimSubsetNumericDispatchReal<uint32>(colonIndex,(const uint32*) getDataPointer(),
(uint32*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_INT64:
getNDimSubsetNumericDispatchReal<int64>(colonIndex,(const int64*) getDataPointer(),
(int64*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_UINT64:
getNDimSubsetNumericDispatchReal<uint64>(colonIndex,(const uint64*) getDataPointer(),
(uint64*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_STRING:
getNDimSubsetNumericDispatchReal<char>(colonIndex,(const char*) getDataPointer(),
(char*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_FUNCPTR_ARRAY:
getNDimSubsetNumericDispatchReal<FuncPtr>(colonIndex,
(const FuncPtr*) getDataPointer(),
(FuncPtr*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_CELL_ARRAY:
getNDimSubsetNumericDispatchReal<Array>(colonIndex,(const Array*) getDataPointer(),
(Array*) qp,outDimsInt,srcDimsInt,
indx, L);
break;
case FM_STRUCT_ARRAY:
getNDimSubsetNumericDispatchBurst<Array>(colonIndex,(const Array*) getDataPointer(),
(Array*) qp,outDimsInt,srcDimsInt,
indx, L, fieldNames().size());
break;
}
Free(indx);
outDims.simplify();
return Array(dataClass(),outDims,qp,sparse(),fieldNames(),
className());
} catch (Exception &e) {
Free(indx);
Free(qp);
throw e;
}
}
Array Array::getDiagonal(int diagonalOrder) {
if (!is2D())
throw Exception("Cannot take diagonal of N-dimensional array.");
int rows = dimensions().getRows();
int cols = dimensions().getColumns();
int outLen;
Dimensions outDims;
int i;
int srcIndex;
if (diagonalOrder < 0) {
outLen = (rows+diagonalOrder) < cols ? (rows+diagonalOrder) : cols;
outLen = (outLen < 0) ? 0 : outLen;
if (outLen == 0)
return Array::emptyConstructor();
outDims = Dimensions(outLen,1);
void *qp;
if (sparse()) {
qp = GetSparseDiagonal(dataClass(), rows, cols, data(), diagonalOrder);
return Array(dataClass(),outDims,qp,true);
} else {
qp = allocateArray(dataClass(),outLen,fieldNames());
for (i=0;i<outLen;i++) {
srcIndex = -diagonalOrder + i*(rows+1);
copyElements(srcIndex,qp,i,1);
}
return Array(dataClass(),outDims,qp,sparse(),fieldNames(),
className());
}
} else {
outLen = rows < (cols-diagonalOrder) ? rows : (cols-diagonalOrder);
outLen = (outLen < 0) ? 0 : outLen;
if (outLen == 0)
return Array::emptyConstructor();
outDims = Dimensions(outLen,1);
void *qp;
if (sparse()) {
qp = GetSparseDiagonal(dataClass(), rows, cols, data(), diagonalOrder);
return Array(dataClass(),outDims,qp,true);
} else {
qp = allocateArray(dataClass(),outLen,fieldNames());
for (i=0;i<outLen;i++) {
srcIndex = diagonalOrder*rows + i*(rows+1);
copyElements(srcIndex,qp,i,1);
}
}
return Array(dataClass(),outDims,qp,sparse(),fieldNames(),
className());
}
}
Array Array::getField(std::string fieldName) {
if (dataClass() != FM_STRUCT_ARRAY)
throw Exception("Attempt to apply field-indexing to non structure-array object.");
if (sparse())
throw Exception("getField not supported for sparse arrays.");
if (isEmpty()) return Array();
if (!isScalar())
throw Exception("Cannot apply structure indexing to arrays in this context.");
const Array *qp = (const Array*) data();
int N = getLength();
int fieldCount = fieldNames().size();
int ndx = getFieldIndex(fieldName);
if (ndx < 0)
throw Exception(std::string("Reference to non-existent field ") + fieldName);
return qp[ndx];
}
ArrayVector Array::getFieldAsList(std::string fieldName) {
if (dataClass() != FM_STRUCT_ARRAY)
throw Exception("Attempt to apply field-indexing to non structure-array object.");
if (sparse())
throw Exception("getFieldAsList not supported for sparse arrays.");
if (isEmpty()) return ArrayVector();
ArrayVector m;
const Array *qp = (const Array*) data();
int N = getLength();
int fieldCount = fieldNames().size();
int ndx = getFieldIndex(fieldName);
if (ndx < 0)
throw Exception(std::string("Reference to non-existent field ") + fieldName);
int i;
for (i=0;i<N;i++)
m.push_back(qp[i*fieldCount+ndx]);
return m;
}
Array Array::getVectorContents(Array& index, Interpreter* m_eval) {
if (dataClass() != FM_CELL_ARRAY)
throw Exception("Attempt to apply contents-indexing to non cell-array object.");
if (sparse())
throw Exception("getVectorContents not supported for sparse arrays.");
if (index.isEmpty()) return Array();
if (isEmpty()) return Array();
if (isColonOperator(index)) {
if (!isScalar()) throw Exception("Attempt to apply {:} to nonscalar argument in scalar context!");
// Get a pointer to our data
const Array* qp = (const Array*) data();
return qp[0];
}
if (!index.isScalar()) throw Exception("Cannot use A{ndx} when ndx is not a scalar in this context (too many right hand sides)");
index.toOrdinalType(m_eval);
// Get a pointer to the index data set
constIndexPtr index_p = (constIndexPtr) index.data();
if (index_p[0] > getLength())
throw Exception("Array index exceeds bounds of cell-array");
// Get a pointer to our data
const Array* qp = (const Array*) data();
return qp[index_p[0]-1];
}
/**
* Return a subset of a cell array as a list.
*/
ArrayVector Array::getVectorContentsAsList(Array& index, Interpreter* m_eval) {
ArrayVector m;
if (dataClass() != FM_CELL_ARRAY)
throw Exception("Attempt to apply contents-indexing to non cell-array object.");
if (sparse())
throw Exception("getVectorContentsAsList not supported for sparse arrays.");
if (index.isEmpty()) return ArrayVector();
if (isEmpty()) return ArrayVector();
if (isColonOperator(index)) {
int cnt = getLength();
// Get a pointer to our data
const Array* qp = (const Array*) data();
for (int i=0;i<cnt;i++)
m.push_back(qp[i]);
return m;
}
index.toOrdinalType(m_eval);
// Get the maximum index
indexType max_index = index.getMaxAsIndex();
// Get our length
int bound = getLength();
if (max_index > bound) throw Exception("Array index exceeds bounds of cell-array");
// Get the length of the index object
int index_length = index.getLength();
// Get a pointer to the index data set
constIndexPtr index_p = (constIndexPtr) index.data();
// Get a pointer to our data
const Array* qp = (const Array*) data();
// Now we copy data from dp to m
for (int i=0;i<index_length;i++)
m.push_back(qp[index_p[i]-1]);
return m;
}
/**
* Return the contents of an cell array as a list.
*/
Array Array::getNDimContents(ArrayVector& index, Interpreter* m_eval) {
if (dataClass() != FM_CELL_ARRAY)
throw Exception("Attempt to apply contents-indexing to non cell-array object.");
if (sparse())
throw Exception("getNDimContentsAsList not supported for sparse arrays.");
// Store the return value here
ArrayVector m;
// Get the number of indexing dimensions
int L = index.size();
Dimensions myDims(dimensions());
bool anyEmpty;
int colonIndex;
Dimensions outDims(L);
constIndexPtr* indx = ProcessNDimIndexes(false,myDims,index,anyEmpty,colonIndex,outDims,true,m_eval);
if (anyEmpty) {
Free(indx);
return Array::emptyConstructor();
}
if (!outDims.isScalar())
throw Exception("Too many values returned by A{N,M,...} expression");
Dimensions currentIndex(L);
const Array *qp = (const Array*) data();
for (int i=0;i<L;i++)
currentIndex.set(i,(int) indx[i][0] - 1);
Free(indx);
return qp[dimensions().mapPoint(currentIndex)];
}
/**
* Return the contents of an cell array as a list.
*/
ArrayVector Array::getNDimContentsAsList(ArrayVector& index, Interpreter* m_eval) {
if (dataClass() != FM_CELL_ARRAY)
throw Exception("Attempt to apply contents-indexing to non cell-array object.");
if (sparse())
throw Exception("getNDimContentsAsList not supported for sparse arrays.");
// Store the return value here
ArrayVector m;
// Get the number of indexing dimensions
int L = index.size();
Dimensions myDims(dimensions());
bool anyEmpty;
int colonIndex;
Dimensions outDims(L);
constIndexPtr* indx = ProcessNDimIndexes(false,myDims,index,anyEmpty,colonIndex,outDims,true,m_eval);
if (anyEmpty) {
Free(indx);
return singleArrayVector(Array::emptyConstructor());
}
Dimensions argPointer(L);
Dimensions currentIndex(L);
const Array *qp = (const Array*) data();
int srcindex;
while (argPointer.inside(outDims)) {
for (int i=0;i<L;i++)
currentIndex.set(i,(int) indx[i][argPointer.get(i)] - 1);
srcindex = dimensions().mapPoint(currentIndex);
m.push_back(qp[srcindex]);
argPointer.incrementModulo(outDims,0);
}
Free(indx);
return m;
}
/********************************************************************************
* Set functions *
********************************************************************************/
void Array::setValue(const Array &x) {
dp = x.dp;
}
/**
*
* This is the vector version of the multidimensional replacement function.
*
* This requires the following steps:
* 1. Compute the maximum along each dimension
* 2. Check that data is either scalar or the right size.
*/
void Array::setVectorSubset(Array& index, Array& rdata, Interpreter* m_eval) {
if (index.isEmpty())
return;
// Check the right-hand-side - if it is empty, then
// we have a delete command in disguise.
if (rdata.isEmpty()) {
deleteVectorSubset(index,m_eval);
return;
}
if (isColonOperator(index)) {
if (rdata.isScalar()) {
index = Array::int32RangeConstructor(1,1,getLength(),true);
} else {
Dimensions myDims;
if (!isEmpty())
myDims = dimensions();
else
myDims = rdata.dimensions();
if (myDims.getElementCount() != rdata.getLength())
throw Exception("assignment A(:) = B requires A and B to be the same size");
dp = rdata.dp;
reshape(myDims);
return;
}
}
// Make sure the index is an ordinal type
index.toOrdinalType(m_eval);
int index_length = index.getLength();
if (index_length == 0) return;
// Get a pointer to the index data set
constIndexPtr index_p = (constIndexPtr) index.data();
int advance;
// Set the right hand side advance pointer to
// - 0 if the rhs is a scalar
// - 1 else
if (rdata.sparse())
rdata.makeDense();
if (rdata.isScalar())
advance = 0;
else if (rdata.getLength() == index_length)
advance = 1;
else
throw Exception("Size mismatch in assignment A(I) = B.\n");
// Compute the maximum index;
indexType max_index = index.getMaxAsIndex();
// If the RHS type is superior to ours, we
// force our type to agree with the inserted data.
// Also, if we are empty, we promote ourselves (regardless of
// our type).
if (!isEmpty() &&
(rdata.dataClass() == FM_STRUCT_ARRAY) &&
(dataClass() == FM_STRUCT_ARRAY)) {
if (rdata.fieldNames().size() > fieldNames().size())
promoteType(FM_STRUCT_ARRAY,rdata.fieldNames());
else
rdata.promoteType(FM_STRUCT_ARRAY,fieldNames());
} else {
if (isEmpty() || rdata.dataClass() > dataClass())
promoteType(rdata.dataClass(),rdata.fieldNames());
// If our type is superior to the RHS, we convert
// the RHS to our type
else if (rdata.dataClass() <= dataClass())
rdata.promoteType(dataClass(),fieldNames());
}
// If the max index is larger than our current length, then
// we have to resize ourselves - but this is only legal if we are
// a vector.
if (sparse()) {
if (dataClass() == FM_LOGICAL)
rdata.promoteType(FM_UINT32);
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
void *qp = SetSparseVectorSubsets(dataClass(),rows,cols,data(),
(const indexType*) index.data(),
index.getDimensionLength(0),
index.getDimensionLength(1),
rdata.getDataPointer(),
advance);
setData(dataClass(),Dimensions(rows,cols),qp,true);
return;
}
vectorResize(max_index);
// Get a writable data pointer
void *qp = getReadWriteDataPointer();
// Now, we copy data from the RHS to our real part,
// computing indices along the way.
indexType srcIndex = 0;
indexType j;
for (int i=0;i<index_length;i++) {
j = index_p[i] - 1;
rdata.copyElements(srcIndex,qp,j,1);
srcIndex += advance;
}
}
/**
* Take the contents of data, and insert this data.
*
* This requires the following steps:
* 1. Compute the maximum along each dimension
* 2. Compute the dimensions of the right hand side
* 3. Check that data is either a scalar or the right size
* 4. If necessary, zero-extend the variable.
* 5. Copy in the result.
*
* This is true for integer arguments - not for logical ones.
* Logical indices need to be converted into integer lists
* before they can be used.
*/
void Array::setNDimSubset(ArrayVector& index, Array& rdata, Interpreter* m_eval) {
constIndexPtr* indx = NULL;
// If the RHS is empty, then we really want to do a delete...
if (rdata.isEmpty()) {
deleteNDimSubset(index,m_eval);
return;
}
// Check for all-scalar version
if (allScalars(index)) {
try {
setNDimSubsetScalars(index,rdata,m_eval);
return;
} catch (Exception &e) {
}
}
// If we are empty, then fill in the colon dimensions with the corresponding
// sizes of data - the problem is that which dimension to take isn't obvious.
// Need some more clarity here...
//
// If the rhs is a vector, we let the first colon operator map to the length
// of the vector. Otherwise, we map each colon operator to the corresponding
// dimension of the RHS.
if (isEmpty()) {
if (rdata.isVector()) {
bool firstcolon = true;
for (int i=0;i<index.size();i++)
if (isColonOperator(index[i]))
if (firstcolon) {
index[i] = Array::int32RangeConstructor(1,1,rdata.getLength(),true);
firstcolon = false;
} else
index[i] = Array::int32RangeConstructor(1,1,1,true);
} else {
// Correction -
// In the assignment of the form:
// g(2,:,4,:) = fs;
// we fill in the colons with the first and second dimensional sizes of fs
int colonDim = 0;
for (int i=0;i<index.size();i++)
if (isColonOperator(index[i]))
index[i] = Array::int32RangeConstructor(1,1,rdata.getDimensionLength(colonDim++),true);
if ((colonDim > 0) && (colonDim < rdata.dimensions().getLength()))
throw Exception("Assignment A(...) = b is not legal");
// Test for colonDim < rdata.get
}
}
try {
int L = index.size();
Dimensions myDims(dimensions());
Dimensions outDims;
bool anyEmpty;
int colonIndex;
indx = ProcessNDimIndexes(true,myDims,index,anyEmpty,colonIndex,outDims,false,m_eval);
if (anyEmpty) {
Free(indx);
return;
}
Dimensions a(L);
// First, we compute the maximum along each dimension.
int dataCount = 1;
for (int i=0;i<L;i++) {
if (isColonOperator(index[i])) {
a.set(i,myDims.get(i));
dataCount *= myDims.get(i);
} else {
a.set(i,index[i].getMaxAsIndex());
dataCount *= index[i].getLength();
}
}
// Next, we compute the dimensions of the right hand side
indexType advance;
if (rdata.sparse())
rdata.makeDense();
if (rdata.isScalar()) {
advance = 0;
} else if (!isEmpty() && (rdata.getLength() == dataCount)) {
advance = 1;
} else if (!isEmpty())
throw Exception("Size mismatch in assignment A(I1,I2,...,In) = B.\n");
else
advance = 1;
// If the RHS type is superior to ours, we
// force our type to agree with the inserted rdata.
if (!isEmpty() &&
(rdata.dataClass() == FM_STRUCT_ARRAY) &&
(dataClass() == FM_STRUCT_ARRAY)) {
if (rdata.fieldNames().size() > fieldNames().size())
promoteType(FM_STRUCT_ARRAY,rdata.fieldNames());
else
rdata.promoteType(FM_STRUCT_ARRAY,fieldNames());
} else {
if (isEmpty() || rdata.dataClass() > dataClass())
promoteType(rdata.dataClass(),rdata.fieldNames());
// If our type is superior to the RHS, we convert
// the RHS to our type
else if (rdata.dataClass() <= dataClass())
rdata.promoteType(dataClass(),fieldNames());
}
// Now, resize us to fit this data
resize(a);
if (sparse()) {
if (dataClass() == FM_LOGICAL)
rdata.promoteType(FM_UINT32);
if (L > 2)
throw Exception("multidimensional indexing (more than 2 dimensions) not legal for sparse arrays in assignment A(I1,I2,...,IN) = B");
SetSparseNDimSubsets(dataClass(),
getDimensionLength(0),
dp->getWriteableData(),
(const indexType*) indx[0],
outDims.get(0),
(const indexType*) indx[1],
outDims.get(1),
rdata.getDataPointer(),advance);
return;
}
myDims = dimensions();
// Get a writable data pointer
void *qp = getReadWriteDataPointer();
int outDimsInt[maxDims];
int srcDimsInt[maxDims];
for (int i=0;i<L;i++) {
outDimsInt[i] = outDims.get(i);
srcDimsInt[i] = myDims.get(i);
}
outDims.simplify();
switch (dataClass()) {
case FM_COMPLEX:
setNDimSubsetNumericDispatchBurst<float>(colonIndex,(float*) qp,
(const float*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L, 2,advance);
break;
case FM_DCOMPLEX:
setNDimSubsetNumericDispatchBurst<double>(colonIndex,(double*) qp,
(const double*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L, 2,advance);
break;
case FM_LOGICAL:
setNDimSubsetNumericDispatchReal<logical>(colonIndex,(logical*) qp,
(const logical*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_FLOAT:
setNDimSubsetNumericDispatchReal<float>(colonIndex,(float*) qp,
(const float*) rdata.getDataPointer()
,outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_DOUBLE:
setNDimSubsetNumericDispatchReal<double>(colonIndex,(double*) qp,
(const double*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_INT8:
setNDimSubsetNumericDispatchReal<int8>(colonIndex,(int8*) qp,
(const int8*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_UINT8:
setNDimSubsetNumericDispatchReal<uint8>(colonIndex,(uint8*) qp,
(const uint8*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_INT16:
setNDimSubsetNumericDispatchReal<int16>(colonIndex,(int16*) qp,
(const int16*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_UINT16:
setNDimSubsetNumericDispatchReal<uint16>(colonIndex,(uint16*) qp,
(const uint16*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_INT32:
setNDimSubsetNumericDispatchReal<int32>(colonIndex,(int32*) qp,
(const int32*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_UINT32:
setNDimSubsetNumericDispatchReal<uint32>(colonIndex,(uint32*) qp,
(const uint32*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_INT64:
setNDimSubsetNumericDispatchReal<int64>(colonIndex,(int64*) qp,
(const int64*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_UINT64:
setNDimSubsetNumericDispatchReal<uint64>(colonIndex,(uint64*) qp,
(const uint64*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_STRING:
setNDimSubsetNumericDispatchReal<char>(colonIndex,(char*) qp,
(const char*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_FUNCPTR_ARRAY:
setNDimSubsetNumericDispatchReal<FuncPtr>(colonIndex,
(FuncPtr*) qp,
(const FuncPtr*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_CELL_ARRAY:
setNDimSubsetNumericDispatchReal<Array>(colonIndex,(Array*) qp,
(const Array*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L,advance);
break;
case FM_STRUCT_ARRAY:
setNDimSubsetNumericDispatchBurst<Array>(colonIndex,(Array*) qp,
(const Array*) rdata.getDataPointer(),
outDimsInt,srcDimsInt,
indx, L, fieldNames().size(),advance);
break;
}
Free(indx);
Dimensions q(dimensions()); q.simplify(); dp->setDimensions(q);
} catch (Exception &e) {
Free(indx);
throw e;
}
}
/**
* This is the vector version of the multidimensional cell-replacement function.
* This is for content-based indexing (curly brackets). Two points that make
* this function different than replaceData are
* 1. If the index is larger than the size, we resize to a vector of sufficient
* length.
* 2. Deletions do not occur.
*/
void Array::setVectorContentsAsList(Array& index, ArrayVector& rdata, Interpreter* m_eval) {
if (sparse())
throw Exception("setVectorContentsAsList not supported for sparse arrays.");
promoteType(FM_CELL_ARRAY);
if (isColonOperator(index)) {
if (getLength() > rdata.size())
throw Exception("Not enough right hand side values to satisy left hand side expression.");
Array *qp = (Array*) getReadWriteDataPointer();
for (int i=0;i<getLength();i++) {
qp[i] = rdata.front();
rdata.pop_front();
}
Dimensions q(dimensions()); q.simplify(); dp->setDimensions(q);
return;
}
index.toOrdinalType(m_eval);
if (rdata.size() < index.getLength())
throw Exception("Not enough right hand side values to satisy left hand side expression.");
// Get the maximum index
indexType max_index = index.getMaxAsIndex();
// Resize us as necessary.
vectorResize(max_index);
// Get a pointer to the dataset
Array *qp = (Array*) getReadWriteDataPointer();
// Get a pointer to the index data set
constIndexPtr index_p = (constIndexPtr) index.data();
// Get the length of the index object
int index_length = index.getLength();
// Copy in the data
for (int i=0;i<index_length;i++) {
int ndx = index_p[i]-1;
qp[ndx] = rdata.front();
rdata.pop_front();
}
Dimensions q(dimensions()); q.simplify(); dp->setDimensions(q);
}
/**
* This is the multidimensional cell-replacement function.
* This is for content-based indexing (curly brackets).
*/
void Array::setNDimContentsAsList(ArrayVector& index, ArrayVector& rdata, Interpreter* m_eval) {
if (sparse())
throw Exception("setNDimContentsAsList not supported for sparse arrays.");
promoteType(FM_CELL_ARRAY);
Dimensions myDims(dimensions());
Dimensions outDims;
bool anyEmpty;
int colonIndex;
// Convert the indexing variables into an ordinal type.
constIndexPtr* indx = ProcessNDimIndexes(false,myDims,index,anyEmpty,colonIndex,outDims,false,m_eval);
int L = index.size();
try {
Dimensions a(L);
// First, we compute the maximum along each dimension.
// We also get pointers to each of the index pointers.
int i;
for (i=0;i<L;i++)
if (isColonOperator(index[i]))
a.set(i,myDims.get(i));
else
a.set(i,index[i].getMaxAsIndex());
// Next, we compute the number of entries in each component.
Dimensions argLengths(L);
Dimensions argPointer(L);
int dataCount = 1;
for (i=0;i<L;i++) {
if (isColonOperator(index[i])) {
argLengths.set(i,myDims.get(i));
dataCount *= myDims.get(i);
} else {
argLengths.set(i,index[i].getLength());
dataCount *= argLengths.get(i);
}
}
if (rdata.size() < dataCount)
throw Exception("Not enough right hand side values to satisfy left hand side expression");
// Resize us as necessary
resize(a);
// Get a writable data pointer
Array *qp = (Array*) getReadWriteDataPointer();
// Now, we copy data from dp to our real part,
// computing indices along the way.
Dimensions currentIndex(dimensions().getLength());
indexType j;
while (argPointer.inside(argLengths)) {
for (i=0;i<L;i++)
currentIndex.set(i,(indexType) indx[i][argPointer.get(i)] - 1);
j = dimensions().mapPoint(currentIndex);
qp[j] = rdata.front();
rdata.pop_front();
argPointer.incrementModulo(argLengths,0);
}
Free(indx);
indx = NULL;
Dimensions q(dimensions()); q.simplify(); dp->setDimensions(q);
} catch (Exception &e) {
Free(indx);
throw e;
}
}
/**
* Set the contents of a field in a structure.
*/
void Array::setFieldAsList(std::string fieldName, ArrayVector& rdata) {
Array *rp = NULL;
if (isEmpty()) {
rvstring names(fieldNames());
names.push_back(fieldName);
promoteType(FM_STRUCT_ARRAY,names);
Dimensions a(1,1);
resize(a);
// setData(FM_STRUCT_ARRAY,dimensions(),NULL,names);
// return;
}
if (sparse())
throw Exception("setFieldAsList not supported for sparse arrays.");
if (dataClass() != FM_STRUCT_ARRAY)
throw Exception("Cannot apply A.field_name = B to non struct-array object A.");
if (rdata.size() < getLength())
throw Exception("Not enough right hand values to satisfy left hand side expression.");
int indexLength = getLength();
int field_ndx = getFieldIndex(fieldName);
if (field_ndx == -1)
field_ndx = insertFieldName(fieldName);
int fieldCount = fieldNames().size();
Array *qp = (Array*) getReadWriteDataPointer();
for (int i=0;i<indexLength;i++) {
qp[i*fieldCount+field_ndx] = rdata.front();
rdata.pop_front();
}
Dimensions q(dimensions()); q.simplify(); dp->setDimensions(q);
}
/**
* Add another fieldname to our structure array.
*/
int Array::insertFieldName(std::string fieldName) {
if (sparse())
throw Exception("insertFieldName not supported for sparse arrays.");
rvstring names(fieldNames());
names.push_back(fieldName);
const Array* qp = (const Array*) data();
Array *rp = (Array*) allocateArray(dataClass(),getLength(),names);
int fN = names.size();
for (int i=0;i<fN-1;i++)
rp[i] = qp[i];
setData(FM_STRUCT_ARRAY,dimensions(),rp,false,names,
className());
return (fN-1);
}
/********************************************************************************
* Delete functions *
********************************************************************************/
/**
* Delete a vector subset of a variable.
*/
void Array::deleteVectorSubset(Array& arg, Interpreter* m_eval) {
void *qp = NULL;
bool *deletionMap = NULL;
try {
// First convert arg to an ordinal type.
if (isColonOperator(arg)) {
setData(dataClass(),Dimensions(0,0),NULL,false,fieldNames(),className());
return;
}
arg.toOrdinalType(m_eval);
if (sparse()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
void *cp = DeleteSparseMatrixVectorSubset(dataClass(),rows,cols,
data(),
(const indexType *)
arg.getDataPointer(),
arg.getLength());
setData(dataClass(),Dimensions(rows,cols),cp,true);
return;
}
// Next, build a deletion map.
int N = getLength();
int i;
deletionMap = arg.getBinaryMap(N);
// Now, we count up the number of elements that remain after deletion.
int newSize = 0;
for (i=0;i<N;i++)
if (!deletionMap[i]) newSize++;
// Special case - if newSize==getLength, the delete is a no-op
if (newSize == getLength()) return;
// Allocate a new space to hold the rdata.
qp = allocateArray(dataClass(),newSize,fieldNames());
// Loop through the indices - copy elements in that
// have not been deleted.
int dstIndex = 0;
for (i=0;i<N;i++)
if (!deletionMap[i])
copyElements(i,qp,dstIndex++,1);
Free(deletionMap);
deletionMap = NULL;
Dimensions newDim;
if (dimensions().isScalar()) {
newDim.reset();
newDim = Dimensions(1,newSize);
} else if (dimensions().isVector()) {
newDim = dimensions();
if (dimensions().get(0) != 1)
newDim.set(0,newSize);
else
newDim.set(1,newSize);
} else {
newDim = Dimensions(1,newSize);
}
setData(dataClass(),newDim,qp,sparse(),
fieldNames(),className());
} catch (Exception &e) {
Free(qp);
Free(deletionMap);
throw e;
}
}
void Array::makeSparse() {
if (!is2D())
throw Exception("Cannot make n-dimensional arrays sparse.");
if (isEmpty())
return;
if (isReferenceType() || isString())
throw Exception("Cannot make strings or reference types sparse.");
if (sparse()) return;
if ((dataClass() != FM_LOGICAL) && (dataClass() < FM_INT32))
promoteType(FM_INT32);
setData(dataClass(),dimensions(),
makeSparseArray(dataClass(),
dimensions().get(0),
dimensions().get(1),
data()),
true,
fieldNames(),
className());
}
int Array::getNonzeros() const {
if (!sparse())
return (dimensions().getElementCount());
if (isEmpty())
return 0;
return CountNonzeros(dataClass(),
dimensions().get(0),
dimensions().get(1),
data());
}
void Array::makeDense() {
if (!sparse())
return;
if (isEmpty()) {
dp->setSparse(false);
return;
}
setData(dataClass(),dimensions(),
makeDenseArray(dataClass(),
dimensions().get(0),
dimensions().get(1),
data()),
false,
fieldNames(),
className());
}
/**
* Delete a subset of a variable.
*/
void Array::deleteNDimSubset(ArrayVector& args, Interpreter* m_eval) {
int singletonReferences = 0;
int singletonDimension = 0;
int i;
Array qp;
bool *indxCovered = NULL;
bool *deletionMap = NULL;
void *cp = NULL;
if (isEmpty()) return;
try {
// Our strategy is as follows. To make the deletion, we need
// one piece of information: the dimension to delete.
// To do so, we first make a pass through the set of arguments,
// checking each one to see if it "covers" its index set.
//
// However, to simplify the testing of
// conditions later on, we must make sure that the length of
// the index list matches our number of dimensions. We extend
// it using 1 references, and throw an exception if there are
// more indices than our dimension set.
for (i=0;i<args.size();i++) {
if (isColonOperator(args[i]))
args[i] = Array::int32RangeConstructor(1,1,dimensions().get(i),true);
args[i].toOrdinalType(m_eval);
}
// First, add enough "1" singleton references to pad the
// index set out to the size of our variable.
if (args.size() < dimensions().getLength())
for (i = args.size();i<dimensions().getLength();i++)
args.push_back(Array::uint32Constructor(1));
// Now cycle through indices one at a time. Count
// the number of non-covering indices. Also track the
// location of the last-occurring non-covering index.
for (i=0;i<args.size();i++) {
qp = args[i];
// Get a binary representation of each index over the range [0,dimensions[i]-1]
indxCovered = qp.getBinaryMap(dimensions().get(i));
// Scan the array, and make sure all elements are true. If not,
// then this is the "singleton" dimension. Kick the singleton
// reference counter, and record the current dimension.
bool allCovered = true;
for (int k=0;allCovered && (k<dimensions().get(i));k++)
allCovered = allCovered && indxCovered[k];
Free(indxCovered);
indxCovered = NULL;
if (!allCovered) {
singletonReferences++;
singletonDimension = i;
}
}
// Now, we check the number of singleton references we
// encountered. There are three cases to check:
// Case 1. No singleton references. This is OK - it
// means we wish to delete the entire variable.
// We set a flag, and proceed to validate the covering
// of each dimension by its corresponding index set.
// Case 2. One singleton reference. This is OK - it
// means we wish to delete a single plane of
// rdata. Retrieve the index of the plane, and store
// it in singletonIndex.
// Case 3. Two or more singleton references. Can't do it -
// throw an error.
if (singletonReferences > 1)
throw Exception("Statement A(...) = [] can only contain one non-colon index.\n");
if (singletonReferences == 0)
singletonDimension = -1;
// If we got this far, the user either entered an expression like
// A(:,:,...,:,s,:,...,:) = [], or something numerically equivalent,
// or the user entered something like A(:,...,:) = [].
// In the latter case (indicated by singletonReferences = 0), we simply
// delete the entire variable, and make it an empty type.
// In the former case, we will have more work to do...
if (singletonReferences != 0) {
// We have to rescan our (now-identified) singleton
// dimension to build a deletion map. The map is
// marked true for each plane we wish to delete.
// The map is the size of the _data_'s dimension.
int M = dimensions().get(singletonDimension);
deletionMap = args[singletonDimension].getBinaryMap(M);
// We can now calculate the new size of the variable in the singletonDimension
// by counting the number of "false" entries in deletionMap.
int newSize = 0;
for (i=0;i<M;i++)
if (!deletionMap[i]) newSize++;
int rowCount = dimensions().get(0);
Dimensions retDims;
// Copy our current dimensions to the output dimensions.
retDims = dimensions();
// Update the singleton dimension to the new size.
retDims.set(singletonDimension,newSize);
// For sparse matrices, we branch here to call the sparse matrix deletion code
if (sparse()) {
int rows = getDimensionLength(0);
int cols = getDimensionLength(1);
if (singletonDimension == 0)
setData(dataClass(),retDims,
DeleteSparseMatrixRows(dataClass(),rows,cols,
data(),deletionMap),true);
else if (singletonDimension == 1)
setData(dataClass(),retDims,
DeleteSparseMatrixCols(dataClass(),cols,
data(),deletionMap),true);
else
throw Exception("sparse matrices do not support deleting n-dimensional planes - they are only 2-D");
Free(deletionMap);
Free(indxCovered);
return;
}
// Allocate space for the return objects data
cp = allocateArray(dataClass(),retDims.getElementCount(),fieldNames());
// Track our offset into the original data & our offset into
// the truncated rdata.
int srcIndex = 0;
int dstIndex = 0;
// Inintialize an ND pointer to the first element in the
// current data structure.
int L = dimensions().getLength();
Dimensions curPos(L);
// Loop until we have exhausted the original data.
while (curPos.inside(dimensions())) {
// Check to see if this column is to be skipped
if (!deletionMap[curPos.get(singletonDimension)]) {
// Copy the data from our original data structure to the
// new data structure, starting from srcIndex, and
// copying to dstIndex.
copyElements(srcIndex,cp,dstIndex,1);
// Advance the destination pointer. - we only do this on a copy
dstIndex ++;
}
// Advance the source pointer - we always do this
srcIndex ++;
curPos.incrementModulo(dimensions(),0);
}
Free(deletionMap);
retDims.simplify();
setData(dataClass(),retDims,cp,sparse(),
fieldNames(),className());
} else {
Dimensions newDims;
setData(dataClass(),newDims,NULL,false,fieldNames(),
className());
}
} catch (Exception &e) {
Free(deletionMap);
Free(cp);
Free(indxCovered);
throw e;
}
}
string Array::getContentsAsString() const {
if (dataClass() != FM_STRING)
throw Exception("Unable to convert supplied object to a string!\n");
int M = getLength();
return string((const char*) data(),M);
}
string Array::getContentsAsStringUpper() const {
string ret(getContentsAsString());
std::transform(ret.begin(),ret.end(),ret.begin(),(int(*)(int))toupper);
return ret;
}
int32 Array::getContentsAsIntegerScalar() {
if (dataClass() == FM_STRING) {
return atoi(getContentsAsString().c_str());
}
int32 *qp;
if (getLength() != 1)
throw Exception("Expected a scalar!\n");
promoteType(FM_INT32);
qp = (int32*) data();
return (*qp);
}
double Array::getContentsAsDoubleScalar() {
if (dataClass() == FM_STRING) {
return atof(getContentsAsString().c_str());
}
double *qp;
if (isComplex() || isReferenceType() || isString())
throw Exception("Expected a real valued scalar");
promoteType(FM_DOUBLE);
qp = (double*) data();
return (*qp);
}
bool Array::isDataClassReferenceType(Class cls) {
return (cls == FM_CELL_ARRAY || cls == FM_STRUCT_ARRAY ||
cls == FM_FUNCPTR_ARRAY);
}
template <class T>
int32 DoCountNNZReal(const void* dp, int len) {
int32 accum = 0;
const T* cp = (const T*) dp;
for (int i=0;i<len;i++)
if (cp[i]) accum++;
return accum;
}
template <class T>
int32 DoCountNNZComplex(const void* dp, int len) {
int32 accum = 0;
const T* cp = (const T*) dp;
for (int i=0;i<len;i++)
if (cp[2*i] || cp[2*i+1]) accum++;
return accum;
}
int32 Array::nnz() const {
if (isEmpty()) return 0;
if (sparse())
return CountNonzeros(dataClass(),
getDimensionLength(0),
getDimensionLength(1),
data());
// OK - its not sparse... now what?
switch (dataClass()) {
case FM_LOGICAL:
return DoCountNNZReal<logical>(data(),getLength());
case FM_INT8:
return DoCountNNZReal<int8>(data(),getLength());
case FM_UINT8:
case FM_STRING:
return DoCountNNZReal<uint8>(data(),getLength());
case FM_INT16:
return DoCountNNZReal<int16>(data(),getLength());
case FM_UINT16:
return DoCountNNZReal<uint16>(data(),getLength());
case FM_INT32:
return DoCountNNZReal<int32>(data(),getLength());
case FM_UINT32:
return DoCountNNZReal<uint32>(data(),getLength());
case FM_INT64:
return DoCountNNZReal<int64>(data(),getLength());
case FM_UINT64:
return DoCountNNZReal<uint64>(data(),getLength());
case FM_FLOAT:
return DoCountNNZReal<float>(data(),getLength());
case FM_DOUBLE:
return DoCountNNZReal<double>(data(),getLength());
case FM_COMPLEX:
return DoCountNNZComplex<float>(data(),getLength());
case FM_DCOMPLEX:
return DoCountNNZComplex<double>(data(),getLength());
case FM_CELL_ARRAY:
return DoCountNNZReal<void*>(data(),getLength());
case FM_STRUCT_ARRAY:
return DoCountNNZReal<void*>(data(),getLength()*fieldNames().size());
case FM_FUNCPTR_ARRAY:
return DoCountNNZReal<void*>(data(),getLength());
}
}
bool Array::anyNotFinite() {
if (sparse())
return SparseAnyNotFinite(dataClass(),
getDimensionLength(1),
data());
switch(dataClass()) {
case FM_FLOAT:
{
const float *sp = (const float *) data();
int len = getLength();
int i=0;
while (i<len)
if (!IsFinite(sp[i++])) return true;
return false;
}
case FM_DOUBLE:
{
const double *sp = (const double *) data();
int len = getLength();
int i=0;
while (i<len)
if (!IsFinite(sp[i++])) return true;
return false;
}
case FM_COMPLEX:
{
const float *sp = (const float *) data();
int len = getLength();
int i=0;
while (i<2*len)
if (!IsFinite(sp[i++])) return true;
return false;
}
case FM_DCOMPLEX:
{
const double *sp = (const double *) data();
int len = getLength();
int i=0;
while (i<2*len)
if (!IsFinite(sp[i++])) return true;
return false;
}
default:
return false;
}
}
int32 ArrayToInt32(const Array& a) {
Array b(a);
return b.getContentsAsIntegerScalar();
}
Array Array::doubleMatrixConstructor(int rows, int cols) {
Dimensions dim(rows,cols);
double *data = (double*) allocateArray(FM_DOUBLE,rows*cols);
return Array(FM_DOUBLE,dim,data,false);
}
Array Array::dcomplexMatrixConstructor(int rows, int cols) {
Dimensions dim(rows,cols);
double *data = (double*) allocateArray(FM_DCOMPLEX,rows*cols);
return Array(FM_DCOMPLEX,dim,data,false);
}
string ArrayToString(const Array& a) {
Array b(a);
return b.getContentsAsString();
}
Array ToSingleArray(const ArrayVector& a) {
if (a.size()>0)
return a.front();
else
return Array::emptyConstructor();
}
double ArrayToDouble(const Array& a) {
Array b(a);
return b.getContentsAsDoubleScalar();
}
uint32 TypeSize(Class cls) {
switch (cls) {
case FM_LOGICAL:
return sizeof(logical);
case FM_UINT8:
case FM_INT8:
return sizeof(int8);
case FM_UINT16:
case FM_INT16:
return sizeof(int16);
case FM_UINT32:
case FM_INT32:
return sizeof(int32);
case FM_UINT64:
case FM_INT64:
return sizeof(int64);
case FM_FLOAT:
return sizeof(float);
case FM_DOUBLE:
return sizeof(double);
}
throw Exception("Unsupported class as argument to TypeSize");
}
string operator+(string a, int d) {
char buf[256];
sprintf(buf,"%d",d);
return a+string(buf);
}
string operator+(int d, string a) {
char buf[256];
sprintf(buf,"%d",d);
return a+string(buf);
}
stringVector operator+(stringVector a, stringVector b) {
for (int i=0;i<b.size();i++)
a.push_back(b[i]);
return a;
}
Array CellArrayFromQStringList(QStringList t) {
ArrayVector tVec;
for (int i=0;i<t.size();i++)
tVec << Array::stringConstructor(t[i].toStdString());
return Array::cellConstructor(tVec);
}
Array Uint32VectorFromQList(QList<uint32> p) {
Array ret(Array::uint32VectorConstructor(p.size()));
for (int i=0;i<p.size();i++)
((uint32*) ret.getDataPointer())[i] = p[i];
return ret;
}
syntax highlighted by Code2HTML, v. 0.9.1