View Javadoc
1   /*
2    * The baseCode project
3    * 
4    * Copyright (c) 2008-2019 University of British Columbia
5    * 
6    * Licensed under the Apache License, Version 2.0 (the "License");
7    * you may not use this file except in compliance with the License.
8    * You may obtain a copy of the License at
9    *
10   *       http://www.apache.org/licenses/LICENSE-2.0
11   *
12   * Unless required by applicable law or agreed to in writing, software
13   * distributed under the License is distributed on an "AS IS" BASIS,
14   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15   * See the License for the specific language governing permissions and
16   * limitations under the License.
17   *
18   */
19  package ubic.basecode.util.r;
20  
21  import java.io.BufferedReader;
22  import java.io.IOException;
23  import java.io.InputStream;
24  import java.io.InputStreamReader;
25  import java.util.ArrayList;
26  import java.util.Arrays;
27  import java.util.HashMap;
28  import java.util.Iterator;
29  import java.util.List;
30  import java.util.Map;
31  
32  import org.apache.commons.collections4.Transformer;
33  import org.apache.commons.collections4.functors.StringValueTransformer;
34  import org.apache.commons.lang3.ArrayUtils;
35  import org.apache.commons.lang3.RandomStringUtils;
36  import org.apache.commons.lang3.StringUtils;
37  import org.rosuda.REngine.REXP;
38  import org.rosuda.REngine.REXPDouble;
39  import org.rosuda.REngine.REXPFactor;
40  import org.rosuda.REngine.REXPGenericVector;
41  import org.rosuda.REngine.REXPInteger;
42  import org.rosuda.REngine.REXPList;
43  import org.rosuda.REngine.REXPLogical;
44  import org.rosuda.REngine.REXPMismatchException;
45  import org.rosuda.REngine.REXPString;
46  import org.rosuda.REngine.RList;
47  import org.slf4j.Logger;
48  import org.slf4j.LoggerFactory;
49  
50  import ubic.basecode.dataStructure.matrix.DoubleMatrix;
51  import ubic.basecode.dataStructure.matrix.ObjectMatrix;
52  import ubic.basecode.dataStructure.matrix.ObjectMatrixImpl;
53  import ubic.basecode.math.linearmodels.LinearModelSummary;
54  import ubic.basecode.util.Configuration;
55  import ubic.basecode.util.r.type.HTest;
56  import ubic.basecode.util.r.type.OneWayAnovaResult;
57  import ubic.basecode.util.r.type.TwoWayAnovaResult;
58  
59  /**
60   * Base class for RClients
61   * 
62   * @author Paul
63   * 
64   */
65  public abstract class AbstractRClient implements RClient {
66  
67      protected static Logger log = LoggerFactory.getLogger( AbstractRClient.class );
68  
69      /**
70       * @param ob
71       * @return
72       */
73      public static String variableIdentityNumber( Object ob ) {
74          return Integer.toString( Math.abs( ob.hashCode() + 1 ) ) + RandomStringUtils.randomAlphabetic( 6 );
75      }
76  
77      /**
78       * Copy a matrix into an array, so that rows are represented consecutively in the array. (RServe has no interface
79       * for passing a 2-d array).
80       * 
81       * @param matrix
82       * @return
83       */
84      private static double[] unrollMatrix( double[][] matrix ) {
85          int rows = matrix.length;
86          int cols = matrix[0].length;
87          double[] unrolledMatrix = new double[rows * cols];
88  
89          int k = 0;
90          for ( int i = 0; i < rows; i++ ) {
91              for ( int j = 0; j < cols; j++ ) {
92                  unrolledMatrix[k] = matrix[i][j];
93                  k++;
94              }
95          }
96          return unrolledMatrix;
97      }
98  
99      /**
100      * Copy a matrix into an array, so that rows are represented consecutively in the array. (RServe has no interface
101      * for passing a 2-d array).
102      * 
103      * @param matrix
104      * @return array representation of the matrix.
105      */
106     private static double[] unrollMatrix( DoubleMatrix<?, ?> matrix ) {
107         // unroll the matrix into an array Unfortunately this makes a
108         // copy of the data...and R will probably make yet
109         // another copy. If there was a way to get the raw element array from the DoubleMatrixNamed, that would
110         // be better.
111         int rows = matrix.rows();
112         int cols = matrix.columns();
113         double[] unrolledMatrix = new double[rows * cols];
114 
115         int k = 0;
116         for ( int i = 0; i < rows; i++ ) {
117             for ( int j = 0; j < cols; j++ ) {
118                 unrolledMatrix[k] = matrix.get( i, j );
119                 k++;
120             }
121         }
122         return unrolledMatrix;
123     }
124 
125     /*
126      * (non-Javadoc)
127      * 
128      * @see ubic.basecode.util.RClient#assignFactor(java.util.List)
129      */
130     @Override
131     public String assignFactor( List<String> strings ) {
132         String variableName = "f." + variableIdentityNumber( strings );
133         return assignFactor( variableName, strings );
134     }
135 
136     /*
137      * (non-Javadoc)
138      * 
139      * @see ubic.basecode.util.RClient#assignFactor(java.lang.String, java.util.List)
140      */
141     @Override
142     public String assignFactor( String factorName, List<String> list ) {
143         String l = assignStringList( list );
144         this.voidEval( factorName + "<-factor(" + l + ")" );
145         return factorName;
146     }
147 
148     /*
149      * (non-Javadoc)
150      * 
151      * @see ubic.basecode.util.RClient#assignMatrix(double[][])
152      */
153     @Override
154     public String assignMatrix( double[][] matrix ) {
155         String matrixVarName = "Matrix_" + variableIdentityNumber( matrix );
156         log.debug( "Assigning matrix with variable name " + matrixVarName );
157         int rows = matrix.length;
158         int cols = matrix[0].length;
159         if ( rows == 0 || cols == 0 ) throw new IllegalArgumentException( "Empty matrix?" );
160         double[] unrolledMatrix = unrollMatrix( matrix );
161         this.assign( "U" + matrixVarName, unrolledMatrix ); // temporary
162         this.voidEval( matrixVarName + "<-matrix(" + "U" + matrixVarName + ", nrow=" + rows + " , ncol=" + cols
163                 + ", byrow=TRUE)" );
164         this.voidEval( "rm(U" + matrixVarName + ")" ); // maybe this saves memory...
165 
166         return matrixVarName;
167     }
168 
169     /*
170      * (non-Javadoc)
171      * 
172      * @see ubic.basecode.util.RClient#assignMatrix(ubic.basecode.dataStructure.matrix.DoubleMatrixNamed)
173      */
174     @Override
175     public String assignMatrix( DoubleMatrix<?, ?> matrix ) {
176         return assignMatrix( matrix, StringValueTransformer.stringValueTransformer() );
177     }
178 
179     /*
180      * (non-Javadoc)
181      * 
182      * @see ubic.basecode.util.RClient#assignMatrix(ubic.basecode.dataStructure.matrix.DoubleMatrix,
183      * org.apache.commons.collections.Transformer)
184      */
185     @Override
186     public String assignMatrix( DoubleMatrix<?, ?> matrix, Transformer rowNameExtractor ) {
187         String matrixVarName = "Matrix_" + variableIdentityNumber( matrix );
188         log.debug( "Assigning matrix with variable name " + matrixVarName );
189         int rows = matrix.rows();
190         int cols = matrix.columns();
191         if ( rows == 0 || cols == 0 ) throw new IllegalArgumentException( "Empty matrix?" );
192         double[] unrolledMatrix = unrollMatrix( matrix );
193         assert ( unrolledMatrix != null );
194         this.assign( "U" + matrixVarName, unrolledMatrix );
195         this.voidEval( matrixVarName + "<-matrix(" + "U" + matrixVarName + ", nrow=" + rows + ", ncol=" + cols
196                 + ", byrow=TRUE)" );
197         this.voidEval( "rm(U" + matrixVarName + ")" ); // maybe this saves memory...
198 
199         if ( matrix.hasColNames() && matrix.hasRowNames() )
200             assignRowAndColumnNames( matrix, matrixVarName, rowNameExtractor );
201         return matrixVarName;
202     }
203 
204     /*
205      * (non-Javadoc)
206      * 
207      * @see ubic.basecode.util.RClient#assignStringList(java.util.List)
208      */
209     @Override
210     public String assignStringList( List<?> strings ) {
211         String variableName = "stringList." + variableIdentityNumber( strings );
212 
213         Object[] array = strings.toArray();
214         String[] sa = new String[array.length];
215         for ( int i = 0; i < array.length; i++ ) {
216             sa[i] = array[i].toString();
217         }
218 
219         assign( variableName, sa );
220         return variableName;
221     }
222 
223     /*
224      * (non-Javadoc)
225      * 
226      * @see ubic.basecode.util.RClient#booleanDoubleArrayEval(java.lang.String, java.lang.String, double[])
227      */
228     @Override
229     public boolean booleanDoubleArrayEval( String command, String argName, double[] arg ) {
230         this.assign( argName, arg );
231         REXP x = this.eval( command );
232         if ( x.isLogical() ) {
233             try {
234                 REXPLogical b = new REXPLogical( new boolean[1], new REXPList( x.asList() ) );
235                 return b.isTRUE()[0];
236             } catch ( REXPMismatchException e ) {
237                 throw new RuntimeException( e );
238             }
239         }
240         return false;
241     }
242 
243     /*
244      * (non-Javadoc)
245      * 
246      * @see ubic.basecode.util.RClient#dataFrame(ubic.basecode.dataStructure.matrix.ObjectMatrix)
247      */
248     @Override
249     public String dataFrame( ObjectMatrix<String, String, Object> matrix ) {
250 
251         /*
252          * Extract columns, convert
253          */
254 
255         List<String> colNames = matrix.getColNames();
256         List<String> rowNames = matrix.getRowNames();
257 
258         assert colNames.size() == matrix.columns();
259 
260         String colV = assignStringList( colNames );
261         String rowV = assignStringList( rowNames );
262 
263         List<String> terms = new ArrayList<String>();
264         for ( int j = 0; j < colNames.size(); j++ ) {
265 
266             Object[] column;
267 
268             Object v = matrix.getEntry( 0, j );
269 
270             if ( v instanceof Number ) {
271                 column = new Double[matrix.rows()];
272             } else if ( v instanceof Boolean ) {
273                 column = new String[matrix.rows()];
274             } else if ( v instanceof String ) {
275                 column = new String[matrix.rows()];
276             } else {
277                 throw new IllegalArgumentException( "Sorry, can't make data frame from values of type "
278                         + v.getClass().getName() );
279             }
280 
281             for ( int i = 0; i < matrix.rows(); i++ ) {
282 
283                 Object value = matrix.getEntry( i, j );
284 
285                 if ( matrix.isMissing( i, j ) ) {
286                     column[i] = null;
287                 } else if ( value instanceof Number ) {
288                     column[i] = ( ( Number ) value ).doubleValue();
289                 } else if ( value instanceof Boolean ) {
290                     column[i] = ( ( Boolean ) value ) ? "T" : "F";
291                 } else if ( value instanceof String ) {
292                     column[i] = value;
293                 }
294             }
295 
296             if ( v instanceof Number ) {
297                 assign( colNames.get( j ), ArrayUtils.toPrimitive( ( Double[] ) column ) );
298                 terms.add( colNames.get( j ) );
299             } else if ( v instanceof Boolean ) {
300                 assignFactor( colNames.get( j ), Arrays.asList( ( String[] ) column ) );
301                 terms.add( colNames.get( j ) );
302             } else if ( v instanceof String ) {
303                 assignFactor( colNames.get( j ), Arrays.asList( ( String[] ) column ) );
304                 terms.add( colNames.get( j ) );
305             }
306 
307         }
308 
309         String varName = "df." + variableIdentityNumber( matrix );
310 
311         String command = varName + "<-data.frame(" + StringUtils.join( terms, "," ) + ", row.names=" + rowV + " )";
312         log.debug( command );
313         eval( command );
314         eval( "names(" + varName + ")<-" + colV );
315 
316         return varName;
317 
318     }
319 
320     @Override
321     public ObjectMatrix<String, String, Object> dataFrameEval( String command ) {
322 
323         REXP df = eval( command );
324 
325         try {
326 
327             RList dfl = df.asList();
328 
329             if ( !df.getAttribute( "class" ).asString().equals( "data.frame" ) ) {
330                 throw new IllegalArgumentException( "Command did not return a dataframe" );
331             }
332 
333             String[] rowNames = df.getAttribute( "row.names" ).asStrings();
334             String[] colNames = df.getAttribute( "names" ).asStrings();
335 
336             assert dfl.size() == colNames.length;
337 
338             ObjectMatrix<String, String, Object> result = new ObjectMatrixImpl<String, String, Object>(
339                     rowNames.length, colNames.length );
340 
341             result.setRowNames( Arrays.asList( rowNames ) );
342             result.setColumnNames( Arrays.asList( colNames ) );
343 
344             for ( int i = 0; i < dfl.size(); i++ ) {
345                 REXP column = ( REXP ) dfl.get( i );
346 
347                 if ( column.isNumeric() ) {
348                     double[] asDoubles = column.asDoubles();
349 
350                     for ( int j = 0; j < rowNames.length; j++ ) {
351                         result.set( j, i, asDoubles[j] );
352                     }
353 
354                 } else {
355                     String[] asStrings = column.asStrings();
356 
357                     for ( int j = 0; j < rowNames.length; j++ ) {
358                         result.set( j, i, asStrings[j] );
359                     }
360 
361                 }
362 
363             }
364 
365             return result;
366         } catch ( REXPMismatchException e ) {
367 
368             throw new RuntimeException( e );
369 
370         }
371 
372     }
373 
374     public abstract void disconnect();
375 
376     /*
377      * (non-Javadoc)
378      * 
379      * @see ubic.basecode.util.RClient#doubleArrayDoubleArrayEval(java.lang.String, java.lang.String, double[])
380      */
381     @Override
382     public double[] doubleArrayDoubleArrayEval( String command, String argName, double[] arg ) {
383         try {
384             this.assign( argName, arg );
385             RList l = this.eval( command ).asList();
386             return l.at( argName ).asDoubles();
387         } catch ( REXPMismatchException e ) {
388             throw new RuntimeException( e );
389         }
390     }
391 
392     /*
393      * (non-Javadoc)
394      * 
395      * @see ubic.basecode.util.RClient#doubleArrayEval(java.lang.String)
396      */
397     @Override
398     public double[] doubleArrayEval( String command ) {
399         REXP r = this.eval( command );
400         if ( r == null ) {
401             return null;
402         }
403 
404         if ( !r.isNumeric() ) {
405             throw new RuntimeException( "Command did not return numbers: " + command + ", result was: " + r );
406         }
407 
408         try {
409             return r.asDoubles();
410         } catch ( REXPMismatchException e ) {
411             throw new RuntimeException( e );
412         }
413     }
414 
415     /*
416      * (non-Javadoc)
417      * 
418      * @see ubic.basecode.util.RClient#doubleArrayTwoDoubleArrayEval(java.lang.String, java.lang.String, double[],
419      * java.lang.String, double[])
420      */
421     @Override
422     public double[] doubleArrayTwoDoubleArrayEval( String command, String argName, double[] arg, String argName2,
423             double[] arg2 ) {
424         this.assign( argName, arg );
425         this.assign( argName2, arg2 );
426         try {
427             return this.eval( command ).asDoubles();
428         } catch ( REXPMismatchException e ) {
429             throw new RuntimeException( e );
430         }
431     }
432 
433     /*
434      * (non-Javadoc)
435      * 
436      * @see ubic.basecode.util.RClient#doubleTwoDoubleArrayEval(java.lang.String, java.lang.String, double[],
437      * java.lang.String, double[])
438      */
439     @Override
440     public double doubleTwoDoubleArrayEval( String command, String argName, double[] arg, String argName2, double[] arg2 ) {
441         this.assign( argName, arg );
442         this.assign( argName2, arg2 );
443         REXP x = this.eval( command );
444         try {
445             return x.asDouble();
446         } catch ( REXPMismatchException e ) {
447             throw new RuntimeException( e );
448         }
449     }
450 
451     /*
452      * (non-Javadoc)
453      * 
454      * @see ubic.basecode.util.RClient#intArrayEval(java.lang.String)
455      */
456     @Override
457     public int[] intArrayEval( String command ) {
458         try {
459             return eval( command ).asIntegers();
460         } catch ( REXPMismatchException e ) {
461             throw new RuntimeException( e );
462         }
463     }
464 
465     /*
466      * (non-Javadoc)
467      * 
468      * @see ubic.basecode.util.RClient#linearModel(double[], java.util.List)
469      */
470     @Override
471     @SuppressWarnings({ "unchecked" })
472     public LinearModelSummary linearModel( double[] data, Map<String, List<?>> factors ) {
473 
474         String datName = RandomStringUtils.randomAlphabetic( 10 );
475         assign( datName, data );
476 
477         for ( String factorName : factors.keySet() ) {
478             List<?> list = factors.get( factorName );
479             if ( list.iterator().next() instanceof String ) {
480                 assignFactor( factorName, ( List<String> ) list );
481             } else {
482                 // treat is as a numeric covariate
483                 List<Double> d = new ArrayList<Double>();
484                 for ( Object object : list ) {
485                     d.add( ( Double ) object );
486                 }
487 
488                 assign( factorName, ArrayUtils.toPrimitive( d.toArray( new Double[] {} ) ) );
489             }
490         }
491 
492         String modelDeclaration = datName + " ~ " + StringUtils.join( factors.keySet(), "+" );
493 
494         String lmName = RandomStringUtils.randomAlphabetic( 10 );
495         String command = lmName + "<-lm(" + modelDeclaration + ", na.action=na.exclude)";
496         log.debug( command );
497         voidEval( command );
498 
499         REXP lmsum = eval( "summary(" + lmName + ")" );
500         REXP anova = eval( "anova(" + lmName + ")" );
501 
502         return new LinearModelSummary( lmsum, anova, factors.keySet().toArray( new String[] {} ) );
503 
504     }
505 
506     /*
507      * (non-Javadoc)
508      * 
509      * @see ubic.basecode.util.RClient#linearModel(double[], ubic.basecode.dataStructure.matrix.ObjectMatrix)
510      */
511     @Override
512     public LinearModelSummary linearModel( double[] data, ObjectMatrix<String, String, Object> d ) {
513 
514         String datName = RandomStringUtils.randomAlphabetic( 10 );
515         assign( datName, data );
516 
517         String df = dataFrame( d );
518 
519         String varNames = StringUtils.join( d.getColNames(), "+" );
520 
521         String lmName = RandomStringUtils.randomAlphabetic( 10 );
522         String command = lmName + "<-lm(" + datName + " ~ " + varNames + ", data=" + df + ", na.action=na.exclude)";
523         voidEval( command );
524 
525         REXP lmsum = eval( "summary(" + lmName + ")" );
526         REXP anova = eval( "anova(" + lmName + ")" );
527 
528         return new LinearModelSummary( lmsum, anova, d.getColNames().toArray( new String[] {} ) );
529 
530     }
531 
532     /**
533      * FIXME only partly implemented, possibly not going to stay.
534      */
535     @Override
536     public List<?> listEval( Class<?> listEntryType, String command ) {
537 
538         REXP rexp = this.eval( command );
539 
540         List<Object> result = new ArrayList<Object>();
541         try {
542             if ( !rexp.isVector() ) {
543                 throw new IllegalArgumentException( "Command did not return some kind of vector" );
544             }
545 
546             if ( rexp instanceof REXPInteger ) {
547                 log.debug( "integer" );
548                 double[][] asDoubleMatrix = rexp.asDoubleMatrix();
549                 for ( double[] ds : asDoubleMatrix ) {
550                     result.add( ds );
551                 }
552 
553                 if ( rexp instanceof REXPFactor ) {
554                     log.info( "factor" );
555                     // not sure what to do...
556                 }
557             } else if ( rexp instanceof REXPGenericVector ) {
558                 log.debug( "generic" );
559                 REXPGenericVector v = ( REXPGenericVector ) rexp;
560                 List<?> tmp = new ArrayList<Object>( v.asList().values() );
561 
562                 if ( tmp.size() == 0 ) return tmp;
563 
564                 for ( Object t : tmp ) {
565                     String clazz = ( ( REXP ) t ).getAttribute( "class" ).asString();
566                     /*
567                      * FIXME!!!!
568                      */
569                     if ( clazz.equals( "htest" ) ) {
570                         try {
571                             result.add( new HTest( ( ( REXP ) t ).asList() ) );
572                         } catch ( REXPMismatchException e ) {
573                             result.add( new HTest() );
574                         }
575                     } else if ( clazz.equals( "lm" ) ) {
576                         throw new UnsupportedOperationException();
577                     } else {
578                         result.add( new HTest() ); // e.g. failed result or something we don't know about yet
579                     }
580                     /*
581                      * todo: support lm objects, anova objects others? pair.htest?
582                      */
583                 }
584 
585             } else if ( rexp instanceof REXPDouble ) {
586                 log.debug( "double" );
587                 double[][] asDoubleMatrix = rexp.asDoubleMatrix();
588                 for ( double[] ds : asDoubleMatrix ) {
589                     result.add( ds );
590                 }
591 
592             } else if ( rexp instanceof REXPList ) {
593                 log.debug( "list" );
594                 if ( rexp.isPairList() ) {
595                     // log.info( "pairlist" ); always true for REXPList.
596                 }
597 
598                 if ( rexp.isLanguage() ) {
599                     throw new UnsupportedOperationException( "Don't know how to deal with vector type of "
600                             + rexp.getClass().getName() );
601                 }
602                 log.debug( rexp.getClass().getName() );
603                 result = new ArrayList<Object>( rexp.asList().values() );
604             } else {
605                 throw new UnsupportedOperationException( "Don't know how to deal with vector type of "
606                         + rexp.getClass().getName() );
607             }
608 
609             return result;
610         } catch ( REXPMismatchException e ) {
611             throw new RuntimeException( e );
612         }
613 
614     }
615 
616     /*
617      * (non-Javadoc)
618      * 
619      * @see ubic.basecode.util.RClient#loadLibrary(java.lang.String)
620      */
621     @Override
622     public boolean loadLibrary( String libraryName ) {
623         try {
624 
625             String userLibPath = Configuration.getString( "rlibpath" );
626             if ( StringUtils.isNotBlank( userLibPath ) ) {
627                 voidEval( ".libPaths(" + userLibPath + ")" );
628             }
629             // List<String> libPaths = stringListEval( ".libPaths()" );
630 
631             List<String> libraries = stringListEval( "installed.packages()[,1]" );
632             if ( !libraries.contains( libraryName ) ) {
633                 return false;
634             }
635         } catch ( Exception e ) {
636             /*
637              * This can happen with a 'Error in gzfile()', no idea why and it seems harmless if uninformative.
638              */
639             log.warn( "Was unable to check for installed library before attempting to load it." );
640         }
641 
642         try {
643             voidEval( "library(" + libraryName + ")" );
644             return true;
645         } catch ( Exception e ) {
646             return false;
647         }
648     }
649 
650     /*
651      * (non-Javadoc)
652      * 
653      * @see ubic.basecode.util.RClient#oneWayAnova(double[], java.util.List)
654      */
655     @Override
656     public OneWayAnovaResult oneWayAnova( double[] data, List<String> factor ) {
657         String f = assignFactor( factor );
658         StringBuffer command = new StringBuffer();
659 
660         assign( "foo", data );
661 
662         String modelDeclaration;
663 
664         modelDeclaration = "foo  ~ " + f;
665 
666         command.append( "anova(aov(" + modelDeclaration + "))" );
667 
668         REXP eval = eval( command.toString() );
669 
670         return new OneWayAnovaResult( eval );
671     }
672 
673     /*
674      * (non-Javadoc)
675      * 
676      * @see ubic.basecode.util.RClient#oneWayAnovaEval(java.lang.String)
677      */
678     @Override
679     public Map<String, OneWayAnovaResult> oneWayAnovaEval( String command ) {
680         REXP rawResult = this.eval( command );
681 
682         if ( rawResult == null ) {
683             return null;
684         }
685 
686         Map<String, OneWayAnovaResult> result = new HashMap<String, OneWayAnovaResult>();
687         try {
688             RList mainList = rawResult.asList();
689             if ( mainList == null ) {
690                 return null;
691             }
692 
693             log.debug( mainList.size() + " results." );
694 
695             for ( int i = 0; i < mainList.size(); i++ ) {
696 
697                 REXP anovaTable = mainList.at( i );
698 
699                 String elementIdentifier = mainList.keyAt( i );
700 
701                 assert elementIdentifier != null;
702 
703                 if ( log.isDebugEnabled() ) log.debug( "Key: " + elementIdentifier );
704 
705                 if ( !anovaTable.isList() || !anovaTable.hasAttribute( "row.names" ) ) {
706                     log.debug( "No anovaresult for " + elementIdentifier );
707                     result.put( elementIdentifier, new OneWayAnovaResult() );
708                     continue;
709                 }
710 
711                 result.put( elementIdentifier, new OneWayAnovaResult( anovaTable ) );
712 
713             }
714         } catch ( REXPMismatchException e ) {
715             throw new RuntimeException( e );
716         }
717 
718         return result;
719 
720     }
721 
722     /*
723      * (non-Javadoc)
724      * 
725      * @see ubic.basecode.util.RClient#remove(java.lang.String)
726      */
727     @Override
728     public void remove( String variableName ) {
729         this.voidEval( "rm(" + variableName + ")" );
730     }
731 
732     /*
733      * (non-Javadoc)
734      * 
735      * @see ubic.basecode.util.RClient#rowApplyLinearModel(java.lang.String, java.lang.String, java.lang.String[])
736      */
737     @Override
738     public Map<String, LinearModelSummary> rowApplyLinearModel( String dataMatrixVarName, String modelFormula,
739             String[] factorNames ) {
740 
741         String lmres = "lmlist." + RandomStringUtils.randomAlphanumeric( 10 );
742 
743         log.info( "Starting model fitting ..." );
744 
745         loadScript( this.getClass().getResourceAsStream( "/ubic/basecode/util/r/linearModels.R" ) );
746 
747         String command = lmres + "<-rowlm(" + modelFormula + ", data.frame(" + dataMatrixVarName + ") )";
748 
749         log.debug( command );
750         this.voidEval( command );
751 
752         log.info( "Model fits complete, summarizing ..." );
753         REXP rawLmSummaries = this.eval( "lapply(" + lmres + ", function(x){ try(summary(x), silent=T)})" );
754 
755         if ( rawLmSummaries == null ) {
756             log.warn( "No results from apply(... lm)" );
757             return null;
758         }
759 
760         log.info( "Summaries done, doing hypothesis tests ..." );
761         /* this tends to be the slow part. */
762         REXP rawAnova = this.eval( "lapply(" + lmres + ", function(x){ try(anova(x), silent=T)})" );
763 
764         Map<String, LinearModelSummary> result = new HashMap<String, LinearModelSummary>();
765         try {
766             log.info( "Processing the results ..." );
767 
768             RList rawLmList = rawLmSummaries.asList();
769             if ( rawLmList == null ) {
770                 log.warn( "Raw lm summary results were null" );
771                 return null;
772             }
773 
774             RList rawAnovaList = rawAnova.asList();
775             if ( rawAnovaList == null ) {
776                 return null;
777             }
778 
779             log.debug( rawLmList.size() + " results." );
780 
781             assert rawLmList.size() == rawAnovaList.size();
782 
783             for ( int i = 0; i < rawLmList.size(); i++ ) {
784 
785                 REXP lmSummary = rawLmList.at( i );
786                 REXP anova = rawAnovaList.at( i );
787 
788                 String elementIdentifier = rawLmList.keyAt( i );
789 
790                 assert elementIdentifier != null;
791                 assert elementIdentifier.equals( rawAnovaList.keyAt( i ) );
792 
793                 if ( log.isDebugEnabled() ) log.debug( "Key: " + elementIdentifier );
794 
795                 if ( !lmSummary.isList() || !lmSummary.getAttribute( "class" ).asString().equals( "summary.lm" ) ) {
796                     log.debug( "No lm for " + elementIdentifier );
797                     result.put( elementIdentifier, new LinearModelSummary( elementIdentifier ) );
798                 } else {
799                     LinearModelSummarytml#LinearModelSummary">LinearModelSummary linearModelSummary = new LinearModelSummary( lmSummary, anova, factorNames );
800                     result.put( elementIdentifier, linearModelSummary );
801                 }
802 
803             }
804         } catch ( REXPMismatchException e ) {
805             throw new RuntimeException( e );
806         }
807 
808         return result;
809     }
810 
811     /*
812      * (non-Javadoc)
813      * 
814      * @see ubic.basecode.util.RClient#stringEval(java.lang.String)
815      */
816     @Override
817     public String stringEval( String command ) {
818         try {
819             return this.eval( command ).asString();
820         } catch ( REXPMismatchException e ) {
821             throw new RuntimeException( e );
822         }
823     }
824 
825     /*
826      * (non-Javadoc)
827      * 
828      * @see ubic.basecode.util.r.RClient#stringListEval(java.lang.String)
829      */
830     @Override
831     public List<String> stringListEval( String command ) {
832         try {
833             REXP eval = this.eval( command );
834 
835             RList v;
836             List<String> results = new ArrayList<String>();
837             if ( eval instanceof REXPString ) {
838                 String[] strs = ( ( REXPString ) eval ).asStrings();
839                 for ( String string : strs ) {
840                     results.add( string );
841                 }
842             } else {
843                 v = eval.asList();
844                 for ( Iterator<?> it = v.iterator(); it.hasNext(); ) {
845                     results.add( ( ( REXPString ) it.next() ).asString() );
846                 }
847             }
848 
849             return results;
850         } catch ( REXPMismatchException e ) {
851             throw new RuntimeException( e );
852         }
853     }
854 
855     /*
856      * (non-Javadoc)
857      * 
858      * @see ubic.basecode.util.r.RClient#twoWayAnova(double[], java.util.List, java.util.List, boolean)
859      */
860     @Override
861     public TwoWayAnovaResult twoWayAnova( double[] data, List<String> factor1, List<String> factor2,
862             boolean includeInteraction ) {
863 
864         String factorA = assignFactor( factor1 );
865         String factorB = assignFactor( factor2 );
866         StringBuffer command = new StringBuffer();
867 
868         assign( "foo", data );
869 
870         String modelDeclaration;
871 
872         if ( includeInteraction ) {
873             modelDeclaration = "foo  ~ " + factorA + "*" + factorB;
874         } else {
875             modelDeclaration = "foo  ~ " + factorA + "+" + factorB;
876         }
877 
878         command.append( "anova(aov(" + modelDeclaration + "))" );
879 
880         REXP eval = eval( command.toString() );
881 
882         return new TwoWayAnovaResult( eval );
883     }
884 
885     /*
886      * (non-Javadoc)
887      * 
888      * @see ubic.basecode.util.RClient#twoWayAnovaEval(java.lang.String, boolean)
889      */
890     @Override
891     public Map<String, TwoWayAnovaResult> twoWayAnovaEval( String command, boolean withInteractions ) {
892         REXP rawResult = this.eval( command );
893 
894         if ( rawResult == null ) {
895             return null;
896         }
897 
898         Map<String, TwoWayAnovaResult> result = new HashMap<String, TwoWayAnovaResult>();
899         try {
900             RList mainList = rawResult.asList();
901             if ( mainList == null ) {
902                 return null;
903             }
904 
905             log.debug( mainList.size() + " results." );
906 
907             for ( int i = 0; i < mainList.size(); i++ ) {
908 
909                 REXP anovaTable = mainList.at( i );
910 
911                 String elementIdentifier = mainList.keyAt( i );
912 
913                 assert elementIdentifier != null;
914 
915                 if ( log.isDebugEnabled() ) log.debug( "Key: " + elementIdentifier );
916 
917                 if ( !anovaTable.isList() || !anovaTable.hasAttribute( "row.names" ) ) {
918                     log.debug( "No anovaresult for " + elementIdentifier );
919                     result.put( elementIdentifier, new TwoWayAnovaResult() );
920                     continue;
921                 }
922 
923                 result.put( elementIdentifier, new TwoWayAnovaResult( anovaTable ) );
924 
925             }
926         } catch ( REXPMismatchException e ) {
927             throw new RuntimeException( e );
928         }
929 
930         return result;
931 
932     }
933 
934     /**
935      * There is a pretty annoying limitation of this. The file must contain only one statement. You can get around this
936      * by using c(x<-1,x<-2). See testScript.R
937      * 
938      * @param is
939      */
940     protected void loadScript( InputStream is ) {
941         try {
942 
943             BufferedReader reader = new BufferedReader( new InputStreamReader( is ) );
944             String line = null;
945             StringBuilder buf = new StringBuilder();
946             while ( ( line = reader.readLine() ) != null ) {
947                 if ( line.startsWith( "#" ) || StringUtils.isBlank( line ) ) {
948                     continue;
949                 }
950                 buf.append( StringUtils.trim( line ) + "\n" );
951             }
952             is.close();
953             log.debug( buf.toString() );
954             this.voidEval( buf.toString() );
955 
956         } catch ( IOException e ) {
957             throw new RuntimeException( e );
958         }
959     }
960 
961     /**
962      * @param matrix
963      * @param matrixVarName
964      * @param rowNameExtractor e.g. you could use StringValueTransformer.getInstance()
965      * @return
966      */
967     private void assignRowAndColumnNames( DoubleMatrix<?, ?> matrix, String matrixVarName, Transformer rowNameExtractor ) {
968 
969         List<Object> rown = new ArrayList<Object>();
970         for ( Object o : matrix.getRowNames() ) {
971             rown.add( rowNameExtractor.transform( o ) );
972         }
973 
974         String rowNameVar = assignStringList( rown );
975         String colNameVar = assignStringList( matrix.getColNames() );
976 
977         String dimcmd = "dimnames(" + matrixVarName + ")<-list(" + rowNameVar + ", " + colNameVar + ")";
978         this.voidEval( dimcmd );
979     }
980 
981 }