A library for working with phylogenetic data.
v0.25.0
heatmap_colorization.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2021 Lucas Czech
4 
5  This program is free software: you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 
18  Contact:
19  Lucas Czech <lczech@carnegiescience.edu>
20  Department of Plant Biology, Carnegie Institution For Science
21  260 Panama Street, Stanford, CA 94305, USA
22 */
23 
32 
37 
40 
41 #include <algorithm>
42 #include <cassert>
43 #include <cmath>
44 #include <functional>
45 #include <limits>
46 #include <stdexcept>
47 
48 namespace genesis {
49 namespace population {
50 
51 // =================================================================================================
52 // Allele Frequency Heatmap
53 // =================================================================================================
54 
55 std::pair<utils::Matrix<utils::Color>, double> HeatmapColorization::spectrum_to_image(
56  Spectrum const& spectrum
57 ) const {
58  using namespace utils;
59 
60  // Check.
61  if( color_map_.empty() ) {
62  throw std::runtime_error(
63  "ColorMap has to be assigned a palette before using HeatmapColorization."
64  );
65  }
66  if( log_scale_ && diverging_scale_ ) {
67  throw std::runtime_error(
68  "Cannot use log and divergnig scale with HeatmapColorization."
69  );
70  }
71 
72  // Edge case.
73  if( spectrum.values.empty() ) {
74  return { {}, 0 };
75  }
76 
77  // Get the row size that we need. We later also confirm that this is consistent
78  // across all spetra, to make sure that the data is actually a matrix/image.
79  assert( ! spectrum.values.empty() );
80  size_t rows = spectrum.values[0].size();
81 
82  // We need two passes through the data: first, find the max entry and the sum per column, then
83  // convert to scale. While doing the first pass, make sure that the data is actually a matrix.
84  double abs_max = 0.0;
85  for( size_t c = 0; c < spectrum.values.size(); ++c ) {
86  auto& col = spectrum.values[c];
87  if( col.size() != rows ) {
88  throw std::runtime_error(
89  "Invalid allele frequency spectrum with inconsistent number of rows."
90  );
91  }
92  // double col_sum = 0.0;
93  // for( auto const& val : col ) {
94  // // assert( std::isfinite(val) );
95  // col_sum += val;
96  // }
97  // if( std::isfinite( col_sum ) && col_sum > 0.0 ) {
98  // for( auto& val : col ) {
99  // val /= col_sum;
100  // abs_max = std::max( abs_max, val );
101  // }
102  // }
103  for( auto& val : col ) {
104  abs_max = std::max( abs_max, val );
105  }
106  }
107  // LOG_DBG << "abs_max " << abs_max;
108 
109  // Now convert to color values.
110  auto image = Matrix<Color>( rows, spectrum.values.size() );
111  for( size_t c = 0; c < spectrum.values.size(); ++c ) {
112  auto const& col = spectrum.values[c];
113  assert( col.size() == rows );
114 
115  // Get the max value of the current column.
116  double col_max = 0.0;
117  for( auto const& val : col ) {
118  col_max = std::max( col_max, val );
119  }
120 
121  // Get the max value that we want to use for normalization.
122  double const used_max = static_cast<double>( max_per_column_ ? col_max : abs_max );
123 
124  // LOG_DBG << "used_max " << used_max;
125 
126  // Do the actual per-bin convertion to color.
127  for( size_t r = 0; r < col.size(); ++r ) {
128  assert( col[r] <= abs_max );
129  assert( col[r] <= col_max );
130 
131  // Get the row where to write the color to.
132  size_t row_idx = invert_vertically_ ? rows - r - 1 : r;
133  assert( row_idx < image.rows() );
134 
135  // Special case: no bin filled at all in this window. That means, there were no variants
136  // in the whole window. If needed, mark with special "empty" color, which we
137  // stored in the mask color of the color map.
138  if( col_max == 0.0 && use_empty_window_color_ ) {
139  image( row_idx, c ) = color_map_( std::numeric_limits<double>::quiet_NaN() );
140  continue;
141  }
142 
143  // Set the color for the current pixel.
144  if( log_scale_ ) {
145  // Special case for log scaling: If either the pixel value or the total max
146  // for the colum is 1 or below, we cannot use log scaling for them (as we are working
147  // with integer number counts here), so we simply use the min value instead.
148 
149  // TODO deactivated for log of low normalized numbers
150  // if( col[r] <= 1 || used_max <= 1 ) {
151  // image( row_idx, c ) = color_map_( 0.0 );
152  // } else {
153  // double frac = std::log( static_cast<double>( col[r] )) / std::log( used_max );
154  // image( row_idx, c ) = color_map_( frac );
155  // }
156 
157  assert( std::isfinite(col[r]) );
158  // assert( col[r] >= 0 );
159  // double frac = std::log( static_cast<double>( col[r] )) / std::log( used_max );
160  // image( row_idx, c ) = color_map_( frac );
161 
162  // TODO 1.0 is not a good min
163  // if( 1.0 > col[r] || col[r] > used_max )
164  // {
165  // LOG_DBG << "fucked " << col[r];
166  // }
167  // assert( 1.0 <= col[r] );
168  // assert( col[r] <= used_max );
169 
170  if( used_max <= 1.0 ) {
171  image( row_idx, c ) = color_map_( std::numeric_limits<double>::quiet_NaN() );
172  continue;
173  }
174 
175  auto norm = ColorNormalizationLogarithmic( 1.0, used_max );
176  image( row_idx, c ) = color_map_( norm(col[r]) );
177 
178  } else if( diverging_scale_ ) {
179 
180  auto norm = ColorNormalizationDiverging( -used_max, used_max );
181  image( row_idx, c ) = color_map_( norm(col[r]) );
182 
183  } else {
184  double frac = static_cast<double>( col[r] ) / used_max;
185  image( row_idx, c ) = color_map_( frac );
186  }
187  }
188  }
189 
190  // Return the image and the appropriate max value used for the color scaling.
191  return { image, max_per_column_ ? 1 : abs_max };
192 }
193 
194 std::pair<utils::SvgGroup, double> HeatmapColorization::spectrum_to_svg(
195  Spectrum const& spectrum,
196  utils::SvgMatrixSettings settings
197 ) const {
198  // Generate the pixel color image matrix.
199  auto const spec_img_and_max = spectrum_to_image( spectrum );
200 
201  // Return the svg group and the max value here.
202  return { make_svg_matrix( spec_img_and_max.first, settings ), spec_img_and_max.second };
203 }
204 
206  HeatmapColorization::Spectrum const& spectrum,
207  std::shared_ptr<utils::BaseOutputTarget> target
208 ) const {
209 
210  // Generate the pixel color image matrix, and write the image to file.
211  auto const spec_img_and_max = spectrum_to_image( spectrum );
212  utils::BmpWriter().write( spec_img_and_max.first, target );
213 
214  // Return only the max value here.
215  return spec_img_and_max.second;
216 }
217 
218 } // namespace population
219 } // namespace genesis
genesis::population::HeatmapColorization::spectrum_to_image
std::pair< utils::Matrix< utils::Color >, double > spectrum_to_image(Spectrum const &spectrum) const
Definition: heatmap_colorization.cpp:55
genesis::utils::make_svg_matrix
SvgGroup make_svg_matrix(Matrix< Color > const &mat, SvgMatrixSettings settings, std::vector< std::string > const &row_labels, std::vector< std::string > const &col_labels)
Definition: formats/svg/matrix.cpp:50
genesis::population::HeatmapColorization::spectrum_to_bmp_file
double spectrum_to_bmp_file(Spectrum const &spectrum, std::shared_ptr< utils::BaseOutputTarget > target) const
Definition: heatmap_colorization.cpp:205
heatmap_colorization.hpp
genesis::population::HeatmapColorization::spectrum_to_svg
std::pair< utils::SvgGroup, double > spectrum_to_svg(Spectrum const &spectrum, utils::SvgMatrixSettings settings={}) const
Definition: heatmap_colorization.cpp:194
logging.hpp
Provides easy and fast logging functionality.
genesis::utils::BmpWriter::write
void write(Matrix< Color > const &image, std::shared_ptr< utils::BaseOutputTarget > target) const
Write a full 24bit RGB Color image to an output target.
Definition: utils/formats/bmp/writer.cpp:53
genesis::population::HeatmapColorization::Spectrum::values
std::vector< std::vector< double > > values
Definition: heatmap_colorization.hpp:72
genesis::utils::SvgMatrixSettings
Definition: formats/svg/matrix.hpp:56
genesis
Container namespace for all symbols of genesis in order to keep them separate when used as a library.
Definition: placement/formats/edge_color.cpp:42
norm_diverging.hpp
genesis::utils::BmpWriter
Write Bitmap image files.
Definition: utils/formats/bmp/writer.hpp:60
genesis::population::HeatmapColorization::Spectrum
Definition: heatmap_colorization.hpp:64
norm_linear.hpp
writer.hpp
normalization.hpp
norm_logarithmic.hpp