Przykład pracy nad nowym modułem

Jako ilustrację pracy w Julii projektu pod postacią modułu spróbujemy stworzyć moduł zajmujący się masami jądrowymi. Kompilacja danych eksperymentalnych ukazuje się cyklicznie jako Atomic Mass Evaluation. Ostatnia edycja z 2020 roku jest dostępna tutaj, a dane w postaci tekstowej są np. tu. W tablicach podawany jest bardzo wygodny opis mas nuklidów (nuklid = neutralny atom) przez mass excess (deficyt masy, defekt masy). Wynika to z faktu, że liczba masowa nuklidu z bardzo dobrym przybliżeniem opisuje jego masę, i potrzebna jest tylko drobna poprawka:

M(A,Z)=A×u+Δ. M(A, Z) = A\times u + \Delta .

Jednostka masy uu zdefiniowana jest arbitralnie jako 1/12×M(12C)1/12\times M(^{12}C), a więc deficyt masy nie opisuje w bezpośredni sposób intuicyjnych wielkości fizycznej. uu jest równe 931.494 MeV/c2^2, a deficyty mas do kilkudziesięciu MeV/c2^2.

Polskie nazewnictwo jest w tym przypadku wyjątkowo niespójne i w większości miejsc pojęcia deficytu (defektu) masy są utożsamiane z różnicą między masą jądra, a sumą mas składników, ale ta wielkość jest po prostu energią wiązania. Po co zatem stosować trzy nazwy na energię wiązania, a żadnej na wielkość Δ\Delta? Niestety za wikipedią ten błąd jest powtarzany tak powszechnie, że trudno jest znaleźć w internecie źródło twierdzące inaczej.

Korzystanie z tych tabel (format pozycyjny) jest mało wygodne. Wprawdzie dość łatwo jest dane wczytać i przedstawić w postaci np. tablicy DataFrame, ale za każdym razem, gdy potrzebna jest nam informacja o danym nuklidzie, musielibyśmy wyszukiwać odpowiedni wiersz tabeli (df[(df.A .== A) .& (df.Z .== Z), :]]). Często zamiast liczb atomowych używane są nazwy pierwiastków lub symbole chemiczne, np "U-238", a nie każdy pamięta wszystkie liczby atomowe. Musielibyśmy napisać zatem wygodne funkcje wyszukujące po nazwie lub A, Z. Moglibyśmy też pokusić się, dla wygody, o zdefiniowanie funkcji łączącej dwa jądra. Pomysłów można mieć pewnie i więcej... to może od razu lepiej zrobić moduł?

Założenia

  • wczytujemy dane w oryginalnej postaci do jakiejś tabeli (jakiej?)

  • dla konkretnego nuklidu mamy strukturę zawierającą interesujące nas informacje (A, Z, delta)

  • funkcje szukające w danych odpowiedniego jądra i zwracające je w postaci struktury

  • operator dodawania

Zaczynamy od stworzenia modułu poprzez wzorce PkgTemplates

using PkgTemplates
t = Template(;user="kmiernik", dir="~/workspace/Julia")
t("MassExcess")

Uwaga!

Od tej pory proponuję przejść do katalogu i używać ulubionego edytora. Fragmenty kodu pojawiające się poniżej służą tylko ilustracji, a nie uruchamianiu w Jupyterze!

Zaczniemy od utworzenia odpowiedniej struktury opisującej nuklid. W pliku MassExcess.jl:

module MassExcess

"""
    Structure holding information on nuclide:
    A, Z, element name, mass excess `me`, uncertainty (`dme`)
    and flag if it is an experimental value (`exp`)
"""
struct Nuclide
    A::Int64
    Z::Int64
    element::String
    me::Float64
    dme::Float64
end


end

Skoro mamy już naszą strukturę, możemy zacząć pisać testy, testujące nasze wymagania. W pliku test/runtest.jl przetestujemy różne postaci konstruktorów naszej struktury:

  • Nuclide(A, Z)

  • Nuclide("symbol", A)

  • Nuclide("symbol-A")

gdzie "symbol" oznacza chemiczny symbol pierwiastka (np. Au, Pb, He), zapisany w dowolnej kombinacji dużych i małych liter (np. He, HE, he). Z tablic odczytamy interesujące nas wartości i stworzymy bezpośrednio nuklid, a następnie zażądamy, ale konstruktory zwracały dokładnie ten sam obiekt.

@testset "MassExcess.jl" begin

    he4 = Nuclide(4, 2, "He", 2424.91587, 0.00015)

    @test Nuclide("He-4") == he4
    @test Nuclide("HE-4") == he4
    @test Nuclide("he-4") == he4

end

Testowanie zaczniemy od razu. Uruchamiamy julię w katalogu głównym projektu. W REPL-u włączymy moduł Revise, który śledzi na bieżąco zmiany, uruchomimy środowisko pakietu i test

julia> using Revise
(@v1.8) pkg> activate .
(MassExcess) pkg> test

Dostajemy

UndefVarError: Nuclide not defined
(...)
     @ ~/workspace/Julia/MassExcess/test/runtests.jl:6
(...)
Test Summary: | Error  Total  Time
MassExcess.jl |     1      1  1.1s
ERROR: LoadError: Some tests did not pass: 0 passed, 0 failed, 1 errored, 0 broken.
in expression starting at /home/krm/workspace/Julia/MassExcess/test/runtests.jl:4
ERROR: Package MassExcess errored during testing

Oczywiście zapominieliśmy wyeksportować strukturę Nuclide. Dodajemy export Nuclide i uruchamiamy test jeszcze raz.

(...)
Test Summary: | Error  Total  Time
MassExcess.jl |     4      4  2.0s
(...)

Żaden test nie przeszedł, bo oczywiście nie mamy jeszcze odpowiednich konstruktorów. Warto zwrócić uwagę, że w poprzedniej próbie był tylko jeden test (@testset doszedł do momentu pierwszego błędu nie objętego testem jednostkowym), natomiast fakt wystąpienia błędu w teście jednostkowym @test nie powoduje przerwania zestawu testów.

Uzupełnianiamy moduł o brakujące informacje. Zaczniemy od nazw pierwiastków, zakodujemy je w postaci tablicy 

# in MassExcess
elements = ("H" ,  "He",  "Li",  "Be",  "B" ,
            "C" ,   "N",   "O",   "F",  "Ne",
            "Na",  "Mg",  "Al",  "Si",   "P",
            "S" ,  "Cl",  "Ar",   "K",  "Ca",
            "Sc",  "Ti",   "V",  "Cr",  "Mn",
            "Fe",  "Co",  "Ni",  "Cu",  "Zn",
            "Ga",  "Ge",  "As",  "Se",  "Br",
            "Kr",  "Rb",  "Sr",   "Y",  "Zr",
            "Nb",  "Mo",  "Tc",  "Ru",  "Rh",
            "Pd",  "Ag",  "Cd",  "In",  "Sn",
            "Sb",  "Te",   "I",  "Xe",  "Cs",
            "Ba",  "La",  "Ce",  "Pr",  "Nd",
            "Pm",  "Sm",  "Eu",  "Gd",  "Tb",
            "Dy",  "Ho",  "Er",  "Tm",  "Yb",
            "Lu",  "Hf",  "Ta",   "W",  "Re",
            "Os",  "Ir",  "Pt",  "Au",  "Hg",
            "Tl",  "Pb",  "Bi",  "Po",  "At",
            "Rn",  "Fr",  "Ra",  "Ac",  "Th",
            "Pa",   "U",  "Np",  "Pu",  "Am",
            "Cm",  "Bk",  "Cf",  "Es",  "Fm",
            "Md",  "No",  "Lr",  "Rf",  "Db",
            "Sg",  "Bh",  "Hs",  "Mt",  "Ds",
            "Rg",  "Cn",  "Nh",  "Fl",  "Mc",
            "Lv",  "Ts",  "Og")

Teraz pierwiastek o liczbie atomowej Z będzie zakodowany w Z-tym elemencie tablicy. Następnie brakuje nam tablicy z danymi. Zapiszemy plik w katalogu data/ame2020.dat i dopiszemy funkcję ładującą tablicę. Aby odpocząć nieco od DataFrames tym razem wczytamy go do struktury NamedTuple, gdzie każdy (nazwany) element będzie wektorem odpowiednich danych:

(Z=[0, 1, 1, 1, 2, ...], A = [1, 1, 2, 3, 3, ...], 
 me = [8071.31806, 7288.971064, 13135.722895, 14949.81090, 14931.21888, ...], ...)

czyli dane z i-tego wiersza tabeli będą zawarte w i-tych elementach odpowiednich wektorów.

Pewną komplikacją jest fakt, iż niektóre dane z tabeli nie są wynikami eksperymentalnymi, a ekstrapolacjami danych. Są one zaznaczone znakiem "#". Ponieważ interesują nas dane ściśle eksperymentalne, te punkty pominiemy przy ładowaniu.

W tabeli informacje o Z jądra znajdziemy w kolumnach 10-14, liczbie masowej A - 15-19, defekcie masy - 29-42, a niepewności w 44-54.

function parseame20(amefile)
    skipto = 37
    iline = 0
    Zs = Int64[]
    As = Int64[]
    ms = Float64[]
    ds = Float64[]
    for line in readlines(amefile)
        iline += 1
        if iline < skipto
            continue
        end
        d = line[29:42]
        if occursin("#", d)
            continue
        end
        push!(Zs, parse(Int64, line[10:14]))
        push!(As, parse(Int64, line[15:19]))
        push!(ms, parse(Float64, d))
        push!(ds, parse(Float64, line[44:54]))
    end    
    (Z=Zs, A=As, me=ms, dme=ds) 
end

Dane chcemy wczytać teraz przy ładowaniu struktury

data = parseame20(pkgdir(MassExcess, "data", "ame2020.dat"))

Funkcja pkgdir(module, "path", ...) zwraca ścieżkę w której zainstalowany jest w danym systemie moduł oraz względne ścieżki do katalogów lub plików.

Nie mamy jeszcze odpowiednich konstruktorów, ale możemy dopisać test wczytywania danych. Pierwszym jądrem powinien być neutron (Z=0, A=1, me=8071.31806), a ostatnim bez znaczka ekstrapolacji Ds-270 (Z=110, me=134681.584). Dopiszemy to do testów

@test MassExcess.data.A[1] == 1
    @test MassExcess.data.Z[1] == 0
    @test MassExcess.data.me[1] == 8071.31806
    @test MassExcess.data.A[end] == 270
    @test MassExcess.data.Z[end] == 110
    @test MassExcess.data.me[end] == 134681.584
Test Summary: | Pass  Error  Total  Time
MassExcess.jl |    6      4     10  1.9s
ERROR: LoadError: Some tests did not pass: 6 passed, 0 failed, 4 errored, 0 broken

Nowe testy (6) przechodzi, 4 błędne to konstruktory, którymi najwyższy czas się zająć. Zaczniemy od postaci Nuclide(A, Z). Musimy znaleźć wiersz, w którym mamy jądro A, Z; jeśli go nie znajdziemy to nadamy mu wartości defektu Inf. Wiersz znajdziemy poleceniem findfirst()

function Nuclide(A::Int64, Z::Int64)
    id = findfirst(i->data.A[i] == A && data.Z[i] == Z, eachindex(data.A))
    if isnothing(id)
        Nuclide(A, Z, elements[Z], Inf, Inf)
    else
        Nuclide(A, Z, elements[Z], data.me[id[1]], data.dme[id[1]])
    end
end

Wszystko jest dobrze, zaliczamy dodatkowy test, ale co jeśli ktoś stworzy jądro Nuclide(208, 300), albo Nuclide(300, 128), albo Nuclide(-1, -1)? Niektóre z nich teoretycznie mają sens (choć nie są znane - (300, 128)), inne w ogóle są bez sensu (Z > A, Z < 0). Wypadałoby odpowiednio potraktować te przypadki. Na przykład nieznanym (jeszcze?) pierwiastkom możemy nadać nazwę zgodną z liczbą Z (np. "120"), a wyrzucić odpowiedni błąd dla złych danych wejściowych.

Ubierzemy to w formę konstruktora (A, Z, me, dme), z którego będą mogły korzystać inne konstruktory. Poprawiamy nasz konstruktor (A, Z).

function Nuclide(A, Z, me, dme)
    @assert A >= Z
    if 1 <= Z <= length(elements)
        name = elements[Z]
    elseif Z == 0
        name = "n"
    elseif Z > length(elements)
        name = "$Z"
    else
        throw(ArgumentError("Element Z = $Z is not valid"))
    end
    Nuclide(A, Z, name, me, dme)
end

function Nuclide(A::Int64, Z::Int64)
    id = findfirst(i->data.A[i] == A && data.Z[i] == Z, eachindex(data.A))
    if isnothing(id)
        Nuclide(A, Z, Inf, Inf)
    else
        Nuclide(A, Z, data.me[id[1]], data.dme[id[1]])
    end
end

Zostały nam jeszcze dwa konstruktory, które przerabiają symbol pierwiastka chemicznego na liczbę Z. Zaczniemy od wersji Nuclide("symbol", A). Potrzebna nam będzie funkcja wyszukująca numer pierwiastka, przy czym musimy pamiętać, że chcieliśmy wersje ("HE", "He", "he"), a tabela zawiera pierwszą literę zawsze wielką.

function getZ(name::String)
    frmname = uppercase(name[1]) * lowercase(name[2:end])
    @assert !isnothing(findfirst(x->x==frmname, elements)) "Element name $name is not known"
    findfirst(x->x==frmname, elements)
end

Konstruktor teraz zmieści się w jednej linijce

Nuclide(A::Int64, element) = Nuclide(A, getZ(element))

Ostatni konstruktor będzie akceptował ciągi w postaci "symbol-A", np. "He-4", "HE-4" lub "he-4". Podzielimy zatem podany ciąg znaków w miejscu myślnika i spróbujemy wywołać poprzednią wersję. Drobna uwaga - split zwraca podciąg (SubString) czyli jedynie widok (view) oryginalnego ciągu (dzięki temu nie tworzy kopii). 

function Nuclide(name::String)
    @assert length(split(name, "-")) == 2 "Nuclide name $name format is not supported"
    e, A = split(name, "-")
    A = parse(Int64, A)
    Nuclide(A, String(e))
end

Ostatnim elementem, o który na razie prosiliśmy jest operator dodawania. Stworzymy najpierw test, w którym trzykrotnie dodamy He-4 i powinniśmy dostać C-12, które jest wyjątkowe, bo ma defekt masy z definicji równy zero.

c12 = he4 + he4 + he4

    @test c12.A == 12
    @test c12.Z == 6
    @test c12.me == 0.0

Aby zdefiniować operator + musimy stworzyć nową metodę funkcji Base.+. W naszym przypadku, jeżeli ktoś poda dwa nuklidy, dostanie trzeci, który stworzymy jako sumę A i Z poszczególnych składników. 

Base.:+(x::Nuclide, y::Nuclide) = Nuclide(x.A + y.A, x.Z + y.Z)

Jako ostatni element możemy jeszcze stworzyć funkcję definującą w jaki sposób ma być wypisywana nasz obiekt Nuclide. Np. może to być Symbol-A (C-12)

Dodawanie modułu do Julii

Nasz moduł jest gotowy. Możemy zapisać teraz zmiany w repozytorium git (git add ., git commit -a). Aby dodać lokalny moduł do Julii:

(@v1.8)] add "/home/user/workspace/MassExcess"

Lub wysłać repozytorium na server (GitHub, GitLab, ...) i dodać je do Julii

(@v1.8)] add https://github.com/kmiernik/MassExcess

Lub zarejestrować moduł w Julia General Registry (tu jest szczegółowa instrukcja) i po zaakceptowaniu dodać go jak każdy inny

(@v1.8)] add MassExcess

Zależności

Jeżeli nasz moduł zależy od innego (powiedzmy DataFrames) to musimy normalnie powiedzieć, że chcemy zaimportować dany moduł

module MassExcess

using DataFrames
(...)

Ale także dodać DataFrames będąc w środowisku naszego pakietu

pkg> activate .
(MassExcess) pkg> add DataFrames

Lub

julia> Pkg.activate(".")
julia> Pkg.add("DataFrames")

Do pliku Manifest.toml dodane zostanie DataFrames i wszystkie potrzebne zależności w aktualnej wersji.

CC BY-SA 4.0 Krzysztof Miernik. Last modified: December 04, 2023. Website built with Franklin.jl and the Julia programming language.