main.jl 16.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

10

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


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


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


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


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


48
49
50
51
52
53
54
55
56
57
58
59
60
61
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


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

66

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

70
    # separate contours into layers
71
    layer_heights = sort(collect(Set([c.pos[end,3] for c in cons])))
72
73
74
75
76
77
78
79
80
81
82
83
84
    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

85
    # loop through contours from previous layer and compare waypoints
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    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
107

108
    return contourdata(cons, G, layers, Dict(), layer_heights[1])
109
110
111
112
113
114
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
115
        return 1 # clockwise
116
    elseif val < 0
117
        return 2 # anticlockwise
118
    else
119
        return 0 # colinear
120
121
122
123
    end
end


124
125
126
127
128
129
130
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


131
132
133
134
135
136
137
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)

138
139
140
141
142
    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)
143
144
145
146
147
148
        return true
    end

    return false
end

149

150
151
152
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)
153

154

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

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

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

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

176
    for cid in cdata.layers[l]
177

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

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

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

188
189
190
191
            # 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 ||
192
                c.pos[i,2] > voxy1 && c.pos[i-1, 2] > voxy1
193

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

201
202
203
204
            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
205

206
207
208
209
210
211
                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
212
                        println("Whole segment inside but something wrong")
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
                    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])
232
                elseif cross_side2 || cross_side4
233
234
235
236
237
238
239
240
241
                    # 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
242
243
244
                        # contour end on the first segment
                        t_start = 0
                        seglen_sofar = 0
245
246
247
248
249
250
251
                    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
252
                else
253
254
255
256
                    # start new contour
                    t_start = t_i
                    seglen_sofar = dist(p_i, c.pos[i, :])
                    seg_now = true
257
                end
258
                continue
259

260
261
262
263
264
            elseif sum([cross_side1, cross_side2, cross_side3, cross_side4]) >= 2
                # intersects twice
                p_is = []
                t_is = []
                if cross_side1
265
266
267
268
269
                    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
270
271
                end
                if cross_side2
272
273
274
275
276
                    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
277
278
                end
                if cross_side3
279
280
281
282
283
                    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
284
285
                end
                if cross_side4
286
287
288
289
290
                    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
291
                end
292

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

297
298
299
300
301
302
303
304
                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
305
                else
306
307
308
309
310
311
312
313
                    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))
314
                end
315
316
317
            end
        end

318
319
320
321
322
323
324
325
326
327
328
        # 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

329
330
        # for those contours find exact segments
    end
331
    return voxmap(seglen, segoffset, segcontours)
332
333
end

334

335
336
337
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
338
    println("Assumed width ", w)
339
340
341
342
343
344
345
346
347
    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)
348
    done_contours = Set{Int}()
349
350
    avail_contours = Set(cdata.layers[1])
    todo_contours = Set(1:length(cdata.contours))
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
    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
366
            elseif length(inneighbors(cdata.G, i)) == 0
367
368
369
370
371
                push!(avail_contours, i)
                continue
            end

            add = true
372
            for j in inneighbors(cdata.G, i)
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
                if !(j in done_contours)
                    add = false
                    break
                end
            end

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

    return rollout
end


389
function valid_swap(rollout::Vector{Int}, i::Int, j::Int, cdata::contourdata)
390
391
392
    # 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
393
    # TODO, leave for now, use check_validity to double check at the end
394
395
396
397
398
399
400
401
402

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

    c1 = rollout[i]
    c2 = rollout[j]
403
    c2_dependson = inneighbors(cdata.G, c2)
404
405
406
407
408

    if c1 in c2_dependson
        return false
    end

409
    c1_dependents = outneighbors(cdata.G, c1)
410
411
412
413
414
415
416
417
418
419
420
421
    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


422
function check_validity(rollout::Vector{Int}, cdata::contourdata)
423
424
425
426
    # make sure a given rollout is valid
    done_contours = Set{Int}()

    for c in rollout
427
        c_dependson = inneighbors(cdata.G, c)
428
429
430
431
432
433
434
435
436
437
438

        if !issubset(c_dependson, done_contours)
            return false
        end

        push!(done_contours, c)
    end
    return true
end


439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
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]
460
461
    cdata = contourdata(contours, 1, 1)
    vm = voxmap(vox, vox_d, cdata)
462
    return vm
463
464
465
end


466
function rollout2time(rollout::Vector{Int}, cdata::contourdata)
467
    # start time of each contour, assuming no travel time
468
    return cumsum([cdata.contours[c].time[end] for c in rollout])
469
470
end

471
472
473
474
# Temperature function
T0 = 215 # extrusion temp
Tc = 25 # room temp
Tcutoff = 100 # temperature above which strain isn't happening
475
k = 0.08 #value unknown
476
477
478
479
480
Temp(t::Number) = Tc + (T0-Tc)*^(-k*t)
# see shape of temp function
#x = 1:100
#y = Temp.(x)
#plot(x,y)
481

482
function calc_cost(rollout::Vector{Int}, cdata::contourdata, vd::voxdata, mat::material)
483
484

    # go from rollout to timestart for each contour
485
    timestart = rollout2time(rollout, cdata)
486
487

    # calculate time at each voxel
488
    voxtimes = [sum(v.seglen.*(v.segoffset + timestart[v.segcontours]))/sum(v.seglen) for v in vd.maps]
489
490

    # calculate temp difference from voxel below it
491
492
493
    considered_voxels = (1:length(vd.below))[(vd.below.!=0) .& (.!isnan.(voxtimes))] # cannot calculate cost is no voxel underneath
    considered_voxels = considered_voxels[.!isnan.(voxtimes[vd.below[considered_voxels]])]
    Δt = voxtimes[considered_voxels] - voxtimes[vd.below[considered_voxels]]
494
    ΔT = Tcutoff .- Temp.(Δt)
495
496

    # calculate stresses at each voxel
497
498
499
500
501
    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)
    rel_voxels = vd.voxels[considered_voxels,:]
502
503

    # calculate cost func at each voxel
504
505
506
507
508
509
510
511
512
513
514
    σ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 = F * (σ11 - σ22).^2 +
        G * ((σ33 - σ11).^2 + (σ33 - σ22).^2) +
        2 * L * (σ12).^2 +
        2 * M * (σ23 + σ31).^2 .- 1.0
    return sum(cost)
515
516
517
end


518
function construct_cost(cdata::contourdata)
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    a = 1
    b = :($a + 5) # using $ 'interpolates' literal expression into the 'quoted' expression
    c = :(5*5 + a)

    d = :($b + $c) #interpolate other expressions to form larger expressions

    ex1 = :(i*2) # construct this according to voxels and contours
    ex2 = quote
        function cost_f(i::Int)
            return $ex1
        end
    end
    return eval(ex2) # return the contructed function
end

Jayant Khatkar's avatar
Jayant Khatkar committed
534
535
536
537
538
539
540
541
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

542
543
544
545
546
547
548
549
550
551
552
# ABS material properties
α_ABS = 100 # 78 - 108
E_ABS = 2.5e9 # 1.19e9 - 2.9e9
σ̄ = 3e7
σ̄_ABS = [
    [σ̄     σ̄*0.5 σ̄*0.5]
    [σ̄*0.5 σ̄     σ̄*0.5]
    [σ̄*0.5 σ̄*0.5 σ̄*0.75]
    ]
ABS = material(α_ABS, E_ABS, σ̄_ABS)

553

Jayant Khatkar's avatar
Jayant Khatkar committed
554
555
obj = "/Users/jayant/phd/tempaware/" * "M1"
contours = clean_contour.(contour.(JSON.parse(open(obj * "contours.json"))))
556
cdata = contourdata(contours, 5, 5) # contour data
Jayant Khatkar's avatar
Jayant Khatkar committed
557
vd = voxdata(obj * "_voxels.csv", cdata) # TODO Need to debug
558
rl = random_rollout(cdata)
Jayant Khatkar's avatar
Jayant Khatkar committed
559
560
@benchmark rl = random_rollout(cdata) # 20/second
@benchmark valid_swap(rl, rand(1:length(contours)), rand(1:length(contours)), cdata) # 1mil/second
561
calc_cost(rl, cdata, vd, ABS) # 500/second