Main Page | Namespace List | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | Examples

pdb_distance_matrix.cc

00001 /* Copyright 2004
00002    Stanford University
00003 
00004    This file is part of the DSR PDB Library.
00005 
00006    The DSR PDB Library is free software; you can redistribute it and/or modify
00007    it under the terms of the GNU Lesser General Public License as published by
00008    the Free Software Foundation; either version 2.1 of the License, or (at your
00009    option) any later version.
00010 
00011    The DSR PDB Library is distributed in the hope that it will be useful, but
00012    WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00013    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00014    License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public License
00017    along with the DSR PDB Library; see the file COPYING.LIB.  If not, write to
00018    the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00019    MA 02111-1307, USA. */
00020 
00021 #include "dsrpdb/PDB.h"
00022 #include <cstdlib>
00023 #include <iostream>
00024 #include <fstream>
00025 #include <dsrpdb/distance.h>
00026 
00027 #ifdef PDB_USE_BOOST_PROGRAM_OPTIONS
00028 #include <boost/program_options.hpp>
00029 #endif
00030 
00031 #ifdef PDB_USE_MAGICK
00032 #include <Magick++.h>
00033 
00034 
00041 int main(int argc, char *argv[]){
00042   std::string input_file, output_file;
00043   bool print_help=false;
00044   bool verbose=false;
00045   bool interactive=false;
00046   bool all_atoms=false;
00047   double contact_map_threshold=std::numeric_limits<double>::infinity();
00048 
00049 #ifdef PDB_USE_BOOST_PROGRAM_OPTIONS
00050   {
00051     boost::program_options::options_description o("Allowed options"), po, ao;
00052     o.add_options()
00053       ("help", boost::program_options::bool_switch(&print_help), "produce help message")
00054       ("verbose,v", boost::program_options::bool_switch(&verbose), "print verbose error messages")
00055       ("all-atoms,a", boost::program_options::bool_switch(&all_atoms),
00056        "Output the distances between all atoms, not just the C_alphas.")
00057       ("image,i",boost::program_options::value<std::string>(&output_file), "Output image file.")
00058       ("contact-threshold,c",boost::program_options::value<double>(&contact_map_threshold), "Output a contact map with the given threshold. (assuming the image argument is passed).")
00059       ("query,q", boost::program_options::bool_switch(&interactive), "Enter a mode to do interactive queries about edge lengths.");
00060     po.add_options()
00061       ("input-pdb", boost::program_options::value< std::string>(&input_file), "input file");
00062 
00063     ao.add(o).add(po);
00064     
00065     boost::program_options::positional_options_description p;
00066     p.add("input-pdb", 1);
00067     p.add("output-image", 1);
00068     
00069     boost::program_options::variables_map vm;
00070     boost::program_options::store(boost::program_options::command_line_parser(argc, argv).
00071                                   options(ao).positional(p).run(), vm);
00072     boost::program_options::notify(vm);
00073     
00074     
00075     
00076     //boost::program_options::store(boost::program_options::parse_command_line(argc, argv, desc), vm);
00077     //boost::program_options::notify(vm);    
00078     
00079     if (input_file.empty() || print_help) {
00080       std::cout << "This program writes the distance matrix to an image file.\n";
00081       std::cout << "usage: " << argv[0] << " input-pdb\n\n";
00082       std::cout << o << "\n";
00083       return EXIT_FAILURE;
00084     }
00085   }
00086 
00087 #else
00088   if (argc != 3){
00089     std::cerr << "The program is built without boost/program_options.hpp.\n";
00090     std::cerr << "useage: " << argv[0] << " [-c] file.pdb template%d[%c].pdb" << std::endl;
00091     return EXIT_FAILURE;
00092   }
00093   input_file = argv[1];
00094   output_file = argv[2];
00095 #endif
00096  
00097   
00098   
00099   
00100   std::ifstream in(input_file.c_str());
00101   if (!in) {
00102     std::cerr << "Error opening input file " << input_file << std::endl;
00103     return EXIT_FAILURE;
00104   }
00105 
00106   dsrpdb::PDB pdb(in, verbose);
00107 
00108   if (verbose) std::cout << "Input PDB has " << pdb.number_of_models() << " models." << std::endl;
00109   if (pdb.number_of_models() != 1){
00110     std::cout << "Attempting to write multiple image files: assuming the output argument has a %d in it.\n";
00111   }
00112  
00113   for (unsigned int i=0;i< pdb.number_of_models(); ++i){
00114     dsrpdb::Model &m= pdb.model(i);
00115     dsrpdb::Matrix arr;
00116     if (all_atoms) {
00117       arr= distance_matrix(m);
00118     } else {
00119       arr= ca_distance_matrix(m);
00120     }
00121 
00122     //std::cout << arr.dim1() << " " << arr.dim2() << std::endl;
00123 
00124     if (interactive) {
00125       std::cout << "Entering interactive mode for model " << i << std::endl;
00126       std::cout << "There are " << arr.dim1() << " atoms." << std::endl;
00127       while (true) {
00128         std::cout << "> " << std::flush;
00129         char buf[100];
00130         std::cin.getline(buf, 100);
00131         if (strcmp(buf, "quit") ==0  || strcmp(buf,"exit") ==0){
00132           break;
00133         } else {
00134           int a,b;
00135           if (sscanf(buf, "%d %d", &a, &b)==2){
00136             if (a > arr.dim1() || b > arr.dim1() || a<0 || b <0){
00137               std::cerr << "Index " << a << " or " << b << " out of range." << std::endl;         
00138             } else {
00139               std::cout << arr[a][b] << std::endl;
00140             }
00141           } else if (sscanf(buf, "r %d %d", &a, &b) ==2) {
00142             if (a > 0 && b > 0) {
00143               dsrpdb::Residue::Index ia(a), ib(b);
00144               if (m.chain(0).has_residue(ia) && m.chain(0).has_residue(ib)) {
00145                 dsrpdb::Point pa= m.chain(0).residue(ia).atom(dsrpdb::Residue::AL_CA).cartesian_coords();
00146                 dsrpdb::Point pb= m.chain(0).residue(ib).atom(dsrpdb::Residue::AL_CA).cartesian_coords();
00147                 dsrpdb::Squared_distance sd;
00148                 double d= std::sqrt(sd(pa,pb));
00149                 std::cout << "The Ca distance is " << d << std::endl;
00150               } else {
00151                 std::cerr << "Invalid residues picked.\n";
00152               }
00153             } else {
00154               std::cerr << "Invalid residues indices picked (must be greater than 0).\n";
00155             }
00156           } else {
00157             std::cerr << "Error parsing input. Please enter a pair of indices or 'exit'." << std::endl;
00158           }
00159          
00160         }
00161       }
00162     }
00163 
00164     if (!output_file.empty()) {
00165 
00166       char buf[100000];
00167       if (pdb.number_of_models() != 1){
00168         sprintf(buf, output_file.c_str(), i);
00169       } else {
00170         sprintf(buf, output_file.c_str());
00171       }
00172       Magick::Geometry geom(arr.dim1(),arr.dim2());
00173 
00174 
00175       Magick::Image im(geom, "red");
00176 
00177       if (contact_map_threshold != std::numeric_limits<double>::infinity()){
00178         for (int j=0; j< arr.dim1(); ++j){
00179           for (int k=0; k< arr.dim2(); ++k){
00180             if (arr[j][k] < contact_map_threshold){
00181               im.pixelColor(j, k, Magick::ColorGray(1));
00182             } else {
00183               im.pixelColor(j,k, Magick::ColorGray(0));
00184             }
00185           }
00186         }
00187       } else {
00188         double max=0;
00189         for (int i=0; i< arr.dim1(); ++i){
00190           for (int j=0; j< arr.dim2(); ++j){
00191             max=std::max(max, arr[i][j]);
00192             //min=std::min(min, arr[i][j]);
00193           }
00194         }
00195         //std::cout << "Maximum distance is " << max << std::endl;
00196         for (int i=0; i< arr.dim1(); ++i){
00197           for (int j=0; j< arr.dim2(); ++j){
00198             im.pixelColor(i,j, Magick::ColorGray(arr[i][j]/max));
00199           }
00200         }
00201       }
00202         
00203       im.write(buf);
00204     }
00205   }
00206   return EXIT_SUCCESS;
00207   }
00208 
00209 #else
00210 
00211 int main(int, char*[]){
00212   bool this_program_requires_image_magick;
00213   std::cerr << "This program requires Image Magick++\n";
00214   return EXIT_FAILURE;
00215 }
00216 
00217 #endif