1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
## Copyright (C) 1996-2015 Piotr Held
##
## This file is part of Octave.
##
## Octave 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;
## either version 3 of the License, or (at your option) any
## later version.
##
## Octave 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.
##
## You should have received a copy of the GNU General Public
## License along with Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn{Function File} {@var{output} =} c2g (@var{d2_out})
##
## This program calculates the Gaussian kernel correlation integral
## and its logarithmic derivatice from correlation sums calculated by d2
## (the 'c2' field of the d2 output).
##
## It uses the following formula to calculate the Gaussian kernel correlation
## integral:
##
## @iftex
## @tex
## $$C_G(r)=\frac{1}{r^2}\int_{0}^{\infty}dx\, e^{-\frac{x^2}{2r^2}}C(x)$$
## @end tex
## @end iftex
## @ifnottex
## @example
## /00 2
## 1 | / x \
## C (r) = --- | dx exp |- ----- | x C(x)
## G 2 | \ 2 /
## r /0 2r
## @end example
## @end ifnottex
##
## And the logarithmic derivative is calculated using:
##
## @iftex
## @tex
## $$D_G(r) = \frac{d}{d\log r}\log{C_G(r)}$$
## @end tex
## @end iftex
## @ifnottex
## @example
## d
## D (r) = ------- log C (r)
## G d log r G
## @end example
## @end ifnottex
##
## @strong{Input}
##
## The input needs to be the output of d2.
##
## @strong{Output}
##
## The output is a struct array of the same length as the input.
## It contains the following fiels:
##
## @table @var
## @item dim
## The dimension for each matrix @var{g}.
## @item g
## Matrix with three columns. The first contains epsilon (the first column
## of field 'c2' from the d2 output), the second is the Gaussian kernel
## correlation integral and the third its logarithmic derivative.
## @end table
##
## @seealso{demo c2g, d2, c2t, av_d2}
##
## @strong{Algorithms}
##
## The algorithms for this functions have been taken from the TISEAN package.
## @end deftypefn
## Author: Piotr Held <pjheld@gmail.com>.
## This function is based on c2g of TISEAN 3.0.1
## https://github.com/heggus/Tisean"
function output = c2g (d2_out)
if (nargin != 1)
print_usage;
endif
if ((!isfield (d2_out, "dim")) || (!isfield (d2_out, "c2")) ...
|| (!isfield (d2_out, "d2")) || (!isfield (d2_out, "h2")))
error ('Octave:invalid-input-arg', "d2_out must be the output of d2");
endif
if (! exist ("__c2g__"))
error ('Octave:tisean', ["__c2g__ was not found in path ",...
"- if package was installed properly the ",...
"compiler might not support C++11 standard"])
endif
# Create the cell that will become the 'g' field of the output
g_out = cell (length (d2_out),1);
# Calculate output for each struct in the input struct array
for i = 1:size(d2_out,1)
tmp = d2_out(i);
# Limit to only the first positive correlation sums
# (do not calculate output for any past first negative sum)
idx_lt0 = min (find (tmp.c2(:,2) <= 0));
if (!isempty (idx_lt0))
tmp.c2 = tmp.c2(1:idx_lt0-1,:);
endif
[s,idx] = sort (tmp.c2(:,1));
emat = log (tmp.c2(idx,1));
cmat = log (tmp.c2(idx,2));
# Create column vectors instead of using loop
k_id = 1:length(tmp.c2)-1;
f = exp((emat(k_id+1).*cmat(k_id).-emat(k_id).*cmat(k_id+1))
./(emat(k_id+1).-emat(k_id)));
d = (cmat(k_id+1).-cmat(k_id))./(emat(k_id+1).-emat(k_id));
a = emat(k_id);
b = emat(k_id+1);
# Create row vector instead of using loop
h = exp(emat).';
# Create output
# __c2g__ performs integration on the function func and funcd
# (depending on whether 'false' or 'true' is passed as the last argument)
# from original TISEAN on vectors (matrices) of parameters
g = sum (__c2g__ (h, f, d, a, b, false));
gd = sum (__c2g__ (h, f, d, a, b, true));
de = emat(end);
cgauss = g ./ (h.^2) + exp (-exp (2.*de) ./ (2.*h.^2));
cgd = gd ./ (h.^4) + (2 + exp (2*de) ./ h.^2) ...
.* exp (-exp (2*de) ./ (2.*h.^2));
g_out{i} = [h.', cgauss.', (-2+cgd./cgauss).'];
endfor
output = struct ("dim", {d2_out.dim}.', "g", g_out);
endfunction
%!demo
%! vals = d2 (henon (5000), 'd', 1, 't', 50);
%! kernel = c2g (vals);
%! do_plot_slope = @(x) semilogx (x{1}(:,1),x{1}(:,3),'r');
%! hold on
%! arrayfun (do_plot_slope, {kernel.g});
%! hold off
%! axis tight
%! xlabel ("Epsilon")
%! ylabel ("logarithmic derivative")
%! title ("Gaussian kernel");
%!###############################################################
%% test against the result obtain from program c2g of TISEAN 3.0.1
%!test
%! c2g_res = [0.00255674845539033 0.00262371218 1.56393838;0.00274151703342795 0.00291572767 1.46379328;0.00293963821604848 0.00321952766 1.37964845;0.00315207871608436 0.00353588979 1.30904245;0.00337987160310149 0.00386583689 1.2498455;0.00362412468530238 0.00421059737 1.20023394;0.00388602935709059 0.00457157521 1.1586616;0.00416686292737722 0.00495032407 1.12383032;0.00446798978373408 0.00534852408 1.09465814;0.00479088025167584 0.00576799177 1.07024431;0.00513710314407945 0.00621065963 1.04983807;0.00550834601745009 0.00667859521 1.03281331;0.00590642029419541 0.00717399921 1.01864672;0.00633325986564159 0.00769920321 1.00690293;0.00679094903171062 0.00825670082 0.997214317;0.00728171039372683 0.00884913001 0.989262819;0.00780794164165854 0.00947930664 0.982763529;0.00837219785898924 0.0101501858 0.977451086;0.00897723622620106 0.0108648958 0.973080397;0.00962599366903305 0.011626686 0.969439983];
%! val = d2 (henon (1000), 'd', 1, 't', 50, 'm', 5);
%! res = c2g (val);
%! assert (cell2mat ({res.g}.')(1:rows(c2g_res),:), c2g_res, -1e-5);
%% input validation
%!error <output> c2g (1)
|
