# Code Question 2.1 (and 2.2 if you set k = 2 and uncomment code)
function ISOMAP(X,k=3)
    n = size(X,1)

    # Compute all distances

    D = geodisic_distance(X,k)

    # Initialize low-dimensional representation with PCA
    model = PCA(X,2)
    Z = model.compress(X)

    funObj(z) = stress(z,D)

    Z[:] = findMin(funObj,Z[:])

    return Z
end

function geodisic_distance(X, k)
    n = size(X,1)
    D = distancesSquared(X,X)
    D = sqrt.(abs.(D))
    Knn = fill(Inf, (n,n))

    for i in 1:n
        x = sortperm(D[i,:])
        kNeighbours = x[2:k+1]

        for idx in kNeighbours
            Knn[i,idx] = D[i,idx]
            Knn[idx, i] = D[idx, i]
            @assert(D[i,idx] == D[idx,i])
        end
    end

    toReturn = zeros(n,n)

    for i in 1:n
        for j in 1:n
            toReturn[i,j] = dijkstra(Knn,i,j)
            end
        end
    end

    # Set Inf distances to max distance in toReturn for Q 2.2
    # Uncomment out this section for 2.2
    #max = 0
    #for i in 1:n
    #    for j in 1:n
    #        val = toReturn[i,j]
    #        if (val > max && val != Inf)
    #            max = val
    #        end
    #    end
    #end

    #for i in 1:n
    #    for j in 1:n
    #        if (toReturn[i,j] == Inf)
    #            toReturn[i,j] = max
    #        end
    #    end
    #end

    return toReturn
end