main.jl 19.8 KB
Newer Older
1
2
using DataFrames
using CSV
3
using JSON
4
using LightGraphs
5
using NearestNeighbors
6
using Statistics
7
using BenchmarkTools
8
using Plots
9
using LinearAlgebra
10

11

12
13
14
15
16
17
18
struct material
    α::Number
    E::Number
    σ̄::AbstractArray{Number, 2} # Yield stress \sigma\bar
end


Jayant Khatkar's avatar
Jayant Khatkar committed
19
mutable struct contour
20
21
22
23
24
    pos
    time
end


25
struct contourdata
26
27
28
    contours::Vector{contour}
    G::SimpleDiGraph
    layers::Vector
29
    travel_dists::Dict
30
    layer_height::Number
31
32
33
end


34
35
36
struct voxmap
    seglen::Vector{Float64}
    segcontours::Vector{Int}
37
    c::Number
38
39
40
41
42
43
44
45
46
47
48
end


struct voxdata
    voxels::DataFrame
    maps::Vector{voxmap}
    below::Vector{Int}
    width::Number
end


49
50
51
52
53
54
55
56
57
58
59
60
61
62
function vecvec_to_matrix(vecvec)
    # convert vector of vectors int a matrix
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(Float32, dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end


63
function contour(d::Dict)
64
    return contour(vecvec_to_matrix(d["pos"]), d["time"])
65
end
66

67

68
function contourdata(cons::Vector{contour}, max_layers::Int, min_dist::Number)
69
    G = LightGraphs.SimpleDiGraph(0)
70

71
    # separate contours into layers
72
    layer_heights = sort(collect(Set([c.pos[end,3] for c in cons])))
73
74
75
76
77
78
79
80
81
82
83
84
85
    layers = [[] for i in 1:length(layer_heights)]
    clayeri = []
    contour_trees = []

    # place contours in layers and construct KDTree for each contour
    for i in 1:length(cons)
        l = searchsorted(layer_heights, cons[i].pos[1,3])[1]
        push!(layers[l], i)
        push!(clayeri, l)
        add_vertex!(G)
        push!(contour_trees, KDTree(transpose(cons[i].pos)))
    end

86
    # loop through contours from previous layer and compare waypoints
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    for i in 1:length(cons)
        l = clayeri[i]

        # add contours from max_layers below
        if l > max_layers
            for c in layers[l-max_layers]
                add_edge!(G, c, i)
            end
        end

        if l == 1 || max_layers == 1
            continue
        end

        for c in layers[l-1]
            # if any points in contour i within min_dist of any points in contour c
            if any([length(b) > 0 for b in inrange(contour_trees[c], transpose(cons[i].pos), min_dist)])
                add_edge!(G, c, i) # mark i dependent on c
            end
        end
    end
108

109
    return contourdata(cons, G, layers, Dict(), layer_heights[1])
110
111
112
113
114
115
end


function seg_helper_orientation(p,q,r)
    val = (q[2]-p[2]) * (r[1]-q[1]) - (q[1]-p[1]) * (r[2]-q[2])
    if val > 0
116
        return 1 # clockwise
117
    elseif val < 0
118
        return 2 # anticlockwise
119
    else
120
        return 0 # colinear
121
122
123
124
    end
end


125
126
127
128
129
130
131
function onseg(p,q,r)
    # check if q lies on segment pr assuming 3 points are colinear
    return ((q[1] <= max(p[1], r[1])) && (q[1] >= min(p[1], r[1])) &&
    (q[2] <= max(p[2], r[2])) && (q[2] >= min(p[2], r[2])))
end


132
133
134
135
136
137
138
function seg_intersect(p1,q1,p2,q2)

    o1 = seg_helper_orientation(p1, q1, p2)
    o2 = seg_helper_orientation(p1, q1, q2)
    o3 = seg_helper_orientation(p2, q2, p1)
    o4 = seg_helper_orientation(p2, q2, q1)

139
140
141
142
143
    if (o1  o2) && (o3  o4) ||
        o1==0 && onseg(p1, p2, q1) ||
        o2==0 && onseg(p1, q2, q1) ||
        o3==0 && onseg(p2, p1, q2) ||
        o4==0 && onseg(p2, q1, q2)
144
145
146
147
148
149
        return true
    end

    return false
end

150

151
152
153
dist(p1, p2) = √sum((p1 -p2).^2)
interpolate(p1, p2, xi, axis) = p1 + (p2-p1)*(xi-p1[axis])/(p2[axis]-p1[axis])
interpolate(p1, p2, x1, xi, x2) = p1 + (p2-p1)*(xi-x1)/(x2-x1)
154

155

156
function voxmap(vox::Vector{Float64}, vox_d::Number, cdata::contourdata)
157
158
159

    # for one vox, get all contours which pass through it
    # only need to search contours in its layer
160
    l = Int(round((vox[3] + cdata.layer_height/2)/cdata.layer_height))
161
162
163
164
    voxx1 = vox[1] + vox_d/2
    voxx2 = vox[1] - vox_d/2
    voxy1 = vox[2] + vox_d/2
    voxy2 = vox[2] - vox_d/2
165

166
    seg_now = false
167
168
169
170
171
172
    seglen = Vector{Number}()
    segoffset = Vector{Number}()
    segcontours = Vector{Int}()
    seglen_sofar = 0
    t_start = 0

Jayant Khatkar's avatar
Jayant Khatkar committed
173
    if l > length(cdata.layers)
174
        return voxmap(seglen, segcontours, 0)
Jayant Khatkar's avatar
Jayant Khatkar committed
175
176
    end

177
    for cid in cdata.layers[l]
178

179
        c = cdata.contours[cid]
180
181
182
183

        # check if contour passes thorough this vox
        for i in 2:size(c.pos)[1]

184
185
186
187
188
            # make sure it is a line segment, not a point
            if c.pos[i-1,1:2] == c.pos[i,1:2]
                continue
            end

189
190
191
192
            # is this line segment completely outside vox?
            if c.pos[i, 1] > voxx1 && c.pos[i-1, 1] > voxx1 ||
                c.pos[i,1] < voxx2 && c.pos[i-1, 1] < voxx2 ||
                c.pos[i,2] < voxy2 && c.pos[i-1, 2] < voxy2 ||
193
                c.pos[i,2] > voxy1 && c.pos[i-1, 2] > voxy1
194

195
                # segment outside vox entirely
196
                if seg_now
197
                    println("Something's gone wrong: segment entirely outside voxel, but last segment inside")
198
199
                end
                continue
200
            end
201

202
203
204
205
            p1inside = c.pos[i-1, 1] < voxx1 && c.pos[i-1, 1] > voxx2 && c.pos[i-1, 2] > voxy2 && c.pos[i-1, 2] < voxy1
            p2inside = c.pos[i, 1] < voxx1 && c.pos[i,1] > voxx2 && c.pos[i,2] > voxy2 && c.pos[i,2] < voxy1
            # is this line segment completely inside vox?
            if p1inside && p2inside
206

207
208
209
210
211
212
                seglen_sofar += dist(c.pos[i], c.pos[i-1]) # append to existing contour
                if !seg_now # start new seg
                    t_start = 0 #  0 bc contour must be starting for this case
                    seg_now = true

                    if i!=2
213
                        println("Whole segment inside but something wrong")
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
                    end
                    continue
                end
            end

            cross_side1 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx1, voxy1], [voxx1, voxy2])
            cross_side2 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx1, voxy1], [voxx2, voxy1])
            cross_side3 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx2, voxy1], [voxx2, voxy2])
            cross_side4 = seg_intersect(c.pos[i-1,:], c.pos[i,:], [voxx2, voxy2], [voxx1, voxy2])

            # does this line segment intersect with vox only once
            if p1inside  p2inside

                # find intersection point
                if cross_side1 || cross_side3
                    # intersection with x
                    xi = [voxx1, voxx2][[cross_side1, cross_side3]][1]
                    p_i = interpolate(c.pos[i-1,:], c.pos[i,:], xi, 1)
                    t_i = interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], xi, c.pos[i,1])
233
                elseif cross_side2 || cross_side4
234
235
236
237
238
239
240
241
242
                    # intersection with y
                    yi = [voxy1, voxy2][[cross_side2, cross_side4]][1]
                    p_i = interpolate(c.pos[i-1,:], c.pos[i,:], yi, 2)
                    t_i = interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], yi, c.pos[i,2])
                end

                if p1inside
                    # end existing segment
                    if !seg_now
243
244
245
                        # contour end on the first segment
                        t_start = 0
                        seglen_sofar = 0
246
247
248
249
250
251
252
                    end
                    seglen_sofar += dist(c.pos[i-1, :], p_i)
                    push!(segcontours, cid)
                    push!(seglen, seglen_sofar)
                    push!(segoffset, (t_i + t_start)/2)
                    seglen_sofar = 0
                    seg_now = false
253
                else
254
255
256
257
                    # start new contour
                    t_start = t_i
                    seglen_sofar = dist(p_i, c.pos[i, :])
                    seg_now = true
258
                end
259
                continue
260

261
262
263
264
265
            elseif sum([cross_side1, cross_side2, cross_side3, cross_side4]) >= 2
                # intersects twice
                p_is = []
                t_is = []
                if cross_side1
266
267
268
269
270
                    p = interpolate(c.pos[i-1,:], c.pos[i,:], voxx1, 1)
                    if !isnan(p[1])
                        push!(p_is, p)
                        push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], voxx1, c.pos[i,1]))
                    end
271
272
                end
                if cross_side2
273
274
275
276
277
                    p = interpolate(c.pos[i-1,:], c.pos[i,:], voxy1, 2)
                    if !isnan(p[1])
                        push!(p_is,p)
                        push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], voxy1, c.pos[i,2]))
                    end
278
279
                end
                if cross_side3
280
281
282
283
284
                    p = interpolate(c.pos[i-1,:], c.pos[i,:], voxx2, 1)
                    if !isnan(p[1])
                        push!(p_is,p)
                        push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,1], voxx2, c.pos[i,1]))
                    end
285
286
                end
                if cross_side4
287
288
289
290
291
                    p = interpolate(c.pos[i-1,:], c.pos[i,:], voxy2, 2)
                    if !isnan(p[1])
                        push!(p_is, p)
                        push!(t_is, interpolate(c.time[i-1], c.time[i], c.pos[i-1,2], voxy2, c.pos[i,2]))
                    end
292
                end
293

294
295
296
297
                if seg_now
                    print("Something's wrong")
                end

298
299
300
301
302
303
304
305
                if length(p_is) >= 2
                    push!(segoffset, mean(t_is))
                    push!(segcontours, cid)
                    if length(p_is) == 2
                        push!(seglen, dist(p_is[1], p_is[2]))
                    else
                        push!(seglen, dist(p_is[1], p_is[3]))
                    end
306
                else
307
308
309
310
311
312
313
314
                    p1inside = c.pos[i-1, 1] <= voxx1 && c.pos[i-1, 1] >= voxx2 &&
                            c.pos[i-1, 2] >= voxy2 && c.pos[i-1, 2] <= voxy1
                    p = p1inside ? c.pos[i-1,:] : c.pos[i,:]
                    t = p1inside ? c.time[i-1] : c.time[i]

                    push!(segcontours, cid)
                    push!(segoffset, (t + t_is[1])/2)
                    push!(seglen, dist(p_is[1], p))
315
                end
316
317
318
            end
        end

319
320
321
322
323
324
325
326
327
328
329
        # if contour ends inside the voxel
        if seg_now
            # end segment
            push!(segcontours, cid)
            push!(seglen, seglen_sofar)
            push!(segoffset, (t_start + last(c.time))/2)
            seglen_sofar = 0
            t_start = 0
            seg_now = false
        end

330
331
        # for those contours find exact segments
    end
332
333
334
335
336
337
338
339
340
341
    c = Float64(segoffset  seglen)/sum(seglen) # constant used for cost calc
    new_seglen = Vector{Float64}()
    new_segcontours = Vector{Int64}()
    for i in 1:length(segcontours)
        if !(segcontours[i] in new_segcontours)
            push!(new_segcontours, segcontours[i])
            push!(new_seglen, sum(seglen[segcontours.==segcontours[i]]))
        end
    end
    return voxmap(new_seglen./sum(new_seglen), new_segcontours, c)
342
343
end

344

345
346
347
function voxdata(fname::String, cdata::contourdata)
    voxels = DataFrames.DataFrame(CSV.File(fname))
    w = dist(Vector(voxels[1, ["x","y","z"]]), Vector(voxels[2, ["x","y","z"]]))
Jayant Khatkar's avatar
Jayant Khatkar committed
348
    println("Assumed width ", w)
349
350
351
352
353
354
355
356
357
    vpos = [[v.x, v.y, v.z] for v in eachrow(voxels)]
    voxms = [voxmap(v, w, cdata) for v in vpos]
    below = indexin([v - [0,0,cdata.layer_height] for v in vpos], vpos)
    replace!(below, nothing=>0)
    return voxdata(voxels, voxms, below, w)
end


function random_rollout(cdata::contourdata)
358
    done_contours = Set{Int}()
359
360
    avail_contours = Set(cdata.layers[1])
    todo_contours = Set(1:length(cdata.contours))
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    rollout = Vector{Int}()

    while length(avail_contours) > 0
        c = rand(avail_contours)
        push!(rollout, c)

        # remove selected contour from todo and avail, add to done
        delete!(avail_contours, c)
        delete!(todo_contours, c)
        push!(done_contours, c)

        # update available contours
        for i in todo_contours
            if i in avail_contours
                continue
376
            elseif length(inneighbors(cdata.G, i)) == 0
377
378
379
380
381
                push!(avail_contours, i)
                continue
            end

            add = true
382
            for j in inneighbors(cdata.G, i)
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
                if !(j in done_contours)
                    add = false
                    break
                end
            end

            if add
                push!(avail_contours, i)
            end
        end
    end

    return rollout
end


399
function valid_swap(rollout::Vector{Int}, i::Int, j::Int, cdata::contourdata)
400
401
402
    # would swapping indices i and j in rollout result in another valid rollout?
    # NOTE THIS FUNCTION DOESNT WORK
    # IT ONLY CHECKS DEPENDENCIES TO A DEPTH OF 1
403
    # TODO, leave for now, use check_validity to double check at the end
404
405
406
407
408
409
410
411
412

    if i>j
        i,j = j,i
    elseif i==j
        return true
    end

    c1 = rollout[i]
    c2 = rollout[j]
413
    c2_dependson = inneighbors(cdata.G, c2)
414
415
416
417
418

    if c1 in c2_dependson
        return false
    end

419
    c1_dependents = outneighbors(cdata.G, c1)
420
421
422
423
424
425
426
427
428
429
430
431
    c_between = rollout[i+1:j-1]

    for c in c_between
        if c in c1_dependents || c in c2_dependson
            return false
        end
    end

    return true
end


432
function check_validity(rollout::Vector{Int}, cdata::contourdata)
433
434
435
436
    # make sure a given rollout is valid
    done_contours = Set{Int}()

    for c in rollout
437
        c_dependson = inneighbors(cdata.G, c)
438
439
440
441
442
443
444
445
446
447
448

        if !issubset(c_dependson, done_contours)
            return false
        end

        push!(done_contours, c)
    end
    return true
end


449
450
451
452
453
454
function swap!(rollout::Vector{Int}, i::Int, j::Int)
    # swap values at ith and jth indices
    rollout[i], rollout[j] = rollout[j], rollout[i]
end


455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
function test_voxmap()
    # create vox
    vox = [0,0,0.5]
    vox_d = 2
    pos1 = [[0.5, -0.5] ones(2)*0.5 ones(2)]
    time1 = [0,1]
    pos2 = [[1.5, 0.5, -0.5, -1.5] ones(4)*0.5 ones(4)]
    time2 = Vector(0:3)
    pos3 = [[-0.5, -0.5] [2, -2] ones(2)]
    time3 = [0, 2.5]
    pos4 = [[0.5, 2.5] [-1.5, -0.5] ones(2)]
    time4 = [0,1]
    pos5 = [[-2,2] [-2,2] ones(2)]
    time5 = [0,1]
    contour1 = contour(pos1, time1)
    contour2 = contour(pos2,time2)
    contour3 = contour(pos3,time3)
    contour4 = contour(pos4,time4)
    contour5 = contour(pos5,time5)

    contours = [contour1, contour2, contour3, contour4, contour5]
476
477
    cdata = contourdata(contours, 1, 1)
    vm = voxmap(vox, vox_d, cdata)
478
    return vm
479
480
481
end


482
483
484
485
# Temperature function
T0 = 215 # extrusion temp
Tc = 25 # room temp
Tcutoff = 100 # temperature above which strain isn't happening
486
k = 0.08 #value unknown
487
488
489
490
491
Temp(t::Number) = Tc + (T0-Tc)*^(-k*t)
# see shape of temp function
#x = 1:100
#y = Temp.(x)
#plot(x,y)
492
493


Jayant Khatkar's avatar
Jayant Khatkar committed
494
495
496
497
498
499
500
501
function clean_contour(c::contour)
    # remove first element of array if second element is the same
    while c.pos[1,:] == c.pos[2,:]
        c.pos = c.pos[2:end, :]
    end
    return c
end

502
503
504
505
506
507
508
509
510
511
512
513

function stress_multiplier!(a::DataFrame, mul::Number)
    a.Sx = a.Sx*mul
    a.Sy = a.Sy*mul
    a.Sz = a.Sz*mul
    a.Txy = a.Txy*mul
    a.Tyz = a.Tyz*mul
    a.Txz = a.Txz*mul
    return
end


514
# ABS material properties
515
α_ABS = 63e-6
516
517
E_ABS = 2.5e9    # 1.19e9 - 2.9e9
σ̄ = 6e7          # mega
518
σ̄_ABS = [
519
520
521
    [σ̄       σ̄*0.5   σ̄*0.375]
    [σ̄*0.5   σ̄       σ̄*0.375]
    [σ̄*0.375 σ̄*0.375 σ̄*0.75 ]
522
523
524
    ]
ABS = material(α_ABS, E_ABS, σ̄_ABS)

525

526
### LOAD IN DATA
Jayant Khatkar's avatar
Jayant Khatkar committed
527
528
obj = "/Users/jayant/phd/tempaware/" * "M1"
contours = clean_contour.(contour.(JSON.parse(open(obj * "contours.json"))))
529
cdata = contourdata(contours, 5, 5) # contour data
530
@time vd = voxdata(obj * "_voxels.csv", cdata)
531
rl = random_rollout(cdata)
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556

#####################################
###### CONSTRUCT COST FUNCTION
mat = ABS
# store contour time lengths inside function
contour_times = [cdata.contours[c].time[end] for c in 1:length(cdata.contours)]

# considered voxels
not_empty_voxels = length.([m.seglen for m in vd.maps]) .>0
valid_voxels = (1:length(vd.below))[(vd.below.!=0) .& not_empty_voxels]
valid_voxels = valid_voxels[not_empty_voxels[vd.below[valid_voxels]]]

#stress_size = .√( vd.voxels.Sx.^2 + vd.voxels.Sy.^2 + vd.voxels.Sz.^2)
#top_voxes = Set(sortperm(stress_size)[1:Int(length(stress_size)/4)]) # top stressors only
#top_valid_voxels = sort(collect(Set(valid_voxels) ∩ top_voxes))
#rel_voxels = vd.voxels[top_valid_voxels,:]

#relmaps = sort(collect(Set(vd.below[top_valid_voxels]) ∪ Set(top_valid_voxels)))
#considered_voxels = indexin(top_valid_voxels, relmaps)
# relbelows = indexin(vd.below[top_valid_voxels], relmaps)
considered_voxels = valid_voxels
relbelows = vd.below[valid_voxels]
rel_voxels = vd.voxels[considered_voxels,:]
relmaps = vd.maps

557
558
559
560
561
562
563
564
565
566
567
568
# voxtimes vectorize
max_contours_per_voxel = maximum([length(r.seglen) for r in relmaps])
vox_contour_id = ones(Int64, length(relmaps), max_contours_per_voxel)
for i in 1:length(relmaps)
    vox_contour_id[i, 1:length(relmaps[i].segcontours)] = relmaps[i].segcontours
end
vox_c = [Float64(v.c) for v in relmaps]
vox_seglen = zeros(length(relmaps), max_contours_per_voxel)
for i in 1:length(relmaps)
    vox_seglen[i, 1:length(relmaps[i].seglen)] = relmaps[i].seglen
end

569
570
571
572
573
574
575
576
# calculate stresses at each voxel
F = (2/mat.σ̄[1,1]^2 - 1/mat.σ̄[3,3]^2)/2
G = 1/(2*mat.σ̄[3,3]^2)
L = 1/(2*mat.σ̄[1,2]^2)
M = 1/(2*mat.σ̄[1,3]^2)

a = quote
    function cost_f(rl::Vector{Int})
577
        # println("Cumulative Sum") # 800k
578
579
        timestart = cumsum([$contour_times[c] for c in rl])

580
581
        #println("Vox times") # 18k
        voxtimes = sum(vox_seglen .* timestart[vox_contour_id], dims=2) .+ vox_c
582

583
        # println("time diff") # 100k
584
585
        Δt = voxtimes[considered_voxels] - voxtimes[relbelows]

586
587
        #println("Temp diff") # 30k
        ΔT = Tcutoff .- ($Tc .+ $(T0-Tc).*.^(-k.*Δt))
588

589
        # println("Temp diff clean") # 500k
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
        replace!(x-> x<0 ? 0 : x, ΔT)

        # println("Stresses") # 6k
        σ11 = rel_voxels.Sx + $(mat.E*mat.α)*ΔT
        σ22 = rel_voxels.Sy + $(mat.E*mat.α)*ΔT
        σ33 = rel_voxels.Sz
        σ12 = rel_voxels.Txy
        σ23 = rel_voxels.Tyz + $((cdata.layer_height/vd.width)*mat.E*mat.α)*ΔT
        σ31 = rel_voxels.Txy + $((cdata.layer_height/vd.width)*mat.E*mat.α)*ΔT
        cost = sum($F * (σ11 - σ22).^2 +
            $G * ((σ33 - σ11).^2 + (σ33 - σ22).^2) +
            $(2 * L) * (σ12).^2 +
            $(2 * M) * (σ23 + σ31).^2)
        return cost
    end
end
eval(a)
cost_f(rl)
608
609

# local search loop
610
611
k = 5 # TODO k is hardcoded in function below (uses global variable)
function best_neighbor!(rl::Vector{Int}, current_cost::Number)
612
    costs = Dict()
613
    for i in 1:1:length(rl)-1
614
        for j in i+1:min(i+k,length(rl)) # TODO length rl is constant
615
616
            if valid_swap(rl, i, j, cdata)
                swap!(rl, i, j)
617
                costs[i,j] = cost_f(rl) # TODO cost function global vairable being used here
618
619
620
621
622
                swap!(rl, i, j)
            end
        end
    end
    v, (i,j) = findmin(costs)
623
624
625
    if abs(v - current_cost) < current_cost/1e6
        return 0
    elseif v < current_cost
626
        swap!(rl, i, j)
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
        return v
    end
    return 0
end


function local_search!(rl::Vector{Int}, max_iter::Int)
    cost_val = Inf
    for l in 1:max_iter
        c = best_neighbor!(rl, cost_val)
        if c  0
            cost_val = c
        else
            break
        end
642
    end
643
644
645
646
647
648
649
650
651
652
653
654
655
    return cost_val
end


rl = random_rollout(cdata)
@time local_search!(rl, 5)

n_starts=10
my_costs = zeros(n_starts)
threads = zeros(n_starts)
@time Threads.@threads for i in 1:n_starts
    my_costs[i] = local_search!(random_rollout(cdata), 5)
    threads[i] = Threads.threadid()
656
end # 12 seconds
657
658
659
660

@time for i in 1:n_starts
    my_costs[i] = local_search!(random_rollout(cdata), 5)
    threads[i] = Threads.threadid()
661
end # 22 seconds
662
663
664


### Visualise distribution of stresses
665
666
667
668
#data = .√( vd.voxels.Sx.^2 + vd.voxels.Sy.^2 + vd.voxels.Sz.^2)
#data = .√(vd.voxels.Txy.^2 + vd.voxels.Tyz.^2 + vd.voxels.Txz.^2)
#histogram(data) # vast majority of voxels near 0 stress - can ignore
#histogram!(sort(data, rev=true)[1:4000]) # vast majority of voxels near 0 stress - can ignore