#!/usr/local/bin/perl


###############################################################################
# Basic SBEAMS setup
###############################################################################
use strict;
use FindBin;

use lib "$FindBin::Bin/../../lib/perl";

use vars qw ($q $sbeams $sbeamsMOD $PROG_NAME
             $current_username $modification_helper );

use SBEAMS::Connection qw($q $log);
use SBEAMS::Connection::Settings;
use SBEAMS::Connection::Tables;

use SBEAMS::PeptideAtlas;
use SBEAMS::PeptideAtlas::Settings;
use SBEAMS::PeptideAtlas::Tables;
use SBEAMS::PeptideAtlas::ModificationHelper;
use SBEAMS::PeptideAtlas::ConsensusSpectrum;
use SBEAMS::Proteomics::SpecViewer;


use File::Basename;

my $denormalize = 0;
my $curveType = 'function';

my $error_state = 0;
$sbeams = new SBEAMS::Connection;
$sbeamsMOD = new SBEAMS::PeptideAtlas;
$sbeamsMOD->setSBEAMS($sbeams);

my $consensus = new SBEAMS::PeptideAtlas::ConsensusSpectrum;
$consensus->setSBEAMS($sbeams);

$PROG_NAME="ShowOneSpectrum";

###############################################################################
# Define global variables if any and execute main()
###############################################################################

my $ymax = 10;

main();


###############################################################################
# Main Program:
#
# If $sbeams->Authenticate() succeeds, print header, process the CGI request,
# print the footer, and end.
###############################################################################
sub main {


  #### Do the SBEAMS authentication and exit if a username is not returned
  exit unless ($current_username = $sbeams->Authenticate(
    permitted_work_groups_ref=>['PeptideAtlas_user','PeptideAtlas_admin',
      'PeptideAtlas_readonly', 'PeptideAtlas_exec'],
    #connect_read_only=>1,
    allow_anonymous_access=>1,
  ));

  #### Print the header, figure and do what the user want, and print footer
  my %parameters;
  $sbeams->parse_input_parameters(q=>$q,parameters_ref=>\%parameters);
  $sbeams->processStandardParameters(parameters_ref=>\%parameters);

  $denormalize = 0 if $parameters{normalize};
  $curveType = $parameters{curveType} if $parameters{curveType};
 
  my %results;
  my $ce_vals;
  my $has_median = 0;
  my $max_intensity = 0;
  # First, grab da peaks
  use Data::Dumper;
  my %spectrum_plot = ( peak_array => [],
                        charge => 0,
                        seq => '' );
  for my $ce ( qw( low mlow medium mhigh high avg  ) ) {

    my $ce_results = mk_plot( %parameters, ce => $parameters{$ce} );
    if ( $ce eq 'medium' ) {
#      die Dumper( $ce_results );
      $spectrum_plot{mz} = $ce_results->{mz};
      $spectrum_plot{charge} = $ce_results->{charge};
      $spectrum_plot{sequence} = $ce_results->{modified_sequence};

      my $cnt = 0;
      for my $mz ( @{$ce_results->{masses}} ) {
        my $intensity = $ce_results->{intensities}->[$cnt];
        push @{$spectrum_plot{peak_array}}, [ $mz, $intensity ];
        $cnt++;
      }
    }


    next unless $ce_results->{max_intensity};

    # Going to base peaks to watch on this guy...
    $has_median++ if $ce eq 'medium';

    $ce_vals ||= $sbeamsMOD->get_Agilent_ce( %{$ce_results} );

    $ce_results->{ce} = $ce_vals->{$ce};
    $ce_results->{level} = ucfirst($ce);

    $parameters{modified_sequence} = $ce_results->{modified_sequence};
    $parameters{charge} = $ce_results->{charge};
    $max_intensity = $ce_results->{max_intensity} if $ce_results->{max_intensity} > $max_intensity; 

    $results{spectra} ||= {};
    $results{spectra}->{$ce} = $ce_results;
  }

  # This has to go into the header
  my $peak_plot = '';
  if ( $has_median ) {
#    die Dumper( %results );
#    my $diestmt = '';
#>    my $cnt = 0;
#    die Dumper( $results{spectra}->medium}->{intensities} );
#    for my $int ( sort { $b <=> $a } @{$results{spectra}->{medium}->{intensities}} ) {
#      $diestmt .= "$results{spectra}->{medium}->{labels}->[$cnt] = $int\n";
#      $cnt++;
#    }
#    die $diestmt;

    my $num_ions = get_num_ions( $parameters{modified_sequence} );
    my $chart_data = get_top_peaks( \%results, $num_ions );
    
    # Trying to get CE and spectra plots on same scale... 
    $parameters{ymax} = $ymax * 1.075;

    $peak_plot = get_chart( data => $chart_data,seq => $parameters{modified_sequence}, chg => $parameters{charge} );
  }


  my $pagename = ( $parameters{pagename} ) ? "$parameters{pagename}" : '';
  my $content = "<H2><CENTER>$pagename</CENTER></H2><BR>\n<TABLE>";

  if ( $ymax < 600 ) {
    $content .= "<TR><TD><H3><FONT color=red><CENTER>This spectrum has low signal, CE plot may not be accurate</CENTER></H3></TD></TR>";
  }
  $content .= "<TR><TD ALIGN=CENTER> <div id='chart_div'></div> </TD></TR>\n";
  if ( $ymax < 600 ) {
    $content .= "<TR><TD><H3><FONT color=red><CENTER>This spectrum has low signal, CE plot may not be accurate</CENTER></H3></TD></TR>";
  }

  if ( $spectrum_plot{charge} ) {
    my $start_scan = 1;
    my $spectrum_name = $spectrum_plot{sequence} . '+' . $spectrum_plot{charge};

    #### Lorikeet me! ##################################################
    my $lorikeet = new SBEAMS::Proteomics::SpecViewer;
    my $lorispectrum = $lorikeet->generateSpectrum(
              charge            => $spectrum_plot{charge},
              modified_sequence => $spectrum_plot{sequence},
              precursor_mass    => $spectrum_plot{mz},
              scan              => $start_scan,
              name              => $spectrum_name,
              a_ions            => $parameters{'ShowA'},
              b_ions            => $parameters{'ShowB'},
              c_ions            => $parameters{'ShowC'},
              x_ions            => $parameters{'ShowX'},
              y_ions            => $parameters{'ShowY'},
              z_ions            => $parameters{'ShowZ'},
              spectrum          => $spectrum_plot{peak_array} 
              );

#    $content .= "<TR><TD>$lorispectrum</TD></TR>";

  }

  $content .= "<TR><TD ALIGN=CENTER><H2>Spectra for $parameters{modified_sequence} <SUP>+$parameters{charge} </SUP></H2></TD></TR>\n";
  $content .= "</TABLE>\n";

	my $url_base = "ShowConsensusSpectrum?plot_only=true;consensus_library_spectrum_id=";
  for my $ce ( qw( low mlow medium mhigh high  ) ) {
	  $content .= "<BR><IFRAME ID=$ce SRC=$url_base" . $parameters{$ce} . "></IFRAME><BR>\n";
    $content .= qq~ 
      <SCRIPT LANGUANGE=javascript>
//        parent.document.getElementById('$ce').height = document['body'].offsetHeight;
//        parent.document.getElementById('$ce').width = document['body'].offsetWidth;
        parent.document.getElementById('$ce').height = 800;
        parent.document.getElementById('$ce').width = 1200;
      </SCRIPT>";
    ~;
	}

  $sbeamsMOD->display_page_header( header_info => $peak_plot );
  print $content;
  $sbeamsMOD->display_page_footer();

  $sbeams->display_page_footer(close_tables=>'YES',
    separator_bar=>'YES',display_footer=>'NO');

} # end main


sub mk_plot {

    my %parameters = @_;
	  $parameters{consensus_library_spectrum_id} = $parameters{ce};

    my %results = ( image => '', lib_desc => '' );
    return \%results unless $parameters{ce};

    #### Define some general variables
    my ($i,$element,$key,$value,$sql);


    my $apply_action  = $q->param('apply_action');

    my @charge;
    push (@charge, 1);
    push (@charge, 2);
    push (@charge, 3);

    $parameters{'ionlab'} = "Horizontal" unless $parameters{'ionlab'};
    my ($labangle,$fjust);
    if ($parameters{'ionlab'} eq "Vertical") {
      $labangle = 90;
      $fjust = 0;
    } else {
      $labangle = 0;
      $fjust = 0.5;
    }

    my $sql = qq~
       SELECT sequence,charge,mz_exact,protein_name,
       protein_name_alt,modified_sequence, consensus_library_name, file_path, entry_idx
       FROM $TBAT_CONSENSUS_LIBRARY_SPECTRUM CLS
       JOIN $TBAT_CONSENSUS_LIBRARY CL ON CL.consensus_library_id = CLS.consensus_library_id
       WHERE consensus_library_spectrum_id = '$parameters{consensus_library_spectrum_id}'
    ~;

    my @rows = $sbeams->selectSeveralColumns($sql);

    foreach my $row (@rows) {
      my ($seq, $chg, $mz, $pn, $pn_alt, $mod_seq, $lib_name, $file_path, $entry_idx ) = @{$row};
      $parameters{'peptide'} = $seq;
      $parameters{'charge'} = $chg;
      $parameters{'precursor_mass'} = $mz;
      $parameters{'protein_name'} = $pn;
      $parameters{'protein_name_alt'} = $pn_alt;
      $parameters{'modified_sequence'} = $mod_seq;
      $parameters{'library_name'} = $lib_name;
      $parameters{'file_path'} = $file_path;
      $parameters{'entry_idx'} = $entry_idx;
    }

    $results{modified_sequence} = $parameters{modified_sequence};
    $results{charge} = $parameters{charge};
    $results{mz} = $parameters{'precursor_mass'};

    $results{lib_desc} = "$parameters{library_name}";
    $results{lib_desc} =~ s/\.sptxt//;

    my $peptide = $parameters{peptide};
    $peptide =~ s/^.\.//;
    $peptide =~ s/\..$//;
    my $precursor_mass = $parameters{'precursor_mass'};
    my $charge = $parameters{'charge'};
    my $protein_name = $parameters{'protein_name'};
    my $protein_name_alt = $parameters{'protein_name_alt'};
    my $modified_sequence = $parameters{'modified_sequence'};

    $parameters{assumed_charge} = $charge;

    ## if no modified_sequence, use unmodified sequence
    unless ($modified_sequence) {
        $modified_sequence = $parameters{peptide};
    }

    #### Calculate peptide mass information
    my $masstype = $parameters{masstype} || 0;

    $modification_helper = new SBEAMS::PeptideAtlas::ModificationHelper();

    if ( !$parameters{spectrum} ) {

      my $peaks = $consensus->get_spectrum_peaks( %parameters, denormalize => 0, strip_unknowns => 1 );
      $peaks->{mz} = $results{mz};
      $peaks->{modified_sequence} = $results{modified_sequence};
      $peaks->{charge} = $results{charge};
#      die Dumper( $peaks );
      return $peaks;
    }

} # end printEntryForm

sub get_top_peaks {
  my $results = shift;
  my $num_peaks = shift || 7;
  return '' unless $results->{spectra}->{medium};
  my %top_peaks;
  $denormalize = 0;
#  die Dumper( $results->{spectra}->{medium} );

  my @medium = $consensus->get_top_n_peaks( n_peaks => $num_peaks,
                                        denormalize => $denormalize,
                                precursor_exclusion => 5,
                                            spectra => $results->{spectra}->{medium} );

#  my @srm_atlas = get_SRMAtlas_peaks();
#  die Dumper( $medium[1] );
                                                     
#  die Dumper( @medium );
  $top_peaks{medium} = $medium[0];
  my $peak_list = $medium[1];
  $ymax = $top_peaks{medium}->{$peak_list->[0]};





  for my $ce ( qw( low mlow mhigh high ) ) {
    $top_peaks{$ce} = $consensus->get_top_n_peaks( n_peaks => $num_peaks,
                                               denormalize => $denormalize,
                                                   spectra => $results->{spectra}->{$ce},
                                                 peak_list => $peak_list );
    
    next unless $top_peaks{$ce};
    for my $ion ( keys( %{$top_peaks{$ce}} ) ) {
      next unless $ion;
      $ymax = $top_peaks{$ce}->{$ion} if $top_peaks{$ce}->{$ion} > $ymax;
    }
  }

  my $fpmax = 1;
  push @{$peak_list}, 'Fp';
  for my $ce ( qw( low mlow medium mhigh high ) ) {
    next unless $results->{spectra}->{$ce};
    my $fp = $results->{spectra}->{$ce}->{postfrag_precursor_intensity};
    $top_peaks{$ce}->{Fp} = $fp; 
    $fpmax = ( $fpmax < $fp ) ? $fp : $fpmax;
  }
  for my $ce ( qw( low mlow medium mhigh high ) ) {
    next unless $results->{spectra}->{$ce};
    $top_peaks{$ce}->{Fp} *= ( (1.2 * $ymax)/ $fpmax ); 
  }


  my $columns .= " data.addColumn( 'string', 'CE' )\n";
  for my $ion ( @{$peak_list} ) {
    $ion = 'Precursor (unfrag, scaled)' if $ion =~ /F/;
#    $ion .= ' (scaled)' if $ion =~ /F/;
    $columns .= " data.addColumn( 'number', '$ion' )\n";
  }
#        data.addRows(4);
#        data.setValue(0, 0, '2004');
#        data.setValue(0, 1, 1000);
#        data.setValue(0, 2, 400);

  my $plot_data = '';
  my %plot_data;
  my $num_rows = 0;
  my $row = 0;
  for my $ce ( qw( low mlow medium mhigh high ) ) {
    next unless $results->{spectra}->{$ce};
    my $col = 0;
    $num_rows++;
    my $ce_label = ( $ce eq 'medium' ) ? "$results->{spectra}->{$ce}->{ce}*" : $results->{spectra}->{$ce}->{ce};
    $plot_data .= "data.setValue( $row, $col, '$ce_label' );\n";
    $col++;
    for my $ion ( @{$peak_list} ) {
      $ion = 'Fp' if $ion =~ /Precursor/;
      $top_peaks{$ce}->{$ion} = sprintf( "%i", $top_peaks{$ce}->{$ion} );

      $plot_data .= "data.setValue( $row, $col, $top_peaks{$ce}->{$ion} );\n";
      $col++;
    }
    $row++;
  }
  my $chart_data = qq~
    $columns
    data.addRows($num_rows);
    $plot_data
  ~;
  return $chart_data;
}

sub get_chart {

  my %args = @_;

  my $data = $args{data} || return '';
  my $ion_info = "$args{seq}+$args{chg}" || '';
  
  my $chart = qq~
    <script type="text/javascript" src="https://www.google.com/jsapi"></script>
    <script type="text/javascript">
      google.load("visualization", "1", {packages:["corechart"]});
      google.setOnLoadCallback(drawChart);
      function drawChart() {
        var data = new google.visualization.DataTable();
        $data
        var chart = new google.visualization.LineChart(document.getElementById('chart_div'));

//         chart.draw(data, {width: 800, height: 600, title: 'Intensity vs CE'});

        chart.draw(data, {width: 800, height: 600, title: 'Intensity vs CE, $ion_info',
                          hAxis: {title: 'CE value (Volts)             * is vendor optimum', showTextEvery: 1,
                                     textStyle: {color: 'black', fontSize: 16, fontWeight: 'bold' },
												          },
                          vAxis: {title: 'Intensity', minValue: 0, maxValue: $ymax },
                          pointSize: 7, curveType: '$curveType'
                         }
            );
      }
    </script>
  ~;
  return $chart;

#  <body>
#    <div id="chart_div"></div>
#  </body>


}

###############################################################################
# get_library_spectrum
###############################################################################
sub get_library_spectrum 
{
  my %args = @_;

  my $inputfile = $args{'inputfile'} || "";
  my $verbose = $args{'verbose'} || "";
  my $consensus_library_spectrum_id = $args{'consensus_library_spectrum_id'} || "";
	
  # For backwards compatibility with NIST tables (TIQAM). 
  $consensus_library_spectrum_id ||= $args{NIST_library_spectrum_id};


  #### Define some general variables
  my ($i,$element,$key,$value,$sql);
  my (@rows,$nrows);


  #### Define the data hash
  my %spectrum;
  my @mass_intensities;


  #### If we have a msms_spectrum_id, get the data from the database
  if ($consensus_library_spectrum_id) {

    #### Extract the actual mass,intensity pairs from database
    $sql = "SELECT mz, relative_intensity ".
           "  FROM $TBAT_CONSENSUS_LIBRARY_SPECTRUM_PEAK ".
           " WHERE consensus_library_spectrum_id = '$consensus_library_spectrum_id'";

    # For backwards compatibility with NIST tables (TIQAM). 
		$sql = $sbeams->evalSQL( $sql );			 
		if ( $args{NIST_library_spectrum_id} ) {
			$sql =~ s/consensus/NIST/gm;
 		  $log->info( "Translated to NIST tables: $sql" );
		}

    my @mass_intensities = $sbeams->selectSeveralColumns($sql);


    #### If we still have no spectrum data, then bail out
    unless (@mass_intensities) {
      print "\nERROR: Unable to get m/z,intensity pairs for ".
	"consensus_library_spectrum_id '$consensus_library_spectrum_id'.\n\n";
      $error_state++;
      return;
    }

    #### Extract rows into two arrays of masses and intensities
    my (@masses,@intensities);
    for ($i=0; $i<=$#mass_intensities; $i++) {
      push(@masses,$mass_intensities[$i]->[0]);
      push(@intensities,$mass_intensities[$i]->[1]);
    }

    $spectrum{n_peaks} = $#mass_intensities + 1;


    #### Put data into hash and return
    $spectrum{masses} = \@masses;
    $spectrum{intensities} = \@intensities;

    return %spectrum;

  #### Otherwise complain and return
  } else {
    print "\nERROR: Unable to determine which consensus_library_spectrum_id to load.\n\n";
    $error_state++;
    return;
  }

}

sub get_num_ions {
  my $sequence = shift;
  $log->info( "initial seq is $sequence" );
  $sequence =~ s/\[\d+\]//g;
  my $length = length $sequence;
  my $n_ions = int( $length/2 + 3.1 ) ;
  $n_ions = 12 if $n_ions > 12;
  $log->info( "mod seq is $sequence, length is $length, n_ions is $n_ions" );
  return $n_ions;
}


###############################################################################
# Fragment
###############################################################################
sub Fragment {
    my $peptide = shift;
    my $length = length($peptide);
    my @residues = ();
    my $i;

    for ($i=0; $i<$length; $i++) {
      if (substr($peptide,$i+1,1) eq '[') {
        push (@residues, substr($peptide,$i,6));
        $i = $i + 5;
      } elsif (substr($peptide,$i+1,1) =~ /\W/) {
        push (@residues, substr($peptide,$i,2));
        $i = $i + 1;
      } else {
        push (@residues, substr($peptide,$i,1));
      }
    }

    #### Return residue array
    return @residues;
}

###############################################################################
# CalcIons -- calculate theoretical ions (including modified masses)
###############################################################################
sub CalcIons {
    my %args = @_;
    my $i;

    my $residues_ref = $args{'Residues'};
    my @residues = @$residues_ref;
    my $charge = $args{'Charge'} || 1;
    my $length = scalar(@residues);
   
    my $modified_sequence = $args{'modified_sequence'};

    my @masses = $modification_helper->getMasses($modified_sequence);

    my $Nterm = 1.0078;
    my $Bion = 0.;
    my $Yion  = 19.0184;  ## H_2 + O

    my @Bcolor = (14) x $length;
    my @Ycolor = (14) x $length;

    my %masslist;
    my (@aminoacids, @indices, @rev_indices, @Bions, @Yions);


    #### Compute the ion masses
    for ($i = 0; $i<$length; $i++) {
      $Bion += $masses[$i];

      #### B index & Y index
      $indices[$i] = $i;
      $rev_indices[$i] = $length-$i;
      $Yion += $masses[ $rev_indices[$i] ]  if ($i > 0);

      #### B ion mass & Y ion mass
      $Bions[$i] = ($Bion + $charge*$Nterm)/$charge;
      $Yions[$i] = ($Yion + $charge*$Nterm)/$charge;
    }

    $masslist{residues} = \@residues;
    $masslist{indices} = \@indices;
    $masslist{Bions} = \@Bions;
    $masslist{Yions} = \@Yions;
    $masslist{rev_indices} = \@rev_indices;

    #### Return reference to a hash of array references
    return (\%masslist);
}

###############################################################################
# printIons
###############################################################################
sub printIons {
  my %args = @_;

  my $masslist_ref = $args{'masslist_ref'};
  my $color = $args{'color'} || 0;
  my $html = $args{'html'} || 0;
  my $charge = $args{'charge'};
  my $length = $args{'length'};
  my $theoretical_spectrum_ref = $args{'theoretical_spectrum_ref'};

# print "SEQ  #     B       Y    +$charge\n";
# print "--- --  ------  ------  --\n";

# my ($bcolbegin, $bcolend, $ycolbegin, $ycolend);

# my (%colors);
# $colors{2} = "#FF0000";
# $colors{4} = "#0000FF";
# $colors{3} = "#218D21";
# $colors{6} = "#F18080";
# $colors{11} = "#00088";
# $colors{10} = "#8FBE8F";


# #### Printing stuff
  for (my $i=0; $i < $length; $i++) {

#   #### If the output is in HTML, define the colorizing tags
#   if ($html) {

#     #### If a color for this B ion mass, set color tags
#     if ($masslist_ref->{Bcolor}->[$i] >= 2) {
#       $bcolbegin = "<FONT COLOR = $colors{$masslist_ref->{Bcolor}->[$i]}>";
#       $bcolend = "</FONT>";
#     #### else no color (default black)
#     } else {
#       $bcolbegin = "";
#       $bcolend = "";
#     }

#     #### If a color for this Y ion mass, set color tags
#     if ($masslist_ref->{Ycolor}->[$i] >= 2) {
#       $ycolbegin = "<FONT COLOR = $colors{$masslist_ref->{Ycolor}->[$i]}>";
#       $ycolend = "</FONT>";
#     #### else no color (default black)
#     } else {
#       $ycolbegin = "";
#       $ycolend = "";
#     }

#   }


    #### Define the m/z columns formats and values
    my $B_format = '%7.1f';
    my $Y_format = '%7.1f';
    my $B_value = $masslist_ref->{Bions}->[$i];
    my $Y_value = $masslist_ref->{Yions}->[$i];


    #### Special case --'s for first row
    if ($i == 0) {
      $Y_format = '%7s';
      $Y_value = '--  ';

    #### Special case --'s for last row
    } elsif ($i == ($length-1)) {
      $B_format = '%7s';
      $B_value = '--  ';
    }

#   #### Print out the data
#   printf "%3s %2d $bcolbegin$B_format$bcolend ".
#          "$ycolbegin$Y_format$ycolend %3d\n",
#     $masslist_ref->{residues}->[$i], $i+1,
#     $B_value, $Y_value, $length-$i;

    #### Fill the theoretical spectrum data in a different format
    #### (Residue,Index,Ion,Charge,m/z)
    $theoretical_spectrum_ref->[2*$length*($charge-1)+$i] =
      [$masslist_ref->{residues}->[$i],$i+1,'B',$charge,$B_value];
    $theoretical_spectrum_ref->[2*$length*($charge-1)+2*$length-1-$i] =
      [$masslist_ref->{residues}->[$i],$length-$i,'Y',
      $charge,$Y_value];


  } # end for

# print "\n";

} # end printIons