#!/usr/bin/perl -w

use Benchmark qw (:all);
use Statistics::Basic qw(vector correlation);
use Inline ( 'C', <<END_C );

double _correlation ( int size, char *inlist1, char *inlist2 ) {
    double *x  = (double *) inlist1;
	double *y  = (double *) inlist2;
    /* use the "arrays" just like normal for C code: */
    int i;
	
	/* compute sum for X and Y */
	double sum_x;
	double sum_y;
	double sum_x2;
	double sum_y2;
	double sum_xy;
    for( i= 0; i < size; ++i ) {
		sum_x += x[i];
		sum_y += y[i];
		sum_x2 += x[i]*x[i];
		sum_y2 += y[i]*y[i];
		sum_xy += x[i]*y[i];
    }
	
	double corr = ( (size*sum_xy) - (sum_x*sum_y) ) / ( sqrt ( (size*sum_x2 - sum_x*sum_x) * (size*sum_y2 - sum_y*sum_y) ) );
	return (corr);
}

END_C

sub corr_c {
    my( $l1, $l2 ) = @_;
    my $size = 0+@$l1;
    my $out = "\0" x ( $size * length pack "d", 0.0 );
	
    my $pack1 = pack "d$size", @$l1;
    my $pack2 = pack "d$size", @$l2;
    my $value = _correlation ( $size, $pack1, $pack2 );
	return ($value);
}

sub corr_bas ( $ $ ) {
	my ($_a, $_b) = @_;
	my $va = vector (@$_a);
	my $vb = vector (@$_b);
	my $cor = correlation($va,$vb);
	return ($cor);
}

sub corr_pdl ( $ $ ) {
	use PDL;
	my ($g1,$g2) = @_;
	my $p1 = pdl ( [@$g1] );
	my $p2 = pdl ( [@$g2] );
	my $X = cat $p1,$p2;
	my $mean = average ($X);
    my $Diff = $X - ($mean->transpose );
    my $Sigma = ( 1 / ( $X->getdim(0) - 1 ) ) * $Diff x transpose( $Diff );
	my $srt = $Sigma->flat;
	my @Value = list ($srt);
	my @T = grep { $_ < 0.999990 } @Value;
	my $R = shift (@T);	
    return $R;
}

my @a = qw (-0.0376361466042004	0.758017549087203	-1.62484700462194	0.803161581184831	0.321057899448435	0.726459611140626	-0.53600636747713	0.0947793616011517	-0.805247085947525	0.920175636968681	0.525199578056304	0.192137486831259	0.422047789696847	-0.328500649401554	0.745565058816148	0.320647462826354	0.899202003461104	-0.109414717253763	-0.70390808641615	0.328927586449054	0.627894949222558	0.348590749317348	0.235350408401976	0.193536895797966	0.191675074564699	0.679184155770257	1.08306449923048	-0.202693305028221	0.612333682431079	0.0887116487088852	0.37199684966577	0.181488076822368	0.161797370192285	0.156345372742766	0.402023366181265	-0.506739741037208	-0.518939419459781	0.343733844610533	0.202770249473228	0.663698999790544	0.180389813745114	1.1057158291283	0.572342167306971	-0.809528898597075	0.176317723754987	0.176357142407161	-0.801069565808923	0.479849408245944	0.783497299206675	0.457877750688472	-0.446369160298885	-1.42071780474148	0.320941321248691	-2.45202681487162	1.49833984673529	0.423977973660329	0.621303685805681	1.09026314733929	0.963671521695143	2.35215280098413	0.0939984688209615	-0.265936037471406	2.12062466617647	0.24066291075711	-0.948515210870928	-1.54202850225898	-0.411877835603484	1.14383111689151	1.68896709463923	1.46245464981868	0.350868553507604	0.369804719406304	0.709866297769061	0.14906150778658	-0.142406818576134	-0.544186625863647	-1.00568410271454	0.747259859392496	2.54506364971975	-0.506705853019587	3.57218365784964	0.361113340976292	-0.0551222653117865	0.819352336358948	-3.28502805900413	-0.930514239119135	-0.459041193379399	-1.13656017332393	0.615325787990152	0.0441904584155815	-0.886753008765268	-1.63958138211096	0.217021819177625	-0.0204092302940296	-0.946430834741571	0.260517818833869	-0.241027620008572	0.592570112258262	-0.324778383352236	-0.845870414378938	-0.695525804916582	-0.454637678698499	-1.05437208956519	-0.172909815072374	0.0649263202382453	-0.364239607880737	-1.00272524349104	-0.26155832222133	-0.473840302378365	0.847863264789135	-0.343213876548939	0.00354436451711844	-0.935585339230541	0.299739932872955	-0.30938396356084	-0.166407619629374	-0.611358420259907	-0.600669118256573	0.21186718176169	0.432690730644746	-0.409016413835568	0.139848096767294	0.00265693517864149	0.284059893927238	-0.122112322319606	0.255825077948685	-0.15616170400785	-0.535732626350526	-0.734949440835133	-0.0563624859191143	0.451449493476319	-0.43468889681053	-4.33121707400563	0.778254284359244	-3.879229271259	-0.476510036126708	0.0255000349818576	-1.5234275667778	0.840400929238361);
my @b = qw (-0.566862899236129	-0.362919907362601	-0.618921815792126	1.49084392126942	-0.194153310741768	-0.865235416095197	-0.131037216514284	0.238975041975721	-0.313114770133457	-0.96878299839573	0.708275856532125	-1.64915919176884	0.207747235677399	-2.19214384576881	-0.915314760498393	-0.818478427289636	-0.726638733553772	-0.382720108345335	-0.0115556049630281	-0.505146564616344	0.379028685879969	-0.109012005265629	-0.66082978597232	-0.829220485089294	0.24539772029878	0.208308130185455	-0.501746299388318	0.208507575558713	0.35883047233513	0.0536180514444205	-0.946270274886411	0.33219566380142	0.41713968232112	0.262430955696001	-1.09358195063564	-2.08059833954678	1.5028855181657	0.444720007869125	0.228601474561978	-0.134142729921479	-1.35358880112752	-0.785396827359728	0.744541163124476	-0.874158528844567	-0.00717483733466419	0.192687219255756	0.403249558201557	0.0792479188838008	0.476091619245819	-0.00229228728602043	0.412496687668963	0.494580760782353	-2.00715447001614	0.677225718476562	1.48848588365331	1.09556988840844	-0.548217420249446	-1.02968364163283	-0.884739562629946	4.10665701672223	3.34166869622101	-0.606170138612289	0.767381440066104	2.28280116898804	-0.0559556902445701	-1.13949257513291	-0.410310241774616	0.531480668688413	-0.486375423048787	2.02118920607841	1.31696296944845	1.13375923755163	0.83180807233043	1.66011905918714	-0.0439823599622661	-0.977913224150234	1.04503680098851	0.939853160478464	0.660071531696586	0.620162660183511	0.555789829696648	1.30097551675608	-1.28726252633491	0.306974332775711	2.06340357849559	1.32248285524116	0.0590831842717602	1.15619836034684	0.501569120399076	1.14832748046601	0.122311876088041	1.00985191967496	-0.7061542396621	-1.83285722269174	-0.770846974431533	0.410151052258053	-0.466981499561544	0.528248932802423	-2.12756443204757	-1.11240365073655	-0.865540824938905	-0.243998917113422	-1.78028235054911	0.615175207272723	-0.373165079766787	-0.329288051422455	1.44847702370934	-0.281650768195238	-0.43473528328612	0.570985708609602	0.00548798159692245	-0.358415736884061	0.528799141043767	-0.39420864780079	0.516283026982441	-0.283726535758159	-0.945465439132822	0.933813873899502	-0.355908943164826	-1.14403703912539	0.352545667157736	0.944707905422389	0.735830929849865	-0.102189714029639	-0.598336726262387	-1.24503408342729	-0.413820107348765	0.341301714445127	-0.131074680050781	-0.739487746857388	0.0523985695101994	-0.202219110444159	-0.295720017815558	0.371370934164971	-0.394445846803775	-1.60220549715511	0.330135241278079	-1.53542079166354	-0.698673108467241);


my $r_cov = timethese( -5, {
	C	=> sub {corr_c(\@a,\@b)},
	pdl	=> sub {corr_pdl(\@a,\@b)},
	bas	=> sub {corr_bas(\@a,\@b)},
} );
cmpthese $r_cov;

#my $output = corr_c (\@a,\@b);
#print "result = $output\n";