Option Explicit

' =====================================================================
' Leitor de instancias TSPLIB (TSP / ATSP) em VBA
' Devolve a matriz de distancias D(0 To n-1, 0 To n-1) do tipo Long
' (a cidade 1 da instancia fica no indice 0, a cidade 2 no indice 1, etc.)
'
' Tipos de peso suportados:
'   EUC_2D, CEIL_2D, ATT, GEO, MAN_2D, MAX_2D  (a partir das coordenadas)
'   EXPLICIT com os formatos: FULL_MATRIX, UPPER_ROW, LOWER_ROW,
'   UPPER_DIAG_ROW, LOWER_DIAG_ROW, UPPER_COL, LOWER_COL,
'   UPPER_DIAG_COL, LOWER_DIAG_COL
' =====================================================================

Private Const PI_TSPLIB As Double = 3.141592   ' valor definido na especificacao TSPLIB95
Private Const RRR As Double = 6378.388         ' raio da Terra usado no tipo GEO

' ---------------------------------------------------------------------
' Funcao principal: le o arquivo e devolve a matriz de distancias
' ---------------------------------------------------------------------
Public Function LerTSPLIB(ByVal caminho As String) As Long()
    Dim f As Integer
    Dim linha As String, chave As String, valor As String, secao As String
    Dim p As Long, i As Long, j As Long
    Dim n As Long
    Dim tipoPeso As String, formato As String
    Dim x() As Double, y() As Double, nCoord As Long
    Dim pesos() As Double, nPesos As Long
    Dim tok() As String
    Dim D() As Long
    Dim conteudo As String, linhas() As String, idx As Long

    ' Le o arquivo inteiro e separa as linhas manualmente. O Line Input do VBA
    ' nao reconhece arquivos com quebra de linha so LF (comum nos .tsp baixados
    ' da internet), o que fazia o arquivo inteiro virar uma unica linha.
    f = FreeFile
    Open caminho For Binary Access Read As #f
    conteudo = Space$(LOF(f))
    Get #f, , conteudo
    Close #f

    conteudo = Replace(conteudo, vbCrLf, vbLf)
    conteudo = Replace(conteudo, vbCr, vbLf)
    linhas = Split(conteudo, vbLf)

    idx = 0
    Do While idx <= UBound(linhas)
        linha = linhas(idx)
        idx = idx + 1
        linha = Trim$(Replace(linha, vbTab, " "))

        If Len(linha) > 0 Then
            If linha Like "[A-Za-z]*" Then
                ' ---- linha de cabecalho ou inicio de secao ----
                p = InStr(linha, ":")
                If p > 0 Then
                    secao = ""
                    chave = UCase$(Trim$(Left$(linha, p - 1)))
                    valor = UCase$(Trim$(Mid$(linha, p + 1)))
                    Select Case chave
                        Case "DIMENSION": n = CLng(Val(valor))
                        Case "EDGE_WEIGHT_TYPE": tipoPeso = valor
                        Case "EDGE_WEIGHT_FORMAT": formato = valor
                    End Select
                Else
                    secao = UCase$(linha)
                    If secao = "EOF" Then Exit Do
                    If n <= 0 Then
                        Err.Raise vbObjectError + 1, "LerTSPLIB", "DIMENSION nao encontrada antes de " & secao
                    End If
                    Select Case secao
                        Case "NODE_COORD_SECTION"
                            ReDim x(1 To n): ReDim y(1 To n): nCoord = 0
                        Case "EDGE_WEIGHT_SECTION"
                            ReDim pesos(1 To n * n): nPesos = 0
                    End Select
                End If
            Else
                ' ---- linha de dados (comeca com numero) ----
                Select Case secao
                    Case "NODE_COORD_SECTION"
                        tok = Split(Colapsa(linha), " ")
                        If UBound(tok) < 2 Then
                            Err.Raise vbObjectError + 2, "LerTSPLIB", "Linha de coordenadas invalida: " & linha
                        End If
                        nCoord = nCoord + 1
                        If nCoord > n Then
                            Err.Raise vbObjectError + 3, "LerTSPLIB", "Mais coordenadas do que DIMENSION"
                        End If
                        x(nCoord) = Val(tok(1))
                        y(nCoord) = Val(tok(2))

                    Case "EDGE_WEIGHT_SECTION"
                        tok = Split(Colapsa(linha), " ")
                        For i = 0 To UBound(tok)
                            nPesos = nPesos + 1
                            If nPesos > n * n Then
                                Err.Raise vbObjectError + 4, "LerTSPLIB", "Pesos demais na EDGE_WEIGHT_SECTION"
                            End If
                            pesos(nPesos) = Val(tok(i))
                        Next i
                End Select
            End If
        End If
    Loop

    If n <= 0 Then Err.Raise vbObjectError + 1, "LerTSPLIB", "DIMENSION nao encontrada"

    ReDim D(0 To n - 1, 0 To n - 1)

    Select Case tipoPeso
        Case "EXPLICIT"
            PreencheExplicita D, pesos, nPesos, n, formato

        Case "EUC_2D", "CEIL_2D", "ATT", "GEO", "MAN_2D", "MAX_2D"
            If nCoord <> n Then
                Err.Raise vbObjectError + 5, "LerTSPLIB", "Numero de coordenadas (" & nCoord & ") diferente de DIMENSION (" & n & ")"
            End If
            For i = 1 To n - 1
                For j = i + 1 To n
                    D(i - 1, j - 1) = CLng(DistCoord(tipoPeso, x(i), y(i), x(j), y(j)))
                    D(j - 1, i - 1) = D(i - 1, j - 1)
                Next j
            Next i

        Case Else
            Err.Raise vbObjectError + 6, "LerTSPLIB", "EDGE_WEIGHT_TYPE nao suportado: " & tipoPeso
    End Select

    LerTSPLIB = D
End Function

' ---------------------------------------------------------------------
' Distancia entre dois pontos conforme o tipo de peso
' ---------------------------------------------------------------------
Private Function DistCoord(ByVal tipo As String, ByVal xi As Double, ByVal yi As Double, _
                           ByVal xj As Double, ByVal yj As Double) As Double
    Dim dx As Double, dy As Double, r As Double, t As Double
    Dim latI As Double, lonI As Double, latJ As Double, lonJ As Double
    Dim q1 As Double, q2 As Double, q3 As Double

    dx = xi - xj
    dy = yi - yj

    Select Case tipo
        Case "EUC_2D"
            DistCoord = Nint(Sqr(dx * dx + dy * dy))

        Case "CEIL_2D"
            DistCoord = Teto(Sqr(dx * dx + dy * dy))

        Case "MAN_2D"
            DistCoord = Nint(Abs(dx) + Abs(dy))

        Case "MAX_2D"
            DistCoord = Application_Max(Nint(Abs(dx)), Nint(Abs(dy)))

        Case "ATT"
            r = Sqr((dx * dx + dy * dy) / 10#)
            t = Nint(r)
            If t < r Then DistCoord = t + 1 Else DistCoord = t

        Case "GEO"
            latI = GeoRad(xi): lonI = GeoRad(yi)
            latJ = GeoRad(xj): lonJ = GeoRad(yj)
            q1 = Cos(lonI - lonJ)
            q2 = Cos(latI - latJ)
            q3 = Cos(latI + latJ)
            DistCoord = Fix(RRR * Acos(0.5 * ((1# + q1) * q2 - (1# - q1) * q3)) + 1#)
    End Select
End Function

' ---------------------------------------------------------------------
' Preenche a matriz a partir da EDGE_WEIGHT_SECTION
' ---------------------------------------------------------------------
Private Sub PreencheExplicita(ByRef D() As Long, ByRef pesos() As Double, _
                              ByVal nPesos As Long, ByVal n As Long, ByVal formato As String)
    Dim i As Long, j As Long, k As Long

    Select Case formato
        Case "FULL_MATRIX"
            If nPesos < n * n Then GoTo Faltam
            For i = 1 To n
                For j = 1 To n
                    k = k + 1
                    D(i - 1, j - 1) = CLng(pesos(k))
                Next j
            Next i

        Case "UPPER_ROW", "LOWER_COL"          ' triangulo superior, sem diagonal
            If nPesos < n * (n - 1) \ 2 Then GoTo Faltam
            For i = 1 To n - 1
                For j = i + 1 To n
                    k = k + 1
                    D(i - 1, j - 1) = CLng(pesos(k)): D(j - 1, i - 1) = D(i - 1, j - 1)
                Next j
            Next i

        Case "LOWER_ROW", "UPPER_COL"          ' triangulo inferior, sem diagonal
            If nPesos < n * (n - 1) \ 2 Then GoTo Faltam
            For i = 2 To n
                For j = 1 To i - 1
                    k = k + 1
                    D(i - 1, j - 1) = CLng(pesos(k)): D(j - 1, i - 1) = D(i - 1, j - 1)
                Next j
            Next i

        Case "UPPER_DIAG_ROW", "LOWER_DIAG_COL" ' triangulo superior, com diagonal
            If nPesos < n * (n + 1) \ 2 Then GoTo Faltam
            For i = 1 To n
                For j = i To n
                    k = k + 1
                    D(i - 1, j - 1) = CLng(pesos(k)): D(j - 1, i - 1) = D(i - 1, j - 1)
                Next j
            Next i

        Case "LOWER_DIAG_ROW", "UPPER_DIAG_COL" ' triangulo inferior, com diagonal
            If nPesos < n * (n + 1) \ 2 Then GoTo Faltam
            For i = 1 To n
                For j = 1 To i
                    k = k + 1
                    D(i - 1, j - 1) = CLng(pesos(k)): D(j - 1, i - 1) = D(i - 1, j - 1)
                Next j
            Next i

        Case Else
            Err.Raise vbObjectError + 7, "LerTSPLIB", "EDGE_WEIGHT_FORMAT nao suportado: " & formato
    End Select
    Exit Sub

Faltam:
    Err.Raise vbObjectError + 8, "LerTSPLIB", "EDGE_WEIGHT_SECTION com menos valores do que o esperado"
End Sub

' ---------------------------------------------------------------------
' Auxiliares
' ---------------------------------------------------------------------
Private Function Colapsa(ByVal s As String) As String
    Do While InStr(s, "  ") > 0
        s = Replace(s, "  ", " ")
    Loop
    Colapsa = s
End Function

' Arredondamento "nint" da TSPLIB (o Round do VBA usa arredondamento bancario)
Private Function Nint(ByVal v As Double) As Double
    Nint = Int(v + 0.5)
End Function

Private Function Teto(ByVal v As Double) As Double
    Dim c As Double
    c = Int(v)
    If c < v Then c = c + 1
    Teto = c
End Function

Private Function Application_Max(ByVal a As Double, ByVal b As Double) As Double
    If a > b Then Application_Max = a Else Application_Max = b
End Function

' Converte coordenada DDD.MM para radianos (formula da especificacao TSPLIB)
Private Function GeoRad(ByVal v As Double) As Double
    Dim graus As Double, minutos As Double
    graus = Fix(v)
    minutos = v - graus
    GeoRad = PI_TSPLIB * (graus + 5# * minutos / 3#) / 180#
End Function

' VBA nao tem Acos nativo
Private Function Acos(ByVal v As Double) As Double
    If v >= 1# Then
        Acos = 0#
    ElseIf v <= -1# Then
        Acos = 4# * Atn(1)
    Else
        Acos = Atn(-v / Sqr(-v * v + 1#)) + 2# * Atn(1)
    End If
End Function

' ---------------------------------------------------------------------
' Exemplo de uso (Excel): escolhe o arquivo e despeja a matriz na planilha
' ---------------------------------------------------------------------
Public Sub ExemploLeitura()
    Dim arq As Variant
    Dim D() As Long
    Dim n As Long

    arq = Application.GetOpenFilename("Instancias TSPLIB (*.tsp;*.atsp),*.tsp;*.atsp")
    If VarType(arq) = vbBoolean Then Exit Sub   ' usuario cancelou

    D = LerTSPLIB(CStr(arq))
    n = UBound(D, 1) + 1   ' matriz indexada de 0 a n-1

    With ActiveSheet
        .Cells.ClearContents
        .Range("A1").Resize(n, n).Value = D
    End With

    MsgBox "Instancia lida: " & n & " cidades. D(0,1) = " & D(0, 1)
End Sub