testing updates

build/make.linux64.lib.py

@@ -15,7 +15,7 @@ builddir = "./build_linux64"
 
 doinstall = True   #copies the build_output to the install dir when finished
 cc = "g++" #compiler
-cflags = "-fPIC"
+cflags = "-fPIC -O3"
 libraries = "-l{}".format(libname)
 libdirs = "-L{} -L{}/lib".format(builddir,commondir)
 linkerflags = "-static -static-libgcc -Wl,-rpath=."

include/amscppnarray/amscppnarray.hpp

@@ -10,6 +10,7 @@
 #include <new>
 #include <initializer_list>
 #include <thread>
+#include <algorithm>
 
 #include <amsmath/amsmath.hpp>
 #include <amsutil/amsutil.hpp>
@@ -88,8 +89,17 @@ template<typename T> class narray
     narray<T> subarray(const narray_size_t ind1, const narray_size_t ind2) const;
 
     //permutations and sorting
+    
 
-
+    //C++11 iterators
+    typedef T* iterator;
+    typedef const T* const_iterator;
+    iterator begin() {return data;}
+    iterator end() {return data + length;}
+    const_iterator begin() const {return data;}
+    const_iterator end() const {return data + length;}
+    const_iterator cbegin() const { return data; }
+    const_iterator cend() const { return data+length;}
 
 };
 
@@ -101,6 +111,8 @@ template<typename T> narray<T> ones(const narray_size_t N);
 narray<float> linspacef(const float low, const float high, narray_size_t N);
 narray<double> linspace(const double low, const double high, narray_size_t N);
 
+//cast/conversion
+
 
 /////////
 //Tests//
@@ -113,6 +125,8 @@ void test_narray3();
 }; //end namespace narray
 }; //end namespace ams
 
+#include <amscppnarray/amscppnarray_llbops.hpp>
+
 #include <amscppnarray/amscppnarray_impl.hpp>
 #include <amscppnarray/amscppnarray_imploper.hpp>
 #include <amscppnarray/amscppnarray_implsubarray.hpp>
@@ -121,4 +135,6 @@ void test_narray3();
 #include <amscppnarray/amscppnarray_random.hpp>
 #include <amscppnarray/amscppnarray_math.hpp>
 
+#include <amscppnarray/amscppnarray_binrw.hpp>
+
 #endif

include/amscppnarray/amscppnarray_binrw.hpp

@@ -0,0 +1,21 @@
+#ifndef __AMSCPPNARRAY_BINRW_HPP__
+#define __AMSCPPNARRAY_BINRW_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+//binary array read/write operations
+
+template<typename T> int fwrite_ndarray(FILE *fp, narray<T> *data, narray<int32_t> *shape);
+
+template<typename T> int fread_ndarray(FILE *fp, narray<T> *data, narray<int32_t> *shape);
+
+
+};
+};
+
+#include <amscppnarray/amscppnarray_binrwimpl.hpp>
+
+#endif

include/amscppnarray/amscppnarray_binrwimpl.hpp

@@ -0,0 +1,153 @@
+#ifndef __AMSCPPNARRAY_BINRWIMPL_HPP__
+#define __AMSCPPNARRAY_BINRWIMPL_HPP__
+
+
+
+namespace ams
+{
+namespace narray
+{
+
+//binary array read/write operations
+
+template<typename T> int fwrite_ndarray(FILE *fp, narray<T> *data, narray<int32_t> *shape)
+{
+    int ret = narray_success;
+    narray_size_t piprod;
+    narray_size_t I;
+    int32_t Nd;
+    narray<int32_t> *lshape=NULL;
+
+    if(fp==NULL)
+    {
+        ret = narray_failure;
+        printf("fwrite_ndarray: fp=NULL\n");
+        return ret;
+    }
+
+    if(shape==NULL)
+    {
+        lshape = new(std::nothrow) narray<int32_t>();
+        lshape->resize(1);
+        lshape->data[0] = data->length;
+    }
+    else
+    {
+        lshape = shape;
+    }
+
+    piprod = 1;
+    for(I=0;I<lshape->length;I++)
+    {
+        if(lshape->data[I]>0)
+        {
+            piprod = piprod*lshape->data[I];
+        }
+        else
+        {
+            piprod = 0; break;
+        }
+    }
+
+    if(piprod!=data->length)
+    {
+        ret = narray_failure;
+        printf("fwrite_ndarray: buffer is size %ld, while shape is size %ld\n",(long) data->length,(long) piprod);
+        return ret;
+    }
+
+    Nd = lshape->length;
+
+    fwrite(&Nd,sizeof(int32_t),1,fp);
+    if(Nd>0)
+    {
+        fwrite(lshape->data,sizeof(int32_t),Nd,fp);
+        if(piprod>0)
+        {
+            fwrite(data->data,sizeof(T),data->length,fp);
+        }
+    }
+
+    if(shape==NULL)
+    {
+        if(lshape!=NULL) {delete lshape; lshape=NULL;}
+    }
+    
+
+    return ret;
+}
+
+template<typename T> int fread_ndarray(FILE *fp, narray<T> *data, narray<int32_t> *shape)
+{
+    int ret = narray_success;
+    narray_size_t piprod;
+    narray_size_t I;
+    int32_t Nd;
+    narray<int32_t> *lshape=NULL;
+
+    if(fp==NULL)
+    {
+        ret = narray_failure;
+        printf("fread_ndarray: fp=NULL\n");
+        return ret;
+    }
+
+    if(shape==NULL)
+    {
+        lshape = new(std::nothrow) narray<int32_t>();
+    }
+    else
+    {
+        lshape = shape;
+    }
+
+
+
+    fread(&Nd,sizeof(int32_t),1,fp);
+    if(Nd<=0)
+    {
+        data->resize(0);
+        lshape->resize(0);
+        ret = narray_failure;
+        return ret;
+    }
+    lshape->resize(Nd);
+    fread(lshape->data,sizeof(int32_t),Nd,fp);
+
+    piprod = 1;
+    for(I=0;I<lshape->length;I++)
+    {
+        if(lshape->data[I]>0)
+        {
+            piprod = piprod*lshape->data[I];
+        }
+        else
+        {
+            piprod = 0; break;
+        }
+    }
+
+    if(piprod<=0)
+    {
+        data->resize(0);
+        lshape->resize(0);
+        ret = narray_failure;
+        return ret;
+    }
+
+    data->resize(piprod);
+    fread(data->data,sizeof(T),data->length,fp);
+
+    if(shape==NULL)
+    {
+        if(lshape!=NULL) {delete lshape; lshape=NULL;}
+    }
+    
+    return ret;
+}
+
+
+};
+};
+
+#endif

include/amscppnarray/amscppnarray_impl.hpp

@@ -21,6 +21,7 @@ template<typename T> int narray<T>::resize(const narray_size_t newlength)
     int ret = narray_failure;
     narray_size_t I;
     T *newdata = NULL;
+    narray_size_t ml;
 
     if(newlength<=0)
     {
@@ -37,14 +38,19 @@ template<typename T> int narray<T>::resize(const narray_size_t newlength)
         return ret;
     }
 
-    for(I=0;I<length&I<newlength;I++)
-    {
-        newdata[I] = data[I];
-    }
-    for(I=length;I<newlength;I++)
-    {
-        newdata[I] = T();
-    }
+    //threaded initialization for speed
+    ml = (length<newlength) ? length : newlength;
+    buffer_copycast<T,T>(newdata,data,ml);
+    buffer_set<T>(newdata,T(),length,newlength);
+
+    // for(I=0;I<length&I<newlength;I++)
+    // {
+    //     newdata[I] = data[I];
+    // }
+    // for(I=length;I<newlength;I++)
+    // {
+    //     newdata[I] = T();
+    // }
 
     if(data!=NULL) {delete[] data; data=NULL;}
     data = newdata;
@@ -60,6 +66,7 @@ template<typename T> narray<T>::narray(const narray<T>& other)
     int res;
     length = 0;
     data = NULL;
+    narray_size_t ml;
 
     //printf("debug: copy constructor called for\n\t%p\n\t%p\n", this, &other);
 
@@ -72,10 +79,13 @@ template<typename T> narray<T>::narray(const narray<T>& other)
         }
         else
         {
-            for(I=0;I<length && I<other.length;I++)
-            {
-                data[I] = other.data[I];
-            }
+            //threaded copy for speed
+            ml = (length<other.length) ? length : other.length;
+            buffer_copycast<T,T>(data,other.data,ml);
+            // for(I=0;I<length && I<other.length;I++)
+            // {
+            //     data[I] = other.data[I];
+            // }
         }
     }
 }
@@ -119,6 +129,7 @@ template<typename T> narray<T>& narray<T>::operator=(const narray<T>& other)
 {
     narray_size_t I;
     int res;
+    narray_size_t ml;
 
     //printf("debug: operator= called for\n\t%p\n\t%p\n", this, &other);
 
@@ -131,10 +142,13 @@ template<typename T> narray<T>& narray<T>::operator=(const narray<T>& other)
         }
         else
         {
-            for(I=0;I<length && I<other.length;I++)
-            {
-                data[I] = other.data[I];
-            }
+            //threaded copy for speed
+            ml = (length<other.length) ? length : other.length;
+            buffer_copycast<T,T>(data,other.data,ml);
+            // for(I=0;I<length && I<other.length;I++)
+            // {
+            //     data[I] = other.data[I];
+            // }
         }
     }
     return *this;

include/amscppnarray/amscppnarray_llbops.hpp

@@ -0,0 +1,23 @@
+#ifndef __AMSCPPNARRAY_LLBOPS_HPP__
+#define __AMSCPPNARRAY_LLBOPS_HPP__
+
+//Low level threaded buffer operations
+
+namespace ams
+{
+namespace narray
+{
+
+    template<typename T1, typename T2> void buffer_copycast(T1 *bufferto, const T2 *bufferfrom, const narray_size_t N);
+    
+    
+    template<typename T> void buffer_set(T *buffer, const T val, const narray_size_t ind1, const narray_size_t ind2);
+
+
+
+};
+};
+
+#include <amscppnarray/amscppnarray_llbopsimpl.hpp>
+
+#endif

include/amscppnarray/amscppnarray_llbopsimpl.hpp

@@ -0,0 +1,166 @@
+#ifndef __AMSCPPNARRAY_LLBOPSIMPL_HPP__
+#define __AMSCPPNARRAY_LLBOPSIMPL_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+template<typename T1, typename T2> void buffer_copycast_tf(
+    T1 *bufferto,
+    const T2 *bufferfrom,
+    const narray_size_t N,
+    int threadnum, int nthreads
+)
+{
+    int I0,I1,Is,I;
+    
+    Is = N/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1))? N: Is*(threadnum+1);
+    for(I=I0;I<I1;I++)
+    {
+        bufferto[I] = (T1) bufferfrom[I];
+    }
+
+    return;
+}
+
+template<typename T1, typename T2> void buffer_copycast(T1 *bufferto, const T2 *bufferfrom, const narray_size_t N)
+{
+    narray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+
+    if(N<=0) return;
+
+    if(N<narray_thread_sz)
+    {
+        for(I=0;I<N;I++)
+        {
+            bufferto[I] = (T1) bufferfrom[I];
+        }
+    }
+    else
+    {
+        //threaded operation
+        
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+        threads.resize(nthreads);
+
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &buffer_copycast_tf<T1,T2>,
+                bufferto,
+                bufferfrom,
+                N,
+                J,nthreads
+            );
+            if(threads[J]==NULL)
+            {
+                printf("buffer_copy_cast: warning, thread %d failed to launch.\n",J);
+            }
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+    return;
+}
+
+template<typename T> void buffer_set_tf(
+    T *bufferto,
+    const T val,
+    const narray_size_t ind1,
+    const narray_size_t ind2,
+    const narray_size_t N,
+    int threadnum, int nthreads
+)
+{
+    int I0,I1,Is,I;
+    
+    Is = N/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1))? N: Is*(threadnum+1);
+
+    I0 = I0 + ind1;
+    I1 = I1 + ind1;
+
+    for(I=I0;I<I1;I++)
+    {
+        bufferto[I] = val;
+    }
+
+    return;
+}
+
+template<typename T> void buffer_set(T *buffer, const T val, const narray_size_t ind1, const narray_size_t ind2)
+{
+    narray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+
+    narray_size_t N = ind2-ind1;
+    if(N<=0) return;
+
+    if(N<narray_thread_sz)
+    {
+        for(I=ind1;I<ind2;I++)
+        {
+            buffer[I] = val;
+        }
+    }
+    else
+    {
+        //threaded operation
+        
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+        threads.resize(nthreads);
+
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &buffer_set_tf<T>,
+                buffer,
+                val,
+                ind1,ind2,N,
+                J,nthreads
+            );
+            if(threads[J]==NULL)
+            {
+                printf("buffer_copy_cast: warning, thread %d failed to launch.\n",J);
+            }
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+    return;
+}
+
+
+};
+};
+
+#endif

include/amscppnarray/amscppnarray_math.hpp

@@ -19,6 +19,7 @@ template<typename T> narray<T> narray_abs(narray<T> &q);
 
 void narray_testmath1();
 void narray_testmath2();
+void narray_testmath3();
 
 }; //end namespaces
 };

include/amscppnarray/amscppnarray_mathimpl.hpp

@@ -101,8 +101,11 @@ namespace narray
     template<typename T> double narray_mean(const narray<T> &q)
     {
         double ret = 0.0;
-        T sum = narray_sum(q);
-        ret = (double)sum/((double)q.length);
+        if(q.length>0)
+        {
+            T sum = narray_sum(q);
+            ret = (double)sum/((double)q.length);
+        }
         return ret;
     }
 
@@ -197,8 +200,12 @@ namespace narray
 
     template<typename T> double narray_stdev(const narray<T> &q)
     {
-        double ret = narray_var(q);
-        ret = ::sqrt(ret);
+        double ret = 0.0;
+        if(q.length>0)
+        {
+            ret = narray_var(q);
+            ret = ::sqrt(ret);
+        }
         return ret;
     }
 

include/amscppnarray/amscppnarray_random.hpp

@@ -18,6 +18,8 @@ namespace rand
 
     static const narray_rands_t dpr32_mod = ( ((narray_rands_t)1) << ((narray_rands_t)30) ) - (narray_rands_t)1;
     static const narray_rands_t dpr32_mult1 = ( (narray_rands_t) 1201633 );
+    static const narray_rands_t dpr32_add1 = ( (narray_rands_t) 293482 );
+    
     extern narray_rands_t dpr32_rseed; // global random number seed for numeric array class
     
     narray_rands_t dpr32_nextseed(narray_rands_t seed);

scripts/amsbinarrayrw.py

@@ -0,0 +1,81 @@
+#!/usr/bin/python3
+import os,sys,math
+import numpy as np
+
+#Write double ndarray
+def fwrite_dndarray(fp,array):
+
+    array = np.array(array,dtype=np.float64,order='F')
+
+    Ndim = np.int32(len(array.shape))
+    shp = np.int32(array.shape)
+
+    q = Ndim.tobytes()
+    fp.write(q)
+    if(Ndim>0):
+        q = shp.tobytes()
+        fp.write(q)
+
+        q = array.tobytes()
+        fp.write(q)
+
+    return
+
+#Read double ndarray
+def fread_dndarray(fp):
+
+    array = np.zeros((),dtype=np.float64,order='F')
+
+    q = fp.read(4)
+    Ndim = np.frombuffer(q,dtype=np.int32,count=1)[0]
+    
+    if(Ndim>0):
+        q = fp.read(4*Ndim)
+        shp = np.frombuffer(q,dtype=np.int32,count=Ndim)
+        piprod = 1
+        for i in range(0,len(shp)):
+            piprod = piprod*shp[i]
+        q = fp.read(8*piprod)
+        array = np.frombuffer(q,dtype=np.float64,count=piprod)
+        array = array.reshape(shp,order='F')
+        
+    return array
+
+#Write float ndarray
+def fwrite_fndarray(fp,array):
+
+    array = np.array(array,dtype=np.float32,order='F')
+
+    Ndim = np.int32(len(array.shape))
+    shp = np.int32(array.shape)
+
+    q = Ndim.tobytes()
+    fp.write(q)
+    if(Ndim>0):
+        q = shp.tobytes()
+        fp.write(q)
+
+        q = array.tobytes()
+        fp.write(q)
+
+    return
+
+#Read float ndarray
+def fread_fndarray(fp):
+
+    array = np.zeros((),dtype=np.float32,order='F')
+
+    q = fp.read(4)
+    Ndim = np.frombuffer(q,dtype=np.int32,count=1)[0]
+    
+    if(Ndim>0):
+        q = fp.read(4*Ndim)
+        shp = np.frombuffer(q,dtype=np.int32,count=Ndim)
+        piprod = 1
+        for i in range(0,len(shp)):
+            piprod = piprod*shp[i]
+        q = fp.read(4*piprod)
+        array = np.frombuffer(q,dtype=np.float32,count=piprod)
+        array = array.reshape(shp,order='F')
+        
+    return array

scripts/checkrandomness.py

@@ -0,0 +1,31 @@
+#!/usr/bin/bash
+
+import os,sys,math
+import numpy as np
+import matplotlib.pyplot as plt
+
+from amsbinarrayrw import *
+
+def plot_random1():
+
+    fp = open("./build_linux64/tmp_testrand.bin","r+b")
+    data = fread_dndarray(fp)
+    N1 = data.shape[0]
+    N2 = int(math.sqrt(N1))
+    N3 = N2*N2
+    data = data[0:N3]
+    data = data.reshape([N2,N2])
+    q = np.zeros([N2,N2,3])
+    q[:,:,0] = data[:,:]
+    q[:,:,1] = data[:,:]
+    q[:,:,2] = data[:,:]
+    
+    hf = plt.figure(1)
+    plt.imshow(q)
+    plt.show()
+    fp.close()
+
+if(__name__=="__main__"):
+    plot_random1()
+
+    pass

src/amscppnarray/amscppnarray_math.cpp

@@ -57,15 +57,13 @@ namespace narray
 
     void narray_testmath2()
     {
-        rand::narray_rands_t rs = 1943981;
+        //rand::narray_rands_t rs = 1943981;
         rand::set_rseed_withtimer();
         printf("rseed: %d\n",(int)rand::get_rseed());
         printf("rseed: %d\n",(int)rand::dpr32_rseed); //< -- this is different than the line above
 
-        //rand::set_rseed(25);
-        printf("rseed: %d\n",(int)rand::get_rseed());
         
-        narray<double> unif = ams::narray::rand::narray_rand(5,&rs);
+        narray<double> unif = ams::narray::rand::narray_rand(100);
         
 
         printf("rseed: %d\n",(int)rand::get_rseed());
@@ -94,5 +92,19 @@ namespace narray
         return;
     }
 
+    void narray_testmath3()
+    {
+        narray<double> rb;
+        const char *fout = "./tmp_testrand.bin";
+        FILE *fp = NULL;
+        fp = fopen(fout,"w+b");
+
+        rb = rand::narray_rand(500*500);
+        fwrite_ndarray<double>(fp,&rb,NULL);
+
+        fclose(fp);
+
+    }
+
 };
 };

src/amscppnarray/amscppnarray_random.cpp

@@ -13,16 +13,16 @@ namespace rand
     narray_rands_t dpr32_nextseed(narray_rands_t seed)
     {
         narray_rands_t sret = seed;
-        sret = ams::mod(sret*dpr32_mult1+1,dpr32_mod);
+        sret = ams::mod(sret*dpr32_mult1+dpr32_add1,dpr32_mod);
         return sret;
     }
     
     double dpr32_randd(narray_rands_t *seed)
     {
         double ret;
-        printf("debug: seed: %d\n",(int)*seed);
+        //printf("debug: seed: %d\n",(int)*seed);
         *seed = ams::narray::rand::dpr32_nextseed(*seed);
-        printf("debug: seed: %d\n",(int)*seed);
+        //printf("debug: seed: %d\n",(int)*seed);
         ret = (double)*seed/(double)(dpr32_mod-1);
         return ret;
     }

src/main.cpp

@@ -4,10 +4,12 @@ int main(int argc, char* argv[])
 {
     int ret = 0;
     printf("ams c++ numeric array library tests.\n");
-    //ams::narray::test_narray1();
-    //ams::narray::test_narray2();
-    //ams::narray::test_narray3();
+    // ams::narray::test_narray1();
+    // ams::narray::test_narray2();
+    // ams::narray::test_narray3();
+    // ams::narray::narray_testmath1();
     ams::narray::narray_testmath2();
+    ams::narray::narray_testmath3();
 
     return ret;
 }