plot histograms of cost function for better understanding

main.jl

@@ -112,4 +112,26 @@ plot(rl_l, cdata)
 @time plot_animate(rl_g, cdata, rate=5)
 
 # plot results
-plot(results)
\ No newline at end of file
+plot(results)
+
+
+# plot before and after histograms
+import Plots: histogram
+i = 8
+obj = "M" * string(i)
+obj_add =  add * obj
+contours = clean_contour.(contour.(JSON.parse(open(obj_add * "contours.json"))))
+cdata = contourdata(contours, 40, 20) # contour data
+vd = voxdata(obj_add * "_voxels.csv", cdata)
+construct_cost_hist(cdata, vd, ABS, td, :cost_hist)
+
+histogram(
+    cost_hist(Vector(1:length(cdata.contours))),
+    bins=50,
+    label="Default " * string(sum(cost_hist(Vector(1:length(cdata.contours))))),
+    opacity=0.5)
+histogram!(
+    cost_hist(Vector{Int}(results[obj]["best_rollout"])),
+    bins=50,
+    label="Opt " * string(sum(cost_hist(Vector{Int}(results[obj]["best_rollout"])))),
+    opacity=0.5)
\ No newline at end of file

utils.jl

@@ -719,6 +719,65 @@ function construct_cost(cdata::contourdata, vd::voxdata, mat::material, td::temp
 end
 
 
+function construct_cost_hist(cdata::contourdata, vd::voxdata, mat::material, td::tempdecay, fname::Symbol=:cost_hist)
+    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]]]
+
+    considered_voxels = valid_voxels
+    relbelows = vd.below[valid_voxels]
+    rel_voxels = vd.voxels[considered_voxels,:]
+    relmaps = vd.maps
+
+    # 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
+
+    # precompute constant values
+    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 $fname(rl::Vector{Int})
+            # voxel times
+            timestart = cumsum([$contour_times[c] for c in rl])
+            voxtimes = sum($vox_seglen .* timestart[$vox_contour_id], dims=2) .+ $vox_c
+
+            # voxel temps
+            Δt = voxtimes[$considered_voxels] - voxtimes[$relbelows]
+            ΔT = $(mat.T_cutoff.-td.ambient) .- $(td.extrusion-td.ambient).*ℯ.^(-$td.decay_rate.*Δt)
+            replace!(x-> x<0 ? 0 : x, ΔT)
+
+            # voxel stresses
+            rel_v= $rel_voxels
+            σ11 = rel_v.Sx + $(mat.E*mat.α)*ΔT
+            σ22 = rel_v.Sy + $(mat.E*mat.α)*ΔT
+            σ33 = rel_v.Sz
+            σ12 = rel_v.Txy
+            σ23 = rel_v.Tyz + $((cdata.layer_height/vd.width)*mat.E*mat.α)*ΔT
+            σ31 = rel_v.Txy + $((cdata.layer_height/vd.width)*mat.E*mat.α)*ΔT
+            return $F * (σ11 - σ22).^2 +
+                $G * ((σ33 - σ11).^2 + (σ33 - σ22).^2) +
+                $(2 * L) * (σ12).^2 +
+                $(2 * M) * (σ23 + σ31).^2
+        end
+    end
+    return eval(a)
+end
+
 function construct_best_neighbor(cdata::contourdata, cost_f::Function, k::Int=50)
     n_contours=length(cdata.contours)
     bn_func = quote
@@ -746,7 +805,7 @@ function construct_best_neighbor(cdata::contourdata, cost_f::Function, k::Int=50
     return eval(bn_func)
 end
 
-
+# TODO - COMPLETE THIS
 function construct_best_neighbor_multi(cdata::contourdata,
         cost_f::Function,
         k::Int=10,