grid-wrap.pl

#!/usr/bin/env perl
use 5.16.0;
use strict;
use Sort::Key qw(nsort);
use rob_misc;
use binary_grid;
use threads;
use threads::shared;
use Fcntl ':flock'; # import LOCK_* constants
use subs qw(logprint logprintf);
use Carp qw(confess);
use File::Path qw(make_path);
$|=1;
############################################################

my $outdir='/tmp/gridtest'; 
my $nmachines=20;

say "Run across (modulo=) ",scalar(@machines)," machines";
my $modulo=scalar(@machines);

my @pids;
my @datafiles;
my $offset=0;
my @spectral_type_strings=('O','B','A','F','G','K','M');
my %spectral_type;  revspechash();
my @stellar_type_strings;
my @short_stellar_type_strings;
my %short_stellar_type_strings;
my $period_binwidth=0.1; # width in logP of histograms
my $log_fp;
defaults();

foreach my $machine (0..$nmachines)
{
    my $datafile="$outdir/results.$offset";
    my $cmd="./src/perl/scripts-flexigrid/grid-canberra.pl modulo=$modulo offset=$offset results_hash_dumpfile=$datafile";
    
    say "Run $cmd";

    my $forkr;

    # run the models : should outsource this to binary_grid::condor module
    if(1){     
	$forkr=fork();
	if($forkr==0)
	{
	    print `$cmd 2>$outdir/grid.$offset.err > $outdir/grid.$offset.out`;
	    exit;
	}
    }

    push(@pids,$forkr);
    push(@datafiles,$datafile);
    
    $offset++;
}

# wait for children to finish
say "Waiting for children (pids: @pids)";
map { waitpid $_,0; }@pids;

say "Children finished\n";


# put all the data together
my $results=merge_results_hash_dumps(@datafiles);

# save ouput
output();

if(0)
{
    use Data::Dumper;
    open($log_fp,">$outdir/datadump")||confess("cannot open $outdir/datadump"); 
    print {$log_fp} Dumper($results);
    close $log_fp;
}

say "Done\n";
exit;

############################################################


sub output
{
    # open log file
    
    # check for an output directory, otherwise assume a default in /tmp
    
    # output binary_grid::grid_options/bse_options
    open($log_fp,">$outdir/grid_options")||confess("cannot open $outdir/grid_options for output");
    print $log_fp "grid_options\n\n";
    map
    {
	print {$log_fp} $_.' : '.$binary_grid::grid_options{$_}."\n";
    }keys %binary_grid::grid_options;
    close $log_fp;
    open($log_fp,">$outdir/bse_options")||confess("cannot open $outdir/bse_options for output");
    print $log_fp "bse_options\n\n";
    map
    {
	print {$log_fp} $_.' : '.$binary_grid::bse_options{$_}."\n";
    }keys %binary_grid::bse_options;
    close $log_fp;
    
    # output grid results
    open ($log_fp,">$outdir/log") || confess ("cannot open $outdir/log for output");
    logprint "\n\n\nCANBERRA Grid results\n============================================================\n\n";
    logprint "SVN Revision ",`svn info|grep Revision|gawk \"\{print \\\$2\}\"`;
    logprint "Run at ",scalar(localtime),"\n";
    
    logprint "\n\nParameters:\n";
    logprint "Initial binary fraction : Raghavan 2010\n";
    foreach my $p ('alpha_ce','lambda_ce','lambda_ionisation','z','wd_sigma')
    {
	my $opt = ($p eq 'lambda_ce' && $binary_grid::bse_options{$p}==-1) 
	    ? 'Fitted to detailed stellar models' 
	    : 
	    sprintf "%g",$binary_grid::bse_options{$p};

	logprint sprintf "% -30s = %s\n",$p eq'z'?'Metallicity':$p,$opt;
    }

    
    my $dtp_wd=1e-30; # total number count of WDs (single or binary)
    my %dtp_wd;

    {
	my $k='number count';
	# get total WD number count
	foreach my $duplicity (keys %{$results->{$k}})
	{ 
	    foreach my $wd_type (keys %{$results->{$k}->{$duplicity}})
	    {
		$dtp_wd += $$results{$k}->{$duplicity}->{$wd_type};
		$dtp_wd{$duplicity} += $$results{$k}->{$duplicity}->{$wd_type};
	    } 
	}
    }

    my $dtp_wdms=0.0;
    {
	# get total WD-MS number count
	my $k='MSWD spectral types';
	foreach my $spectral_type (nsort keys %{$results->{$k}})
	{
	    $dtp_wdms += $$results{$k}{$spectral_type};
	}    
    }


    foreach my $k (sort keys %$results)
    {
	next if($k eq 'thread_number');
	if($k eq 'number count')
	{
	    logprint sprintf "WD binary fraction %3.2f %%\n",$dtp_wd{'Binary'}/$dtp_wd*100.0;
	    foreach my $duplicity (keys %{$results->{$k}})
	    {
		foreach my $wd_type (nsort keys %{$results->{$k}->{$duplicity}})
		{
		    logprint sprintf "% -6s % -30s %3.2f %%\n",
		    $duplicity,
		    $stellar_type_strings[$wd_type], 
		    $$results{$k}->{$duplicity}->{$wd_type}/$dtp_wd*100.0;
		} 
	    }
	}
	elsif($k eq 'companion types')
	{
	    my $test_tot=0.0;
	    logprint "\n\nWhite Dwarf - ? pairs (as a fraction of all White Dwarf Binaries)\n";
	    foreach my $wd1 (nsort keys %{$results->{$k}})
	    {
		foreach my $wd2 (nsort keys %{$results->{$k}->{$wd1}})
		{
		    logprint sprintf "% -40s - % -40s     % 5.2f %%\n",
		    $stellar_type_strings[$wd1],
		    $stellar_type_strings[$wd2],
		    $$results{$k}{$wd1}{$wd2}
		    /$dtp_wd{'Binary'}*100.0;
		}
	    }
	    logprint "\n";
	}
	elsif($k eq 'test')
	{
	    # test initial binary fraction
	    my $test_tot=0.0;
	    logprint "\n\nTest numbers\n\n";
	    foreach my $duplicity (keys %{$results->{$k}})
	    {
		$test_tot += $$results{$k}->{$duplicity};
	    }
	    foreach my $duplicity (nsort keys %{$results->{$k}})
	    {
		logprint sprintf "MS $duplicity %3.2f %%\n", 100.0*$$results{$k}->{$duplicity}/$test_tot;
	    }
	}
	elsif($k eq 'test dist spectypes')
	{
	    print "ZAMS Binary fraction by spectral type : ",keys %{$results->{$k}->{'Binary'}},"\n\n";
	    foreach my $spec_type (@spectral_type_strings)
	    {
		my $t=$$results{$k}->{'Single'}->{$spec_type}+$$results{$k}->{'Binary'}->{$spec_type};
		logprint sprintf "%s %5.2f %%\n",$spec_type,$$results{$k}->{'Binary'}->{$spec_type}/$t*100.0 if($t>0.0);
	    }
	}
	elsif($k eq 'test dist')
	{
	    # test initial distributions
	    foreach my $duplicity (keys %{$results->{$k}})
	    {
		foreach my $dist (keys %{$results->{$k}->{$duplicity}})
		{
		    my $f="$outdir/init_".$duplicity.'_'.$dist.'.dat';
		    logprint "Output distribution to $f\n";
		    open(FP,">$f")||confess("can't open $f");
		    foreach my $x (nsort keys %{$results->{$k}->{$duplicity}->{$dist}})
		    {
			printf FP "%g %g\n",$x,$$results{$k}{$duplicity}{$dist}{$x};
		    }
		    close FP;
		}
	    }
	}

	elsif($k eq 'MSWD spectral types')
	{
	    my $f="$outdir/MSWD_spectral_types.dat";
	    open(FP,'>'.$f)||confess("cannot open $f for output");
	    logprint "MSWD binaries: spectral types (see datafile $f)\n";
	    pad_datahash($results->{$k},1);
	    foreach my $spectral_type (nsort keys %{$results->{$k}})
	    {
		printf FP "%d %g\n",$spectral_type,$$results{$k}{$spectral_type}/$dtp_wdms;
		logprint sprintf "MSWD spectral type %d = %s = %g\n",
		$spectral_type,$spectral_type_strings[$spectral_type],$$results{$k}{$spectral_type}/$dtp_wdms;
	    } 
	    logprint "\n";
	    close FP;
	}
	# misc period distributions
	elsif($k =~/(\S+) period distribution/)
	{
	    my $distname=$1;
	    my $f="$outdir/$distname".'_period_distribution.dat';
	    #print "Pad datahash $results->{$k}\n";
	    pad_datahash($results->{$k},$period_binwidth);
	    open(FP,'>'.$f)||confess("cannot open $f for output");
	    logprint "$1 binaries: period distribution (see datafile $f)\n";
	    foreach my $logP (nsort keys %{$results->{$k}})
	    {
		printf FP "%g %g\n",$logP,$$results{$k}{$logP};
	    } 
	    close FP;
	}
    }

    # compare to Holberg 2009
    {
	my %obs=(
	    # Holberg 2009
	    'Single WDs'=>0.70, 
	    'WD +dM'=>0.12,
	    'Sirius-Like'=>0.08,
	    'WD + WD'=>0.10,
	    'Total WDs'=>1.0,
	    );
	map{$obs{$_}*=100.0}keys %obs;

	my $format="\% -30s \% 5.2f \%\% (% 5.2f)\n";
	logprint "\n=============\nHolberg 2009 comparison\n\n";

	# Single WD fraction
	my $single_WD = $dtp_wd{'Single'}/$dtp_wd*100.0;
	logprintf $format,'Single WDs',$single_WD,$obs{'Single WDs'};

	# WD + dM fraction
	my $WD_dM=$$results{'MSWD spectral types'}{$spectral_type{'M'}}/$dtp_wd*100.0;
	logprintf $format,'WD +dM',$WD_dM,$obs{'WD +dM'};
	
	# Sirius-like fraction
	my $Sirius=0.0;
	map{$Sirius+=$$results{'MSWD spectral types'}{$spectral_type{$_}}}('O','B','A','F','G','K');
	$Sirius = $Sirius/$dtp_wd * 100.0;
	logprintf $format,'Sirius-Like',$Sirius,$obs{'Sirius-Like'};

	# WD + WD fraction
	my $WDWD=0.0;
	foreach my $wd1 (10,11,12)
	{
	    foreach my $wd2 (10,11,12)
	    {
		$WDWD += $$results{'companion types'}{$wd1}{$wd2};
	    }
	}
	$WDWD = $WDWD/$dtp_wd*100.0;
        logprintf $format,'WD + WD',$WDWD,$obs{'WD + WD'};
	
	# total (should be 100%!)
	logprintf $format,'Total WDs',($single_WD + $WD_dM + $Sirius + $WDWD),$obs{'Total WDs'};
    }

    logprint "\n============================================================\n\n\n";

    close $log_fp;
}

sub logprint
{
    print {$log_fp} "@_";
    print "@_";
}
sub logprintf
{
    my $format=shift @_;
    printf {$log_fp} $format,@_;
    printf $format,@_;
}

sub revspechash
{
    for(my $i=0;$i<=$#spectral_type_strings;$i++)
    {
	$spectral_type{$spectral_type_strings[$i]}=$i;
    }
}

sub spec_to_float
{
    my $s=$_[0]; # spectral type to check
    if($s=~/([A-Z])(\d)/)
    {
	return $spectral_type{$1}+0.1*$2;
    }
    else
    {
	print "Failed to convert spec type $_[1] to float\n";
	exit;
    }
}

sub spec_type_sorter
{
    return spec_to_float($a) <=> spec_to_float($b);
}

sub pad_datahash
{
    # pad datahash $h so it has zeros where there's no data
    my $h=shift;
    my $d=shift;
    my @k=nsort keys %$h;
    my $xmax=$k[$#k]; # max value
    my $x=$k[0]; # start (min) value
    while($x<$xmax)
    {
	$$h{$x}=0.0 if(!defined($$h{$x}));
	$x=sprintf '%g',$x+$d;
    }
}

sub defaults
{
    grid_defaults();
    @stellar_type_strings=@{(stellar_type_strings())[2]};
    @short_stellar_type_strings=@{(stellar_type_strings())[3]};
    %short_stellar_type_strings=%{(stellar_type_strings())[4]};
}