mirror of
https://asciireactor.com/otho/psdlag-agn.git
synced 2024-11-21 21:35:07 +00:00
strongly modifying program to use the new clag
This commit is contained in:
caes
parent
142f013d80
commit
f3a88ca3fa
@ -3,3 +3,31 @@
|
|||||||
# This script is meant to encapsulate all functions in proper order to produce
|
# This script is meant to encapsulate all functions in proper order to produce
|
||||||
# the analyses. Parameters that this script should control: Δt, σ type, input
|
# the analyses. Parameters that this script should control: Δt, σ type, input
|
||||||
# dataset(lightcurve directory), more?
|
# dataset(lightcurve directory), more?
|
||||||
|
|
||||||
|
mkdir -p analyses/tables
|
||||||
|
|
||||||
|
#for analysis in analyses/*
|
||||||
|
for echo_curve in $lc_dir/*
|
||||||
|
do
|
||||||
|
# Grab and determine labels of analyses, skip if over the same band.
|
||||||
|
ref_band=$(basename $analysis|sed 's|\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}|\1|')
|
||||||
|
echo_band=$(basename $analysis|sed 's|[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}|\1|')
|
||||||
|
if [[ $ref_band == $echo_band ]]; then continue; fi
|
||||||
|
err_str=$(basename $analysis|sed 's|[^≺]*[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*;\(σ∊[CLM][MFC]\)}|\1|')
|
||||||
|
|
||||||
|
echo "Tabling PSD and time lags referred to ${ref_band} for $echo_band{${err_str}}."
|
||||||
|
|
||||||
|
# Propagate tables into analyses/tables
|
||||||
|
echoPSD_tabfile=analyses/tables/PSD_${echo_band}_\{${err_str}\}.tab
|
||||||
|
refPSD_tabfile=analyses/tables/PSD_${ref_band}_\{${err_str}\}.tab
|
||||||
|
timelag_tabfile=analyses/tables/timelag_${ref_band}_≺_${echo_band}_\{${err_str}\}.tab
|
||||||
|
|
||||||
|
# Output curves to temporary files using perl script, move tables to
|
||||||
|
# permanent location. This just assumes there are no conflicts.
|
||||||
|
scripts/extract_tables.pl $analysis > /dev/null
|
||||||
|
mv tmp.echoPSD $echoPSD_tabfile
|
||||||
|
mv tmp.refPSD $refPSD_tabfile
|
||||||
|
mv tmp.timelag $timelag_tabfile
|
||||||
|
done
|
||||||
|
|
||||||
|
rm tmp.*
|
||||||
@ -7,10 +7,6 @@ import getopt
|
|||||||
sys.path.insert(1,"/usr/local/science/clag/")
|
sys.path.insert(1,"/usr/local/science/clag/")
|
||||||
import clag
|
import clag
|
||||||
|
|
||||||
# psd1_output = open("psd1.tmp")
|
|
||||||
# psd2_output = open("psd2.tmp")
|
|
||||||
# lag_output = open("lag.tmp")
|
|
||||||
|
|
||||||
# For jupyter notebook
|
# For jupyter notebook
|
||||||
# %pylab inline
|
# %pylab inline
|
||||||
|
|
||||||
@ -76,7 +72,7 @@ P1.logLikelihood(inpars)
|
|||||||
|
|
||||||
|
|
||||||
## Now do the fitting and find the best fit psd values at the given frequency bins ##
|
## Now do the fitting and find the best fit psd values at the given frequency bins ##
|
||||||
psd1, psd1e = clag.optimize(P1, inpars)
|
ref_psd, ref_psd_err = clag.optimize(P1, inpars)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -85,7 +81,7 @@ psd1, psd1e = clag.optimize(P1, inpars)
|
|||||||
fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.)
|
fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.)
|
||||||
|
|
||||||
#loglog(fqd, 0.1*fqd**(-1.5), label='input psd')
|
#loglog(fqd, 0.1*fqd**(-1.5), label='input psd')
|
||||||
#errorbar(fqd[1:-1], psd1[1:-1], yerr=psd1e[1:-1], fmt='o', ms=10, label='fit')
|
#errorbar(fqd[1:-1], ref_psd[1:-1], yerr=ref_psd_err[1:-1], fmt='o', ms=10, label='fit')
|
||||||
|
|
||||||
# load second lightcurve
|
# load second lightcurve
|
||||||
|
|
||||||
@ -110,7 +106,7 @@ lc2_strength_err = [lc2_table[i, 2] for i in index]
|
|||||||
|
|
||||||
P2 = clag.clag('psd', lc2_time, lc2_strength, lc2_strength_err, dt, fqL)
|
P2 = clag.clag('psd', lc2_time, lc2_strength, lc2_strength_err, dt, fqL)
|
||||||
|
|
||||||
psd2, psd2e = clag.optimize(P2, inpars)
|
echo_psd, echo_psd_err = clag.optimize(P2, inpars)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -123,26 +119,18 @@ Cx = clag.clag('cxd',
|
|||||||
[list(i) for i in zip(lc1_time,lc1_time)],
|
[list(i) for i in zip(lc1_time,lc1_time)],
|
||||||
[list(i) for i in zip(lc1_strength,lc2_strength)],
|
[list(i) for i in zip(lc1_strength,lc2_strength)],
|
||||||
[list(i) for i in zip(lc1_strength_err,lc2_strength_err)],
|
[list(i) for i in zip(lc1_strength_err,lc2_strength_err)],
|
||||||
dt, fqL, psd1, psd2)
|
dt, fqL, ref_psd, echo_psd)
|
||||||
|
|
||||||
inpars = np.concatenate( (0.3*(psd1*psd2)**0.5, psd1*0+1.) )
|
inpars = np.concatenate( (0.3*(ref_psd*echo_psd)**0.5, ref_psd*0+1.) )
|
||||||
p, pe = clag.optimize(Cx, inpars)
|
p, pe = clag.optimize(Cx, inpars)
|
||||||
phi, phie = p[nfq:], pe[nfq:]
|
phi, phie = p[nfq:], pe[nfq:]
|
||||||
lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd)
|
lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd)
|
||||||
cx, cxe = p[:nfq], pe[:nfq]
|
cx, cxe = p[:nfq], pe[:nfq]
|
||||||
|
|
||||||
np.savetxt("lag.out",lag)
|
np.savetxt("freq.out",fqL.reshape((-1,len(fqL))))
|
||||||
|
np.savetxt("ref_psd.out",[ref_psd,ref_psd_err])
|
||||||
|
np.savetxt("echo_psd.out",[echo_psd,echo_psd_err])
|
||||||
|
np.savetxt("crsspctrm.out",[cx,cxe])
|
||||||
|
np.savetxt("timelag.out",[lag,lage])
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
## plot ##
|
|
||||||
|
|
||||||
#semilogx(fqd, fqd*0+1.0, label='input phase lag')
|
|
||||||
#ylim([0.8, 1.2])
|
|
||||||
#errorbar(fqd[1:-1], phi[1:-1], yerr=phie[1:-1], fmt='o', ms=10, label='fit')
|
|
||||||
|
|
||||||
|
|||||||
@ -152,7 +152,8 @@ do
|
|||||||
sed "s@%ROOTDIR%@$(pwd)@g"|
|
sed "s@%ROOTDIR%@$(pwd)@g"|
|
||||||
sed "s@%OUTPUTFILE%@${outputfile_noUTF}@g" > $submitscript
|
sed "s@%OUTPUTFILE%@${outputfile_noUTF}@g" > $submitscript
|
||||||
cd thor
|
cd thor
|
||||||
qsub $(basename $submitscript) >> submissions
|
# qsub $(basename $submitscript) >> submissions
|
||||||
|
echo CURRENTLY NOT IMPLEMENTED!!!
|
||||||
cd ..
|
cd ..
|
||||||
fi
|
fi
|
||||||
else
|
else
|
||||||
@ -160,7 +161,10 @@ do
|
|||||||
then
|
then
|
||||||
echo -n " Analysis already exists; skipping."
|
echo -n " Analysis already exists; skipping."
|
||||||
else
|
else
|
||||||
time bin/psdlag tmp.psdlagargs > $outputfile
|
#time bin/psdlag tmp.psdlagargs > $outputfile
|
||||||
|
time python scripts/analyze_lightcurve.py "$ref_curve $echo_curve"
|
||||||
|
scripts/extract_tables.pl
|
||||||
|
|
||||||
fi
|
fi
|
||||||
fi
|
fi
|
||||||
echo ""
|
echo ""
|
||||||
|
|||||||
@ -23,54 +23,40 @@ if (!${^UTF8LOCALE}) {
|
|||||||
say encode($charset,"You are not using UTF-8 encoding. :(");
|
say encode($charset,"You are not using UTF-8 encoding. :(");
|
||||||
}
|
}
|
||||||
my $charset=$ENV{LANG};
|
my $charset=$ENV{LANG};
|
||||||
our $outputfilename=decode($charset,$ARGV[0]);
|
|
||||||
open $outputfile,'<',$outputfilename or die $!;
|
|
||||||
|
|
||||||
my $star_linenum_1=0;
|
|
||||||
my $star_linenum_2=0;
|
|
||||||
my $star_bytenum_1=0;
|
|
||||||
my $star_bytenum_2=0;
|
|
||||||
|
|
||||||
my $linenum=0;
|
# This program attempts to fill function_bin with the tabulated PSDs
|
||||||
while(<$outputfile>) {
|
# and time lags.
|
||||||
$linenum++;
|
our %function_bin = ();
|
||||||
if ($_ =~ /^\*+$/) {
|
|
||||||
$star_linenum_1=$star_linenum_2;
|
open $freq_file,'<',"freq.out" or die $!;
|
||||||
$star_linenum_2 = $linenum;
|
open $ref_psd_file,'<',"ref_psd.out" or die $!;
|
||||||
$star_bytenum_1=$star_bytenum_2;
|
open $echo_psd_file,'<',"echo_psd.out" or die $!;
|
||||||
#$star_bytenum_2 = $outputfile.tell();
|
open $crsspctrm_file,'<',"crsspctrm.out" or die $!;
|
||||||
$star_bytenum_2 = tell($outputfile);
|
open $timelag_file,'<',"timelag.out" or die $!;
|
||||||
}
|
|
||||||
}
|
=pod
|
||||||
if ($debug) {
|
|
||||||
say encode($charset,"Final set found between lines "),
|
This section collects the various quantities.
|
||||||
$star_linenum_1,
|
|
||||||
" and ",
|
=cut
|
||||||
$star_linenum_2;
|
|
||||||
|
@bin_bounds = split /\s/,<$freq_file>;
|
||||||
|
@ref_psd = split /\s/,<$ref_psd_file>;
|
||||||
|
@ref_psd_σ = split /\s/,<$ref_psd_file>;
|
||||||
|
@echo_psd = split /\s/,<$echo_psd_file>;
|
||||||
|
@echo_psd_σ = split /\s/,<$echo_psd_file>;
|
||||||
|
@crsspctrm_psd = split /\s/,<$crsspctrm_file>;
|
||||||
|
@crsspctrm_psd_σ = split /\s/,<$crsspctrm_file>;
|
||||||
|
@timelag = split /\s/,<$timelag_file>;
|
||||||
|
@timelag_σ = split /\s/,<$timelag_file>;
|
||||||
|
|
||||||
say encode($charset,"Final set found between bytes "),
|
|
||||||
$star_bytenum_1,
|
|
||||||
" and ",
|
|
||||||
$star_bytenum_2;
|
|
||||||
}
|
|
||||||
seek($outputfile,$star_bytenum_1,0);
|
|
||||||
<$outputfile>;
|
|
||||||
$bin_bounds_line = <$outputfile>;
|
|
||||||
#say $bin_bounds_line;
|
|
||||||
$bin_bounds_line =~ s/^#\s*(.*)$/$1/;
|
|
||||||
#$line =~ /
|
|
||||||
@bin_bounds = split /\s/,$bin_bounds_line;
|
|
||||||
if ($debug) {
|
if ($debug) {
|
||||||
print encode($charset,"Found bin boundaries: ");
|
print encode($charset,"Found bin boundaries: ");
|
||||||
foreach (@bin_bounds) {print encode($charset,"$_ ");}
|
foreach (@bin_bounds) {print encode($charset,"$_ ");}
|
||||||
say encode($charset," ");
|
say encode($charset," ");
|
||||||
}
|
}
|
||||||
|
my $count = 0;
|
||||||
|
|
||||||
|
|
||||||
# This captures the frequency bins.
|
|
||||||
our %function_bin = ();
|
|
||||||
|
|
||||||
my $upper_bound = 0;
|
my $upper_bound = 0;
|
||||||
my $lower_bound = 0;
|
my $lower_bound = 0;
|
||||||
foreach (@bin_bounds) {
|
foreach (@bin_bounds) {
|
||||||
@ -81,10 +67,29 @@ foreach (@bin_bounds) {
|
|||||||
my $Δ = ($upper_bound - $lower_bound)/2;
|
my $Δ = ($upper_bound - $lower_bound)/2;
|
||||||
#say ($μ,":",$Δ);
|
#say ($μ,":",$Δ);
|
||||||
# push(@freq_coords_mean,$μ);
|
# push(@freq_coords_mean,$μ);
|
||||||
# push(@freq_coords_err,$Δ);
|
# push(@freq_coords_σ,$Δ);
|
||||||
$function_bin{$μ} = {"Δ" => $Δ};
|
$function_bin{$μ} = {"Δ" => $Δ};
|
||||||
|
$function_bin{$μ}{"ref_PSD_μ"} = $ref_psd[$count];
|
||||||
|
$function_bin{$μ}{"ref_PSD_σ"} = $ref_psd_σ[$count];
|
||||||
|
$function_bin{$μ}{"echo_PSD_μ"} = $echo_psd[$count];
|
||||||
|
$function_bin{$μ}{"echo_PSD_σ"} = $echo_psd_σ[$count];
|
||||||
|
$function_bin{$μ}{"crsspctrm_μ"} = $crosssp_psd[$count];
|
||||||
|
$function_bin{$μ}{"crsspctrm_σ"} = $crosssp_psd_σ[$count];
|
||||||
|
# $function_bin{$μ}{"φdiff_μ"} = $μ;
|
||||||
|
# $function_bin{$μ}{"φdiff_σ"} = $σ;
|
||||||
|
# $μ = $μ/(2*PI*$μ);
|
||||||
|
# $σ = $σ/(2*PI*$μ);
|
||||||
|
$function_bin{$μ}{"timelag_μ"} = $timelag[$count];
|
||||||
|
$function_bin{$μ}{"timelag_σ"} = $timelag_σ[$count];
|
||||||
|
$count = $count + 1;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
close $freq_file;
|
||||||
|
close $ref_psd_file;
|
||||||
|
close $echo_psd_file;
|
||||||
|
close $crsspctrm_file;
|
||||||
|
close $timelag_file;
|
||||||
|
|
||||||
$numbins = keys %function_bin;
|
$numbins = keys %function_bin;
|
||||||
say encode($charset,"$numbins frequency bins captured in output.");
|
say encode($charset,"$numbins frequency bins captured in output.");
|
||||||
|
|
||||||
@ -97,140 +102,6 @@ if($verbose) {
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
=pod
|
|
||||||
|
|
||||||
This section collects the various quantities. The mode counter
|
|
||||||
increments to designate the data being captured.
|
|
||||||
|
|
||||||
=cut
|
|
||||||
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,"");
|
|
||||||
say encode($charset,
|
|
||||||
" New Curve -- Reference Curve PSD");
|
|
||||||
say encode($charset,
|
|
||||||
"──────────────────────────────────────────────────");
|
|
||||||
}
|
|
||||||
|
|
||||||
foreach ( sort { $a <=> $b } keys %function_bin ) {
|
|
||||||
my $μ = my $σ = <$outputfile>;
|
|
||||||
$μ =~ s/^([\-\+e0-9\.]+)+\s+[\-\+e0-9\.]+\s*$/$1/;
|
|
||||||
$σ =~ s/^[\-\+e0-9\.]+\s+([\-\+e0-9\.]+)\s*$/$1/;
|
|
||||||
# if ($μ<0) {
|
|
||||||
# $μ = abs($μ);
|
|
||||||
# $function_bin{$_}{"φdiff_μ"} = PI;
|
|
||||||
# }
|
|
||||||
$function_bin{$_}{"ref_PSD_μ"} = $μ;
|
|
||||||
$function_bin{$_}{"ref_PSD_σ"} = $σ;
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,
|
|
||||||
"freq = " .
|
|
||||||
sprintf('%.3f',$_) .
|
|
||||||
": μ = " .
|
|
||||||
sprintf('%10.3e',$μ) .
|
|
||||||
"; σ = " .
|
|
||||||
sprintf('%10.3e',$σ));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,"");
|
|
||||||
say encode($charset,
|
|
||||||
" New Curve -- Reverberating Curve PSD");
|
|
||||||
say encode($charset,
|
|
||||||
"──────────────────────────────────────────────────");
|
|
||||||
}
|
|
||||||
foreach ( sort { $a <=> $b } keys %function_bin ) {
|
|
||||||
my $μ = my $σ = <$outputfile>;
|
|
||||||
$μ =~ s/^([\-\+e0-9\.]+)+\s+[\-\+e0-9\.]+\s*$/$1/;
|
|
||||||
$σ =~ s/^[\-\+e0-9\.]+\s+([\-\+e0-9\.]+)\s*$/$1/;
|
|
||||||
# if ($μ<0) {
|
|
||||||
# $μ = abs($μ);
|
|
||||||
# $function_bin{$_}{"φdiff_μ"} = PI;
|
|
||||||
# }
|
|
||||||
$function_bin{$_}{"echo_PSD_μ"} = $μ;
|
|
||||||
$function_bin{$_}{"echo_PSD_σ"} = $σ;
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,
|
|
||||||
"freq = " .
|
|
||||||
sprintf('%.3f',$_) .
|
|
||||||
": μ = " .
|
|
||||||
sprintf('%10.3e',$μ) .
|
|
||||||
"; σ = " .
|
|
||||||
sprintf('%10.3e',$σ));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,"");
|
|
||||||
say encode($charset,
|
|
||||||
" New Curve -- Cross Spectrum PSD");
|
|
||||||
say encode($charset,
|
|
||||||
"──────────────────────────────────────────────────");
|
|
||||||
}
|
|
||||||
foreach ( sort { $a <=> $b } keys %function_bin ) {
|
|
||||||
my $μ = my $σ = <$outputfile>;
|
|
||||||
$μ =~ s/^([\-\+e0-9\.]+)+\s+[\-\+e0-9\.]+\s*$/$1/;
|
|
||||||
$σ =~ s/^[\-\+e0-9\.]+\s+([\-\+e0-9\.]+)\s*$/$1/;
|
|
||||||
# if ($μ<0) {
|
|
||||||
# $μ = abs($μ);
|
|
||||||
# $function_bin{$_}{"φdiff_μ"} = PI;
|
|
||||||
# }
|
|
||||||
$function_bin{$_}{"crsspctrm_PSD_μ"} = $μ;
|
|
||||||
$function_bin{$_}{"crsspctrm_PSD_σ"} = $σ;
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,
|
|
||||||
"freq = " .
|
|
||||||
sprintf('%.3f',$_) .
|
|
||||||
": μ = " .
|
|
||||||
sprintf('%10.3e',$μ) .
|
|
||||||
"; σ = " .
|
|
||||||
sprintf('%10.3e',$σ));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,"");
|
|
||||||
say encode($charset,
|
|
||||||
" New Curve -- Phase Difference φ");
|
|
||||||
say encode($charset,
|
|
||||||
"──────────────────────────────────────────────────");
|
|
||||||
}
|
|
||||||
foreach ( sort { $a <=> $b } keys %function_bin ) {
|
|
||||||
my $μ = my $σ = <$outputfile>;
|
|
||||||
$μ =~ s/^([\-\+e0-9\.]+)+\s+[\-\+e0-9\.]+\s*$/$1/;
|
|
||||||
$σ =~ s/^[\-\+e0-9\.]+\s+([\-\+e0-9\.]+)\s*$/$1/;
|
|
||||||
if (exists $function_bin{$_}{"φdiff_μ"}) {
|
|
||||||
$μ = $function_bin{$_}{"φdiff_μ"} + $μ;
|
|
||||||
}
|
|
||||||
$function_bin{$_}{"φdiff_μ"} = $μ;
|
|
||||||
$function_bin{$_}{"φdiff_σ"} = $σ;
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,
|
|
||||||
"freq = " .
|
|
||||||
sprintf('%.3f',$_) .
|
|
||||||
": μ = " .
|
|
||||||
sprintf('%10.3e',$μ) .
|
|
||||||
"; σ = " .
|
|
||||||
sprintf('%10.3e',$σ));
|
|
||||||
}
|
|
||||||
$μ = $μ/(2*PI*$_);
|
|
||||||
$σ = $σ/(2*PI*$_);
|
|
||||||
$function_bin{$_}{"timelag_μ"} = $μ;
|
|
||||||
$function_bin{$_}{"timelag_σ"} = $σ;
|
|
||||||
if ($verbose) {
|
|
||||||
say encode($charset,
|
|
||||||
" timelag: μ = " .
|
|
||||||
sprintf('%10.3e',$μ) .
|
|
||||||
"; σ = " .
|
|
||||||
sprintf('%10.3e',$σ)) .
|
|
||||||
"\n";
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
say encode($charset,"");
|
say encode($charset,"");
|
||||||
|
|
||||||
if($debug) {
|
if($debug) {
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -1,769 +0,0 @@
|
|||||||
/*************************************************************************
|
|
||||||
Copyright (c) Sergey Bochkanov (ALGLIB project).
|
|
||||||
|
|
||||||
>>> SOURCE LICENSE >>>
|
|
||||||
This program is free software; you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation (www.fsf.org); either version 2 of the
|
|
||||||
License, or (at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
A copy of the GNU General Public License is available at
|
|
||||||
http://www.fsf.org/licensing/licenses
|
|
||||||
>>> END OF LICENSE >>>
|
|
||||||
*************************************************************************/
|
|
||||||
#ifndef _alglibmisc_pkg_h
|
|
||||||
#define _alglibmisc_pkg_h
|
|
||||||
#include "ap.h"
|
|
||||||
#include "alglibinternal.h"
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
ae_int_t s1;
|
|
||||||
ae_int_t s2;
|
|
||||||
ae_int_t magicv;
|
|
||||||
} hqrndstate;
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
ae_int_t n;
|
|
||||||
ae_int_t nx;
|
|
||||||
ae_int_t ny;
|
|
||||||
ae_int_t normtype;
|
|
||||||
ae_matrix xy;
|
|
||||||
ae_vector tags;
|
|
||||||
ae_vector boxmin;
|
|
||||||
ae_vector boxmax;
|
|
||||||
ae_vector nodes;
|
|
||||||
ae_vector splits;
|
|
||||||
ae_vector x;
|
|
||||||
ae_int_t kneeded;
|
|
||||||
double rneeded;
|
|
||||||
ae_bool selfmatch;
|
|
||||||
double approxf;
|
|
||||||
ae_int_t kcur;
|
|
||||||
ae_vector idx;
|
|
||||||
ae_vector r;
|
|
||||||
ae_vector buf;
|
|
||||||
ae_vector curboxmin;
|
|
||||||
ae_vector curboxmax;
|
|
||||||
double curdist;
|
|
||||||
ae_int_t debugcounter;
|
|
||||||
} kdtree;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS C++ INTERFACE
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib
|
|
||||||
{
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Portable high quality random number generator state.
|
|
||||||
Initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
Fields:
|
|
||||||
S1, S2 - seed values
|
|
||||||
V - precomputed value
|
|
||||||
MagicV - 'magic' value used to determine whether State structure
|
|
||||||
was correctly initialized.
|
|
||||||
*************************************************************************/
|
|
||||||
class _hqrndstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_hqrndstate_owner();
|
|
||||||
_hqrndstate_owner(const _hqrndstate_owner &rhs);
|
|
||||||
_hqrndstate_owner& operator=(const _hqrndstate_owner &rhs);
|
|
||||||
virtual ~_hqrndstate_owner();
|
|
||||||
alglib_impl::hqrndstate* c_ptr();
|
|
||||||
alglib_impl::hqrndstate* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::hqrndstate *p_struct;
|
|
||||||
};
|
|
||||||
class hqrndstate : public _hqrndstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
hqrndstate();
|
|
||||||
hqrndstate(const hqrndstate &rhs);
|
|
||||||
hqrndstate& operator=(const hqrndstate &rhs);
|
|
||||||
virtual ~hqrndstate();
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
|
|
||||||
*************************************************************************/
|
|
||||||
class _kdtree_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_kdtree_owner();
|
|
||||||
_kdtree_owner(const _kdtree_owner &rhs);
|
|
||||||
_kdtree_owner& operator=(const _kdtree_owner &rhs);
|
|
||||||
virtual ~_kdtree_owner();
|
|
||||||
alglib_impl::kdtree* c_ptr();
|
|
||||||
alglib_impl::kdtree* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::kdtree *p_struct;
|
|
||||||
};
|
|
||||||
class kdtree : public _kdtree_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
kdtree();
|
|
||||||
kdtree(const kdtree &rhs);
|
|
||||||
kdtree& operator=(const kdtree &rhs);
|
|
||||||
virtual ~kdtree();
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
HQRNDState initialization with random values which come from standard
|
|
||||||
RNG.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void hqrndrandomize(hqrndstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
HQRNDState initialization with seed values
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void hqrndseed(const ae_int_t s1, const ae_int_t s2, hqrndstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function generates random real number in (0,1),
|
|
||||||
not including interval boundaries
|
|
||||||
|
|
||||||
State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
double hqrnduniformr(const hqrndstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function generates random integer number in [0, N)
|
|
||||||
|
|
||||||
1. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()
|
|
||||||
2. N can be any positive number except for very large numbers:
|
|
||||||
* close to 2^31 on 32-bit systems
|
|
||||||
* close to 2^62 on 64-bit systems
|
|
||||||
An exception will be generated if N is too large.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
ae_int_t hqrnduniformi(const hqrndstate &state, const ae_int_t n);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Random number generator: normal numbers
|
|
||||||
|
|
||||||
This function generates one random number from normal distribution.
|
|
||||||
Its performance is equal to that of HQRNDNormal2()
|
|
||||||
|
|
||||||
State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
double hqrndnormal(const hqrndstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Random number generator: random X and Y such that X^2+Y^2=1
|
|
||||||
|
|
||||||
State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void hqrndunit2(const hqrndstate &state, double &x, double &y);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Random number generator: normal numbers
|
|
||||||
|
|
||||||
This function generates two independent random numbers from normal
|
|
||||||
distribution. Its performance is equal to that of HQRNDNormal()
|
|
||||||
|
|
||||||
State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 02.12.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void hqrndnormal2(const hqrndstate &state, double &x1, double &x2);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Random number generator: exponential distribution
|
|
||||||
|
|
||||||
State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 11.08.2007 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
double hqrndexponential(const hqrndstate &state, const double lambdav);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function generates random number from discrete distribution given by
|
|
||||||
finite sample X.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
State - high quality random number generator, must be
|
|
||||||
initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
X - finite sample
|
|
||||||
N - number of elements to use, N>=1
|
|
||||||
|
|
||||||
RESULT
|
|
||||||
this function returns one of the X[i] for random i=0..N-1
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 08.11.2011 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
double hqrnddiscrete(const hqrndstate &state, const real_1d_array &x, const ae_int_t n);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function generates random number from continuous distribution given
|
|
||||||
by finite sample X.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
State - high quality random number generator, must be
|
|
||||||
initialized with HQRNDRandomize() or HQRNDSeed().
|
|
||||||
X - finite sample, array[N] (can be larger, in this case only
|
|
||||||
leading N elements are used). THIS ARRAY MUST BE SORTED BY
|
|
||||||
ASCENDING.
|
|
||||||
N - number of elements to use, N>=1
|
|
||||||
|
|
||||||
RESULT
|
|
||||||
this function returns random number from continuous distribution which
|
|
||||||
tries to approximate X as mush as possible. min(X)<=Result<=max(X).
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 08.11.2011 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
double hqrndcontinuous(const hqrndstate &state, const real_1d_array &x, const ae_int_t n);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function serializes data structure to string.
|
|
||||||
|
|
||||||
Important properties of s_out:
|
|
||||||
* it contains alphanumeric characters, dots, underscores, minus signs
|
|
||||||
* these symbols are grouped into words, which are separated by spaces
|
|
||||||
and Windows-style (CR+LF) newlines
|
|
||||||
* although serializer uses spaces and CR+LF as separators, you can
|
|
||||||
replace any separator character by arbitrary combination of spaces,
|
|
||||||
tabs, Windows or Unix newlines. It allows flexible reformatting of
|
|
||||||
the string in case you want to include it into text or XML file.
|
|
||||||
But you should not insert separators into the middle of the "words"
|
|
||||||
nor you should change case of letters.
|
|
||||||
* s_out can be freely moved between 32-bit and 64-bit systems, little
|
|
||||||
and big endian machines, and so on. You can serialize structure on
|
|
||||||
32-bit machine and unserialize it on 64-bit one (or vice versa), or
|
|
||||||
serialize it on SPARC and unserialize on x86. You can also
|
|
||||||
serialize it in C++ version of ALGLIB and unserialize in C# one,
|
|
||||||
and vice versa.
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreeserialize(kdtree &obj, std::string &s_out);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function unserializes data structure from string.
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreeunserialize(std::string &s_in, kdtree &obj);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
KD-tree creation
|
|
||||||
|
|
||||||
This subroutine creates KD-tree from set of X-values and optional Y-values
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
XY - dataset, array[0..N-1,0..NX+NY-1].
|
|
||||||
one row corresponds to one point.
|
|
||||||
first NX columns contain X-values, next NY (NY may be zero)
|
|
||||||
columns may contain associated Y-values
|
|
||||||
N - number of points, N>=0.
|
|
||||||
NX - space dimension, NX>=1.
|
|
||||||
NY - number of optional Y-values, NY>=0.
|
|
||||||
NormType- norm type:
|
|
||||||
* 0 denotes infinity-norm
|
|
||||||
* 1 denotes 1-norm
|
|
||||||
* 2 denotes 2-norm (Euclidean norm)
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
|
|
||||||
1. KD-tree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory
|
|
||||||
requirements.
|
|
||||||
2. Although KD-trees may be used with any combination of N and NX, they
|
|
||||||
are more efficient than brute-force search only when N >> 4^NX. So they
|
|
||||||
are most useful in low-dimensional tasks (NX=2, NX=3). NX=1 is another
|
|
||||||
inefficient case, because simple binary search (without additional
|
|
||||||
structures) is much more efficient in such tasks than KD-trees.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreebuild(const real_2d_array &xy, const ae_int_t n, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
|
|
||||||
void kdtreebuild(const real_2d_array &xy, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
KD-tree creation
|
|
||||||
|
|
||||||
This subroutine creates KD-tree from set of X-values, integer tags and
|
|
||||||
optional Y-values
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
XY - dataset, array[0..N-1,0..NX+NY-1].
|
|
||||||
one row corresponds to one point.
|
|
||||||
first NX columns contain X-values, next NY (NY may be zero)
|
|
||||||
columns may contain associated Y-values
|
|
||||||
Tags - tags, array[0..N-1], contains integer tags associated
|
|
||||||
with points.
|
|
||||||
N - number of points, N>=0
|
|
||||||
NX - space dimension, NX>=1.
|
|
||||||
NY - number of optional Y-values, NY>=0.
|
|
||||||
NormType- norm type:
|
|
||||||
* 0 denotes infinity-norm
|
|
||||||
* 1 denotes 1-norm
|
|
||||||
* 2 denotes 2-norm (Euclidean norm)
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
|
|
||||||
1. KD-tree creation have O(N*logN) complexity and O(N*(2*NX+NY)) memory
|
|
||||||
requirements.
|
|
||||||
2. Although KD-trees may be used with any combination of N and NX, they
|
|
||||||
are more efficient than brute-force search only when N >> 4^NX. So they
|
|
||||||
are most useful in low-dimensional tasks (NX=2, NX=3). NX=1 is another
|
|
||||||
inefficient case, because simple binary search (without additional
|
|
||||||
structures) is much more efficient in such tasks than KD-trees.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreebuildtagged(const real_2d_array &xy, const integer_1d_array &tags, const ae_int_t n, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
|
|
||||||
void kdtreebuildtagged(const real_2d_array &xy, const integer_1d_array &tags, const ae_int_t nx, const ae_int_t ny, const ae_int_t normtype, kdtree &kdt);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
K-NN query: K nearest neighbors
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
X - point, array[0..NX-1].
|
|
||||||
K - number of neighbors to return, K>=1
|
|
||||||
SelfMatch - whether self-matches are allowed:
|
|
||||||
* if True, nearest neighbor may be the point itself
|
|
||||||
(if it exists in original dataset)
|
|
||||||
* if False, then only points with non-zero distance
|
|
||||||
are returned
|
|
||||||
* if not given, considered True
|
|
||||||
|
|
||||||
RESULT
|
|
||||||
number of actual neighbors found (either K or N, if K>N).
|
|
||||||
|
|
||||||
This subroutine performs query and stores its result in the internal
|
|
||||||
structures of the KD-tree. You can use following subroutines to obtain
|
|
||||||
these results:
|
|
||||||
* KDTreeQueryResultsX() to get X-values
|
|
||||||
* KDTreeQueryResultsXY() to get X- and Y-values
|
|
||||||
* KDTreeQueryResultsTags() to get tag values
|
|
||||||
* KDTreeQueryResultsDistances() to get distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
ae_int_t kdtreequeryknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const bool selfmatch);
|
|
||||||
ae_int_t kdtreequeryknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
R-NN query: all points within R-sphere centered at X
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
X - point, array[0..NX-1].
|
|
||||||
R - radius of sphere (in corresponding norm), R>0
|
|
||||||
SelfMatch - whether self-matches are allowed:
|
|
||||||
* if True, nearest neighbor may be the point itself
|
|
||||||
(if it exists in original dataset)
|
|
||||||
* if False, then only points with non-zero distance
|
|
||||||
are returned
|
|
||||||
* if not given, considered True
|
|
||||||
|
|
||||||
RESULT
|
|
||||||
number of neighbors found, >=0
|
|
||||||
|
|
||||||
This subroutine performs query and stores its result in the internal
|
|
||||||
structures of the KD-tree. You can use following subroutines to obtain
|
|
||||||
actual results:
|
|
||||||
* KDTreeQueryResultsX() to get X-values
|
|
||||||
* KDTreeQueryResultsXY() to get X- and Y-values
|
|
||||||
* KDTreeQueryResultsTags() to get tag values
|
|
||||||
* KDTreeQueryResultsDistances() to get distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double r, const bool selfmatch);
|
|
||||||
ae_int_t kdtreequeryrnn(const kdtree &kdt, const real_1d_array &x, const double r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
K-NN query: approximate K nearest neighbors
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
X - point, array[0..NX-1].
|
|
||||||
K - number of neighbors to return, K>=1
|
|
||||||
SelfMatch - whether self-matches are allowed:
|
|
||||||
* if True, nearest neighbor may be the point itself
|
|
||||||
(if it exists in original dataset)
|
|
||||||
* if False, then only points with non-zero distance
|
|
||||||
are returned
|
|
||||||
* if not given, considered True
|
|
||||||
Eps - approximation factor, Eps>=0. eps-approximate nearest
|
|
||||||
neighbor is a neighbor whose distance from X is at
|
|
||||||
most (1+eps) times distance of true nearest neighbor.
|
|
||||||
|
|
||||||
RESULT
|
|
||||||
number of actual neighbors found (either K or N, if K>N).
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
significant performance gain may be achieved only when Eps is is on
|
|
||||||
the order of magnitude of 1 or larger.
|
|
||||||
|
|
||||||
This subroutine performs query and stores its result in the internal
|
|
||||||
structures of the KD-tree. You can use following subroutines to obtain
|
|
||||||
these results:
|
|
||||||
* KDTreeQueryResultsX() to get X-values
|
|
||||||
* KDTreeQueryResultsXY() to get X- and Y-values
|
|
||||||
* KDTreeQueryResultsTags() to get tag values
|
|
||||||
* KDTreeQueryResultsDistances() to get distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
ae_int_t kdtreequeryaknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const bool selfmatch, const double eps);
|
|
||||||
ae_int_t kdtreequeryaknn(const kdtree &kdt, const real_1d_array &x, const ae_int_t k, const double eps);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
X-values from last query
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
X - possibly pre-allocated buffer. If X is too small to store
|
|
||||||
result, it is resized. If size(X) is enough to store
|
|
||||||
result, it is left unchanged.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
X - rows are filled with X-values
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
1. points are ordered by distance from the query point (first = closest)
|
|
||||||
2. if XY is larger than required to store result, only leading part will
|
|
||||||
be overwritten; trailing part will be left unchanged. So if on input
|
|
||||||
XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get
|
|
||||||
XY = [[1,2],[C,D]]. This is done purposely to increase performance; if
|
|
||||||
you want function to resize array according to result size, use
|
|
||||||
function with same name and suffix 'I'.
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
* KDTreeQueryResultsXY() X- and Y-values
|
|
||||||
* KDTreeQueryResultsTags() tag values
|
|
||||||
* KDTreeQueryResultsDistances() distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsx(const kdtree &kdt, real_2d_array &x);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
X- and Y-values from last query
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
XY - possibly pre-allocated buffer. If XY is too small to store
|
|
||||||
result, it is resized. If size(XY) is enough to store
|
|
||||||
result, it is left unchanged.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
XY - rows are filled with points: first NX columns with
|
|
||||||
X-values, next NY columns - with Y-values.
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
1. points are ordered by distance from the query point (first = closest)
|
|
||||||
2. if XY is larger than required to store result, only leading part will
|
|
||||||
be overwritten; trailing part will be left unchanged. So if on input
|
|
||||||
XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get
|
|
||||||
XY = [[1,2],[C,D]]. This is done purposely to increase performance; if
|
|
||||||
you want function to resize array according to result size, use
|
|
||||||
function with same name and suffix 'I'.
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
* KDTreeQueryResultsX() X-values
|
|
||||||
* KDTreeQueryResultsTags() tag values
|
|
||||||
* KDTreeQueryResultsDistances() distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsxy(const kdtree &kdt, real_2d_array &xy);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Tags from last query
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
Tags - possibly pre-allocated buffer. If X is too small to store
|
|
||||||
result, it is resized. If size(X) is enough to store
|
|
||||||
result, it is left unchanged.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
Tags - filled with tags associated with points,
|
|
||||||
or, when no tags were supplied, with zeros
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
1. points are ordered by distance from the query point (first = closest)
|
|
||||||
2. if XY is larger than required to store result, only leading part will
|
|
||||||
be overwritten; trailing part will be left unchanged. So if on input
|
|
||||||
XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get
|
|
||||||
XY = [[1,2],[C,D]]. This is done purposely to increase performance; if
|
|
||||||
you want function to resize array according to result size, use
|
|
||||||
function with same name and suffix 'I'.
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
* KDTreeQueryResultsX() X-values
|
|
||||||
* KDTreeQueryResultsXY() X- and Y-values
|
|
||||||
* KDTreeQueryResultsDistances() distances
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultstags(const kdtree &kdt, integer_1d_array &tags);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Distances from last query
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
KDT - KD-tree
|
|
||||||
R - possibly pre-allocated buffer. If X is too small to store
|
|
||||||
result, it is resized. If size(X) is enough to store
|
|
||||||
result, it is left unchanged.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - filled with distances (in corresponding norm)
|
|
||||||
|
|
||||||
NOTES
|
|
||||||
1. points are ordered by distance from the query point (first = closest)
|
|
||||||
2. if XY is larger than required to store result, only leading part will
|
|
||||||
be overwritten; trailing part will be left unchanged. So if on input
|
|
||||||
XY = [[A,B],[C,D]], and result is [1,2], then on exit we will get
|
|
||||||
XY = [[1,2],[C,D]]. This is done purposely to increase performance; if
|
|
||||||
you want function to resize array according to result size, use
|
|
||||||
function with same name and suffix 'I'.
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
* KDTreeQueryResultsX() X-values
|
|
||||||
* KDTreeQueryResultsXY() X- and Y-values
|
|
||||||
* KDTreeQueryResultsTags() tag values
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsdistances(const kdtree &kdt, real_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
X-values from last query; 'interactive' variant for languages like Python
|
|
||||||
which support constructs like "X = KDTreeQueryResultsXI(KDT)" and
|
|
||||||
interactive mode of interpreter.
|
|
||||||
|
|
||||||
This function allocates new array on each call, so it is significantly
|
|
||||||
slower than its 'non-interactive' counterpart, but it is more convenient
|
|
||||||
when you call it from command line.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsxi(const kdtree &kdt, real_2d_array &x);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
XY-values from last query; 'interactive' variant for languages like Python
|
|
||||||
which support constructs like "XY = KDTreeQueryResultsXYI(KDT)" and
|
|
||||||
interactive mode of interpreter.
|
|
||||||
|
|
||||||
This function allocates new array on each call, so it is significantly
|
|
||||||
slower than its 'non-interactive' counterpart, but it is more convenient
|
|
||||||
when you call it from command line.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsxyi(const kdtree &kdt, real_2d_array &xy);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Tags from last query; 'interactive' variant for languages like Python
|
|
||||||
which support constructs like "Tags = KDTreeQueryResultsTagsI(KDT)" and
|
|
||||||
interactive mode of interpreter.
|
|
||||||
|
|
||||||
This function allocates new array on each call, so it is significantly
|
|
||||||
slower than its 'non-interactive' counterpart, but it is more convenient
|
|
||||||
when you call it from command line.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultstagsi(const kdtree &kdt, integer_1d_array &tags);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Distances from last query; 'interactive' variant for languages like Python
|
|
||||||
which support constructs like "R = KDTreeQueryResultsDistancesI(KDT)"
|
|
||||||
and interactive mode of interpreter.
|
|
||||||
|
|
||||||
This function allocates new array on each call, so it is significantly
|
|
||||||
slower than its 'non-interactive' counterpart, but it is more convenient
|
|
||||||
when you call it from command line.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 28.02.2010 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void kdtreequeryresultsdistancesi(const kdtree &kdt, real_1d_array &r);
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
void hqrndrandomize(hqrndstate* state, ae_state *_state);
|
|
||||||
void hqrndseed(ae_int_t s1,
|
|
||||||
ae_int_t s2,
|
|
||||||
hqrndstate* state,
|
|
||||||
ae_state *_state);
|
|
||||||
double hqrnduniformr(hqrndstate* state, ae_state *_state);
|
|
||||||
ae_int_t hqrnduniformi(hqrndstate* state, ae_int_t n, ae_state *_state);
|
|
||||||
double hqrndnormal(hqrndstate* state, ae_state *_state);
|
|
||||||
void hqrndunit2(hqrndstate* state, double* x, double* y, ae_state *_state);
|
|
||||||
void hqrndnormal2(hqrndstate* state,
|
|
||||||
double* x1,
|
|
||||||
double* x2,
|
|
||||||
ae_state *_state);
|
|
||||||
double hqrndexponential(hqrndstate* state,
|
|
||||||
double lambdav,
|
|
||||||
ae_state *_state);
|
|
||||||
double hqrnddiscrete(hqrndstate* state,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_state *_state);
|
|
||||||
double hqrndcontinuous(hqrndstate* state,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_bool _hqrndstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _hqrndstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _hqrndstate_clear(void* _p);
|
|
||||||
void _hqrndstate_destroy(void* _p);
|
|
||||||
void kdtreebuild(/* Real */ ae_matrix* xy,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t nx,
|
|
||||||
ae_int_t ny,
|
|
||||||
ae_int_t normtype,
|
|
||||||
kdtree* kdt,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreebuildtagged(/* Real */ ae_matrix* xy,
|
|
||||||
/* Integer */ ae_vector* tags,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t nx,
|
|
||||||
ae_int_t ny,
|
|
||||||
ae_int_t normtype,
|
|
||||||
kdtree* kdt,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_int_t kdtreequeryknn(kdtree* kdt,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
ae_int_t k,
|
|
||||||
ae_bool selfmatch,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_int_t kdtreequeryrnn(kdtree* kdt,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
double r,
|
|
||||||
ae_bool selfmatch,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_int_t kdtreequeryaknn(kdtree* kdt,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
ae_int_t k,
|
|
||||||
ae_bool selfmatch,
|
|
||||||
double eps,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsx(kdtree* kdt,
|
|
||||||
/* Real */ ae_matrix* x,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsxy(kdtree* kdt,
|
|
||||||
/* Real */ ae_matrix* xy,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultstags(kdtree* kdt,
|
|
||||||
/* Integer */ ae_vector* tags,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsdistances(kdtree* kdt,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsxi(kdtree* kdt,
|
|
||||||
/* Real */ ae_matrix* x,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsxyi(kdtree* kdt,
|
|
||||||
/* Real */ ae_matrix* xy,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultstagsi(kdtree* kdt,
|
|
||||||
/* Integer */ ae_vector* tags,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreequeryresultsdistancesi(kdtree* kdt,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void kdtreealloc(ae_serializer* s, kdtree* tree, ae_state *_state);
|
|
||||||
void kdtreeserialize(ae_serializer* s, kdtree* tree, ae_state *_state);
|
|
||||||
void kdtreeunserialize(ae_serializer* s, kdtree* tree, ae_state *_state);
|
|
||||||
ae_bool _kdtree_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _kdtree_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _kdtree_clear(void* _p);
|
|
||||||
void _kdtree_destroy(void* _p);
|
|
||||||
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
10661
src/inc/alglib/ap.cpp
10661
src/inc/alglib/ap.cpp
File diff suppressed because it is too large
Load Diff
1575
src/inc/alglib/ap.h
1575
src/inc/alglib/ap.h
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -1,267 +0,0 @@
|
|||||||
/*************************************************************************
|
|
||||||
Copyright (c) Sergey Bochkanov (ALGLIB project).
|
|
||||||
|
|
||||||
>>> SOURCE LICENSE >>>
|
|
||||||
This program is free software; you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation (www.fsf.org); either version 2 of the
|
|
||||||
License, or (at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
A copy of the GNU General Public License is available at
|
|
||||||
http://www.fsf.org/licensing/licenses
|
|
||||||
>>> END OF LICENSE >>>
|
|
||||||
*************************************************************************/
|
|
||||||
#ifndef _diffequations_pkg_h
|
|
||||||
#define _diffequations_pkg_h
|
|
||||||
#include "ap.h"
|
|
||||||
#include "alglibinternal.h"
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
ae_int_t n;
|
|
||||||
ae_int_t m;
|
|
||||||
double xscale;
|
|
||||||
double h;
|
|
||||||
double eps;
|
|
||||||
ae_bool fraceps;
|
|
||||||
ae_vector yc;
|
|
||||||
ae_vector escale;
|
|
||||||
ae_vector xg;
|
|
||||||
ae_int_t solvertype;
|
|
||||||
ae_bool needdy;
|
|
||||||
double x;
|
|
||||||
ae_vector y;
|
|
||||||
ae_vector dy;
|
|
||||||
ae_matrix ytbl;
|
|
||||||
ae_int_t repterminationtype;
|
|
||||||
ae_int_t repnfev;
|
|
||||||
ae_vector yn;
|
|
||||||
ae_vector yns;
|
|
||||||
ae_vector rka;
|
|
||||||
ae_vector rkc;
|
|
||||||
ae_vector rkcs;
|
|
||||||
ae_matrix rkb;
|
|
||||||
ae_matrix rkk;
|
|
||||||
rcommstate rstate;
|
|
||||||
} odesolverstate;
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
ae_int_t nfev;
|
|
||||||
ae_int_t terminationtype;
|
|
||||||
} odesolverreport;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS C++ INTERFACE
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib
|
|
||||||
{
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
|
|
||||||
*************************************************************************/
|
|
||||||
class _odesolverstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_odesolverstate_owner();
|
|
||||||
_odesolverstate_owner(const _odesolverstate_owner &rhs);
|
|
||||||
_odesolverstate_owner& operator=(const _odesolverstate_owner &rhs);
|
|
||||||
virtual ~_odesolverstate_owner();
|
|
||||||
alglib_impl::odesolverstate* c_ptr();
|
|
||||||
alglib_impl::odesolverstate* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::odesolverstate *p_struct;
|
|
||||||
};
|
|
||||||
class odesolverstate : public _odesolverstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
odesolverstate();
|
|
||||||
odesolverstate(const odesolverstate &rhs);
|
|
||||||
odesolverstate& operator=(const odesolverstate &rhs);
|
|
||||||
virtual ~odesolverstate();
|
|
||||||
ae_bool &needdy;
|
|
||||||
real_1d_array y;
|
|
||||||
real_1d_array dy;
|
|
||||||
double &x;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
|
|
||||||
*************************************************************************/
|
|
||||||
class _odesolverreport_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_odesolverreport_owner();
|
|
||||||
_odesolverreport_owner(const _odesolverreport_owner &rhs);
|
|
||||||
_odesolverreport_owner& operator=(const _odesolverreport_owner &rhs);
|
|
||||||
virtual ~_odesolverreport_owner();
|
|
||||||
alglib_impl::odesolverreport* c_ptr();
|
|
||||||
alglib_impl::odesolverreport* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::odesolverreport *p_struct;
|
|
||||||
};
|
|
||||||
class odesolverreport : public _odesolverreport_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
odesolverreport();
|
|
||||||
odesolverreport(const odesolverreport &rhs);
|
|
||||||
odesolverreport& operator=(const odesolverreport &rhs);
|
|
||||||
virtual ~odesolverreport();
|
|
||||||
ae_int_t &nfev;
|
|
||||||
ae_int_t &terminationtype;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Cash-Karp adaptive ODE solver.
|
|
||||||
|
|
||||||
This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys
|
|
||||||
(here Y may be single variable or vector of N variables).
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
Y - initial conditions, array[0..N-1].
|
|
||||||
contains values of Y[] at X[0]
|
|
||||||
N - system size
|
|
||||||
X - points at which Y should be tabulated, array[0..M-1]
|
|
||||||
integrations starts at X[0], ends at X[M-1], intermediate
|
|
||||||
values at X[i] are returned too.
|
|
||||||
SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!
|
|
||||||
M - number of intermediate points + first point + last point:
|
|
||||||
* M>2 means that you need both Y(X[M-1]) and M-2 values at
|
|
||||||
intermediate points
|
|
||||||
* M=2 means that you want just to integrate from X[0] to
|
|
||||||
X[1] and don't interested in intermediate values.
|
|
||||||
* M=1 means that you don't want to integrate :)
|
|
||||||
it is degenerate case, but it will be handled correctly.
|
|
||||||
* M<1 means error
|
|
||||||
Eps - tolerance (absolute/relative error on each step will be
|
|
||||||
less than Eps). When passing:
|
|
||||||
* Eps>0, it means desired ABSOLUTE error
|
|
||||||
* Eps<0, it means desired RELATIVE error. Relative errors
|
|
||||||
are calculated with respect to maximum values of Y seen
|
|
||||||
so far. Be careful to use this criterion when starting
|
|
||||||
from Y[] that are close to zero.
|
|
||||||
H - initial step lenth, it will be adjusted automatically
|
|
||||||
after the first step. If H=0, step will be selected
|
|
||||||
automatically (usualy it will be equal to 0.001 of
|
|
||||||
min(x[i]-x[j])).
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
State - structure which stores algorithm state between subsequent
|
|
||||||
calls of OdeSolverIteration. Used for reverse communication.
|
|
||||||
This structure should be passed to the OdeSolverIteration
|
|
||||||
subroutine.
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 01.09.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void odesolverrkck(const real_1d_array &y, const ae_int_t n, const real_1d_array &x, const ae_int_t m, const double eps, const double h, odesolverstate &state);
|
|
||||||
void odesolverrkck(const real_1d_array &y, const real_1d_array &x, const double eps, const double h, odesolverstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function provides reverse communication interface
|
|
||||||
Reverse communication interface is not documented or recommended to use.
|
|
||||||
See below for functions which provide better documented API
|
|
||||||
*************************************************************************/
|
|
||||||
bool odesolveriteration(const odesolverstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function is used to launcn iterations of ODE solver
|
|
||||||
|
|
||||||
It accepts following parameters:
|
|
||||||
diff - callback which calculates dy/dx for given y and x
|
|
||||||
ptr - optional pointer which is passed to diff; can be NULL
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 01.09.2009 by Bochkanov Sergey
|
|
||||||
|
|
||||||
*************************************************************************/
|
|
||||||
void odesolversolve(odesolverstate &state,
|
|
||||||
void (*diff)(const real_1d_array &y, double x, real_1d_array &dy, void *ptr),
|
|
||||||
void *ptr = NULL);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
ODE solver results
|
|
||||||
|
|
||||||
Called after OdeSolverIteration returned False.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
State - algorithm state (used by OdeSolverIteration).
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
M - number of tabulated values, M>=1
|
|
||||||
XTbl - array[0..M-1], values of X
|
|
||||||
YTbl - array[0..M-1,0..N-1], values of Y in X[i]
|
|
||||||
Rep - solver report:
|
|
||||||
* Rep.TerminationType completetion code:
|
|
||||||
* -2 X is not ordered by ascending/descending or
|
|
||||||
there are non-distinct X[], i.e. X[i]=X[i+1]
|
|
||||||
* -1 incorrect parameters were specified
|
|
||||||
* 1 task has been solved
|
|
||||||
* Rep.NFEV contains number of function calculations
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 01.09.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void odesolverresults(const odesolverstate &state, ae_int_t &m, real_1d_array &xtbl, real_2d_array &ytbl, odesolverreport &rep);
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
void odesolverrkck(/* Real */ ae_vector* y,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
ae_int_t m,
|
|
||||||
double eps,
|
|
||||||
double h,
|
|
||||||
odesolverstate* state,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_bool odesolveriteration(odesolverstate* state, ae_state *_state);
|
|
||||||
void odesolverresults(odesolverstate* state,
|
|
||||||
ae_int_t* m,
|
|
||||||
/* Real */ ae_vector* xtbl,
|
|
||||||
/* Real */ ae_matrix* ytbl,
|
|
||||||
odesolverreport* rep,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_bool _odesolverstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _odesolverstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _odesolverstate_clear(void* _p);
|
|
||||||
void _odesolverstate_destroy(void* _p);
|
|
||||||
ae_bool _odesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _odesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _odesolverreport_clear(void* _p);
|
|
||||||
void _odesolverreport_destroy(void* _p);
|
|
||||||
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
@ -1,691 +0,0 @@
|
|||||||
/*************************************************************************
|
|
||||||
Copyright (c) Sergey Bochkanov (ALGLIB project).
|
|
||||||
|
|
||||||
>>> SOURCE LICENSE >>>
|
|
||||||
This program is free software; you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation (www.fsf.org); either version 2 of the
|
|
||||||
License, or (at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
A copy of the GNU General Public License is available at
|
|
||||||
http://www.fsf.org/licensing/licenses
|
|
||||||
>>> END OF LICENSE >>>
|
|
||||||
*************************************************************************/
|
|
||||||
#ifndef _fasttransforms_pkg_h
|
|
||||||
#define _fasttransforms_pkg_h
|
|
||||||
#include "ap.h"
|
|
||||||
#include "alglibinternal.h"
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS C++ INTERFACE
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib
|
|
||||||
{
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex FFT.
|
|
||||||
|
|
||||||
Array size N may be arbitrary number (composite or prime). Composite N's
|
|
||||||
are handled with cache-oblivious variation of a Cooley-Tukey algorithm.
|
|
||||||
Small prime-factors are transformed using hard coded codelets (similar to
|
|
||||||
FFTW codelets, but without low-level optimization), large prime-factors
|
|
||||||
are handled with Bluestein's algorithm.
|
|
||||||
|
|
||||||
Fastests transforms are for smooth N's (prime factors are 2, 3, 5 only),
|
|
||||||
most fast for powers of 2. When N have prime factors larger than these,
|
|
||||||
but orders of magnitude smaller than N, computations will be about 4 times
|
|
||||||
slower than for nearby highly composite N's. When N itself is prime, speed
|
|
||||||
will be 6 times lower.
|
|
||||||
|
|
||||||
Algorithm has O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..N-1] - complex function to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
A - DFT of a input array, array[0..N-1]
|
|
||||||
A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 29.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fftc1d(complex_1d_array &a, const ae_int_t n);
|
|
||||||
void fftc1d(complex_1d_array &a);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex inverse FFT.
|
|
||||||
|
|
||||||
Array size N may be arbitrary number (composite or prime). Algorithm has
|
|
||||||
O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
See FFTC1D() description for more information about algorithm performance.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..N-1] - complex array to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
A - inverse DFT of a input array, array[0..N-1]
|
|
||||||
A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 29.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fftc1dinv(complex_1d_array &a, const ae_int_t n);
|
|
||||||
void fftc1dinv(complex_1d_array &a);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional real FFT.
|
|
||||||
|
|
||||||
Algorithm has O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..N-1] - real function to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
F - DFT of a input array, array[0..N-1]
|
|
||||||
F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
F[] satisfies symmetry property F[k] = conj(F[N-k]), so just one half
|
|
||||||
of array is usually needed. But for convinience subroutine returns full
|
|
||||||
complex array (with frequencies above N/2), so its result may be used by
|
|
||||||
other FFT-related subroutines.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 01.06.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fftr1d(const real_1d_array &a, const ae_int_t n, complex_1d_array &f);
|
|
||||||
void fftr1d(const real_1d_array &a, complex_1d_array &f);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional real inverse FFT.
|
|
||||||
|
|
||||||
Algorithm has O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
F - array[0..floor(N/2)] - frequencies from forward real FFT
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
A - inverse DFT of a input array, array[0..N-1]
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just one
|
|
||||||
half of frequencies array is needed - elements from 0 to floor(N/2). F[0]
|
|
||||||
is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd, then
|
|
||||||
F[floor(N/2)] has no special properties.
|
|
||||||
|
|
||||||
Relying on properties noted above, FFTR1DInv subroutine uses only elements
|
|
||||||
from 0th to floor(N/2)-th. It ignores imaginary part of F[0], and in case
|
|
||||||
N is even it ignores imaginary part of F[floor(N/2)] too.
|
|
||||||
|
|
||||||
When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
|
|
||||||
- you can pass either either frequencies array with N elements or reduced
|
|
||||||
array with roughly N/2 elements - subroutine will successfully transform
|
|
||||||
both.
|
|
||||||
|
|
||||||
If you call this function using reduced arguments list - "FFTR1DInv(F,A)"
|
|
||||||
- you must pass FULL array with N elements (although higher N/2 are still
|
|
||||||
not used) because array size is used to automatically determine FFT length
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 01.06.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fftr1dinv(const complex_1d_array &f, const ae_int_t n, real_1d_array &a);
|
|
||||||
void fftr1dinv(const complex_1d_array &f, real_1d_array &a);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex convolution.
|
|
||||||
|
|
||||||
For given A/B returns conv(A,B) (non-circular). Subroutine can automatically
|
|
||||||
choose between three implementations: straightforward O(M*N) formula for
|
|
||||||
very small N (or M), overlap-add algorithm for cases where max(M,N) is
|
|
||||||
significantly larger than min(M,N), but O(M*N) algorithm is too slow, and
|
|
||||||
general FFT-based formula for cases where two previois algorithms are too
|
|
||||||
slow.
|
|
||||||
|
|
||||||
Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - complex function to be transformed
|
|
||||||
M - problem size
|
|
||||||
B - array[0..N-1] - complex function to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..N+M-2].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that A is zero at T<0, B is zero too. If one or both
|
|
||||||
functions have non-zero values at negative T's, you can still use this
|
|
||||||
subroutine - just shift its result correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convc1d(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex non-circular deconvolution (inverse of ConvC1D()).
|
|
||||||
|
|
||||||
Algorithm has M*log(M)) complexity for any M (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - convolved signal, A = conv(R, B)
|
|
||||||
M - convolved signal length
|
|
||||||
B - array[0..N-1] - response
|
|
||||||
N - response length, N<=M
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - deconvolved signal. array[0..M-N].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
deconvolution is unstable process and may result in division by zero
|
|
||||||
(if your response function is degenerate, i.e. has zero Fourier coefficient).
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that A is zero at T<0, B is zero too. If one or both
|
|
||||||
functions have non-zero values at negative T's, you can still use this
|
|
||||||
subroutine - just shift its result correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convc1dinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional circular complex convolution.
|
|
||||||
|
|
||||||
For given S/R returns conv(S,R) (circular). Algorithm has linearithmic
|
|
||||||
complexity for any M/N.
|
|
||||||
|
|
||||||
IMPORTANT: normal convolution is commutative, i.e. it is symmetric -
|
|
||||||
conv(A,B)=conv(B,A). Cyclic convolution IS NOT. One function - S - is a
|
|
||||||
signal, periodic function, and another - R - is a response, non-periodic
|
|
||||||
function with limited length.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
S - array[0..M-1] - complex periodic signal
|
|
||||||
M - problem size
|
|
||||||
B - array[0..N-1] - complex non-periodic response
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..M-1].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that B is zero at T<0. If it has non-zero values at
|
|
||||||
negative T's, you can still use this subroutine - just shift its result
|
|
||||||
correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convc1dcircular(const complex_1d_array &s, const ae_int_t m, const complex_1d_array &r, const ae_int_t n, complex_1d_array &c);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional circular complex deconvolution (inverse of ConvC1DCircular()).
|
|
||||||
|
|
||||||
Algorithm has M*log(M)) complexity for any M (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - convolved periodic signal, A = conv(R, B)
|
|
||||||
M - convolved signal length
|
|
||||||
B - array[0..N-1] - non-periodic response
|
|
||||||
N - response length
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - deconvolved signal. array[0..M-1].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
deconvolution is unstable process and may result in division by zero
|
|
||||||
(if your response function is degenerate, i.e. has zero Fourier coefficient).
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that B is zero at T<0. If it has non-zero values at
|
|
||||||
negative T's, you can still use this subroutine - just shift its result
|
|
||||||
correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convc1dcircularinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional real convolution.
|
|
||||||
|
|
||||||
Analogous to ConvC1D(), see ConvC1D() comments for more details.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - real function to be transformed
|
|
||||||
M - problem size
|
|
||||||
B - array[0..N-1] - real function to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..N+M-2].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that A is zero at T<0, B is zero too. If one or both
|
|
||||||
functions have non-zero values at negative T's, you can still use this
|
|
||||||
subroutine - just shift its result correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convr1d(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional real deconvolution (inverse of ConvC1D()).
|
|
||||||
|
|
||||||
Algorithm has M*log(M)) complexity for any M (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - convolved signal, A = conv(R, B)
|
|
||||||
M - convolved signal length
|
|
||||||
B - array[0..N-1] - response
|
|
||||||
N - response length, N<=M
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - deconvolved signal. array[0..M-N].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
deconvolution is unstable process and may result in division by zero
|
|
||||||
(if your response function is degenerate, i.e. has zero Fourier coefficient).
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that A is zero at T<0, B is zero too. If one or both
|
|
||||||
functions have non-zero values at negative T's, you can still use this
|
|
||||||
subroutine - just shift its result correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convr1dinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional circular real convolution.
|
|
||||||
|
|
||||||
Analogous to ConvC1DCircular(), see ConvC1DCircular() comments for more details.
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
S - array[0..M-1] - real signal
|
|
||||||
M - problem size
|
|
||||||
B - array[0..N-1] - real response
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..M-1].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that B is zero at T<0. If it has non-zero values at
|
|
||||||
negative T's, you can still use this subroutine - just shift its result
|
|
||||||
correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convr1dcircular(const real_1d_array &s, const ae_int_t m, const real_1d_array &r, const ae_int_t n, real_1d_array &c);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex deconvolution (inverse of ConvC1D()).
|
|
||||||
|
|
||||||
Algorithm has M*log(M)) complexity for any M (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..M-1] - convolved signal, A = conv(R, B)
|
|
||||||
M - convolved signal length
|
|
||||||
B - array[0..N-1] - response
|
|
||||||
N - response length
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - deconvolved signal. array[0..M-N].
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
deconvolution is unstable process and may result in division by zero
|
|
||||||
(if your response function is degenerate, i.e. has zero Fourier coefficient).
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that B is zero at T<0. If it has non-zero values at
|
|
||||||
negative T's, you can still use this subroutine - just shift its result
|
|
||||||
correspondingly.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void convr1dcircularinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional complex cross-correlation.
|
|
||||||
|
|
||||||
For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
|
|
||||||
|
|
||||||
Correlation is calculated using reduction to convolution. Algorithm with
|
|
||||||
max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
|
|
||||||
about performance).
|
|
||||||
|
|
||||||
IMPORTANT:
|
|
||||||
for historical reasons subroutine accepts its parameters in reversed
|
|
||||||
order: CorrC1D(Signal, Pattern) = Pattern x Signal (using traditional
|
|
||||||
definition of cross-correlation, denoting cross-correlation as "x").
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
Signal - array[0..N-1] - complex function to be transformed,
|
|
||||||
signal containing pattern
|
|
||||||
N - problem size
|
|
||||||
Pattern - array[0..M-1] - complex function to be transformed,
|
|
||||||
pattern to search withing signal
|
|
||||||
M - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - cross-correlation, array[0..N+M-2]:
|
|
||||||
* positive lags are stored in R[0..N-1],
|
|
||||||
R[i] = sum(conj(pattern[j])*signal[i+j]
|
|
||||||
* negative lags are stored in R[N..N+M-2],
|
|
||||||
R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j]
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
|
|
||||||
on [-K..M-1], you can still use this subroutine, just shift result by K.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void corrc1d(const complex_1d_array &signal, const ae_int_t n, const complex_1d_array &pattern, const ae_int_t m, complex_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional circular complex cross-correlation.
|
|
||||||
|
|
||||||
For given Pattern/Signal returns corr(Pattern,Signal) (circular).
|
|
||||||
Algorithm has linearithmic complexity for any M/N.
|
|
||||||
|
|
||||||
IMPORTANT:
|
|
||||||
for historical reasons subroutine accepts its parameters in reversed
|
|
||||||
order: CorrC1DCircular(Signal, Pattern) = Pattern x Signal (using
|
|
||||||
traditional definition of cross-correlation, denoting cross-correlation
|
|
||||||
as "x").
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
Signal - array[0..N-1] - complex function to be transformed,
|
|
||||||
periodic signal containing pattern
|
|
||||||
N - problem size
|
|
||||||
Pattern - array[0..M-1] - complex function to be transformed,
|
|
||||||
non-periodic pattern to search withing signal
|
|
||||||
M - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..M-1].
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void corrc1dcircular(const complex_1d_array &signal, const ae_int_t m, const complex_1d_array &pattern, const ae_int_t n, complex_1d_array &c);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional real cross-correlation.
|
|
||||||
|
|
||||||
For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
|
|
||||||
|
|
||||||
Correlation is calculated using reduction to convolution. Algorithm with
|
|
||||||
max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
|
|
||||||
about performance).
|
|
||||||
|
|
||||||
IMPORTANT:
|
|
||||||
for historical reasons subroutine accepts its parameters in reversed
|
|
||||||
order: CorrR1D(Signal, Pattern) = Pattern x Signal (using traditional
|
|
||||||
definition of cross-correlation, denoting cross-correlation as "x").
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
Signal - array[0..N-1] - real function to be transformed,
|
|
||||||
signal containing pattern
|
|
||||||
N - problem size
|
|
||||||
Pattern - array[0..M-1] - real function to be transformed,
|
|
||||||
pattern to search withing signal
|
|
||||||
M - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - cross-correlation, array[0..N+M-2]:
|
|
||||||
* positive lags are stored in R[0..N-1],
|
|
||||||
R[i] = sum(pattern[j]*signal[i+j]
|
|
||||||
* negative lags are stored in R[N..N+M-2],
|
|
||||||
R[N+M-1-i] = sum(pattern[j]*signal[-i+j]
|
|
||||||
|
|
||||||
NOTE:
|
|
||||||
It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
|
|
||||||
on [-K..M-1], you can still use this subroutine, just shift result by K.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void corrr1d(const real_1d_array &signal, const ae_int_t n, const real_1d_array &pattern, const ae_int_t m, real_1d_array &r);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional circular real cross-correlation.
|
|
||||||
|
|
||||||
For given Pattern/Signal returns corr(Pattern,Signal) (circular).
|
|
||||||
Algorithm has linearithmic complexity for any M/N.
|
|
||||||
|
|
||||||
IMPORTANT:
|
|
||||||
for historical reasons subroutine accepts its parameters in reversed
|
|
||||||
order: CorrR1DCircular(Signal, Pattern) = Pattern x Signal (using
|
|
||||||
traditional definition of cross-correlation, denoting cross-correlation
|
|
||||||
as "x").
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
Signal - array[0..N-1] - real function to be transformed,
|
|
||||||
periodic signal containing pattern
|
|
||||||
N - problem size
|
|
||||||
Pattern - array[0..M-1] - real function to be transformed,
|
|
||||||
non-periodic pattern to search withing signal
|
|
||||||
M - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
R - convolution: A*B. array[0..M-1].
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 21.07.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void corrr1dcircular(const real_1d_array &signal, const ae_int_t m, const real_1d_array &pattern, const ae_int_t n, real_1d_array &c);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional Fast Hartley Transform.
|
|
||||||
|
|
||||||
Algorithm has O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..N-1] - real function to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
A - FHT of a input array, array[0..N-1],
|
|
||||||
A_out[k] = sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)), j=0..N-1)
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 04.06.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fhtr1d(real_1d_array &a, const ae_int_t n);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
1-dimensional inverse FHT.
|
|
||||||
|
|
||||||
Algorithm has O(N*logN) complexity for any N (composite or prime).
|
|
||||||
|
|
||||||
INPUT PARAMETERS
|
|
||||||
A - array[0..N-1] - complex array to be transformed
|
|
||||||
N - problem size
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
A - inverse FHT of a input array, array[0..N-1]
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 29.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void fhtr1dinv(real_1d_array &a, const ae_int_t n);
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
void fftc1d(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state);
|
|
||||||
void fftc1dinv(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state);
|
|
||||||
void fftr1d(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* f,
|
|
||||||
ae_state *_state);
|
|
||||||
void fftr1dinv(/* Complex */ ae_vector* f,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* a,
|
|
||||||
ae_state *_state);
|
|
||||||
void fftr1dinternaleven(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* buf,
|
|
||||||
fasttransformplan* plan,
|
|
||||||
ae_state *_state);
|
|
||||||
void fftr1dinvinternaleven(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* buf,
|
|
||||||
fasttransformplan* plan,
|
|
||||||
ae_state *_state);
|
|
||||||
void convc1d(/* Complex */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convc1dinv(/* Complex */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convc1dcircular(/* Complex */ ae_vector* s,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* c,
|
|
||||||
ae_state *_state);
|
|
||||||
void convc1dcircularinv(/* Complex */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convr1d(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convr1dinv(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convr1dcircular(/* Real */ ae_vector* s,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* c,
|
|
||||||
ae_state *_state);
|
|
||||||
void convr1dcircularinv(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convc1dx(/* Complex */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_bool circular,
|
|
||||||
ae_int_t alg,
|
|
||||||
ae_int_t q,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void convr1dx(/* Real */ ae_vector* a,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* b,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_bool circular,
|
|
||||||
ae_int_t alg,
|
|
||||||
ae_int_t q,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void corrc1d(/* Complex */ ae_vector* signal,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* pattern,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void corrc1dcircular(/* Complex */ ae_vector* signal,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Complex */ ae_vector* pattern,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Complex */ ae_vector* c,
|
|
||||||
ae_state *_state);
|
|
||||||
void corrr1d(/* Real */ ae_vector* signal,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* pattern,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* r,
|
|
||||||
ae_state *_state);
|
|
||||||
void corrr1dcircular(/* Real */ ae_vector* signal,
|
|
||||||
ae_int_t m,
|
|
||||||
/* Real */ ae_vector* pattern,
|
|
||||||
ae_int_t n,
|
|
||||||
/* Real */ ae_vector* c,
|
|
||||||
ae_state *_state);
|
|
||||||
void fhtr1d(/* Real */ ae_vector* a, ae_int_t n, ae_state *_state);
|
|
||||||
void fhtr1dinv(/* Real */ ae_vector* a, ae_int_t n, ae_state *_state);
|
|
||||||
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
@ -1,837 +0,0 @@
|
|||||||
/*************************************************************************
|
|
||||||
Copyright (c) Sergey Bochkanov (ALGLIB project).
|
|
||||||
|
|
||||||
>>> SOURCE LICENSE >>>
|
|
||||||
This program is free software; you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation (www.fsf.org); either version 2 of the
|
|
||||||
License, or (at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
A copy of the GNU General Public License is available at
|
|
||||||
http://www.fsf.org/licensing/licenses
|
|
||||||
>>> END OF LICENSE >>>
|
|
||||||
*************************************************************************/
|
|
||||||
#ifndef _integration_pkg_h
|
|
||||||
#define _integration_pkg_h
|
|
||||||
#include "ap.h"
|
|
||||||
#include "alglibinternal.h"
|
|
||||||
#include "linalg.h"
|
|
||||||
#include "specialfunctions.h"
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
ae_int_t terminationtype;
|
|
||||||
ae_int_t nfev;
|
|
||||||
ae_int_t nintervals;
|
|
||||||
} autogkreport;
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
double a;
|
|
||||||
double b;
|
|
||||||
double eps;
|
|
||||||
double xwidth;
|
|
||||||
double x;
|
|
||||||
double f;
|
|
||||||
ae_int_t info;
|
|
||||||
double r;
|
|
||||||
ae_matrix heap;
|
|
||||||
ae_int_t heapsize;
|
|
||||||
ae_int_t heapwidth;
|
|
||||||
ae_int_t heapused;
|
|
||||||
double sumerr;
|
|
||||||
double sumabs;
|
|
||||||
ae_vector qn;
|
|
||||||
ae_vector wg;
|
|
||||||
ae_vector wk;
|
|
||||||
ae_vector wr;
|
|
||||||
ae_int_t n;
|
|
||||||
rcommstate rstate;
|
|
||||||
} autogkinternalstate;
|
|
||||||
typedef struct
|
|
||||||
{
|
|
||||||
double a;
|
|
||||||
double b;
|
|
||||||
double alpha;
|
|
||||||
double beta;
|
|
||||||
double xwidth;
|
|
||||||
double x;
|
|
||||||
double xminusa;
|
|
||||||
double bminusx;
|
|
||||||
ae_bool needf;
|
|
||||||
double f;
|
|
||||||
ae_int_t wrappermode;
|
|
||||||
autogkinternalstate internalstate;
|
|
||||||
rcommstate rstate;
|
|
||||||
double v;
|
|
||||||
ae_int_t terminationtype;
|
|
||||||
ae_int_t nfev;
|
|
||||||
ae_int_t nintervals;
|
|
||||||
} autogkstate;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS C++ INTERFACE
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib
|
|
||||||
{
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Integration report:
|
|
||||||
* TerminationType = completetion code:
|
|
||||||
* -5 non-convergence of Gauss-Kronrod nodes
|
|
||||||
calculation subroutine.
|
|
||||||
* -1 incorrect parameters were specified
|
|
||||||
* 1 OK
|
|
||||||
* Rep.NFEV countains number of function calculations
|
|
||||||
* Rep.NIntervals contains number of intervals [a,b]
|
|
||||||
was partitioned into.
|
|
||||||
*************************************************************************/
|
|
||||||
class _autogkreport_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_autogkreport_owner();
|
|
||||||
_autogkreport_owner(const _autogkreport_owner &rhs);
|
|
||||||
_autogkreport_owner& operator=(const _autogkreport_owner &rhs);
|
|
||||||
virtual ~_autogkreport_owner();
|
|
||||||
alglib_impl::autogkreport* c_ptr();
|
|
||||||
alglib_impl::autogkreport* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::autogkreport *p_struct;
|
|
||||||
};
|
|
||||||
class autogkreport : public _autogkreport_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
autogkreport();
|
|
||||||
autogkreport(const autogkreport &rhs);
|
|
||||||
autogkreport& operator=(const autogkreport &rhs);
|
|
||||||
virtual ~autogkreport();
|
|
||||||
ae_int_t &terminationtype;
|
|
||||||
ae_int_t &nfev;
|
|
||||||
ae_int_t &nintervals;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This structure stores state of the integration algorithm.
|
|
||||||
|
|
||||||
Although this class has public fields, they are not intended for external
|
|
||||||
use. You should use ALGLIB functions to work with this class:
|
|
||||||
* autogksmooth()/AutoGKSmoothW()/... to create objects
|
|
||||||
* autogkintegrate() to begin integration
|
|
||||||
* autogkresults() to get results
|
|
||||||
*************************************************************************/
|
|
||||||
class _autogkstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
_autogkstate_owner();
|
|
||||||
_autogkstate_owner(const _autogkstate_owner &rhs);
|
|
||||||
_autogkstate_owner& operator=(const _autogkstate_owner &rhs);
|
|
||||||
virtual ~_autogkstate_owner();
|
|
||||||
alglib_impl::autogkstate* c_ptr();
|
|
||||||
alglib_impl::autogkstate* c_ptr() const;
|
|
||||||
protected:
|
|
||||||
alglib_impl::autogkstate *p_struct;
|
|
||||||
};
|
|
||||||
class autogkstate : public _autogkstate_owner
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
autogkstate();
|
|
||||||
autogkstate(const autogkstate &rhs);
|
|
||||||
autogkstate& operator=(const autogkstate &rhs);
|
|
||||||
virtual ~autogkstate();
|
|
||||||
ae_bool &needf;
|
|
||||||
double &x;
|
|
||||||
double &xminusa;
|
|
||||||
double &bminusx;
|
|
||||||
double &f;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Computation of nodes and weights for a Gauss quadrature formula
|
|
||||||
|
|
||||||
The algorithm generates the N-point Gauss quadrature formula with weight
|
|
||||||
function given by coefficients alpha and beta of a recurrence relation
|
|
||||||
which generates a system of orthogonal polynomials:
|
|
||||||
|
|
||||||
P-1(x) = 0
|
|
||||||
P0(x) = 1
|
|
||||||
Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
|
|
||||||
|
|
||||||
and zeroth moment Mu0
|
|
||||||
|
|
||||||
Mu0 = integral(W(x)dx,a,b)
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
Alpha – array[0..N-1], alpha coefficients
|
|
||||||
Beta – array[0..N-1], beta coefficients
|
|
||||||
Zero-indexed element is not used and may be arbitrary.
|
|
||||||
Beta[I]>0.
|
|
||||||
Mu0 – zeroth moment of the weight function.
|
|
||||||
N – number of nodes of the quadrature formula, N>=1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -2 Beta[i]<=0
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* 1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 2005-2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgeneraterec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Computation of nodes and weights for a Gauss-Lobatto quadrature formula
|
|
||||||
|
|
||||||
The algorithm generates the N-point Gauss-Lobatto quadrature formula with
|
|
||||||
weight function given by coefficients alpha and beta of a recurrence which
|
|
||||||
generates a system of orthogonal polynomials.
|
|
||||||
|
|
||||||
P-1(x) = 0
|
|
||||||
P0(x) = 1
|
|
||||||
Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
|
|
||||||
|
|
||||||
and zeroth moment Mu0
|
|
||||||
|
|
||||||
Mu0 = integral(W(x)dx,a,b)
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
Alpha – array[0..N-2], alpha coefficients
|
|
||||||
Beta – array[0..N-2], beta coefficients.
|
|
||||||
Zero-indexed element is not used, may be arbitrary.
|
|
||||||
Beta[I]>0
|
|
||||||
Mu0 – zeroth moment of the weighting function.
|
|
||||||
A – left boundary of the integration interval.
|
|
||||||
B – right boundary of the integration interval.
|
|
||||||
N – number of nodes of the quadrature formula, N>=3
|
|
||||||
(including the left and right boundary nodes).
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -2 Beta[i]<=0
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* 1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 2005-2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategausslobattorec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const double b, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Computation of nodes and weights for a Gauss-Radau quadrature formula
|
|
||||||
|
|
||||||
The algorithm generates the N-point Gauss-Radau quadrature formula with
|
|
||||||
weight function given by the coefficients alpha and beta of a recurrence
|
|
||||||
which generates a system of orthogonal polynomials.
|
|
||||||
|
|
||||||
P-1(x) = 0
|
|
||||||
P0(x) = 1
|
|
||||||
Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
|
|
||||||
|
|
||||||
and zeroth moment Mu0
|
|
||||||
|
|
||||||
Mu0 = integral(W(x)dx,a,b)
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
Alpha – array[0..N-2], alpha coefficients.
|
|
||||||
Beta – array[0..N-1], beta coefficients
|
|
||||||
Zero-indexed element is not used.
|
|
||||||
Beta[I]>0
|
|
||||||
Mu0 – zeroth moment of the weighting function.
|
|
||||||
A – left boundary of the integration interval.
|
|
||||||
N – number of nodes of the quadrature formula, N>=2
|
|
||||||
(including the left boundary node).
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -2 Beta[i]<=0
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* 1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 2005-2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategaussradaurec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
|
|
||||||
nodes.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of nodes, >=1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. N is too large to obtain
|
|
||||||
weights/nodes with high enough accuracy.
|
|
||||||
Try to use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] with weight
|
|
||||||
function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of nodes, >=1
|
|
||||||
Alpha - power-law coefficient, Alpha>-1
|
|
||||||
Beta - power-law coefficient, Beta>-1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. Alpha or Beta are too close
|
|
||||||
to -1 to obtain weights/nodes with high enough
|
|
||||||
accuracy, or, may be, N is too large. Try to
|
|
||||||
use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N/Alpha/Beta was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) with
|
|
||||||
weight function W(x)=Power(x,Alpha)*Exp(-x)
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of nodes, >=1
|
|
||||||
Alpha - power-law coefficient, Alpha>-1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. Alpha is too close to -1 to
|
|
||||||
obtain weights/nodes with high enough accuracy
|
|
||||||
or, may be, N is too large. Try to use
|
|
||||||
multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N/Alpha was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategausslaguerre(const ae_int_t n, const double alpha, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns nodes/weights for Gauss-Hermite quadrature on (-inf,+inf) with
|
|
||||||
weight function W(x)=Exp(-x*x)
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of nodes, >=1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. May be, N is too large. Try to
|
|
||||||
use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N/Alpha was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
W - array[0..N-1] - array of quadrature weights.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gqgenerategausshermite(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Computation of nodes and weights of a Gauss-Kronrod quadrature formula
|
|
||||||
|
|
||||||
The algorithm generates the N-point Gauss-Kronrod quadrature formula with
|
|
||||||
weight function given by coefficients alpha and beta of a recurrence
|
|
||||||
relation which generates a system of orthogonal polynomials:
|
|
||||||
|
|
||||||
P-1(x) = 0
|
|
||||||
P0(x) = 1
|
|
||||||
Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
|
|
||||||
|
|
||||||
and zero moment Mu0
|
|
||||||
|
|
||||||
Mu0 = integral(W(x)dx,a,b)
|
|
||||||
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
Alpha – alpha coefficients, array[0..floor(3*K/2)].
|
|
||||||
Beta – beta coefficients, array[0..ceil(3*K/2)].
|
|
||||||
Beta[0] is not used and may be arbitrary.
|
|
||||||
Beta[I]>0.
|
|
||||||
Mu0 – zeroth moment of the weight function.
|
|
||||||
N – number of nodes of the Gauss-Kronrod quadrature formula,
|
|
||||||
N >= 3,
|
|
||||||
N = 2*K+1.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -5 no real and positive Gauss-Kronrod formula can
|
|
||||||
be created for such a weight function with a
|
|
||||||
given number of nodes.
|
|
||||||
* -4 N is too large, task may be ill conditioned -
|
|
||||||
x[i]=x[i+1] found.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -2 Beta[i]<=0
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes,
|
|
||||||
in ascending order.
|
|
||||||
WKronrod - array[0..N-1] - Kronrod weights
|
|
||||||
WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
|
|
||||||
corresponding to extended Kronrod nodes).
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 08.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gkqgeneraterec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre
|
|
||||||
quadrature with N points.
|
|
||||||
|
|
||||||
GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed table) is
|
|
||||||
used depending on machine precision and number of nodes.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of Kronrod nodes, must be odd number, >=3.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. N is too large to obtain
|
|
||||||
weights/nodes with high enough accuracy.
|
|
||||||
Try to use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes, ordered in
|
|
||||||
ascending order.
|
|
||||||
WKronrod - array[0..N-1] - Kronrod weights
|
|
||||||
WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
|
|
||||||
corresponding to extended Kronrod nodes).
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gkqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi
|
|
||||||
quadrature on [-1,1] with weight function
|
|
||||||
|
|
||||||
W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of Kronrod nodes, must be odd number, >=3.
|
|
||||||
Alpha - power-law coefficient, Alpha>-1
|
|
||||||
Beta - power-law coefficient, Beta>-1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -5 no real and positive Gauss-Kronrod formula can
|
|
||||||
be created for such a weight function with a
|
|
||||||
given number of nodes.
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. Alpha or Beta are too close
|
|
||||||
to -1 to obtain weights/nodes with high enough
|
|
||||||
accuracy, or, may be, N is too large. Try to
|
|
||||||
use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* +1 OK
|
|
||||||
* +2 OK, but quadrature rule have exterior nodes,
|
|
||||||
x[0]<-1 or x[n-1]>+1
|
|
||||||
X - array[0..N-1] - array of quadrature nodes, ordered in
|
|
||||||
ascending order.
|
|
||||||
WKronrod - array[0..N-1] - Kronrod weights
|
|
||||||
WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
|
|
||||||
corresponding to extended Kronrod nodes).
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gkqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
|
|
||||||
|
|
||||||
Reduction to tridiagonal eigenproblem is used.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of Kronrod nodes, must be odd number, >=3.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
Info - error code:
|
|
||||||
* -4 an error was detected when calculating
|
|
||||||
weights/nodes. N is too large to obtain
|
|
||||||
weights/nodes with high enough accuracy.
|
|
||||||
Try to use multiple precision version.
|
|
||||||
* -3 internal eigenproblem solver hasn't converged
|
|
||||||
* -1 incorrect N was passed
|
|
||||||
* +1 OK
|
|
||||||
X - array[0..N-1] - array of quadrature nodes, ordered in
|
|
||||||
ascending order.
|
|
||||||
WKronrod - array[0..N-1] - Kronrod weights
|
|
||||||
WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
|
|
||||||
corresponding to extended Kronrod nodes).
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gkqlegendrecalc(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Returns Gauss and Gauss-Kronrod nodes for quadrature with N points using
|
|
||||||
pre-calculated table. Nodes/weights were computed with accuracy up to
|
|
||||||
1.0E-32 (if MPFR version of ALGLIB is used). In standard double precision
|
|
||||||
accuracy reduces to something about 2.0E-16 (depending on your compiler's
|
|
||||||
handling of long floating point constants).
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
N - number of Kronrod nodes.
|
|
||||||
N can be 15, 21, 31, 41, 51, 61.
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS:
|
|
||||||
X - array[0..N-1] - array of quadrature nodes, ordered in
|
|
||||||
ascending order.
|
|
||||||
WKronrod - array[0..N-1] - Kronrod weights
|
|
||||||
WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
|
|
||||||
corresponding to extended Kronrod nodes).
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 12.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void gkqlegendretbl(const ae_int_t n, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, double &eps);
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Integration of a smooth function F(x) on a finite interval [a,b].
|
|
||||||
|
|
||||||
Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
|
|
||||||
is calculated with accuracy close to the machine precision.
|
|
||||||
|
|
||||||
Algorithm works well only with smooth integrands. It may be used with
|
|
||||||
continuous non-smooth integrands, but with less performance.
|
|
||||||
|
|
||||||
It should never be used with integrands which have integrable singularities
|
|
||||||
at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
|
|
||||||
cases.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
A, B - interval boundaries (A<B, A=B or A>B)
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
State - structure which stores algorithm state
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
AutoGKSmoothW, AutoGKSingular, AutoGKResults.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 06.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void autogksmooth(const double a, const double b, autogkstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Integration of a smooth function F(x) on a finite interval [a,b].
|
|
||||||
|
|
||||||
This subroutine is same as AutoGKSmooth(), but it guarantees that interval
|
|
||||||
[a,b] is partitioned into subintervals which have width at most XWidth.
|
|
||||||
|
|
||||||
Subroutine can be used when integrating nearly-constant function with
|
|
||||||
narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
|
|
||||||
subroutine can overlook them.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
A, B - interval boundaries (A<B, A=B or A>B)
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
State - structure which stores algorithm state
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
AutoGKSmooth, AutoGKSingular, AutoGKResults.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 06.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void autogksmoothw(const double a, const double b, const double xwidth, autogkstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Integration on a finite interval [A,B].
|
|
||||||
Integrand have integrable singularities at A/B.
|
|
||||||
|
|
||||||
F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known
|
|
||||||
alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates
|
|
||||||
from below can be used (but these estimates should be greater than -1 too).
|
|
||||||
|
|
||||||
One of alpha/beta variables (or even both alpha/beta) may be equal to 0,
|
|
||||||
which means than function F(x) is non-singular at A/B. Anyway (singular at
|
|
||||||
bounds or not), function F(x) is supposed to be continuous on (A,B).
|
|
||||||
|
|
||||||
Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
|
|
||||||
is calculated with accuracy close to the machine precision.
|
|
||||||
|
|
||||||
INPUT PARAMETERS:
|
|
||||||
A, B - interval boundaries (A<B, A=B or A>B)
|
|
||||||
Alpha - power-law coefficient of the F(x) at A,
|
|
||||||
Alpha>-1
|
|
||||||
Beta - power-law coefficient of the F(x) at B,
|
|
||||||
Beta>-1
|
|
||||||
|
|
||||||
OUTPUT PARAMETERS
|
|
||||||
State - structure which stores algorithm state
|
|
||||||
|
|
||||||
SEE ALSO
|
|
||||||
AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 06.05.2009 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void autogksingular(const double a, const double b, const double alpha, const double beta, autogkstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function provides reverse communication interface
|
|
||||||
Reverse communication interface is not documented or recommended to use.
|
|
||||||
See below for functions which provide better documented API
|
|
||||||
*************************************************************************/
|
|
||||||
bool autogkiteration(const autogkstate &state);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
This function is used to launcn iterations of the 1-dimensional integrator
|
|
||||||
|
|
||||||
It accepts following parameters:
|
|
||||||
func - callback which calculates f(x) for given x
|
|
||||||
ptr - optional pointer which is passed to func; can be NULL
|
|
||||||
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 07.05.2009 by Bochkanov Sergey
|
|
||||||
|
|
||||||
*************************************************************************/
|
|
||||||
void autogkintegrate(autogkstate &state,
|
|
||||||
void (*func)(double x, double xminusa, double bminusx, double &y, void *ptr),
|
|
||||||
void *ptr = NULL);
|
|
||||||
|
|
||||||
|
|
||||||
/*************************************************************************
|
|
||||||
Adaptive integration results
|
|
||||||
|
|
||||||
Called after AutoGKIteration returned False.
|
|
||||||
|
|
||||||
Input parameters:
|
|
||||||
State - algorithm state (used by AutoGKIteration).
|
|
||||||
|
|
||||||
Output parameters:
|
|
||||||
V - integral(f(x)dx,a,b)
|
|
||||||
Rep - optimization report (see AutoGKReport description)
|
|
||||||
|
|
||||||
-- ALGLIB --
|
|
||||||
Copyright 14.11.2007 by Bochkanov Sergey
|
|
||||||
*************************************************************************/
|
|
||||||
void autogkresults(const autogkstate &state, double &v, autogkreport &rep);
|
|
||||||
}
|
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
//
|
|
||||||
// THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS)
|
|
||||||
//
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
|
||||||
namespace alglib_impl
|
|
||||||
{
|
|
||||||
void gqgeneraterec(/* Real */ ae_vector* alpha,
|
|
||||||
/* Real */ ae_vector* beta,
|
|
||||||
double mu0,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategausslobattorec(/* Real */ ae_vector* alpha,
|
|
||||||
/* Real */ ae_vector* beta,
|
|
||||||
double mu0,
|
|
||||||
double a,
|
|
||||||
double b,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategaussradaurec(/* Real */ ae_vector* alpha,
|
|
||||||
/* Real */ ae_vector* beta,
|
|
||||||
double mu0,
|
|
||||||
double a,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategausslegendre(ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategaussjacobi(ae_int_t n,
|
|
||||||
double alpha,
|
|
||||||
double beta,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategausslaguerre(ae_int_t n,
|
|
||||||
double alpha,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gqgenerategausshermite(ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* w,
|
|
||||||
ae_state *_state);
|
|
||||||
void gkqgeneraterec(/* Real */ ae_vector* alpha,
|
|
||||||
/* Real */ ae_vector* beta,
|
|
||||||
double mu0,
|
|
||||||
ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* wkronrod,
|
|
||||||
/* Real */ ae_vector* wgauss,
|
|
||||||
ae_state *_state);
|
|
||||||
void gkqgenerategausslegendre(ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* wkronrod,
|
|
||||||
/* Real */ ae_vector* wgauss,
|
|
||||||
ae_state *_state);
|
|
||||||
void gkqgenerategaussjacobi(ae_int_t n,
|
|
||||||
double alpha,
|
|
||||||
double beta,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* wkronrod,
|
|
||||||
/* Real */ ae_vector* wgauss,
|
|
||||||
ae_state *_state);
|
|
||||||
void gkqlegendrecalc(ae_int_t n,
|
|
||||||
ae_int_t* info,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* wkronrod,
|
|
||||||
/* Real */ ae_vector* wgauss,
|
|
||||||
ae_state *_state);
|
|
||||||
void gkqlegendretbl(ae_int_t n,
|
|
||||||
/* Real */ ae_vector* x,
|
|
||||||
/* Real */ ae_vector* wkronrod,
|
|
||||||
/* Real */ ae_vector* wgauss,
|
|
||||||
double* eps,
|
|
||||||
ae_state *_state);
|
|
||||||
void autogksmooth(double a,
|
|
||||||
double b,
|
|
||||||
autogkstate* state,
|
|
||||||
ae_state *_state);
|
|
||||||
void autogksmoothw(double a,
|
|
||||||
double b,
|
|
||||||
double xwidth,
|
|
||||||
autogkstate* state,
|
|
||||||
ae_state *_state);
|
|
||||||
void autogksingular(double a,
|
|
||||||
double b,
|
|
||||||
double alpha,
|
|
||||||
double beta,
|
|
||||||
autogkstate* state,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_bool autogkiteration(autogkstate* state, ae_state *_state);
|
|
||||||
void autogkresults(autogkstate* state,
|
|
||||||
double* v,
|
|
||||||
autogkreport* rep,
|
|
||||||
ae_state *_state);
|
|
||||||
ae_bool _autogkreport_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _autogkreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _autogkreport_clear(void* _p);
|
|
||||||
void _autogkreport_destroy(void* _p);
|
|
||||||
ae_bool _autogkinternalstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _autogkinternalstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _autogkinternalstate_clear(void* _p);
|
|
||||||
void _autogkinternalstate_destroy(void* _p);
|
|
||||||
ae_bool _autogkstate_init(void* _p, ae_state *_state, ae_bool make_automatic);
|
|
||||||
ae_bool _autogkstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic);
|
|
||||||
void _autogkstate_clear(void* _p);
|
|
||||||
void _autogkstate_destroy(void* _p);
|
|
||||||
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -1,2 +0,0 @@
|
|||||||
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue
Block a user