Numerical differentiation using the finite difference formula implemented with AVX

During another one of my assembly-related escapades, I came across a following problem : using the finite difference formula, create a function that would calculate the values of a given function’s derivative, given the beginning and end point to be inspected, as well as the precision, e.g. how many samples (individual values) to take in that area. The function being inspected was a simple quadratic function, given by the general formula f(x) = ax^2 + bx + c \; (a,b,c \in \mathbb{R}). Of course, it has a very trivial derivative of f'(x) = 2ax + b, but the computer doesn’t necessarily know that.

The finite difference formula can be summarized as follows : if we need to calculate n \; \in \mathbb{Z} values of a function’s derivative on the area of <a, b> \; (a,b \in \mathbb{R}), then the value of the derivative can be approximated by the following formula :

f'(x) = \frac{f(x + \Delta x) - f(x - \Delta x)}{2 \Delta x}

where \Delta x = \frac{b-a}{n}, e.g. the size of a single „step” according to the selected precision – the bigger the precision, the „smaller” the step is and thus the function is calculated with a greater granularity.

Writing some sequential code for this formula is obviously trivial. Parallelization of the computations can be a more difficult task, though. In my approach, I chose to use „masks” (couldn’t come up with a better name for these, unfortunately : they kind of resemble bitmasks, that’s why the name came to me, I guess) which represent the increments (or decrements) of their respective arguments. There are four such „masks” :

  • the „dx mask”, which contains the value of \Delta x multiplied by 8, since there are 8 values calculated in one iteration of the loop
  • the „2dx mask”, simply containing \Delta x multiplied by 2 – this is the value by which the calculated values are divided by
  • the „plus mask”, which contains the values of \Delta x multiplied by 1 to 8. This mask is used to calculate the value of $latex x + \Delta x$ in the loop, which is later used as an argument for the function.
  • the „minus mask”, which contains the values of \Delta x multiplied by -1 to 6. Analogously to the „plus mask”, it is used to calculate the value of $latex x – \Delta x$.

Each of those masks takes up a whole YMM register and is added to the current quadratic function arguments, producing new arguments for which the value of the derivative is calculated.

I do realize that this all may be quite confusing, so let me illustrate this on a diagram. It attempts to show (though poorly) what is done in one iteration of the loop.

YMM5 represents the „current position” within the function, e.g. the current argument. After adding the „plus” and „minus masks” to these values, every respective element of the resulting YMM register contains the argument to the function which corresponds to the argument required by the formula of the theorem. Therefore, after calculating the value of that function, subtracting the respective values from each other and dividing these results by 2 \Delta x, we obtain the „dval”, which is the final value of the derivative of the function at the given point. That value is written to memory. The „current position” is incremented by 8 \Delta x and the loop is executed again (unless the required number of iterations has been reached).

If there are any leftover values, they are calculated by normal sequential code, operating on one element in one iteration. The sequential loop will make at most 7 iterations, so it’s not a bottleneck by any means.

I also came across a very interesting issue related to memory performance and its correlation with the alignment of the address. All the SSE instructions required the destination (or source, in the case of using instructions operating on packed data) memory operand’s address to be 16-byte aligned. If it wasn’t aligned, a general protection fault was reported while executing that instruction. The one of the few SSE instructions that didn’t require the address to be 16-byte aligned was movups, but using non-16-byte aligned addresses resulted in a dramatic performance decrease. I found this not to be true anymore in the case of AVX. The Intel Basic Architecture manual says that

Most VEX-encoded SIMD numeric and data processing instruction semantics withmemory operand have relaxed memory alignment requirements than instructions encoded using SIMD prefixes.
(section 13.1.3)


With the exception of explicitly aligned 16 or 32 byte SIMD load/store instructions, most VEX-encoded, arithmetic and data processing instructions operate in a flexible environment regarding memory address alignment, i.e. VEX-encoded instruction with 32-byte or 16-byte load semantics will support unaligned load operation by default. Memory arguments for most instructions with VEX prefix operate normally without causing #GP(0) on any byte-granularity alignment (unlike Legacy SSE instructions).
(section 13.3)

I was curious to find out what kind of penalty an unaligned memory access actually introduces. In order to do that, I compared the performance of the function when the destination array was allocated via standard malloc (which is guaranteed to return an 8-byte aligned address when using glibc, I don’t know about any other libc implementations) and mmap (which returns addresses located on a page boundary, which is 4KB in case of x86). The results showed that there is a 5% increase in performance when using the address returned by mmap. I thought that this might be improved even more by changing the instruction from movups to movaps, self-modifying the instruction after checking the address’ alignment. The results were not satisfactory, though : differences in performance were negligible, around 0.041%, which most certainly can’t be considered an improvement considering all the extra code that must be added in order to allow the function to self-modify itself. So, the good news are : you can always use movups for memory transfers in AVX, since it doesn’t introduce any kind of performance penalty when the pointer is aligned, and doesn’t fail if it’s unaligned. Just keep the addresses aligned yourself and you’ll be fine.

Here’s the source of the implementation. The assembly is written using the AMD64 ABI, was compiled under gcc 4.6.3 and nasm 2.09.10. The program actually using the function is written like a benchmark : it performs 10000 executions of the function with precision 5000007, randomizing the parameters beforehand, and sums the number of ticks returned by rdtsc, writing the average number of ticks from all the iterations to stdout on exit.

Regarding my last joystick-related project : had to put it on hold until further notice. This is the kind of thing that bites you on the ass once you leave it for some time and then want to come back to it, since coming back to it takes a lot of time. I’d like to continue it, obviously, but it doesn’t seem likely nowadays.

Opublikowano komputer | Otagowano , , , , , , | Dodaj komentarz

Windows 98 DDK (Driver Development Kit) and VXD joystick drivers

In the past few days, I started working on reverse-engineering a joystick driver written for Windows 98. It’s a six-button joystick which uses the ordinary (and, of course, obsolete by today’s standards) gameport. Since the gameport provides only four actual pins for representing buttons, drivers for devices that had more buttons had to use some sort of multiplexing. The main motivation behind this is to find out what kind of multiplexing is employed – and also, hopefully, write a driver for this joystick for Linux.

The first problem related to this was finding a utility that would read the contents of the VXD file containing the driver. The format of this file is the LX/LE „Linear Executable” executable format, which was used in OS/2 and early versions of Windows before the PE format took over. Fortunately, after some googling (and trying to write my own dumper for this format) I found a great utility called DUMPLX. It’s simply amazing, and features a full disassembler of the executable contents, and disassembles VMM/VxD calls as normal instructions, which was just what I was looking for.

The second problem was finding the definitions of the constants and structures used in the program. Microsoft doesn’t distribute old DDKs anymore, but I managed to get a hold of a copy of the Windows 98 DDK, which contains all the include files that I needed to work out what all the magic numbers mean. I put it up if on my mediafire account if anybody needs it : click.

Now, back to disassembling…

Opublikowano Uncategorized | Dodaj komentarz

How to upload a song without a music video to YouTube?

Quite often, you might find yourself in a situation in which you want to upload a song to YouTube which doesn’t have an accompanying music video. YouTube, being – obviously – a video sharing site, doesn’t accept only audio files and will reject them when uploading. I was wondering if there was a simple way to create a video that would contain the album cover shown throughout the whole length of the audio track.

Of course, the most obvious solution to this would be to use a video editing tool like the Windows Movie Maker, Sony Vegas, or Adobe Premiere, but that seemed like an overkill for such a – supposedly – simple task. Fortunately, the ffmpeg project can do it all automagically in one command, which was just the solution that I was looking for.

But let’s have a look at what YouTube wants from us first so that we can upload our content in the highest possible quality. The advanced encoding specifications show that the highest audio bitrate is allowed in 720p video files and equals 384 kbps. The recommended audio codec is AAC-LC, so you can use the Nero AAC Encoder or FAAC in order to encode the audio file into AAC (hopefully, from a lossless source!).

If you’ve got your audio files ready, it’s time to start making them into YouTube-compatible 720p files so that the audio bitrate is preserved and the resulting upload has the same high quality audio track. That’s where the aforementioned ffmpeg comes in. Hold on to your seats, now.

ffmpeg -loop 1 -r 5 -i $IMAGE_FILE -i $AUDIO_FILE -c:v libx264 -preset slow -crf 18 -tune stillimage -c:a copy
-filter:v "[in] scale=-1:720, pad=1280:720:640-iw/2 [out]" -shortest output.mp4

This is the command that will take your $AUDIO_FILE, put it together with $IMAGE_FILE and produce output.mp4, which is a video file with 1280×720 resolution with the $IMAGE_FILE centered in its every frame encoded with x264, and $AUDIO_FILE as the audio track. The length of output.mp4 is equal to the length of $AUDIO_FILE.

Let us have a look at the parameters themselves, though :

  • -loop 1 states that the input files should be looped indefinitely. This is so that the video track consists of repeated frames and not just one single frame, which is what ffmpeg would do by default – and YouTube doesn’t like it, since it requires video and audio tracks to be of equal lengths.
  • -r 5 specifies the frame rate for the video track. Since our consists of just a still image, it doesn’t make sense to have the video track run at a full 25 / 23.976 fps, which would just unnecessarily enlarge the size of the resulting file. Having -r 1 makes the video track much longer than the audio track, for some bizarre reason.
  • -i $IMAGE_FILE and -i $AUDIO_FILE are both input file declarations. The format of the files and whether they’re video or audio files is detected by ffmpeg at runtime.
  • -c:v libx264 chooses libx264 as the codec used to encode the video track. Thus, the resulting data is encoded in YouTube-compatible H.264 format.
  • -preset slow -crf 18 -tune stillimage are options passed to the x264 encoder. If you want to get to the bottom of them, I suggest reading this page.
  • -c:a copy specifies that no processing whatsoever should be done to the audio track and that the track in the resulting file should be an exact copy of the input audio track.
  • -filter:v „[in] scale=-1:720, pad=1280:720:640-iw/2 [out]” is the filter specification for the output video track. Here, we specify that the input image should first be pre-scaled to be 720 pixels high, while preserving the aspect ratio. Then, we pad the scaled image (assuming „normalized” album cover dimensions, we’re now dealing with a 720×720 image) with black borders to the size of 1280×720, and put the actual image at (640-iw/2), which happens to be the middle of the image. iw is the width of the input, e.g. what’s received from the scale filter.
  • -shortest tells ffmpeg to stop encoding when the shortest encoded track ends. In our case this is the audio track, since the image is looped indefinitely.
  • output.mp4 is just the output file name. ffmpeg automagically performs all the necessary muxing and handles output to the MP4 container.

That’s pretty much all there is to it. The resulting videos can be uploaded to YouTube without any further hassle. If you’re curious how the effect looks, here’s one of the videos that I prepared that way : click.

On a side note : it is possible to skip the first step, which consists of preparing the audio files manually with an external AAC encoder, and just have it all done by ffmpeg, thus simplifying the process even more. However, support for encoding AAC in ffmpeg is experimental and the sound quality of the result usually leaves much to be desired. If you would like to try it out, though, then use the following command.

ffmpeg -loop 1 -r 5 -i $IMAGE_FILE -i $AUDIO_FILE -c:v libx264 -preset slow -crf 18 -tune stillimage -c:a aac
-b:a 384k -filter:v "[in] scale=-1:720, pad=1280:720:640-iw/2 [out]" -shortest output.mp4

Hope you found this useful,

Opublikowano komputer | Otagowano , , , , , , | Dodaj komentarz

Hey 2012-02-12 Klub Eter, Wrocław [FLAC/MP3]

Koncert był absolutnie fenomenalny. Wszystkim fanom zespołu polecam ściągnięcie nagrania, a szczególnie posłuchanie przewspaniałej moim skromnym zdaniem wersji „Mimo wszystko”.

Nagranie dostępne w wersji MP3-V2 (na mediafire) oraz FLAC (jako torrent).

Info :

12 February 2012
Klub Eter, Wrocław
Audience Recording
Lineage : Church Audio CAFS > Tascam DR-07 (24bit/48kHz) > Sound Forge 10
(normalization, conversion to 16bit/44.1kHz, track splitting) >
FLAC v1.2.1 (-V -8) > foobar2000 v1.1.11 (tagging, ReplayGain scanning)
Lineup :
Katarzyna Nosowska - vocals
Paweł Krawczyk - guitar, trombone, backup vocals
Marcin Żabiełowicz - guitar
Jacek Chrzanowski - bass
Robert Ligiewicz - drums
Tracklisting :
01 - Fate [3:04]
02 - Muka [3:37]
03 - Umieraj stąd [3:55]
04 - A Ty? [3:47]
05 - That's a Lie [3:58]
06 - [sic!] [3:46]
07 - Cudzoziemka w raju kobiet [3:33]
08 - Piersi ćwierć [5:02]
09 - Miłość! Uwaga! Ratunku! Pomocy! [4:40]
10 - Faza delta [4:37]
11 - Z rejestru strasznych snów [3:23]
12 - Mimo wszystko [4:03]
13 - Stygnę [5:14]
14 - Teksański [2:59]
15 - Luli lali [4:41]
16 - Kto tam? Kto jest w środku? [3:30]
17 - Zazdrość [2:06]
18 - Vanitas [5:17]
19 - encore break [2:11]
20 - Boję się o nas [5:17]
21 - Cisza, ja i czas [3:28]
22 - Schizophrenic Family [2:36]
23 - Moja i Twoja nadzieja [4:27]
Info :
A concert from the 2012 winter tour by Hey. Very good setlist. Unfortunately, I
was a bit too late to get a good place and I was standing in a rather chatty
crowd, which unfortunately shows in the quieter passages and songs. There is
also a cut in "Kto tam? Kto jest w środku?" around 2:53 - sadly, it comes from
the master. The recording is rather good overall, so I hope you can enjoy it.
Please don't sell this recording in any form.

W razie wygaśnięcia linka lub torrenta proszę o kontakt.

Opublikowano życie | Otagowano , , , , , , | Dodaj komentarz

Symulowane wyżarzanie w asemblerze x86 (AMD64 + SSE)

Program, który powstał w sumie jako zakład z jednym z kolegów na kierunku. Celem było napisanie algorytmu symulowanego wyżarzania w asemblerze, bez łączenia z biblioteką C. Naturalnie, podjąłem ten zakład. Pierwszym problemem było opracowanie w jakiś sposób obliczania funkcji e^x, co w sumie sprawiło, że projekt utknął w miejscu na bardzo długi czas. Po bardzo dużej ilości googlowania i przeglądania zestawu instrukcji udało się jednak tego dokonać.

Nie było też do dyspozycji generatora liczb pseudolosowych, przez co byłem zmuszony do zastosowania bardzo prostego (i pewnie mało efektywnego pod względem prawdziwej losowości) algorytmu, którego kiedyś znalazłem na Wikipedii (teraz nie umiem go znaleźć).

Jedyną funkcją udostępnianą jest funkcja annealing() o następującym prototypie :
void annealing(float Tstart,float Tmin,float alpha,uint32_t *sol,int solsize,int(*crit)(uint32_t*,int));

Tstart, Tmin i alpha to parametry samego procesu wyżarzania (opisane w jednym z artykułów), *sol to wskaźnik na tablicę zawierającą solsize elementów 32-bitowych, które są początkowym rozwiązaniem algorytmu. Do oceny rozwiązania wykorzystywana jest funkcja przekazywana za pomocą wskaźnika *crit, przyjmująca wskaźnik na tablicę i jej rozmiar, oraz zwracająca liczbę całkowitą będącą oceną przekazanego do niej rozwiązania. Przed wywołaniem funkcji należy również wypełnić pola dwuelementowej tablicy seed, która jest używana jako ziarno dla generatora liczb pseudolosowych.

Kod funkcji został napisany zgodnie z ABI dla Linuksa AMD64 z rozszerzeniami SSE (jako że każdy procesor obsługujący operacje 64-bitowe ma też SSE). Po modyfikacjach na pewno będzie bezproblemowo działać na zwykłych i386.

Jako przykładowe zastosowanie tej funkcji napisałem krótki program w C, który znajduje rozwiązanie problemu TSP dla instancji gr17. Wątpię, żeby moja implementacja była w jakikolwiek sposób szybsza bądź lepsza od analogicznej w C – pewnie jest nawet gorsza. Program powstał po to, żeby pokazać, że można.

Źródełka tutaj. Do kompilacji wymagany NASM. Public domain.

Opublikowano komputer | Otagowano , , , , | Dodaj komentarz

Anonimowe udostępnianie plików i folderów w Windows XP Pro

Dzisiaj prezentacja rozwiązania problemu, nad którym zastanawiałem się przez dość długi czas, mianowicie – jak sprawić, żeby (oczywiście nie wszystkie, ale chociażby wybrane) udostępniane w Windows XP Professional foldery były dostępne dla osób, które nie mają na tym komputerze konta. Rozwiązanie takie jest bardzo wygodne w małych sieciach lokalnych, jeśli chcemy udostępnić jakiś określony zasób każdemu, kto ma do tej sieci dostęp.

W Windows XP Home Edition zrobienie czegoś takiego jest banalnie proste, co wynika z okrojenia edycji Home z wszelkich opcji zabezpieczeń, powiedzmy, „enterprajsowych” – chociaż tutaj bardziej po prostu chodzi o stosowanie zwykłego schematu wywodzącego się prosto z Windows NT, bez żadnych jego uproszczeń.

Ale do rzeczy. Zaznaczam od razu, że wszelkie opisy tyczą się anglojęzycznej wersji systemu z zainstalowanym SP3. Komputer ten nie jest przyłączony do domeny i pracuje jako część grupy roboczej. Do naszego problemu mamy dwa rozwiązania :

  • przydzielić każdemu użytkownikowi łączącemu się z zasobem dokładnie takie same prawa, bez względu na to, kim on jest,
  • przydzielić inne uprawnienia użytkownikom anonimowym, a inne tym, którzy łączą się podając poprawną kombinację user/pass.

Pierwsze rozwiązanie jest idealną kalką „prostego udostępniania plików” znanego z Home Edition. Możemy je również wprowadzić na Professional, i to na dwa sposoby – pierwszym jest zaznaczenie odpowiedniej opcji w Opcjach folderów w Panelu sterowania.

Drugim sposobem natomiast (który nie spowoduje ograniczenia dostępnych opcji dotyczących zabezpieczeń zasobów współdzielonych oraz zabezpieczeń samych plików i folderów) jest zmodyfikowanie odpowiedniego ustawienia w Local Security Policy, gałąź Local Policies oraz Security Options nazwanego Network access: Sharing and security model for local accounts.

Ustawienie to ma dwie możliwe wartości, których chyba nie trzeba chyba dodatkowo wyjaśniać :

Jeśli chcemy, aby każdy łączący się z komputerem był rozpoznawany jako Guest (i otrzymywał tylko i wyłącznie jego uprawnienia), wybieramy drugą opcję.

Jedna z powyższych ścieżek wystarczy nam, aby zaimplementować „proste udostępnianie plików”. Jeśli szanowny czytelnik chce jednak mieć dokładniejszą kontrolę nad udostępnianymi zasobami, uprasza się o dalsze czytanie.

A więc : myśląc w ogóle o udostępnianiu anonimowym czegokolwiek, musimy mieć włączone konto gościa. Możemy tego dokonać przez aplet User Accounts w panelu sterowania, lub za pomocą ustawienia Accounts: Guest account status w Local Security Policy.

Musimy też odpowiednio zmodyfikować prawa, jakie posiada gość, a mianowicie jedno – prawo do dostępu do komputera z sieci, które domyślnie w XP Pro temu użytkownikowi jest odebrane. Znowu z pomocą przychodzi Local Security Policy :

Jeśli widzimy, że jest tam wspomniany „Guest”, to oczywiście usuwamy go. Teraz jest pora na modyfikację trzech kluczowych ustawień w Local Security Policy, mianowicie :

  • Network access: Let Everyone permissions apply to anonymous users
  • Network access: Named pipes that can be accessed anonymously
  • Network access: Shares that can be accessed anonymously

Pierwsze z nich zmieniamy oczywiście na „Enabled”. W przeciwnym razie połączenia od użytkowników anonimowych nie będą mapowane do uprawnień nadanym grupie „Everyone”, przez co nie będą miały żadnych uprawnień i dla jakiegokolwiek zasobu dostaniemy po prostu „Access denied”.

Drugie ustawienie powinniśmy zmienić – a mianowicie dopisać tam srvsvc tylko, jeśli chcemy, żeby anonimowi użytkownicy mieli możliwość wylistowania wszystkich udostępnianych przez komputer zasobów. Jeśli tego nie zrobimy, użytkownik będzie musiał znać nazwę zasobu – niemożliwe będzie pobranie nazw wszystkich, do których można się ewentualnie dostać.

Trzecie ustawienie musimy zmieniać za każdym razem, gdy tworzymy lub usuwamy nowy anonimowo udostępniany zasób. Jeśli go tutaj nie dodamy, to nawet pomimo nadania uprawnień grupie „Everyone” anonimowi użytkownicy nie będą mieli do niego dostępu!

Local Security Policy możemy już zamknąć. Teraz wystarczy utworzyć zasób, pamiętając o uprawnieniach, jakie musimy nadać. Po wybraniu naszego katalogu wchodzimy więc do „Sharing” w jego właściwościach :

Nadajemy mu żądaną przez nas nazwę i ewentualnie jakiś elokwentny komentarz, a następnie modyfikujemy uprawnienia :

Wychodzimy z okienka zatwierdzając zmiany i wracamy do właściwości naszego folderu. Musimy jednak jeszcze nadać uprawnienia do dostępu do niego z poziomu samego systemu plików, a nie tylko jako udostępnianego zasobu. W tym celu przechodzimy na zakładkę obok we właściwościach.

Widzimy tutaj, że „Everyone” nie ma zupełnie praw dostępu do tego katalogu. Tak więc nawet jeśli anonimowemu użytkownikowi udałoby się uzyskać dostęp do zasobu (a teraz już by się udało), nie mógłby on w ogóle odczytać folderu reprezentowanego przez ten zasób. Dodajemy więc grupę „Everyone”, modyfikujemy uprawnienia, jakie chcemy nadać, i zatwierdzamy zmiany.

Należy w tym momencie zaznaczyć, że uprawnienia przydzielone do zasobu i do folderu to dwie różne rzeczy. Mając pełne uprawnienia do zasobu, lecz read-only dla folderu, nic w nim nie zapiszemy – i na odwrót. Dla przykładu, logując się do tego zasobu jako „Daniel” miałbym tylko uprawnienia do jego odczytu – ponieważ do zasobu tego ma tylko uprawnienia „Everyone”, i są to uprawnienia RO. Nikt nam jednak nie broni dodać Daniela do listy uprawnień zasobu i nadać mu ich więcej.

Nie zapominamy oczywiście do dodaniu naszego nowego zasobu do listy anonimowo dostępnych w Local Security Policy.

To wszystko. :) Zasób taki jest teraz dostępny anonimowo, czyli może się na niego zalogować ktokolwiek, podając zupełnie dowolną nazwę użytkownika, która nie występuje na tym komputerze, i zupełnie dowolne hasło. Jeśli nazwa już występuje – będziemy zobowiązani podać hasło dla tego użytkownika.

Co najciekawsze, Google jak i dokumentacja MS na temat anonimowego udostępniania w XP Pro z grupą roboczą niewiele potrafią powiedzieć. W dojściu do wielu z tych rzeczy o wiele bardziej pomógł mi Linux, a konkretniej mówiąc – smbclient. Wstydź się, Bill.

Opublikowano Uncategorized | Dodaj komentarz

Symulowane wyżarzanie

Do czytelników : artykuł ten jest głównie adresowany do kolegów i koleżanek z kierunku, jednak powinien zainteresować każdego szukającego informacji o tym algorytmie metaheurystycznym.

Jako że wczoraj nie poszedłem na wykład dr. Cabana z systemów operacyjnych i męczą mnie z tego powodu wyrzuty sumienia, postanowiłem swój czas poświęcić na czyn społeczny i napisać krótki, acz konkretny artykuł na temat symulowanego wyżarzania. Przedstawione tutaj informacje będą trochę miksem notatek z wykładu dr. Kozika i moich własnych doświadczeń.

Algorytm symulowanego wyżarzania jest algorytmem wyszukującym przybliżone rozwiązanie konkretnej instancji jakiegoś problemu NP-trudnego. Jak wie każdy mniej więcej ogarnięty, przestrzeń możliwych rozwiązań problemów NP-trudnych jest z reguły na tyle duża, że odpada jej dokładne przebadanie, tj. sprawdzanie każdego możliwego rozwiązania po kolei. Dla problemu komiwojażera na przykład mamy \frac{(n-1)!}{2} rozwiązań, co już dla 20 miast daje liczbę na tyle dużą, że wykonanie takiego algorytmu zajmie niewiarygodnie dużo czasu. Wykorzystując algorytmy przybliżone jesteśmy w stanie znaleźć rozwiązanie, które może i nie będzie najlepszym (optymalnym), jednak będzie „wystarczająco dobrym”, a do tego zostanie ono znalezione w czasie w miarę rozsądnym.

Algorytmy metaheurystyczne zawsze działają inaczej, ponieważ w swoim działaniu opierają się na losowości – najbardziej ogólny algorytm metaheurystyczny będzie wyglądał tak :

while(i < limit_iteracji)
	wylosuj nowe rozwiązanie;
	jeśli jest lepsze od poprzednio wylosowanego, zapisz je;

Taki „naiwny” algorytm jednak działa bardzo kiepsko, ponieważ za każdym razem przeskakuje w zupełnie inny punkt obszaru możliwych rozwiązań, a możliwe jest, że jeśli jesteśmy w jakimś punkcie, to rozwiązanie, które będzie od niego lepsze znajduje się zaraz obok niego – dlatego zdefiniujemy coś takiego jak otoczenie rozwiązania jako rozwiązanie mu sąsiadujące na przykład przez zamieniony jeden element. Czyli sąsiedztwem permutacji { 3, 7, 1, 2, 4, 5, 6 } może być permutacja { 3, 7, 1, 4, 2, 5, 6 }. Pierwszym pomysłem powinno być, by algorytm dostawał „jakieś” rozwiązanie początkowe (chociażby wylosowane) i działał tak długo, aż w sąsiedztwie nie znajdują się żadne lepsze rozwiązania. Jest z tym związany jeden mankament – aby to zobrazować, powiedzmy, że chcemy zminimalizować poniższą funkcję na obszarze ograniczonym wykresem.

Na pierwszy rzut oka widać, że funkcja ma globalne minimum w x=1 i lokalne minimum w okolicach x=-0.8. Jeśli rozwiązaniem początkowym naszego algorytmu będzie na przykład x = -0.25, to nie ma możliwości, by kiedykolwiek dotarł on do globalnego minimum, ponieważ nigdy nie będzie się ono znajdowało w sąsiedztwie badanego rozwiązania.

Modyfikacją takiego „naiwnego” algorytmu jest właśnie symulowane wyżarzanie, w którym zezwalamy na zmianę rozwiązania z pewnym prawdopodobieństwem, nawet wtedy, gdy jest ono gorsze – bo możliwe, że w sąsiedztwie gorszego będą lepsze niż te, które mamy w sąsiedztwie rozwiązania aktualnie rozpatrywanego. Dzięki temu symulowane wyżarzanie potrafi „unikać” problemów minimum lokalnego. Poza tym, cytując dr. Kozika „Własność losowości jest taka, że po pewnym czasie daje dobre wyniki”, co tutaj ma chyba jakiś sens.

Algorytm cechuje się następującymi parametrami :

  • początkową temperaturą T,
  • końcową temperaturą Tmin,
  • funkcją prawdopodobieństwa P,
  • funkcją obniżania temperatury G.

Funkcja prawdopodobieństwa określa, z jakim prawdopodobieństwem będziemy zmieniać rozwiązanie na gorsze. W „oryginalnej” odmianie symulowanego wyżarzania funkcja ta była zależna od oceny rozwiązania nowego, starego i aktualnej temperatury. Bierze się to z analogii tego algorytmu do procesu metalurgicznego – póki temperatura jest wysoka, akceptujemy zmiany na rozwiązania początkowo gorsze, ale potencjalnie lepsze. Zakładamy, że w miarę postępowania algorytmu i obniżania temperatury dojdziemy do rozwiązania „w miarę dobrego”, tak więc „ryzykowne” zmiany na rozwiązania początkowo gorsze nie będą już potrzebne. Z tego wynika też funkcja obniżania temperatury – wykonywana jest ona po każdej iteracji algorytmu i w oryginalnej odmianie była zależna od aktualnej temperatury i aktualnego numeru iteracji.

Po przeczytaniu tego wszystkiego poniższy pseudokod powinien być w miarę oczywisty.

i = 0
wybierz losowo rozwiązanie A
while(T > Tmin)
	stwórz rozwiązanie B przez zamianę dwóch losowych elementów w A
	if(Kryterium(B) < Kryterium(A)) A = B
	else if(random(0,1) < P()) A = B
	T = G()
	i = i+1

Jak widać, sama idea algorytmu jest bardzo prosta, jeszcze prostsze jest to, jak to zastosować w przypadku problemów przez nas rozpatrywanych – zarówno w przypadku komiwojażera, jak i szeregowania wystarczy po prostu zamienić ze sobą dwa elementy aktualnie rozpatrywanej permutacji i przeliczyć dla nich kryterium. I to wszystko.

Jedynym problemem jest dobranie funkcji P i G tak, aby algorytm dawał dobre rozwiązania. W przypadku P odpowiedź jest prosta :

P = e^{-\frac{f(S') - f(S)}{T}}

Zakładając, że e to stała Eulera, f() to funkcja określająca kryterium rozwiązania, S’ to „nowe” rozwiązanie, S to „stare” rozwiązanie, a T to aktualna temperatura. Natomiast funkcja G :

G = \frac{\Gamma}{\log(k + k_0)} lub G = \alpha T_k

Gdzie gamma, alpha i k0 to „jakieś stałe” (twierdzenie, z którego pochodzi ten wzór mówi tylko tyle, że „istnieją”), k to aktualna iteracja algorytmu, a Tk to aktualna temperatura.

Dla zainteresowanych jako „further reading” polecam :

Mam nadzieję, że coś rozjaśniłem. Proszę rzucać pytania, jeśli się pojawią.

Opublikowano komputer | 8 Komentarzy