exp performance

density.m

function [density, pts_PCn, di_min] = eval_density(data, points)
[n,d] = size(data);
n_points = size(points,1);

Datan = (data-ones(n,1)*mean(data))./(ones(n,1)*std(data));
[U,S,V] = svd(Datan,0);
PC = Datan*V;
PCn = PC./(ones(n,1)*std(PC));

Pointsn = (points-ones(n_points,1)*mean(data))./(ones(n_points,1)*std(data));
Points_PC = Pointsn*V;
pts_PCn = Points_PC./(ones(n_points,1)*std(PC));


h_bar = n^(-1/(d+4));
Constant=1/n/(2*pi)^(d/2)/h_bar^d;

density = zeros(n_points,1);
di_min = zeros(n_points,1);

for i = 1:n_points
    Di_2 = (ones(n,1)*pts_PCn(i,:)-PCn).^2*ones(d,1);
    density(i,1) = Constant*ones(1,n)*exp(-0.5*Di_2/h_bar^2);
    Di_2(i,1) = inf;
    di_min(i,1) = sqrt(min(Di_2));
end

end

eval_density_bcast_prof

The following results were obtained by running this julia code:

```
julia> Profile.clear()

julia> Profile.@profile eval_density_bcast(x,x);

julia> Profile.@profile eval_density_bcast(x,x);

julia> save_prof("eval_density_bcast_prof.txt")
```

---

Profile results (note call to myexp is on line 135):


```
1    abstractarray.jl; ndims; line: 60
1    broadcast.jl; broadcast!; line: 220
1    reducedim.jl; _mapreducedim!; line: 186
2    reducedim.jl; reduced_dims; line: 10
4254 task.jl; anonymous; line: 89
 4254 REPL.jl; eval_user_input; line: 60
  4254 profile.jl; anonymous; line: 14
   1    /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 115
    1 /Users/sglyon/Desktop/perf_exp.jl; call; line: 17
     1 broadcast.jl; broadcast!; line: 227
      1 broadcast.jl; _F_; line: 32
   1    /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 120
    1 /Users/sglyon/Desktop/perf_exp.jl; call; line: 38
     1 /Users/sglyon/Desktop/perf_exp.jl; normalize; line: 9
      1 broadcast.jl; broadcast!; line: 227
       1 broadcast.jl; _F_; line: 96
   1    /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 121
    1 array.jl; transpose!; line: 1343
     1 array.jl; transposeblock!; line: 1371
      1 array.jl; transposeblock!; line: 1366
       1 array.jl; transposeblock!; line: 1371
        1 array.jl; transposeblock!; line: 1371
         1 array.jl; transposeblock!; line: 1366
          1 array.jl; transposeblock!; line: 1366
           1 array.jl; transposeblock!; line: 1371
            1 array.jl; transposeblock!; line: 1371
             1 array.jl; transposeblock!; line: 1366
              1 array.jl; transposeblock!; line: 1360
   3535 /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 134
    1374 array.jl; .^; line: 758
     2  array.jl; reshape; line: 137
     67 array.jl; reshape; line: 140
    67   broadcast.jl; broadcast!; line: 220
    1371 broadcast.jl; broadcast!; line: 227
     188 broadcast.jl; _F_; line: 32
     34  broadcast.jl; _F_; line: 88
      10 broadcast.jl; check_broadcast_shape; line: 51
       1 abstractarray.jl; ndims; line: 60
      9  broadcast.jl; check_broadcast_shape; line: 54
      8  broadcast.jl; check_broadcast_shape; line: 55
       2 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; getindex; (unknown line)
      3  broadcast.jl; check_broadcast_shape; line: 56
      2  broadcast.jl; check_broadcast_shape; line: 57
     161 broadcast.jl; _F_; line: 95
     980 broadcast.jl; _F_; line: 96
     1   broadcast.jl; _F_; line: 189
    3    broadcast.jl; broadcast!; line: 643
     2 dict.jl; ht_keyindex2; line: 544
     1 dict.jl; ht_keyindex2; line: 545
    1    broadcast.jl; broadcast!; line: 651
    2    broadcast.jl; broadcast_shape; line: 30
    1    broadcast.jl; broadcast_shape; line: 31
    2    broadcast.jl; broadcast_shape; line: 32
     2 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; fill!; (unknown line)
    1    broadcast.jl; broadcast_shape; line: 33
    24   broadcast.jl; broadcast_shape; line: 34
    2    broadcast.jl; broadcast_shape; line: 35
    3    broadcast.jl; broadcast_shape; line: 36
    1    broadcast.jl; broadcast_shape; line: 37
    3    broadcast.jl; broadcast_shape; line: 38
    2    broadcast.jl; broadcast_shape; line: 39
     1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; !=; (unknown line)
    5    broadcast.jl; broadcast_shape; line: 40
    28   broadcast.jl; broadcast_shape; line: 45
     8 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; append_any; (unknown line)
    4    multidimensional.jl; anonymous; line: 287
     2 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; checkbounds; (unknown line)
    7    multidimensional.jl; anonymous; line: 288
     2 multidimensional.jl; anonymous; line: 263
      2 multidimensional.jl; anonymous; line: 231
    566  reducedim.jl; mapreducedim; line: 238
     1   array.jl; fill!; line: 204
     87  reducedim.jl; _mapreducedim!; line: 65
     45  reducedim.jl; _mapreducedim!; line: 72
     1   reducedim.jl; _mapreducedim!; line: 186
      1 reducedim.jl; check_reducedims; line: 162
     54  reducedim.jl; _mapreducedim!; line: 202
     17  reducedim.jl; _mapreducedim!; line: 203
     57  reducedim.jl; _mapreducedim!; line: 204
     192 reducedim.jl; _mapreducedim!; line: 206
     33  reducedim.jl; _mapreducedim!; line: 208
     44  reducedim.jl; reduced_dims; line: 13
      44 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; getindex; (unknown line)
       1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; call; (unknown line)
     16  reducedim.jl; reduced_dims; line: 17
   642  /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 135
    520 operators.jl; myexp; line: 358
    1   reduce.jl; _mapreduce; line: 145
    121 reduce.jl; _mapreduce; line: 149
     60 reduce.jl; mapreduce_pairwise_impl; line: 109
      34 reduce.jl; mapreduce_pairwise_impl; line: 109
       19 reduce.jl; mapreduce_pairwise_impl; line: 109
        9  reduce.jl; mapreduce_pairwise_impl; line: 109
         9 reduce.jl; mapreduce_pairwise_impl; line: 106
          9 simdloop.jl; mapreduce_seq_impl; line: 171
        10 reduce.jl; mapreduce_pairwise_impl; line: 110
         10 reduce.jl; mapreduce_pairwise_impl; line: 106
          10 simdloop.jl; mapreduce_seq_impl; line: 171
       15 reduce.jl; mapreduce_pairwise_impl; line: 110
        1  reduce.jl; mapreduce_pairwise_impl; line: 106
        11 reduce.jl; mapreduce_pairwise_impl; line: 109
         11 reduce.jl; mapreduce_pairwise_impl; line: 106
          11 simdloop.jl; mapreduce_seq_impl; line: 171
        3  reduce.jl; mapreduce_pairwise_impl; line: 110
         3 reduce.jl; mapreduce_pairwise_impl; line: 106
          3 simdloop.jl; mapreduce_seq_impl; line: 171
      26 reduce.jl; mapreduce_pairwise_impl; line: 110
       13 reduce.jl; mapreduce_pairwise_impl; line: 109
        10 reduce.jl; mapreduce_pairwise_impl; line: 109
         10 reduce.jl; mapreduce_pairwise_impl; line: 106
          10 simdloop.jl; mapreduce_seq_impl; line: 171
        3  reduce.jl; mapreduce_pairwise_impl; line: 110
         3 reduce.jl; mapreduce_pairwise_impl; line: 106
          3 simdloop.jl; mapreduce_seq_impl; line: 171
       13 reduce.jl; mapreduce_pairwise_impl; line: 110
        7 reduce.jl; mapreduce_pairwise_impl; line: 109
         7 reduce.jl; mapreduce_pairwise_impl; line: 106
          7 simdloop.jl; mapreduce_seq_impl; line: 171
        6 reduce.jl; mapreduce_pairwise_impl; line: 110
         6 reduce.jl; mapreduce_pairwise_impl; line: 106
          6 simdloop.jl; mapreduce_seq_impl; line: 171
     61 reduce.jl; mapreduce_pairwise_impl; line: 110
      24 reduce.jl; mapreduce_pairwise_impl; line: 109
       14 reduce.jl; mapreduce_pairwise_impl; line: 109
        9 reduce.jl; mapreduce_pairwise_impl; line: 109
         9 reduce.jl; mapreduce_pairwise_impl; line: 106
          9 simdloop.jl; mapreduce_seq_impl; line: 171
        5 reduce.jl; mapreduce_pairwise_impl; line: 110
         5 reduce.jl; mapreduce_pairwise_impl; line: 106
          5 simdloop.jl; mapreduce_seq_impl; line: 171
       10 reduce.jl; mapreduce_pairwise_impl; line: 110
        4 reduce.jl; mapreduce_pairwise_impl; line: 109
         4 reduce.jl; mapreduce_pairwise_impl; line: 106
          4 simdloop.jl; mapreduce_seq_impl; line: 171
        6 reduce.jl; mapreduce_pairwise_impl; line: 110
         6 reduce.jl; mapreduce_pairwise_impl; line: 106
          6 simdloop.jl; mapreduce_seq_impl; line: 171
      37 reduce.jl; mapreduce_pairwise_impl; line: 110
       21 reduce.jl; mapreduce_pairwise_impl; line: 109
        10 reduce.jl; mapreduce_pairwise_impl; line: 109
         10 reduce.jl; mapreduce_pairwise_impl; line: 106
          10 simdloop.jl; mapreduce_seq_impl; line: 171
        11 reduce.jl; mapreduce_pairwise_impl; line: 110
         11 reduce.jl; mapreduce_pairwise_impl; line: 106
          11 simdloop.jl; mapreduce_seq_impl; line: 171
       16 reduce.jl; mapreduce_pairwise_impl; line: 110
        11 reduce.jl; mapreduce_pairwise_impl; line: 109
         11 reduce.jl; mapreduce_pairwise_impl; line: 106
          11 simdloop.jl; mapreduce_seq_impl; line: 171
        5  reduce.jl; mapreduce_pairwise_impl; line: 110
         5 reduce.jl; mapreduce_pairwise_impl; line: 106
          5 simdloop.jl; mapreduce_seq_impl; line: 171
   74   /Users/sglyon/Desktop/perf_exp.jl; eval_density_bcast; line: 137
    73 reduce.jl; _mapreduce; line: 149
     38 reduce.jl; mapreduce_impl; line: 252
     34 reduce.jl; mapreduce_impl; line: 254
```

eval_density_ones_prof

The following results were obtained by running this julia code:

```
julia> Profile.clear()

julia> Profile.@profile eval_density_ones(x,x);

julia> Profile.@profile eval_density_ones(x,x);

julia> save_prof("eval_density_ones_prof.txt")
```

---

Profile Results (note call to my_exp is on line 174):

```
1    /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; _start; (unknown line)
 1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; run_repl; (unknown line)
  1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; run_frontend; (unknown line)
   1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; run_interface; (unknown line)
    1 REPL.jl; anonymous; line: 602
     1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; transition; (unknown line)
      1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; activate; (unknown line)
       1 LineEdit.jl; refresh_multi_line; line: 212
2    abstractarray.jl; trailingsize; line: 97
1    linalg/matmul.jl; gemm_wrapper!; line: 269
3169 task.jl; anonymous; line: 89
 3169 REPL.jl; eval_user_input; line: 60
  3169 profile.jl; anonymous; line: 14
   1    /Users/sglyon/Desktop/perf_exp.jl; eval_density_ones; line: 155
    1 /Users/sglyon/Desktop/perf_exp.jl; call; line: 38
     1 /Users/sglyon/Desktop/perf_exp.jl; normalize; line: 9
      1 broadcast.jl; broadcast!; line: 220
   2770 /Users/sglyon/Desktop/perf_exp.jl; eval_density_ones; line: 173
    113  array.jl; -; line: 791
    1133 array.jl; -; line: 793
    1389 array.jl; .^; line: 758
     68 array.jl; reshape; line: 140
    83   linalg/matmul.jl; *; line: 76
     1  array.jl; reshape; line: 140
     1  linalg/matmul.jl; gemm_wrapper!; line: 269
      1 linalg/matmul.jl; lapack_size; line: 288
     4  linalg/matmul.jl; gemm_wrapper!; line: 278
     66 linalg/matmul.jl; gemm_wrapper!; line: 281
      65 linalg/blas.jl; gemm!; line: 560
    1    linalg/matmul.jl; gemv!; line: 210
    23   linalg/matmul.jl; gemv!; line: 214
     23 linalg/blas.jl; gemv!; line: 269
    2    multidimensional.jl; anonymous; line: 287
     1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; checkbounds; (unknown line)
    6    multidimensional.jl; anonymous; line: 288
     1 multidimensional.jl; anonymous; line: 263
     1 multidimensional.jl; anonymous; line: 269
   288  /Users/sglyon/Desktop/perf_exp.jl; eval_density_ones; line: 174
    2   linalg/matmul.jl; gemv!; line: 214
     1 /usr/local/Cellar/julia/HEAD/lib/julia/sys.dylib; stride; (unknown line)
     1 linalg/blas.jl; gemv!; line: 269
    210 operators.jl; myexp; line: 357
   110  /Users/sglyon/Desktop/perf_exp.jl; eval_density_ones; line: 176
    108 reduce.jl; _mapreduce; line: 149
     57 reduce.jl; mapreduce_impl; line: 253
     51 reduce.jl; mapreduce_impl; line: 254

```

eval_density_prof.txt

The following results were obtained by running this julia code:

```
julia> Profile.clear()

julia> Profile.@profile eval_density(x,x);

julia> Profile.@profile eval_density(x,x);

julia> save_prof("eval_density_prof.txt")
```

---

Profile results (note call to myexp is on line 96):

757 task.jl; anonymous; line: 89
 757 REPL.jl; eval_user_input; line: 60
  757 profile.jl; anonymous; line: 14
   1   /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 70
    1 /Users/sglyon/Desktop/perf_exp.jl; call; line: 17
     1 broadcast.jl; broadcast!; line: 227
      1 broadcast.jl; _F_; line: 32
   40  /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 91
   488 /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 94
   102 /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 96
   82  /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 97
   43  /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 98
   1   /Users/sglyon/Desktop/perf_exp.jl; eval_density; line: 101

julia_repl_timing

julia> reload("perf_exp.jl")

julia> x = randn(9500, 2);

julia> eval_density(x,x); eval_density_bcast(x,x); eval_density_ones(x,x);

julia> @time eval_density(x,x);
elapsed time: 8.397094688 seconds (1 MB allocated)

julia> @time eval_density_bcast(x,x);
elapsed time: 10.203116128 seconds (4147 MB allocated, 2.90% gc time in 190 pauses with 0 full sweep)

julia> @time eval_density_ones(x,x);
elapsed time: 10.545925151 seconds (5518 MB allocated, 3.71% gc time in 251 pauses with 0 full sweep)

matlab_timing

>> x = randn(9500,2);
>> tic; eval_density(x,x); toc;
Elapsed time is 2.610879 seconds.
>> tic; eval_density(x,x); toc;
Elapsed time is 2.072053 seconds.

perf_exp.jl

# PCA is a type I need in my function below
normalize(x::Matrix, dim::Int=1) = (x .- mean(x, dim)) ./ std(x, dim)

"""
Normalize `x` by subtracting mean of `other` and dividing by the
standard deviation of `other`. The mean and standard deviation are computed
along dimension `dim`, which has a default value of `1`.
"""
normalize(x::Matrix, other::Matrix, dim::Int=1) =
    (x .- mean(other, dim)) ./ std(other, dim)

"Construct a `PCA` object, given data"
function PCA(data::Matrix)
    datan = normalize(data)   # normalize data
    U, S, V = svd(datan)      # compute svd
    PC = datan*V              # compute principal components (PCs)
    PCn = (PC ./ std(PC, 1))  # normalize PCs
    PCA(data, V, PC, PCn)     # return object
end


"""
Given a matrix `data` and a `PCA` representation `other`, compute a
principal components representation of `data` using the `svd` of other
"""
function PCA(data::Matrix, other::PCA)
    datan = normalize(data, other.data)
    PC = datan*other.V
    PCn = (PC ./ std(other.PC, 1))
    return PCA(data, other.V, PC, PCn)
end

"""
Given a matrix `data` and a `PCA` representation `other`, compute a
principal components representation of `data` using the `svd` of other
"""
function PCA(data::Matrix, other::PCA)
    datan = normalize(data, other.data)
    PC = datan*other.V
    PCn = (PC ./ std(other.PC, 1))
    return PCA(data, other.V, PC, PCn)
end

# helper function to save profile results
function save_prof(f_name="profile_resuts.txt")
    s = open(f_name,"w")
    Profile.print(s,cols = 500)
    close(s)
end

# Somehow the apple libm is faster than the one julia calls on OSX
@osx? (
         begin
             myexp(x::Float64) = ccall((:exp, :libm), Float64, (Float64,), x)
         end
       : begin
             myexp(x::Float64) = exp(x)
         end
       )

@vectorize_1arg Number myexp

function eval_density(data::Matrix, points::Matrix)

    n, d = size(data)

    n_points = size(points, 1)

    # transform data
    data_PCA = PCA(data)
    PC = data_PCA.PC
    PCn = data_PCA.PCn'  # do `'` so I can operate on columns at a time

    # transform points
    points_PCA = PCA(points, data_PCA)
    pts_PCn = points_PCA.PCn'  # same comment about `'`

    # Define constants
    hbar = n^(-1.0/(d+4.0))
    hbar2 = hbar^2
    constant = myexp(-0.5/hbar2)/(n*hbar^(d) * (2π)^(d/2))

    # allocate space
    density = Array(Float64, n_points)
    Di_min = Array(Float64, n_points)

    # apply formula (2)
    for i=1:n_points  # loop over all points
        dens_i = 0.0
        min_di2 = Inf
        for j=1:n_points  # loop over all other points
            d_i2_j = 0.0
            for k=1:d  # loop over d
                @inbounds d_i2_j += (pts_PCn[k, i] - PCn[k, j])^2
            end
            dens_i += myexp(d_i2_j)
            if i != j && d_i2_j < min_di2
                min_di2 = d_i2_j
            end
        end
        density[i] = constant * dens_i
        Di_min[i] = sqrt(min_di2)
    end

    return density, points_PCA.PCn, Di_min
end

function eval_density_bcast(data::Matrix, points::Matrix)

    n, d = size(data)

    n_points = size(points, 1)

    # transform data
    data_PCA = PCA(data)
    PC = data_PCA.PC
    PCn = data_PCA.PCn'  # do `'` so I can operate on columns at a time

    # transform points
    points_PCA = PCA(points, data_PCA)
    pts_PCn = points_PCA.PCn'  # same comment about `'`

    # Define constants
    hbar = n^(-1.0/(d+4.0))
    hbar2 = hbar^2
    constant = myexp(-0.5/hbar2)/(n*hbar^(d) * (2π)^(d/2))

    # allocate space
    density = Array(Float64, n_points)
    Di_min = Array(Float64, n_points)

    # apply formula (2)
    for i=1:n_points
        D_i2 = sum((pts_PCn[:, i] .- PCn).^2, 1)
        density[i] = constant * sum(myexp(D_i2))
        D_i2[i] = Inf
        Di_min[i] = sqrt(minimum(D_i2))
    end

    return density, points_PCA.PCn, Di_min
end

function eval_density_ones(data::Matrix, points::Matrix)

    n, d = size(data)

    n_points = size(points, 1)

    # transform data
    data_PCA = PCA(data)
    PC = data_PCA.PC
    PCn = data_PCA.PCn

    # transform points
    points_PCA = PCA(points, data_PCA)
    pts_PCn = points_PCA.PCn

    # Define constants
    hbar = n^(-1.0/(d+4.0))
    constant = myexp(-0.5/(hbar^2))/(n*hbar^(d) * (2π)^(d/2))

    # allocate space
    density = Array(Float64, n_points)
    Di_min = Array(Float64, n_points)

    # pre-allocate vectors of ones so I don't re-create them each loop
    ones_n = ones(n)
    ones_d = ones(d)
    ones_1n = ones(1, n)

    # apply formula (2)
    for i=1:n_points
        D_i2 = (ones_n*pts_PCn[i, :] - PCn).^2 * ones_d
        density[i] = constant * (ones_1n*myexp(D_i2))[1]  # get element of 1×1 matrix
        D_i2[i] = Inf
        Di_min[i] = sqrt(minimum(D_i2))
    end

    return density, points_PCA.PCn, Di_min
end