Type consistency will help. R has the unusual convention that 0 is a
floating-point constant. In Julia 0 is an integer. Use '0.0' for a
floating point value in places like
if R <= r0
F0 = log10(r0/R)
else
F0 = 0
end
which could be rewritten
F0 = R <= r0 ? log10(r0/R) : 0.0
Then profile the code execution using @profile.
On Wednesday, December 24, 2014 10:47:41 AM UTC-6, [email protected] wrote:
>
> 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
>
>
