Hi All,
I have just started coding a few functions in JULIA with the hope I can
gain some speed compared with my existing R code.
However the code I have is running 5 times slower than my R code, so I am
guessing I have written this very poorly!!
The function takes arguments which are arrays, and carries out a simple
calculation. In R my code is vectorised nicely. I have written this JULIA
code using a loop as I believe this is where JULIA excels??
So please can anybody give me some ideas on how to achieve maximum warp
speed with JULIA. I am also interested in how to make the loop parallel for
more speed. Thank you in advance.
If I cannot convince myself of the speed gains with a simple function like
this I may just keep life simple and stay with R!!
Currently on my machine the below in JULIA takes 31.5secs and my R code
takes about 6secs.
Here are my arguments and function code. Thanks for all your help in
advance.
T = 0
Mw = ones(100000000)*7
Rjb = 100
Eps = 0
Magdep = 1
function RFrb( T, Mw, Rjb, Eps, Magdep)
if Magdep == 0
#Table 5 Coefficients (Magnitude Independant, Stress Parameter Model)
Ti = [0, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, 0.16, 0.2, 0.25, 0.31, 0.4,
0.5, 0.63, 0.79, 1, 1.25, 1.59, 2, 2.5, 3.13, 4, 5]
c1 = [-0.0135, 0.8282, 0.4622, 0.2734, 0.0488, -0.2112, -0.5363, -0.9086,
-1.3733, -1.918, -2.5107, -3.1571, -3.8516, -4.5556, -5.2405, -5.8909, -6.4633,
-6.925, -7.296, -7.5053, -7.5569, -7.451, -7.1688, -6.8063]
c2 = [0.6889, 0.5976, 0.6273, 0.6531, 0.6945, 0.7517, 0.8319, 0.93, 1.0572,
1.2094, 1.3755, 1.5549, 1.7429, 1.9258, 2.0926, 2.2357, 2.3419, 2.4037, 2.4189,
2.3805, 2.2933, 2.1598, 1.9738, 1.7848]
c3 = [-0.0488, -0.0418, -0.0391, -0.0397, -0.042, -0.046, -0.0521, -0.0597,
-0.0698, -0.0819, -0.0949, -0.1087, -0.1228, -0.136, -0.1471, -0.1557, -0.1605,
-0.1612, -0.1573, -0.1492, -0.1376, -0.1228, -0.1048, -0.0879]
c4 = [-1.8987, -2.1321, -1.7242, -1.5932, -1.4913, -1.4151, -1.3558, -1.309,
-1.2677, -1.2315, -1.1992, -1.1677, -1.1354, -1.1015, -1.0659, -1.0279,
-0.9895, -0.9545, -0.9247, -0.9128, -0.9285, -0.9872, -1.1274, -1.3324]
c5 = [0.2151, 0.2159, 0.1644, 0.1501, 0.1405, 0.134, 0.1296, 0.1264, 0.1237,
0.1213, 0.1189, 0.116, 0.1126, 0.1084, 0.1035, 0.0981, 0.0925, 0.0879, 0.0848,
0.0855, 0.0915, 0.105, 0.1325, 0.1691]
c6 = [-1.9063, -2.053, -1.6849, -1.5698, -1.4807, -1.413, -1.3579, -1.312,
-1.2684, -1.227, -1.1881, -1.1494, -1.1099, -1.0708, -1.0328, -0.9969, -0.9665,
-0.9462, -0.9421, -0.9658, -1.0264, -1.1349, -1.3132, -1.5158]
c7 = [0.174, 0.1676, 0.127, 0.1161, 0.1084, 0.1027, 0.0985, 0.0948, 0.091,
0.0872, 0.0833, 0.0791, 0.0746, 0.07, 0.0655, 0.0612, 0.0577, 0.0558, 0.0567,
0.0619, 0.0729, 0.0914, 0.1207, 0.1533]
c8 = [-2.0131, -1.5148, -1.4513, -1.535, -1.6563, -1.7821, -1.8953, -1.9863,
-2.0621, -2.1196, -2.1598, -2.1879, -2.2064, -2.2171, -2.222, -2.2229, -2.2211,
-2.2178, -2.2137, -2.211, -2.2108, -2.2141, -2.2224, -2.2374]
c9 = [0.0887, 0.1163, 0.091, 0.0766, 0.0657, 0.0582, 0.052, 0.0475, 0.0434,
0.0396, 0.0361, 0.0328, 0.0294, 0.0261, 0.0229, 0.0197, 0.0167, 0.0139, 0.0111,
0.0086, 0.0067, 0.006, 0.0079, 0.0142]
c10 = [-0.002747, -0.004463, -0.004355, -0.003939, -0.003449, -0.002987,
-0.002569, -0.002234, -0.001944, -0.001708, -0.001522, -0.001369, -0.00124,
-0.001129, -0.001033, -0.000945, -0.000863, -0.000785, -0.000701, -0.000618,
-0.000535, -0.000458, -0.000397, -0.000387]
c11 = [1.5473, 1.1096, 1.1344, 1.1493, 1.2154, 1.2858, 1.3574, 1.426, 1.4925,
1.5582, 1.6049, 1.6232, 1.632, 1.6109, 1.5735, 1.5262, 1.4809, 1.471, 1.5183,
1.6365, 1.8421, 2.1028, 2.4336, 2.6686]
sigT = [0.436, 0.449, 0.445, 0.442, 0.438, 0.433, 0.428, 0.422, 0.416,
0.409, 0.402, 0.395, 0.387, 0.378, 0.369, 0.36, 0.35, 0.341, 0.331,
0.323, 0.315, 0.308, 0.299, 0.291]
sigB = [0.409, 0.417, 0.417, 0.416, 0.414, 0.41, 0.405, 0.399, 0.392,
0.384, 0.376, 0.366, 0.356, 0.345, 0.333, 0.32, 0.307, 0.294, 0.28,
0.267, 0.254, 0.242, 0.227, 0.214]
sigW = [0.153, 0.167, 0.155, 0.149, 0.143, 0.14, 0.138, 0.138, 0.139,
0.141, 0.144, 0.148, 0.152, 0.156, 0.16, 0.164, 0.168, 0.172, 0.177,
0.181, 0.186, 0.19, 0.195, 0.198]
else
#Table 6 Coefficients (Magnitude Dependant, Stress Parameter Model)
Ti = [0, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, 0.16, 0.2, 0.25, 0.31, 0.4,
0.5, 0.63, 0.79, 1, 1.25, 1.59, 2, 2.5, 3.13, 4, 5]
c1 = [-2.6934, -1.9654, -2.3216, -2.4879, -2.6647, -2.8579, -3.092, -3.3595,
-3.6974, -4.1004, -4.5462, -5.0372, -5.565, -6.0933, -6.5914, -7.0402,
-7.4028, -7.6577, -7.8128, -7.8368, -7.7341, -7.4991, -7.1376, -6.7757]
c2 = [1.7682, 1.7265, 1.7592, 1.7771, 1.8009, 1.8325, 1.8771, 1.9336, 2.0096,
2.1032, 2.2068, 2.318, 2.4307, 2.5325, 2.6123, 2.6616, 2.6715, 2.6402, 2.5609,
2.444, 2.2967, 2.1232, 1.9232, 1.7502]
c3 = [-0.1366, -0.1346, -0.1328, -0.1328, -0.1336, -0.1353, -0.1381, -0.1419,
-0.1472, -0.1537, -0.1607, -0.1679, -0.1745, -0.1795, -0.182, -0.1813, -0.1767,
-0.1686, -0.1559, -0.1409, -0.1244, -0.1072, -0.0893, -0.0753]
c4 = [-1.8544, -2.1011, -1.7198, -1.591, -1.4916, -1.4169, -1.3587, -1.3133,
-1.2717, -1.2341, -1.1985, -1.1625, -1.1242, -1.0837, -1.0432, -1.0023,
-0.9634, -0.9299, -0.9021, -0.8896, -0.9012, -0.9638, -1.1238, -1.3603]
c5 = [0.2123, 0.2101, 0.1635, 0.1502, 0.1409, 0.1349, 0.1312, 0.1289, 0.1262,
0.1234, 0.1198, 0.1155, 0.1104, 0.1043, 0.0985, 0.0927, 0.0877, 0.0842, 0.0829,
0.0857, 0.0975, 0.1202, 0.1549, 0.196]
c6 = [-1.8467, -2.0063, -1.6582, -1.5447, -1.4576, -1.3909, -1.3363, -1.2908,
-1.2472, -1.2055, -1.1658, -1.1258, -1.0847, -1.0436, -1.0044, -0.9683,
-0.9394, -0.9226, -0.9242, -0.9546, -1.027, -1.1604, -1.3647, -1.5804]
c7 = [0.159, 0.1555, 0.1192, 0.1087, 0.1013, 0.0959, 0.0919, 0.0884, 0.0847,
0.0809, 0.077, 0.0728, 0.0682, 0.0635, 0.0591, 0.0552, 0.0526, 0.0519, 0.0544,
0.0613, 0.0747, 0.0965, 0.1284, 0.1613]
c8 = [-1.8809, -1.3684, -1.3348, -1.4304, -1.5603, -1.6918, -1.8091, -1.9038,
-1.9844, -2.0474, -2.0935, -2.1276, -2.1519, -2.168, -2.1775, -2.182, -2.1827,
-2.1811, -2.1782, -2.1751, -2.1717, -2.1763, -2.1901, -2.207]
c9 = [0.0681, 0.0946, 0.073, 0.0599, 0.0494, 0.0423, 0.0363, 0.0322, 0.0286,
0.0254, 0.0227, 0.0201, 0.0176, 0.0152, 0.0129, 0.0109, 0.0092, 0.0079, 0.0069,
0.0064, 0.0066, 0.0077, 0.0114, 0.0191]
c10 = [-0.002888, -0.004626, -0.004488, -0.004056, -0.003549, -0.003077,
-0.002651, -0.002311, -0.002015, -0.001772, -0.001579, -0.001419, -0.001282,
-0.001166, -0.001066, -0.000977, -0.000898, -0.000827, -0.000756, -0.000691,
-0.000634, -0.000571, -0.000515, -0.000518]
c11 = [2.1589, 1.3437, 1.3363, 1.3942, 1.4587, 1.5466, 1.6583, 1.7807, 1.854,
1.9055, 1.9052, 1.8732, 1.8142, 1.7238, 1.6524, 1.59, 1.5549, 1.5712, 1.6701,
1.9205, 2.6233, 3.5221, 4.0984, 4.3313]
sigT = [0.335, 0.348, 0.343, 0.339, 0.335, 0.33, 0.325, 0.32, 0.314, 0.307,
0.301, 0.294, 0.288, 0.282, 0.276, 0.272, 0.268, 0.265, 0.262, 0.26, 0.258,
0.256, 0.253, 0.251]
sigB = [0.298, 0.304, 0.305, 0.304, 0.302, 0.299, 0.295, 0.289, 0.282, 0.273,
0.263, 0.253, 0.242, 0.231, 0.22, 0.21, 0.201, 0.193, 0.185, 0.177, 0.169,
0.161, 0.152, 0.144]
sigW = [0.154, 0.169, 0.157, 0.15, 0.144, 0.14, 0.138, 0.137, 0.138, 0.141,
0.145, 0.15, 0.155, 0.161, 0.167, 0.172, 0.177, 0.181, 0.186, 0.191, 0.195,
0.199, 0.203, 0.205]
end
tidLO = find( Ti .<= T )
tidHI = find( Ti .>= T )
tidLO = tidLO[length(tidLO)]
tidHI = tidHI[1]
if tidLO == tidHI # we have a known period
tid = tidHI
Ti = Ti[tid]
c1 = c1[tid]
c2 = c2[tid]
c3 = c3[tid]
c4 = c4[tid]
c5 = c5[tid]
c6 = c6[tid]
c7 = c7[tid]
c8 = c8[tid]
c9 = c9[tid]
c10 = c10[tid]
c11 = c11[tid]
sigT = sigT[tid]
sigB = sigB[tid]
sigW = sigW[tid]
Sa = zeros(length(Mw))
for i = 1:length(Mw)
R = sqrt(Rjb^2 + c11^2)
r0 = 10
r1 = 50
r2 = 200
if R <= r0
F0 = log10(r0/R)
else
F0 = 0
end
if R <= r1
F1 = log10(R/1.0)
else
F1 = log10(r1/1.0)
end
if R <= r2
F2 = 0
else
F2 = log10(R/r2)
end
return Sa[i] = 10^(c1 + c2*Mw[i] + c3*Mw[i]^2 + (c4+c5*Mw[i])*F0 +
(c6+c7*Mw[i])F1 + (c8+c9*Mw[i])*F2 + c10*R + Eps*sigT) / 981
end
