001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.survival.cox;
022
023import org.biojava.nbio.survival.cox.matrix.Matrix;
024
025import java.util.ArrayList;
026import java.util.LinkedHashMap;
027
028//import java.util.Collections;
029
030/**
031 *
032 * @author Scooter Willis <willishf at gmail dot com>
033 */
034public class ResidualsCoxph {
035
036        /**
037         *
038         */
039        public enum Type {
040
041                /**
042                 *
043                 */
044                dfbeta,
045                /**
046                 *
047                 */
048                dfbetas,
049                /**
050                 *
051                 */
052                score;
053        }
054
055        /**
056         *
057         * @param ci
058         * @param type
059         * @param useWeighted
060         * @param cluster
061         * @return
062         * @throws Exception
063         */
064        public static double[][] process(CoxInfo ci, Type type, boolean useWeighted, ArrayList<String> cluster) throws Exception {
065                Type otype = type;
066                if (type == Type.dfbeta || type == Type.dfbetas) {
067                        type = Type.score;
068                        //if missing weighted is a required so never missing
069                } //64 2 625 310
070
071                double[][] rr = null;
072                if (type == Type.score) {
073                        rr = CoxScore.process(ci.method, ci.survivalInfoList, ci, false);
074                }
075
076                //debug
077//        if (false) {
078//            for (int i = 0; i < ci.survivalInfoList.size(); i++) {
079//                SurvivalInfo si = ci.survivalInfoList.get(i);
080//                System.out.print("residuals " + si.getOrder() + " " + si.getClusterValue());
081//                for (int j = 0; j < 2; j++) {
082//                    System.out.print(" " + rr[i][j]);
083//                }
084//                System.out.println();
085//            }
086//        }
087
088
089                double[][] vv = null;
090                if (ci.getNaiveVariance() != null) {
091                        vv = ci.getNaiveVariance();
092                } else {
093                        vv = ci.getVariance();
094                }
095                if (otype == Type.dfbeta) {
096                        //rr <- rr %*% vv
097                        rr = Matrix.multiply(rr, vv);
098                } else if (otype == Type.dfbetas) {
099                        //rr <- (rr %*% vv) %*% diag(sqrt(1/diag(vv)))
100                        double[][] d1 = Matrix.multiply(rr, vv);
101                        double[][] d2 = Matrix.diag(Matrix.sqrt(Matrix.oneDivide(Matrix.diag(vv))));
102                        rr = Matrix.multiply(d1, d2);
103                }
104
105
106
107                if (useWeighted) {
108                        double[] weighted = ci.getWeighted();
109                        rr = Matrix.scale(rr, weighted);
110                }
111                if (cluster != null && cluster.size() > 0) {
112                        rr = rowsum(rr, cluster);
113                }
114
115
116                return rr;
117        }
118
119        /**
120         * From R in residuals.coxph.S rowsum(rr, collapse)
121         *
122         * @param rr
123         * @param sets
124         * @return
125         */
126        private static double[][] rowsum(double[][] rr, ArrayList<String> sets) throws Exception {
127                LinkedHashMap<String, Double> sumMap = new LinkedHashMap<String, Double>();
128                if (rr.length != sets.size()) {
129                        throw new Exception("Cluster value for each sample are not of equal length n=" + rr.length + " cluster length=" + sets.size());
130                }
131                double[][] sum = null;
132                for (int j = 0; j < rr[0].length; j++) {
133                        for (int i = 0; i < sets.size(); i++) {
134                                String s = sets.get(i);
135                                Double v = sumMap.get(s); //get in order
136                                if (v == null) {
137                                        v = 0.0;
138                                }
139                                v = v + rr[i][j];
140                                sumMap.put(s, v);
141
142                        }
143                        if (sum == null) {
144                                sum = new double[sumMap.size()][rr[0].length];
145                        }
146
147                        ArrayList<String> index = new ArrayList<String>(sumMap.keySet());
148                        //sorting does seem to make a difference in test cases at the .0000000001
149           //     ArrayList<Integer> in = new ArrayList<Integer>();
150           //     for (String s : index) {
151           //         in.add(Integer.parseInt(s));
152           //     }
153           //     Collections.sort(index);
154
155                        for (int m = 0; m < index.size(); m++) {
156                                String key = index.get(m);
157                                sum[m][j] = sumMap.get(key);
158                        }
159
160                        sumMap.clear();
161                }
162
163                return sum;
164
165        }
166
167        /**
168         * @param args the command line arguments
169         */
170        public static void main(String[] args) {
171                // TODO code application logic here
172        }
173}