Tuesday, July 25, 2017

Gayrı-Lineer Dinamik - Ders 16

Wednesday, June 28, 2017

Disk Temizliği

Arada sırada Ubuntu / Unix'imizde "bahar temizliği" iyi olabilir; ardı ardına kurulan paketler, işletilen programların yarattığı onbellek (cache) dosyaları arkada kalıp diski bitirebilir. Temizlik için bazı programlar:

En basiti komut satırında du,

sudo du -hx --max-depth=1 /

Dizin bazında kullanım raporu verir. Benzer bir raporu görsel araç üzerinden dosya idarecisi Nemo ile görebiliriz, herhangi bir dizin üzerinde sağ tıklama ve Open With | Disk Usage Analyzer.

Otomatik temizleme icin bleachbit: Kurulum apt-get install bleachbit ile, gksudo ile başlatıp solda temizlenecek şeyler seçilir. APT temizliği iyi oluyor, diğer bir sürü seçenek te var.

Thursday, June 8, 2017

GraphHopper: Pür Java ile Haritada Yol Bulmak

Yol bulmak için OSRM projesinden bahsettik. OSRM çok hızlı çalışır, sonuçları kullanışlı, fakat C++ ile yazılmış, eğer pür Java bazlı bir çözüm istiyorsak (ki böylece kodları mesela Android'e rahatça koyabilelim) o zaman Graphhopper projesi var. GH ismi bir kelime oyunu, grasshopper bilinebileceği gibi çekirge demek, graph (yani çizit) hopper "bir çizitte sağ sola zıplayan" bir görüntü zihinde canlandırıyor, ki bu görüntü pek gerçekten uzak değil, haritada yol bulma algoritmaları hakikaten bilgisayar bilimdeki çizit yapısını temel alıyorlar. Kod şurada,

https://github.com/graphhopper/graphhopper

Kodun derleme sistemi gradle / Maven bazlı, fakat o kadar çetrefil derleme işlemlerine gerek yok, önce GH içinden gerekli java kodlarını çıkartalım, sonra çok baz bir derleme sistemi ile ve ek birkaç jar ile derlemeyi kendimiz dışarıda halledelim. Alternatif bir proje / dizin yaratılır, altında lib, src, resources dizinleri olur. GH indirilir, şimdi alternatif dizin altında

cp -r [GH]/core/src/main/java/* src/
cp -r [GH]/core/src/main/resources/* resources
cp -r [GH]/reader-osm/src/main/java/* src/

ile sadece gerekli Java kodlarını alırız. Şimdi gereken jarları alalım, 


Bu jarları nasıl bulduk? Ya derlerken olmayan jar için gelen hata mesajlara bakıp bu class ismini "vs.vs.Class maven jar" ile Google'da ararız, ya da [GH]/pom.xml ya da hangi pom.xml var ise versiyon sayısını oradan buluruz ve yine Google'da ararız. Mesela jackson-databind lazım, "jackson-databind maven jar" bize gerekli bağlantıyı verir. Zaten görüldüğü gibi bağlantılarda belli bir kalıp da var, class ismi ve versiyon ile bağlantı tahmin de edilebilir.

Jar'lar lib altına gidiyor tabii. Derleme için Ant build.xml

<project name="gh" default="dist" basedir=".">
  <property name="src" location="src"/>
  <property name="build" location="build"/>
  <property name="dist" location="bin"/>
  <target name="init">
    <tstamp/>
    <mkdir dir="${build}"/>
  </target>
  <path id="build.classpath">
    <fileset dir="lib">
      <include name="**/*.jar"/>
    </fileset>
  </path>
  <path id="compile.classpath">
    <pathelement location ="${build}"/>
    <fileset dir="lib">
      <include name="**/*.jar"/>
    </fileset>
  </path>
  <target name="compile" depends="init"
        description="compile the source">
    <javac destdir="${build}">
      <src path="${src}"/>
    <classpath refid="build.classpath"/>
  </javac>
  </target>
  <target name="resources">
    <copy todir="${build}" includeEmptyDirs="no">
      <fileset dir="resources">
        <patternset>
          <include name="**/*.*"/>
        </patternset>
      </fileset>
    </copy>
  </target>
  <target name="dist" depends="compile" description="generate the distribution">
    <mkdir dir="${dist}"/>
    <jar jarfile="${dist}/ghopper.jar" basedir="${build}"/>
  </target>
    <target name="test" depends="compile">
      <java fork="yes" classname="Test" failonerror="true">
        <classpath refid="compile.classpath"/>
      </java>
  </target>
  <target name="clean"
        description="clean up">
    <delete dir="${build}"/>
    <delete dir="${dist}"/>
  </target>
</project>

Bu kadar, ant ile derleriz. Kodların işleyişi için OSM harita dosyaları gerekli, bu dosyalar daha önce de bahsedilen 


sitesinden alınabiliyor. Mesela Türkiye için 


Bu dosyayı bunzip2 ile açıp gzip ile tekrar zipleyelim, çünkü GH .gz dosyası istiyor. Ya da .osm dosyası olduğu gibi bırakılır, GH onu da işler.

GH sürekli direk .gz dosyasını da kullanmaz bu arada, onu işleyip bazı ara dosyalar üretmesi lazım (daha önce bahsedilen çizit dosyaları bunlar, GH OSM formatını alıp daha hızlı işleyiş için çizit yapısına geçiyor). Ara dosyaların üretilmesi ve yol bulma testi çin src altında Test.java yaratalım, 

import com.graphhopper.*;
import com.graphhopper.routing.util.EncodingManager;
import com.graphhopper.reader.osm.GraphHopperOSM;

public class Test {

    static String tmpGraphFile = "/tmp/gh";

    public static void createGraph() {
        String fin = "[DB DIZINI]/turkey-latest.osm.gz";
        GraphHopper gh = new GraphHopperOSM().setStoreOnFlush(true).
            setEncodingManager(new EncodingManager("foot")).setCHEnabled(false).
            setGraphHopperLocation(tmpGraphFile).
            setDataReaderFile(fin);
        gh.importOrLoad();
    }

    public static void testPath() {
        GraphHopper gh = new GraphHopperOSM().
            setEncodingManager(new EncodingManager("foot")).setCHEnabled(false).
            setGraphHopperLocation(tmpGraphFile);
        gh.importOrLoad();
        GHResponse res = gh.route(new GHRequest(40.987659, 29.036428, 
                                                40.992186, 29.039228).
                                 setVehicle("foot"));
        System.out.println(""+res );
    }

    public static void main(String[] args) throws Exception{
        createGraph();
        testPath();
        testPath();
    }
}

Isletmek icin ant test. Araba tarifleri için setEncodingManager ve setVehicle çağrılarında car kullanılır. Kodda createGraph'ın sadece bir kez çağırılması yeterli, bu çağrı /tmp/gh altında gerekli çizit dosyalarını yaratır (çağrı birkaç dakika sürebilir), diğer tüm yol bulma çağrıları bundan sonra çizit dosyalarını kullanıp hızlı şekilde cevabı üretir. Artık elimizde pür Java bazlı, çok hızlı işleyen, mobilde İnternet bağlantısız çalışabilecek harita tarif kodları var. Üstteki kodların kendisi fazla yer tutmuyor zaten (biraz da bunun için daha büyük GH içinden sadece gerekli java ve jar dosyalarını çekip çıkardık, az sayıda java ve 10 tane jar dosyasi ile elimize yol tarifi yapabilen çekirdek bir kod geçti), çizit dosyaları ise bir şehrin büyüklüğüne göre 10-20 MB arasında olur.

Kodun teorik, algoritmik yapısını merak edenler için: kısa yol algoritmaları çoğunlukla Dijkstra algoritmasını kullanırlar (Dijkstra'yı bilgisayar bilim notlarında isledik), OSRM ve Graphhopper da böyle yapıyor, yanlız Dijkstra yaklaşımına bazı ekler yapıyorlar, yaklaşımın adı "kısalan hiyerarşiler (contracting hierarchies)"; yol ağ yapısına bazı noktalarda haritaya dışarıdan kısa yol / zıplamalar ekliyorlar, böylece Dijkstra daha hızlı, verimli işliyor.

Dil, derleme hakkında bazı bilgiler:  

Monday, June 5, 2017

Zaman Serilerinde Tepe Noktası Bulmak (Peak Detection)

Matlab / Octave kullananlar zaman serisi tepe noktası analizinde peakutils adlı bir aracı kullanıyorlar. Bu kod Python'a taşınmış,

https://github.com/atjacobs/PeakUtils

Kod içinden gerekli fonksiyonu çıkarttık,

import numpy as np

def indexes(y, thres, min_dist):
    thres *= np.max(y) - np.min(y)
    dy = np.diff(y)
    peaks = np.where((np.hstack([dy, 0.]) < 0.)
                     & (np.hstack([0., dy]) > 0.)
                     & (y > thres))[0]

    if peaks.size > 1 and min_dist > 1:
        highest = peaks[np.argsort(y[peaks])][::-1]
        rem = np.ones(y.size, dtype=bool)
        rem[peaks] = False

        for peak in highest:
            if not rem[peak]:
                sl = slice(max(0, peak - min_dist), peak + min_dist + 1)
                rem[sl] = True
                rem[peak] = False

        peaks = np.arange(y.size)[~rem]

    return peaks

Parametreler tepe noktalarının arasında en az ne kadar mesafe, ayrıca noktaların en az ne kadar genliğe (y yönünde) sahip olmaları gerektiği. Eldeki bir veri üzerinde indexes(d,thres=0.3,min_dist=3) çağrısı sonrası

Java için benzer bir çağrı, alttaki kütüphaneden,

https://github.com/JorenSix/TarsosDSP

import java.util.*;

public static LinkedList<Integer> indexes(double[] data,
                                          int width,
                                          double threshold,
                                          double decayRate,
                                          boolean isRelative) {

    LinkedList<Integer> peaks = new LinkedList<Integer>();
    int maxp = 0;
    int mid = 0;
    int end = data.length;
    double av = data[0];
    while (mid < end) {
        av = decayRate * av + (1 - decayRate) * data[mid];
        if (av < data[mid])
            av = data[mid];
        int i = mid - width;
        if (i < 0)
            i = 0;
        int stop = mid + width + 1;
        if (stop > data.length)
            stop = data.length;
        maxp = i;
        for (i++; i < stop; i++)
            if (data[i] > data[maxp])
                maxp = i;
        if (maxp == mid) {
            if (overThreshold(data, maxp, width, threshold, isRelative,av)){
                peaks.add(new Integer(maxp));
            }
        }
        mid++;
    }
    return peaks;
}

public static boolean overThreshold(double[] data, int index, int width,
                                    double threshold, boolean isRelative,
                                    double av) {

    int pre = 3;
    int post = 1;

    if (data[index] < av)
        return false;
    if (isRelative) {
        int iStart = index - pre * width;
        if (iStart < 0)
            iStart = 0;
        int iStop = index + post * width;
        if (iStop > data.length)
            iStop = data.length;
        double sum = 0;
        int count = iStop - iStart;
        while (iStart < iStop)
            sum += data[iStart++];
        return (data[index] > sum / count + threshold);
    } else
        return (data[index] > threshold);
}

Güzel. Fakat bu tepe noktası çağrıları toptan şekilde çalışıyor, verinin tamamını alıyorlar, fonksiyon verinin tamamına bakiyor. Soru şu: canlı zamanda, sürekli akan veri üzerinde tepe noktası nasıl buluruz? En son N tane noktayı hatırlamak mümkün olsun diyelim (ama verinin tamamı değil).

İlk akla gelebilecek basit çözüm bir genlik alt limiti ymin'i aşan tepe noktasını seçmek, bu nokta ardından xmin kadar bekleriz, ylim'den üstte yeni bir tepe nokta var ise, bunu da alırız, böyle devam ederiz. Fakat basit çözüm biraz aceleci aslında, belki xmin sonrası bulunan tepe noktasına çok yakın "daha iyi (yani daha yukarıda)" bir nokta daha var, tipik olarak hemen bir önceki noktanın üzerine atlamak istemiyoruz, daha yukarıdaki noktayı istiyoruz.

Bu durum için yine ilk akla gelen çözüm xmin'i arttırmak olabilir, fakat global bir parametreyle bu şekilde oynamak çoğunlukla iyi sonuç vermez, büyük ihtimalle analizin diğer yerlerinde başka yanlışlara yol açar.

Çözüm yaklaşımı değiştirip veriler için bir kayan pencere kullanmak. Tepe noktası "en son N öğenin içinde orta noktanın maksimum olması demek" diye bir tanım getirebiliriz. Kayan pencere  için CircularFifoQueue kullanılabilir, bu bir org.apache.commons sınıfı, bir tür sınırlandırılmış kuyruk (queue). N öğe tutsun diye tanimlarsak N öğe ardından verilen N+1'inci öğe için 1. öğe dışarı atılır. O zaman tepe noktası bulan kod şöyle olabilir

import java.io.*;
import java.util.*;
import java.util.Queue;
import org.apache.commons.collections4.queue.CircularFifoQueue;

public class TestPeak {

    public static class PeakFinder {
        double xmin = 0;
        double ymin = 0;
        int mid = 0;
        CircularFifoQueue idxs = null;
        CircularFifoQueue vals = null;
        // "null" olmayan degeri temsil etmek icin kullaniyoruz
        // elinde process'e verecek id'si olmayan cagrilar icin
        // kullanisli olabilir
        public static int DUMMY = 1; 
        public PeakFinder(int xmin, double ymin){
            this.xmin = xmin;
            this.ymin = ymin;
            idxs = new CircularFifoQueue(xmin);
            vals = new CircularFifoQueue(xmin);
            this.mid = (int)vals.size() / 2;
        }
        // process hem id hem deger aliyor, fakat id aslinda
        // tutulup oldugu gibi geri veriliyor, cagri yapana
        // yardimci olmasi icin bunu yapiyoruz, fonksiyonun ic
        // isleyisi icin onemli degil.
        public int process(double val, int idx) {
            vals.add(new Double(val));
            idxs.add(new Integer(idx));
            if (vals.size() < xmin) {
                return Integer.MIN_VALUE;
            }
            if (vals.get(mid) > this.ymin && vals.get(mid) == Collections.max(vals)) {
                return idxs.get(mid);
            }

           // Integer.MIN_VALUE null degeri yerine kullanildi, 
           // cunku int yerel tip, obje degil
            return Integer.MIN_VALUE;
        }
    }
        
    public static void main(String[] args) throws Exception{
        java.util.Random r = new java.util.Random();
        r.setSeed(0);
        double samples[] = new double[1024];
        for ( int n = 0; n < samples.length; n++ )
            {
                double noise = r.nextGaussian() * Math.sqrt(10);
                samples[n] = Math.sin( n * Math.PI / 18. ) + noise;
            }
            
        // uretilen veriyi dosyaya yazalim, ustteki Python kodu ayni veriyi
        // kullanacak.
        PrintWriter writer = new PrintWriter("out.txt");
        for (int i = 0; i < samples.length; i++) writer.println(samples[i]);
        writer.close();

        ArrayList res = new ArrayList();
        PeakFinder pf = new PeakFinder(6, 6.0);
        for ( int n = 0; n < samples.length; n++ ) {
            int idx = pf.process(samples[n], n);
            if (idx > Integer.MIN_VALUE) res.add(idx);
        }

        System.out.println(""+res );
    }
}

Sonuc,
Tekrarlamak gerekirse: tepe noktası "sağında ve solunda vadiler olan orta noktadır" tanımı üzerinden  sonucu aldık. Pencereyi sürekli kaydırıyoruz (daha doğrusu her yeni veriyi N tane sınırlı öğe içerebilen bir kuyruğa sürekli sondan ekliyoruz), ve bu sırada orta noktanın maksimum olup olmadığına bakıyoruz. 

JAMA - Java ile Matris Islemleri

Python ile prototip hatta servis tarafı sayısal işlem kodları rahatça yazılıyor; fakat bazen lineer cebir yapan kodları Java'ya çevirmek gerekebilir, mesela Android üzerinde işlemesi gereken kodlar. Pür Java ile yazılmış kullanışlı liner cebir kütüphanelerinden biri JAMA:

http://math.nist.gov/javanumerics/jama/#Package

Matrix Class Doc


İhtiyaç olan en önemli özellikler temel işlemler, toplama, çıkartma, ve matris arası çarpım, matris tersi (matrix inversion) ve devriği (transpose) Ek iyi olabilecek özellikler en az kareler (least squares) çözümü, SVD, LU, Cholesky ayrıştırması, vs. Bu özelliklerin hepsi Jama'da var.

Üstteki ilk bağlantıdan jar indirilir, altta kısa bir test,

import Jama.*;

public class Test {

    public static void main(String[] args) throws Exception{

        Matrix A = new Matrix(3, 3);

        A.set(1,0,10.0);
        A.set(1,1,10.0);
        A.set(1,2,10.0);

        System.out.println("A\n"+toString(A));

        System.out.println("A get\n"+A.get(1,1));

        double[][] bvals = {{1.,2.,3},{4.,5.,6.},{7.,8.,10.}};
        Matrix B = new Matrix(bvals);
        System.out.println("B\n"+toString(B));

        Matrix C = A.times(B);

        System.out.println("Carpim\n"+toString(C));

        Matrix D = A.plus(B);

        System.out.println("Toplam\n"+toString(D));

    }

    public static String toString(Matrix m) {
        String s = "";
        for (int i=0;i<m.getRowDimension(); i++){
            for (int j=0;j<m.getColumnDimension(); j++) s += "" + m.get(i,j) + " ";
            s += "\n";
        }
        return s;
    }
}

Sonuclar

A
0.0 0.0 0.0 
10.0 10.0 10.0 
0.0 0.0 0.0 

A get
10.0
B
1.0 2.0 3.0 
4.0 5.0 6.0 
7.0 8.0 10.0 

Carpim
0.0 0.0 0.0 
120.0 150.0 190.0 
0.0 0.0 0.0 

Toplam
1.0 2.0 3.0 
14.0 15.0 16.0 
7.0 8.0 10.0 

Bazı püf noktaları:

Android ortamı her ne kadar tam gömülü (embedded) bir ortam sayılmasa da, servis tarafı kodlanıyormuş gibi davranmaktan kaçınmak iyi olur, mesela telefon ivmeölçerinden gelen verileri hızlı bir şekilde işlemek gerekiyor, her işlem için bir matris çarpımı lazım, bu durumda her yeni veri iletimi için yeni matrisler yaratıp onları çarpmaya gerek yok, her seferinde hafızaya obje ekleme, onun çöp toplayıcısı tarafından silinmesi işleri ağırlaştırır. En iyisi sadece iki Matrix objesi yaratıp ama değerlerini her seferinde değiştirip çarpımı bu aynı objeler üzerinde yapmak.

Ek kodlar, mesela Jama kodlarında çapraz çarpım verilmemiş, bu benzeri ekler,

import Jama.Matrix;
import java.net.*;
import java.util.*;
import java.io.*;
import java.text.*;

public class Util {

    // Helper class hidden here to provide additional matrix methods
    // taken from https://github.com/thorstenwagner/ij-ellipsesplit
    public static class Matrix2  {
        private double[][] A;
        private int m, n;
        public Matrix2(int m, int n) {
            this.m = m;
            this.n = n;
            A = new double[m][n];
        }
        public Matrix2(double[][] A) {
            m = A.length;
            n = A[0].length;
            for (int i = 0; i < m; i++) {
                if (A[i].length != n) {
                    throw new IllegalArgumentException("All rows must have the same length.");
                }
            }
            this.A = A;
        }

        public Matrix2 crossProduct(Matrix2 B) {
            if (m != B.m || n != B.n) {
                throw new IllegalArgumentException("Matrix dimensions must agree");
            }
            Matrix2 X = new Matrix2(m, n);
            if (m == 1 && n == 3) {
                X.A[0][0] = A[0][1] * B.A[0][2] - A[0][2] * B.A[0][1];
                X.A[0][1] = A[0][2] * B.A[0][0] - A[0][0] * B.A[0][2];
                X.A[0][2] = A[0][0] * B.A[0][1] - A[0][1] * B.A[0][0];
            } else if (m == 3 && n == 1) {
                X.A[0][0] = A[1][0] * B.A[2][0] - A[2][0] * B.A[1][0];
                X.A[1][0] = A[2][0] * B.A[0][0] - A[0][0] * B.A[2][0];
                X.A[2][0] = A[0][0] * B.A[1][0] - A[1][0] * B.A[0][0];
            } else {
                throw new IllegalArgumentException(
                                                   "Matrix must be a 3-element row or column vector");
            }
            return X;
        }

        public double[][] getArray() {
            return A;
        }
    }

    public static String toString(Matrix m) {
        String s = "";
        for (int i=0;i<m.getRowDimension(); i++){
            for (int j=0;j<m.getColumnDimension(); j++){
                s += "" + m.get(i,j) + " ";
            }
            s += "\n";
        }
        return s;
    }

    // project u onto plane with normal n
    public static Matrix proj (double[] u, double[] n) {
        Matrix um = new Matrix(u, 3);
        Matrix un = new Matrix(n, 3);
        Matrix tmp1 = um.transpose().times(un);
        Matrix tmp2 = un.transpose().times(un);
        Matrix tmp3 = un.times(tmp1.get(0,0) / tmp2.get(0,0));
        Matrix res = um.minus(tmp3);
        return res;
    }

    public static Matrix cross(Matrix a, Matrix b) {
        Matrix2 aa = new Matrix2(a.getArray());
        Matrix2 bb = new Matrix2(b.getArray());
        Matrix2 res = aa.crossProduct(bb);
        return new Matrix(res.getArray());
    }
}