SPLINT package example
SPLINT is a QCDNUM add-on for converting results computed on the evolution grid to cubic splines. This is convenient for integrating/differentiating these results.
using QCDNUM
using PrintfExample evolution grid
First, we need a QCDNUM evolution grid to work with. For this, we use the QCDNUM example program described in the example notebook.
xmin = Float64.([1.0e-4])
iwt = Int32.([1])
ng = 1
nxin = 100
iosp = 3
nx = 10
qq = Float64.([2e0, 1e4])
wt = Float64.([1e0, 1e0])
nq = 1
nqin = 60
ngq = 2
itype = 1
as0 = 0.364
r20 = 2.0
q2c = 3.0
q2b = 25.0
q0 = 2.0
iqt = 999
def = Float64.([0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,
0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,
0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]);
nfin = 0
x = 1.0e-3
q = 1.0e3
pdf = Array{Float64}(undef, 13)
qmz2 = 8315.25;
function func(ipdf, x)::Float64
i = ipdf[]
xb = x[]
adbar = 0.1939875
f = 0
if (i == 0)
ag = 1.7
f = ag * xb^-0.1 * (1.0-xb)^5.0
end
if (i == 1)
ad = 3.064320
f = ad * xb^0.8 * (1.0-xb)^4.0
end
if (i == 2)
au = 5.107200
f = au * xb^0.8 * (1.0-xb)^3.0
end
if (i == 3)
f = 0.0
end
if (i == 4)
f = adbar * xb^-0.1 * (1.0-xb)^6.0
end
if (i == 5)
f = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
end
if (i == 6)
xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
f = 0.2 * (xdbar + xubar)
end
if (i == 7)
f = 0.0
end
if (i == 8)
f = 0.0
end
if (i == 9)
f = 0.0
end
if (i == 10)
f = 0.0
end
if (i == 11)
f = 0.0
end
if (i == 12)
f = 0.0
end
return f
end
wrapped_func = WrappedPDF(func)FunctionWrappers.FunctionWrapper{Float64, Tuple{Int32, Float64}}(Ptr{Nothing} @0x00007f9ead334a40, Ptr{Nothing} @0x00007f9ee1b60048, Base.RefValue{typeof(Main.func)}(Main.func), typeof(Main.func))Let's set up and run this example evolution.
QCDNUM.qcinit(-6, " ")
nx = QCDNUM.gxmake(xmin, iwt, ng, nxin, iosp)
nq = QCDNUM.gqmake(qq, wt, ngq, nqin)
nw = QCDNUM.fillwt(itype)
QCDNUM.setord(3)
QCDNUM.setalf(as0, r20)
iqc = QCDNUM.iqfrmq(q2c)
iqb = QCDNUM.iqfrmq(q2b)
iqt = QCDNUM.iqfrmq(1e11)
QCDNUM.setcbt(0, iqc, iqb, 0)
iq0 = QCDNUM.iqfrmq(q0)
eps = QCDNUM.evolfg(itype, wrapped_func, def, iq0)
pdf = QCDNUM.allfxq(itype, x, q, 0, 1)
asmz, a, b = QCDNUM.asfunc(qmz2)
csea = 2*pdf[3];
@printf("x, q, CharmSea = %0.4e, %0.4e, %0.4e\n", x, q, csea)
@printf("as(mz2) = %0.4e", asmz) +---------------------------------------------------------------------+
| |
| If you use QCDNUM, please refer to: |
| |
| M. Botje, Comput. Phys. Commun. 182(2011)490, arXiV:1005.1481 |
| |
+---------------------------------------------------------------------+
FILLWT: start unpolarised weight calculations
Subgrids 1 Subgrid points 100
Pij LO
Pij NLO
Pij NNLO
Aij LO
Aij NLO
Aij NNLO
FILLWT: weight calculations completed
x, q, CharmSea = 1.0000e-03, 1.0000e+03, 1.8708e+00
as(mz2) = 1.1807e-01SPLINT
For the SPLINT example, we need a PDF stored and the corresponding iset and ipdf. We will use this to input the function to spline.
iset = itype
ipdf = 1; # component of PDF, -6 to 6
QCDNUM.fvalij(iset, ipdf, 1, 1, 1) # check we get values0.4889149594161791Now that we have run the example program, we can load SPLINT, and write iset and ipdf to the user space.
QCDNUM.ssp_spinit(100)
QCDNUM.ssp_uwrite(1, Float64(iset))
QCDNUM.ssp_uwrite(2, Float64(ipdf))
+---------------------------------------+
| You are using SPLINT version 20220308 |
+---------------------------------------+Make spline object
iasp = QCDNUM.isp_s2make(5, 5)205Define function to read from QCDNUM into spline
function func(ix, iq, first)::Float64
ix = ix[] # deref ptr
iq = iq[]
iset = Int32(QCDNUM.dsp_uread(1))
ipdf = Int32(QCDNUM.dsp_uread(2))
return QCDNUM.fvalij(iset, ipdf, ix, iq, 1)
end
fun = @cfunction(func, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))Ptr{Nothing} @0x00007f9ead337110Fill the spline and set no kinematic limit
QCDNUM.ssp_s2fill(iasp, fun, 0.0)We can use some helper functions to query the spline properties...
QCDNUM.isp_splinetype(iasp)
nu, u1, u2, nv, v1, v2, n = QCDNUM.ssp_splims(iasp)(22, 0.00010000000000000009, 1.0, 14, 2.0, 10000.00000000001, 308)... or copy the nodes locally and print a summary.
xarray = QCDNUM.ssp_unodes(iasp, nu, nu);
qqarray = QCDNUM.ssp_vnodes(iasp, nv, nv);
QCDNUM.ssp_nprint(iasp)
N IX XNODE NQMA IQ QNODE NXMI
1 1 0.10000E-03 14 1 0.20000E+01 1
2 5 0.14454E-03 14 2 0.23106E+01 1
3 10 0.22909E-03 14 7 0.47555E+01 1
4 15 0.36308E-03 14 12 0.97874E+01 1
5 20 0.57544E-03 14 17 0.20144E+02 1
6 25 0.91201E-03 14 22 0.41458E+02 1
7 30 0.14454E-02 14 27 0.85327E+02 1
8 35 0.22909E-02 14 32 0.17561E+03 1
9 40 0.36308E-02 14 37 0.36143E+03 1
10 45 0.57544E-02 14 42 0.74388E+03 1
11 50 0.91201E-02 14 47 0.15310E+04 1
12 55 0.14454E-01 14 52 0.31510E+04 1
13 60 0.22909E-01 14 57 0.64851E+04 1
14 65 0.36308E-01 14 60 0.10000E+05 1
15 70 0.57544E-01 14 - - -
16 75 0.91201E-01 14 - - -
17 80 0.14454E+00 14 - - -
18 85 0.22909E+00 14 - - -
19 90 0.36308E+00 14 - - -
20 95 0.57544E+00 14 - - -
21 100 0.91201E+00 14 - - -
22 101 0.10000E+01 14 - - -There are routines to evaluate the function and its integral at desired x and qq values/ranges.
x = 0.1
q = 100.0
QCDNUM.dsp_funs2(iasp, x, q, 1)
x1 = 0.01
x2 = 0.1
q1 = 10.0
q2 = 100.0
rs = 370.0
np = 4
QCDNUM.dsp_ints2(iasp, x1, x2, q1, q2, rs, np)3.7096736983595426It is also possible to set user nodes, if the automatically chosen ones are not satisfactory.
xarr = Float64.([1e-2, 5e-2, 1e-1, 0.5])
qarr = Float64.([10, 1e2, 1e3, 5e3])
nx = length(xarr)
nq = length(qarr)
iasp = QCDNUM.isp_s2user(xarr, nx, qarr, nq)5562Let's check the update nodes:
QCDNUM.ssp_nprint(iasp)
N IX XNODE NQMA IQ QNODE NXMI
1 51 0.10000E-01 0 12 0.97874E+01 0
2 68 0.47863E-01 0 28 0.98578E+02 0
3 76 0.10000E+00 0 44 0.99286E+03 0
4 93 0.47863E+00 0 55 0.48588E+04 0This page was generated using Literate.jl.