

                Submitted for perusal, comment, improvements, and/or critique. The presentation is in 3 sections:  motivation, code, and comment.

 

Motivation:

 

            As a new-comer to R from matrix oriented Gauss and Mata, I miss the tools for using a vector (and operator) to ‘sweep’ across a matrix.

 

            Here is how these work. If M is I rows by J columns, then one entry corresponding to each row of M is provided any 1 by I  (row) vector, R. This can be swept ‘down’ the columns with a result that is I rows, J columns with entries {M[i,j] op R[i]) } where op can be any atomic binary operator including comparisons, such as  !=, or arithmetic operators, such as  ^.

 

            A J by 1 column vector can be applied analogously to ‘sweep’ through the rows.

 

            It is natural to implement these as R op M and M op C where R, M, and C conform as for matrix multiplication. 

                        

            Following Gauss, an arguably analogous outer product is returned when 2 vectors are equal length are submitted. Scalars are allowed. And so are matrices with identical numbers of  rows and numberrs of columns.

 

            Code: 

            

            dop<-function(ml,mr,op) {#  inspired by  Gauss dot op, more plainly order sensitive

                

                                if (!is.matrix(ml))

                                                print("Warning: in dop(), Left arg will be coerced to matrix")

                                if (!is.matrix(mr))

                                                print("Warning: in dop(), Right arg will be coerced to matrix")

 

                                if (is.vector(ml) & ! is.matrix(ml))ml<-matrix(ml,nrow=1,ncol=length(ml))

                                                #risky business, force conversion of Left-hand vector

                                m.L<-as.matrix(ml);m.R<-as.matrix(mr); #conserve old Left and Right objects for debug

                                I=nrow(m.L);J=ncol(m.L);K=nrow(m.R);L=ncol(m.R)

 

                                if ( I == K & J==L)  # standard element-wise conformity for binary op

                                {return(mapply(op,m.L,m.R)); }      else 

# the elses are just for show since the if’s form a partition. Unless my testing was mistaken, of course

                                 if (!( I==1 | J==1 |K==1|L==1))  # reject if no vector and not same size

                                                {print(" DOT op requires same shapes or vector  argument")

                                                return(NA); }                 else

                                 if ( (I==1 & J==1) | (K==1 & L==1)) # at least one arg is scalar, mapply works

                                                {return(mapply(op,m.L,m.R));}  else

                                 if (I==1 & L==1 & J==K)                     #this and next mimic Gauss convention(s)          

                                              {return((matrix(outer(m.R,m.L,op),nrow=J,ncol=J )));} else # r1.*r2' 

                                                                                # n.b outer returns matrix but matrix() left in for clarity

                                 if (J==1 & K==1 & I==L )

                                                {return(matrix(outer(m.L,m.R,op),nrow=I,ncol=I ));}else # r1'.*r2

                                 if (I==1 & J==K) # m.L is row vec, hence sweep columns 

{return( (matrix(mapply(op,t(m.L),m.R),nrow=nrow(m.R),ncol=ncol(m.R))));}else

                                 if (J==K & L==1) # m.R is col vec, hence sweep rows

                {   return(t(matrix(qq1<-mapply(op,t(m.L),m.R),nrow=ncol(m.L),ncol=nrow(m.L))));}else 

                                 print(c(" conformity failure with I,J,K,L= ",paste(I,J,K,L)));return(NA)

}

 

#a further step is to implement the various binary operators using the %text% convention.

# some construction tools

#set.of.ops<-c(getGroupMembers(Compare),getGroupMembers(Arith))

#names.for.ops<-as.array(c(

                "eq","gt","lt","ne","ge","le","plus","minus","x","exp","mod","dmod","div"))

 

`%eq%`<-function(x,y) dop(x,y,op= "==")   

`%gt%`<-function(x,y) dop(x,y,op=">")    

`%lt%`<-function(x,y) dop(x,y,op="<")    

`%ne%`<-function(x,y) dop(x,y,op="!=")   

`%ge%`<-function(x,y) dop(x,y,op="<=")   

`%le%`<-function(x,y) dop(x,y,op=">=")   

`%plus%`<-function(x,y) dop(x,y,op="+")  

`%minus%`<-function(x,y) dop(x,y,op="-") 

`%x%`<-function(x,y) dop(x,y,op="*")     

`%exp%`<-function(x,y) dop(x,y,op="^")   

`%mod%`<-function(x,y) dop(x,y,op="%%")  

`%dmod%`<-function(x,y) dop(x,y,op="%/%")

`%div%`<-function(x,y) dop(x,y,op="/")   

 

Comments:

 

            I have found these constructs very helpful in the past. One illustration is using logical comparison for a vector of standards applied to some set of variables.

 

            There are many different ways to try to achieve these results in R. I would love to hear if I have re-invented the wheel or over-looked some approach that is more general, more straightforward, or what have you. 



            Could anyone suggest more elegant code that exploits the obvious parallelism(s) of that last mass of  functions?

             

               The core construction is relatively strict in imposing explicit conformability requirements. This implicitly rejects the usual R approach of freely allowing multiple copies of shorter entities. While this will, in principle, return the ‘correct’ result whenever entities of the right size and shape are presented in the right order, it also seems to raise  risk of returning unintended nonsense without adequate warning when presented with entities not in the right order or of inapt sizes. 

 

And I have, on related grounds, imposed restrictions to matrix objects. 

 

I would be interested in sentiments on the pluses and minuses of such decisions.



While one could confine this to a (new?) class, it seems general enough to leave it unconfined. Would wiser folk disagree? 

 

            My present approach is to add this code to my list of files that are source() –d at startup. That bulk is growing. Any thoughts on how to manage this?

            Perhaps for some this would seem to be a wasteful proliferation of operators. Is there perhaps some convenient mechanism so that this collection could be switched on and then off as needed? (I bet my question misses some obvious point which is the point of asking it.)

 

            Those illustrate my questions but I presume others may see other avenues worth pursuing. 

 

            Thanks in advance for any comments and/or reactions.

Steven Rytina
Department of Sociology
McGill University
	[[alternative HTML version deleted]]


