#!/usr/bin/perl -w
# Get experiment data from the database...
# Use -tabulate 1 to get output with 1 or more columns for each sub-experiment
# Use -locus LocusesList or -locus all to get output by locus rather than by spot
#
# When using -locus, use sum. or count. or mean. or sd. to get summaries of the
# relevant variables instead of comma-delimited lists...

use strict;
use lib "../Perl/modules";
use GenomicsUtils;
use Experiment;
use SubExperiment;
use Chip;
use Genome;
use Gene;
use lib "../stat";
use Vector;
use Getopt::Long;

{
    # attributes should be prefixed with SpotAttr. or ArrayDataAttr. if they are to
    # be fetched from those tables
    my $usage = "getExperiments.pl [-host host -debug 10]\n"
	. "   [-taxon T | -condition C1,..,CN | -chipId c | -spotSet ... | -expId exp1,..,expN | -subExpId se1,..,seN]\n"
	. "   [-attributes locusId,log2RatioNormalized] [-failedOK 0] [-dyeCorr .6] [-keepNulls 0]\n"
	. "   [-locus all | -locus o1,...oN [ -spots 1 ]] [-tabulate 1]\n";

    my $host = undef;
    my $debug = 0;
    my $taxon = "";		# taxon to fetch all data for
    my $condition = "";		# conditions to fetch all data for
    my $expId = "";		# list of experiments to fetch
    my $subExpId = "";		# list of subexperiments to fetch
    my $chipId = "";		# chip to fetch data for
    my $locus = "";		# list of locuss to fetch, or "all", or empty if not summarizing by locus
    my $spotSet = "";		# list of spots to fetch
    my $spots = 0;		# 1 if show spots for the given locuss instead of summarizing them

    # both per-spot and per-subexperiment attributes to load
    my $attrs = "";
    my $failedOK = 0;		# 0 if only keep spots with FAILED==0
    my $dyeCorrReq = .6;        # only keep spots with dyeCorr >= this value
    my $keepNulls = 0;		# keep spots in subexps even if missing a requested field?

    my $tabulate = 1;		# all the sub-experiments side by side

    (GetOptions('host=s' => \$host,
		'debug=i' => \$debug,
		'taxon=i' => \$taxon,
		'condition=s' => \$condition,
		'expId=s' => \$expId,
		'subExpId=s' => \$subExpId,
		'chipId=i' => \$chipId,
		'locus=s' => \$locus,
		'spotSet=s' => \$spotSet,
		'attributes=s' => \$attrs,
		'failedOK=s' => \$failedOK,
		'dyeCorrReq=f' => \$dyeCorrReq,
		'keepNulls=i' => \$keepNulls,
		'spots=i' => \$spots,
		'tabulate=i' => \$tabulate)
     && ($taxon || $expId || $subExpId || $chipId)) || die $usage;
    die "Cannot specify both -spotSet and -locus:\n$usage" if $spotSet && $locus;

    if ($attrs eq "") {
	if ($locus && !$spots) {
	    $attrs = "locusId,spotId,mean.log2RatioNormalized";
	} else {
	    $attrs = "spotId,locusId,log2RatioNormalized";
	}
	$attrs = "subExpId,$attrs" if !$tabulate;
    }
    GenomicsUtils::setDebug($debug) if $debug > 0;
    GenomicsUtils::connect($host);

    my @attrs = split /,/, $attrs;
    my @attrsNoFunc = @attrs;
    my @attrsFunc = ();
    my @func = ("sum", "count", "mean", "sd", "first");

    if ($locus && !$spots) {
	$keepNulls = 1; # always keep NULL values for spots, since we summarize over them
	foreach my $attr (@attrsNoFunc) {
	    if ($attr =~ m/^([a-z]+)\.(.*)$/) {
		die "Unknown function $1" if !Vector::within($1, @func);
		$attr = $2;
		push @attrsFunc, $1;
	    } else {
		push @attrsFunc, "";
	    }
	}
    }

    my %attrTypes = %{ &SubExperiment::attrTypes(@attrsNoFunc) };
    my @knownSpotAttrs = qw( spotId chipId sourceId seq mismatch verificationCode
			     contaminated sector sectorRow sectorCol );
    my @knownDataAttrs = qw(subExpId log2Ratio log2RatioNormalized control treatment
			    dyeCorr leftPix rightPix topPix bottomPix );
    map { die "Unknown attribute $_\n"
	      . "Spot attributes: " . join(",",@knownSpotAttrs) . "\n"
	      . "Data attributes: " . join(",", @knownDataAttrs) . "\n"
	      . "Other attributes: locusId\n"
	      . "Attr-table prefixes: SpotAttr. and ArrayDataAttr.\n"
	      if $attrTypes{$_} eq "";
      } (keys %attrTypes);
    my @attrSpotInfo = grep { $attrTypes{$_} eq "SpotInfo" } @attrsNoFunc;
    my @attrArrayData = grep { $attrTypes{$_} eq "ArrayData" } @attrsNoFunc;
    my @attrOrf = grep { $attrTypes{$_} eq "locusId" } @attrsNoFunc;
    my @attrSpotAttr = grep { $attrTypes{$_} eq "SpotAttr" } @attrsNoFunc;
    my @attrArrayDataAttr = grep { $attrTypes{$_} eq "ArrayDataAttr" } @attrsNoFunc;
    die "Must specify at least one ArrayData attribute\n" if scalar @attrArrayData == 0;

    my @subexps = ();
    my @exps = ();
    if ($subExpId ne "") {
	die $usage if $taxon ne "" || $condition ne "" || $expId ne "" || $chipId ne "";
	my @seIds = split /,/, $subExpId;
	foreach my $seId (@seIds) {
	    die "Bad subExpId $seId" if !($seId =~ m/^[0-9]+$/);
	    push @subexps, SubExperiment::new('subExpId' => $seId);
	}
	my @expIds = sort {$a <=> $b} (map { $_->expId() } @subexps);
	@exps = map {Experiment::new('expId'=>$_)} Vector::uniq(@expIds);
    } elsif ($expId ne "") {
	die $usage if $taxon ne "" || $condition ne "" || $chipId ne "";
	my @ids = split /,/, $expId;
	map { die "Bad expId $_" if !($_ =~ m/^[0-9]+$/) } @ids;
	@exps = map { Experiment::new('expId' => $_ ) } @ids;
    } elsif ($chipId ne "") {
	# handle chip&taxon combination by test below
	my $chip = Chip::new('chipId' => $chipId);
	@exps = $chip->experiments();
    } else {
	my @params = ("taxonomyId", $taxon);
	if ($condition ne "") {
	    my @conditions = split /,/, $condition;
	    push @params, "condition", \@conditions;
	}
	@exps = Experiment::getExperiments(@params);
    }
    if (scalar @exps > 0 && scalar @subexps == 0) {
	foreach my $exp (@exps) {
	    push @subexps, $exp->subExperiments()
		if $taxon eq "" || $exp->sampleTaxId() eq $taxon;
	}
    }
    die "No sub-experiments found for the given criteria\n" if scalar @subexps == 0;
    my @taxa = Vector::uniq(sort( map {$_->sampleTaxId()} @exps) );
    die "Getting data for more than one taxon: " . (map {$_->sampleTaxId()} @exps) if scalar @taxa != 1;

    my @subexpids = sort {$a <=> $b} (map { $_->subExpId() } @subexps);
    print STDERR "Getting data for " . scalar @subexpids . " sub-experiments: "
	. join(",", @subexpids), "\n";

    # Now figure out the locuss
    my @locusIds = (); # if specifying specific locuss
    my %locusidToName = ();
    if ($locus && $locus ne "all") {
	my $genome = Genome::new('taxonomyId'=>$taxa[0]);
	my @names = split /,/, $locus;
	my $hash = $genome->getOrfIdsByName(@names);
	foreach my $name (@names) {
	    die "Unknown gene $name" if !exists $hash->{$name};
	    my $list = $hash->{$name};
	    die "locus name $name is not unique\n" if scalar @$list != 1;
	    push @locusIds, $list->[0];
	    $locusidToName{$list->[0]} = $name;
	}
    }

    # First fill the per-spot info
    my @allchipIds = map { $_->chipId() } @exps;
    my @chipIds = Vector::uniq(sort {$a <=> $b} @allchipIds);
    my @chips = map { Chip::new('chipId'=>$_) } @chipIds;
    my %spotInfo = (); # spot -> list
    my %spotAttrs = (); # spot -> list
    my %spotOrf = (); # spot -> output (comma-joined)
    my %chipToOrfSpots = ();
    my %chipToSpotList = (); # comma-separated; for selects if just looking at a few locuss

    foreach my $chip (@chips) {
	my $cid = $chip->chipId();
	if ($locus) {
	    $chipToOrfSpots{$cid} = $chip->getUniqueSpotsForORFs(@locusIds);
	    if (scalar @locusIds > 0) {
		my @spots = ();
		foreach my $list (values %{ $chipToOrfSpots{$cid} }) {
		    push @spots, @$list;
		}
		$chipToSpotList{$cid} = "spotId IN (" . join(",", @spots) . ")";
	    }
	} elsif ($spotSet) {
	    $chipToSpotList{$cid} = "spotId IN ($spotSet)";
	}
	$chipToSpotList{$cid} = "" if !exists $chipToSpotList{$cid};

	if (scalar @attrSpotInfo > 0) {
	    my $result = $chip->getSpotInfo(\@attrSpotInfo, $failedOK, $chipToSpotList{$cid});
	    foreach my $spot (keys %$result) {
		$spotInfo{$spot} = $result->{$spot};
	    }
	}
	if (scalar @attrOrf > 0) {
	    if ($locus) {
		foreach my $locus (keys %{ $chipToOrfSpots{$cid} }) {
		    map { $spotOrf{$_} = $locus } @{ $chipToOrfSpots{$cid}{$locus} };
		}
	    } else {
		my $result = $chip->getSpotORFs($chipToSpotList{$cid} ? "SpotInfo." . $chipToSpotList{$cid} : "");
		foreach my $spot (keys %$result) {
		    $spotOrf{$spot} = join(",", @{ $result->{$spot} });
		}
	    }
	}
	if (scalar @attrSpotAttr > 0) {
	    my @attrShort = @attrSpotAttr;
	    map s/^SpotAttr\.//, @attrShort;
	    my $result = $chip->getSpotAttrs(\@attrShort, $chipToSpotList{$cid});
	    # first set up the keys
	    foreach my $attr (keys %$result) {
		foreach my $spot (keys %{ $result->{$attr} }) {
		    $spotAttrs{$spot} = [];
		}
	    }
	    foreach my $spot (keys %spotAttrs) {
		my $list = $spotAttrs{$spot};
		foreach my $attr (@attrShort) {
		    my $value = "";
		    $value = $result->{$attr}{$spot} if exists $result->{$attr}{$spot};
		    push @$list, $value;
		}
	    }
	}
    }

    # Now get the per-subexperiment info
    my %data = (); # spot -> subexp -> list
    my %dataAttrs = (); # spot -> subexp -> list

    foreach my $subexp (@subexps) {
	my $subexpId = $subexp->subExpId();
	my $cid = Experiment::new('expId'=>$subexp->expId())->chipId();
	if (scalar @attrArrayData > 0) {
	    my $result = $subexp->getSpotData(\@attrArrayData,
					      $failedOK, $dyeCorrReq, $keepNulls,
					      $chipToSpotList{$cid} ? "ArrayData." .$chipToSpotList{$cid} : "" );
	    foreach my $spot (keys %$result) {
		$data{$spot} = {} if !exists $data{$spot};
		$data{$spot}{$subexpId} = $result->{$spot};
	    }
	}

	if (scalar @attrArrayDataAttr > 0) {
	    my @attrShort = @attrArrayDataAttr;
	    map s/^ArrayDataAttr\.//, @attrShort;
	    my $result = $subexp->getSpotAttrs(\@attrShort, $chipToSpotList{$cid});
	    foreach my $attr (keys %$result) {
		foreach my $spot (keys %{ $result->{$attr} }) {
		    $dataAttrs{$spot} = {} if !exists $dataAttrs{$spot};
		    $dataAttrs{$spot}{$subexpId} = [];
		}
	    }
	    foreach my $spot (keys %dataAttrs) {
		my $list = $dataAttrs{$spot}{$subexpId};
		foreach my $attr (@attrShort) {
		    my $value = "";
		    $value = $result->{$attr}{$spot} if exists $result->{$attr}{$spot};
		    push @$list, $value;
		}
	    }
	}
    }

    # Output -- only those spots in %data
    my %seToSpotToOut = ();

    my @attrLineIndices = 0..$#attrs;	# actual attributes per line - complicated if tabulating

    if ($tabulate) {
	$seToSpotToOut{0} = {}; # no subexp index

	my @attrIndicesFront = grep {within($attrTypes{$attrsNoFunc[$_]}, "SpotInfo", "SpotAttr", "locusId")} 0..$#attrs;
	my @attrIndicesSE = grep {!within($attrTypes{$attrsNoFunc[$_]}, "SpotInfo", "SpotAttr", "locusId")} 0..$#attrs;
	@attrLineIndices = @attrIndicesFront;
	map { push @attrLineIndices, @attrIndicesSE } @subexpids;

	my @printFront = map $attrs[$_], @attrIndicesFront;
	map s/^SpotAttr\.//, @printFront;
	my @printPerSE = map $attrs[$_], @attrIndicesSE;
	map s/^ArrayDataAttr\.//, @printPerSE;
	print "locusName\t" if scalar @locusIds > 0 && !$spots; # print names used to find genes
	print join("\t", @printFront);
	foreach my $subexpid (@subexpids) {
	    print "\t" unless $subexpid == $subexpids[0] && scalar @printFront == 0;
	    print join("\t", map { "$_.$subexpid" } @printPerSE);
	}
	print "\n";
	
	foreach my $spot (keys %data) {
	    # per-spot
	    my @out = ();
	    my @outSpotInfo = ();
	    my @outSpotAttrs = ();
	    my $outOrf = "";
	    
	    @outSpotInfo = @{ $spotInfo{$spot} } if exists $spotInfo{$spot};
	    @outSpotAttrs = @{ $spotAttrs{$spot} } if exists $spotAttrs{$spot};
	    $outOrf = $spotOrf{$spot} if exists $spotOrf{$spot};
	    
	    foreach my $attr (map $attrsNoFunc[$_],@attrIndicesFront) {
		if ($attrTypes{$attr} eq "SpotInfo") {
		    push @out, (scalar @outSpotInfo > 0 ? shift @outSpotInfo : "");
		} elsif ($attrTypes{$attr} eq "SpotAttr") {
		    push @out, (scalar @outSpotAttrs > 0 ? shift @outSpotAttrs : "");
		} elsif ($attrTypes{$attr} eq "locusId") {
		    push @out, $outOrf;
		}
	    }
	    # per-subexp
	    foreach my $subexpid (@subexpids) {
		my @outData = ();
		my @outDataAttrs = ();
		@outData = @{ $data{$spot}{$subexpid} } if exists $data{$spot}{$subexpid};
		@outDataAttrs = @{ $dataAttrs{$spot}{$subexpid} } if exists $dataAttrs{$spot}{$subexpid};
		foreach my $attr (map $attrsNoFunc[$_], @attrIndicesSE) {
		    if ($attrTypes{$attr} eq "ArrayData") {
			push @out, (scalar @outData > 0 ? shift @outData : "");
		    } elsif ($attrTypes{$attr} eq "ArrayDataAttr") {
			push @out, (scalar @outDataAttrs > 0 ? shift @outDataAttrs : "");
		    }
		}
	    }
	    $seToSpotToOut{0}{$spot} = \@out;
	}
    } else {
	my @printAttrs = @attrs;
	map { s/^SpotAttr\.//; s/^ArrayDataAttr\.//; } @printAttrs;
	print "locusName\t" if scalar @locusIds > 0 && !$spots;
	print join("\t",@printAttrs),"\n";
	foreach my $subexpid (@subexpids) {
	    $seToSpotToOut{$subexpid} = {};
	    foreach my $spot (keys %data) {
		if (exists $data{$spot}{$subexpid}) {
		    my @out = ();
		    my @outSpotInfo = ();
		    my @outSpotAttrs = ();
		    my $outOrf = "";
		    my @outData = @{ $data{$spot}{$subexpid} };
		    my @outDataAttrs = ();

		    @outSpotInfo = @{ $spotInfo{$spot} } if exists $spotInfo{$spot};
		    @outSpotAttrs = @{ $spotAttrs{$spot} } if exists $spotAttrs{$spot};
		    $outOrf = $spotOrf{$spot} if exists $spotOrf{$spot};
		    @outDataAttrs = @{ $dataAttrs{$spot}{$subexpid} }
		        if exists $dataAttrs{$spot}{$subexpid};
		    foreach my $attr (@attrsNoFunc) {
			if ($attrTypes{$attr} eq "SpotInfo") {
			    push @out, (scalar @outSpotInfo > 0 ? shift @outSpotInfo : "");
			} elsif ($attrTypes{$attr} eq "SpotAttr") {
			    push @out, (scalar @outSpotAttrs > 0 ? shift @outSpotAttrs : "");
			} elsif ($attrTypes{$attr} eq "locusId") {
			    push @out, $outOrf;
			} elsif ($attrTypes{$attr} eq "ArrayData") {
			    push @out, (shift @outData);
			} elsif ($attrTypes{$attr} eq "ArrayDataAttr") {
			    push @out, (scalar @outDataAttrs > 0 ? shift @outDataAttrs : "");
			}
		    }
		    $seToSpotToOut{$subexpid}{$spot} = \@out;
		}
	    }
	}
    }
    my @list = @subexpids;
    @list = (0) if $tabulate;

    my @showOrfs;
    my %locusToSpots = (); # locusid -> spotid -> 1
    if (scalar @locusIds > 0) {
	@showOrfs = @locusIds;
    } elsif ($locus && !$spots) {
	my %locusFound = ();
	# build the list from all chips
	foreach my $chip (keys %chipToOrfSpots) {
	    map { $locusFound{$_} = 1 } (keys %{ $chipToOrfSpots{$chip} });
	}
	@showOrfs = sort {$a <=> $b} keys %locusFound;
    }
    foreach my $locusid (@showOrfs) {
	$locusToSpots{$locusid} = {};
	foreach my $chip (keys %chipToOrfSpots) {
	    if (exists $chipToOrfSpots{$chip}{$locusid}) {
		map { $locusToSpots{$locusid}{$_} = 1 } @{ $chipToOrfSpots{$chip}{$locusid} };
	    }
	}
    }
    foreach my $key (@list) {
	if ($locus && !$spots) {
	    # group spots by locus & do summarization & report name if necessary
	    foreach my $locusid (@showOrfs) {
		my @spots = sort {$a <=> $b} keys %{ $locusToSpots{$locusid} };
		next if scalar @spots == 0;
		my @out = ();
		foreach my $i (0..$#attrLineIndices) {
		    my $iAttr = $attrLineIndices[$i];
		    my $attr = $attrs[$iAttr];
		    my $attrNoFunc = $attrsNoFunc[$iAttr];
		    my $func = $attrsFunc[$iAttr];
		    my $type = $attrTypes{$attrNoFunc};
		    my @values = map { $seToSpotToOut{$key}{$_}->[$i] } @spots;
		    my $value = "";
		    # Handle some ids manually in case all spots were missing
		    if ($attr eq "locusId") {
			$value = $locusid;
		    } elsif ($attr eq "spotId") {
			$value = join(",",@spots);
		    } elsif ($attr eq "subExpId" && !$tabulate) {
			$value = $key;
		    } elsif ($func eq "") {
			$value = join(",", map { defined($_) ? $_ : "" } @values);
		    } elsif ($func eq "mean") {
			my @def = grep {defined($_) && $_ ne ""} @values;
			$value = (scalar @def > 0) ? Vector::mean(@def) : "";
		    } elsif ($func eq "sum") {
			my @def = grep {defined($_) && $_ ne ""} @values;
			$value = (scalar @def > 0) ? Vector::sum(@def) : "";
		    } elsif ($func eq "count") {
			my @def = grep {defined($_) && $_ ne ""} @values;
			$value = scalar @def;
		    } elsif ($func eq "sd") {
			my @def = grep {defined($_) && $_ ne ""} @values;
			my $value = (scalar @def > 1) ? Vector::sd(@def) : "";
		    } elsif ($func eq "first") {
			my @def = grep {defined($_) && $_ ne ""} @values;
			$value = (scalar @def > 0) ? $def[0] : "";
		    }
		    push @out, $value;
		}
		print $locusidToName{$locusid}."\t" if exists $locusidToName{$locusid};
		print join("\t", @out),"\n";
	    }
	} else {
	    foreach my $spot (sort {$a <=> $b} keys %{$seToSpotToOut{$key}}) {
		print join("\t", map {defined($_) ? $_ : ""} @{ $seToSpotToOut{$key}{$spot}}),"\n";
	    }
	}
    }
}
