include("utils.jl")

### ABS material properties
α_ABS = 63e-6
E_ABS = 2.5e9    # 1.19e9 - 2.9e9
σ̄ = 6e7          # mega
Tcutoff = 100 # temperature above which strain isn't happening
σ̄_ABS = [
    [σ̄       σ̄*0.5   σ̄*0.375]
    [σ̄*0.5   σ̄       σ̄*0.375]
    [σ̄*0.375 σ̄*0.375 σ̄*0.75 ]
    ]
ABS = material(α_ABS, E_ABS, σ̄_ABS, Tcutoff)
td = tempdecay(215, 25, 0.04) # extrusion temp, room temp, decay rate
tdslow = tempdecay(215, 25, 0.01)
tdfast = tempdecay(215, 25, 0.08)
# visualise_tempdecay(tdfast)


### LOAD IN DATA
obj = "/Users/jayant/phd/tempaware/models/" * "M5"
contours = clean_contour.(contour.(JSON.parse(open(obj * "contours.json"))))
cdata = contourdata(contours, 40, 5) # contour data
@time vd = voxdata(obj * "_voxels.csv", cdata)
#stress_multiplier!(vd.voxels, 10)
rl = random_rollout(cdata)

### CONSTRUCT COST FUNCTION GIVEN MODEL ADN MATERIAL
construct_cost(cdata, vd, ABS, td, :cost_f)
c = cost_f(rl)

# local search loop
k = 150 # TODO k is hardcoded in function below (uses global variable)
function best_neighbor!(rl::Vector{Int}, current_cost::Number)
    costs = Dict()
    for i in 1:1:length(rl)-1
        for j in i+1:min(i+k,length(rl)) # TODO length rl is constant
            if valid_swap(rl, i, j, cdata)
                swap!(rl, i, j)
                costs[i,j] = cost_f(rl) # TODO cost function global vairable being used here
                swap!(rl, i, j)
            end
        end
    end
    v, (i,j) = findmin(costs)
    if abs(v - current_cost) < current_cost/1e6
        return 0
    elseif v < current_cost
        swap!(rl, i, j)
        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
    end
    return cost_val
end


# single local search
rl = random_rollout(cdata)
cost_f(rl)
max_iterations = 50
@time local_search!(rl, max_iterations)

# run many times
n_starts=10
opt_costs = zeros(n_starts)
start_costs = zeros(n_starts)
threads = zeros(n_starts)
@time Threads.@threads for i in 1:n_starts
    rl = random_rollout(cdata)
    start_costs[i] = cost_f(rl)
    opt_costs[i] = local_search!(rl, max_iterations)
    threads[i] = Threads.threadid()
end

@time for i in 1:n_starts
    my_costs[i] = local_search!(random_rollout(cdata), 5)
    threads[i] = Threads.threadid()
end

### Visualise distribution of stresses
data_stress = .√( vd.voxels.Sx.^2 + vd.voxels.Sy.^2 + vd.voxels.Sz.^2)
data_strain = .√(vd.voxels.Txy.^2 + vd.voxels.Tyz.^2 + vd.voxels.Txz.^2)
data = data_stress
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

### Visualise voxmap
plot(vd, 46, cdata)
plot(Vector(1:length(rl)), cdata)
plot(rl, cdata)

rl = random_rollout(cdata)
rl = Vector(1:length(rl))
@time plot_animate(rl, cdata)