coding-experiments

Evaluating math expressions in F#

2009-10-02T10:34:00.000-07:00

This time we will try to write math expressions evaluator in F#. At first - why we need that ? Well, math expressions evaluator is useful in computer algebra systems - for evaluating some equations. Also such evaluator may be useful for financial institutions where also some kind of formula computation/evaluation is needed and etc. It`s a pitty that such eval function is not implemented in F# language. Possible reasons is that F# is compiled language... Anyway we will try to write such universal eval function which accepts arbitrary math expression and returns result to the user. So what are possible ways of implementing such function in F# ?

- by writing so called domain-specific language - math expressions compiler. This is achieved by using some tools such as fsyacc and fslex. Defining custom tokens and etc.

- by using .NET System.CodeDom.Compiler libraries. This is achieved by generating class definition "on the fly" and with the help of CodeDom - compiling this class in memory.

- another different approach - is try to write F# script which recursively analyzes math expression and evaluates it piece by piece until it finds the answer. This solution does not use tools mentioned above.

So the last approach will be discussed here. Main algorithm of evaluating math expressions is displayed below:

F# code which implements this idea is placed here.
There are also some eval() function tests implemented in this code. So that after executing this script you will get output something like:

Conclusions:

It is possible to generalize eval() function algorithm in about 90~ lines of F# code without using CodeDom model or fsyacc/fslex tools. This makes solution pretty much "homogeneous" which can be useful in some situations. Also this eval() implementation takes as input arbitrary math expression. As for the bad side - solution is very slow (can be hundreds of times slower than plain F# function call), because of joining/spliting lists constantly. So if you need fast implementation of eval - you should optimize/rewrite this code if possible or take a look in other above mentioned solutions. But in systems where speed issues are not critical - this solution may be OK.

Have fun with F# and eval !!

XNA quest for centroid of polygon

2009-09-17T23:07:00.000-07:00

Now we will talk about XNA Framework and How to find geometrical center of polygon (aka. so called centroid). At first - Why do we need centroid at all ? Well... for example - for some graphics applications we may want drag&drop shape by it`s centroid. Or maybe in some game applications we may need to find shape`s center of gravity, for simulating physics of balancing objects around center of gravity. And in cases when shape material density is uniform (the same over shape) then shape geometrical center (centroid) is also center of gravity (or center of mass). So lets begin search of centroid...
For this you will need:

- Visual Studio 2008
- XNA 3.1 Framework

First - some notes about XNA framework. It`s a pitty that such great 2D/3D .NET library doesn`t have 2D primitives drawing routines. Well, actually it has some, but for being able to draw some 2D primitives like points and lines, you should dig deeply in 3D knowledge - for drawing mentioned primitives in MSDN way you should define 3D vectors and 3D viewports !! What the hell ? If user only wants 2D stuff, it means that all 3D stuff should be abstracted out of the view. So that at the end, framework should have separate levels of abstraction into 2D and 3D. But this is not the case about drawing 2D primtives... So one way of drawing point in non-MSDN way is to generate small rectangle texture in run-time and draw this generated texture on screen. Next thing is drawing a simple line in non-MSDN way. Recipe for that - we take generated 2D rectangle texture earlier and draw it as line by scaling and rotating it as needed. After defining this stuff we have all prerequisites for target problem.
Now- How we will find centroid of polygon ? This is not that hard. According to wikipedia centroid formulas can be computed simply, assuming we have polygon vertex coordinates and assuming polygon is not complex.
So centroid computation algorithm in C# is this:

_______________________________________________________


public static Vector2 CalculateCentroid(Vector2[] points, 
                                        int lastPointIndex)
{
    float area = 0.0f;
    float Cx = 0.0f;
    float Cy = 0.0f;
    float tmp = 0.0f;
    int k;

    for (int i = 0; i <= lastPointIndex; i++)
    {
        k = (i + 1) % (lastPointIndex + 1);
        tmp = points[i].X * points[k].Y - 
              points[k].X * points[i].Y;
        area += tmp;
        Cx += (points[i].X + points[k].X) * tmp;
        Cy += (points[i].Y + points[k].Y) * tmp;
    }
    area *= 0.5f;
    Cx *= 1.0f / (6.0f * area);
    Cy *= 1.0f / (6.0f * area);

    return new Vector2(Cx, Cy);
}

___________________________________________

NOTE: this algorithm doesn`t apply for complex polygons.
But How simply check that polygon is not complex, or not self-intersecting ? One way of doing this is brute-force way -> we can check every polygon edge pair for intersection. If there are some edges which intersects, then polygon is complex, otherwise - simple. So by using this brute-force solution, complex polygon check simplifies to lines self-intersection problem. And How do we check if lines self-intersects ? One way is to find-out 2 linear equations of these lines: y = a1*x + b1 ; y = a2*x + b2.
Then self-intersection point is simply where these equations equals:
a1*x + b1 = a2*x + b2 --> x_intersection = (b2-b1)/(a1-a2). Substituting into any of these two equations we get y_intersection = (a1*b2-a2*b1)/(a1-a2).
Next step - we need to check does lines segments intersects, not the whole lines. This is done just comparing does x_intersection is within x_min,x_max of both lines, the similar check apply for y_intersection. So line intersection check simplifies to finding linear equations slope and intercept constants. This is one way of self-intersection test of lines. But there is also other way to do it (if we don`t need exact intersection point) - suppose we have line1 with end-points l1p1,l1p2 and line2 with end-points l2p1,l2p2. Define function DIRECTION(point1,point2,point3) which returns the direction [CLOCKWISE,COUNTER_CLOCKWISE,LINE] in which we traverse path point1->point2->point3. By having such function we can easily define when 2 line segments intersects. Below are one example when lines intersects:

Another example when lines do not intersects:

So lines intersects when :
direction(l1p1,l1p2,l2p1) != direction(l1p1,l1p2,l2p2) AND
direction(l2p1,l2p2,l1p1) != direction(l2p1,l2p2,l1p2)

This type of lines intersection test is implemented in polygon self-intersection test. So you can download centroid XNA 3.1 project and try experimenting by drawing various polygons and seeing where centroid is computed.
Final part - some results from this project (centroid is draw as black pixel):
Some convex polygons:

Now, more interesting -> concave polygons:

From the concave polygon examples (especially from first concave) can be seen that centroid is not the average of (max,min) coordinates in shape. However in some cases it can be that centroid coincides with AVERAGE(Min,Max) of figure (for example for rectangle). But in general - it not. If to talk informally, then centroid is simply the average of ALL points coordinates in shape.
Finally just some examples of complex polygons - as evidence that application detects complex (or self-intersecting) polygons :-) :

Actually we can also compute centroids of complex polygon by using info from above mentioned wikipedia site. That is - to get centroid of complex polygon we need:
1. Decompose complex polygon into N simple polygons.
2. Compute area and centroid of each simple polygon by method already explained in this article.
3. Average computed centroids by some formula.
The main task is number 1. - How effectively split complex polygon into simple ones.
This task is left as exercise to the reader :-)

Have fun with centroids, polygons, XNA and 2D/3D programming in general !!

Langton's ant Or order out of chaos

2009-06-21T12:23:00.000-07:00

This time we will talk about Langton's ant - the simple system of simple rules, but with complex behavior. Simple system of Langton`s ant is defined as :
* At a white square, turn 90° left, flip the color of the square, move forward one unit
* At a black square, turn 90° right, flip the color of the square, move forward one unit
This is so called L/R system. But what if for the one color we will have not the one direction, but the direction chain ? That is- for the white squares we will have some chain of 'LRLRR...' events and for the black squares we will have another chain of 'RRLRRL...' events ?? These events should cycle in a row. Then more interesting patterns should occure. So this experiment tests this idea of multiple direction row for the same state/color. Experiment script is here.
Results are somewhat interesting. At first, the classic Langton`s ant example-

Can we have system with shorter steps to "highway" ? Seems - yes we can. This is it:

But this is not the end. There is system which is almost "highway" from the scratch:

And finally- some exotic general Langton`s ant systems:

Conclusion
Seems that systems L/R and L/LRL is unstable and switches to stable system L/RLL.
Another conclusion is that there are some other systems that gets to the different types of "highway" - the order out of chaos.

Have fun with Langton`s ant !!

Detecting copy-move forgery in images

2009-03-15T10:13:00.000-07:00

This time I`ll talk about digital image forensic and how to detect copy-move forgery in images. I`ve implemented some ad-hoc algorithm for detection of this kind of forgeries. Algorithm is robust, so it can detect forgeries in lossy image formats- such as JPEG. Please keep in mind - this algorithm is experimental toy, if you want more general solutions - you should read this paper or maybe this to get some ideas in the field.

Robust detection algorithm steps
1. Blur image for eliminating image details
2. Convert image to degraded palette
3. Decompose image into small NxN pixel blocks
4. Alphabetically order these blocks by their pixel values
5. Extract only these adjacent blocks which have small absolute color difference
6. Cluster these blocks into clusters by intersection area among blocks
7. Extract only these clusters which are bigger than block size
8. Extract only these clusters which have similar cluster, by using some sort of similarity function (in this case Hausdorff distance between clusters)
9. Draw discovered similar clusters on image

Copy-move detection script usage instructions
1. At first try to execute script with default parameters-
python detect_copymove.py image
2. Then try to lower block color deviation threshold-
python detect_copymove.py image --blcoldev=0.05
3. Finally - run script in manual mode and try to spot similar regions by eyes-
python detect_copymove.py image --blcoldev=0.05 --imauto=0
If by trying all 3 steps - no copy-move tamperings revealed - there is a good chance that really there are no such tamperings in image. BTW, script has more parameters, full list of them -
python detect_copymove.py --help

Experiments - object cloning forgery
This type of forgery tries to deceive viewer that there were more objects than really was.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected cloning of one tank, but did not detected cloning of tank in further row. This can be because of small size of that tank and because that tank is on boundary of image.

Experiments - object hiding forgery
In this type of forgery viewer is deceived that there were less objects than in reality was.
Original image:

Doctored image:

Copy-move script output:

Experiments - mixed type forgery
In this type of forgery some objects was hidden, some - cloned.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected 2 forgeries - cloning of dog, and cloning of background (which was done with intention to hide other dog)

Conclusions
Copy-move forgery detection script works not bad. However it is also not perfect of course (may not detect all tamperings, so called false negative rate is not zero of course). And I`ve not tested very well, but I suppose script also have false positive rate - i.e. can detect areas which actually was not tampered. But all in all - I think this algorithm is pretty obvious and can be simply implemented and used as first line of defense against copy-move forgeries.

Have fun in copy-move forgeries hunting !!

Evolving first lady of the internet

2009-02-21T10:32:00.000-08:00

At last I ported selfish gene algorithm into Python. This algorithm pseudo-code is:
1. Encode some solution into set of genes.
2. Set default probability for each gene.
3. Generate 2 random solutions, according to genes probability.
4. Compare these 2 solutions:
for better solution- increase it`s genes probability by some value
for worse solution- decrease it`s genes probability by some value
5. Repeat everything from 3, until needed.

Now, interesting part. What can we do with selfish gene algorithm ? After seeing this post about genetic programming and evolutionary art, I`ve decided to try something similar. But only by using selfish gene algorithm. I`ve tried 2 experiments - tried to evolve picture of first lady of the internet composed of polygons. And second experiment - lenna picture is evolved as some number of lines.
So experiment idea is to generate 2 random images, composed of random polygons (or lines in other experiment) and to compare these 2 images. For image which is more similar to original lenna picture - we increase polygons probability, for other picture - decrease polygons probability. In the long run - "good polygons" tends to group together.
Below are the results of these experiments. Evolved pictures of lenna are compiled as frames of animated GIF image. N - is the iteration number (starting from zero):

Lena evolved as set of polygons

Original Lena	Lena evolution: 39614 iterations 100 polygons 3.5 hours experiment code

Lena evolved as set of lines

Original Lena	Lena evolution: 69471 iterations 200 lines 2.5 hours experiment code

Conclusions:
Selfish gene algorithm (sort of evolution strategy algorithm) is suitable for solving search and optimization problems, including generation of evolutionary art :-)
Below are final iteration pictures better quality than gif:

Have fun with selfish gene algorithm, genetic algorithms and evolution art !

Project Euler p83 and Dijkstra's algorithm

2008-12-03T04:28:00.000-08:00

Problem definition:
-----------------------------------------------------------
In the 5 by 5 matrix below, the minimal path sum from the top left to the bottom right, by moving left, right, up, and down, is indicated in bold and is equal to 2297.

131 673 234 103 18
201 96 342 965 150
630 803 746 422 111
537 699 497 121 956
805 732 524 37 331

Find the minimal path sum, in matrix.txt , a 31K text file containing a 80 by 80 matrix, from the top left to the bottom right by moving left, right, up, and down.

-----------------------------------------------------------
Problem is solved with the help of Dijkstra's algorithm
Answer is calculated in a 2.5 seconds. Also it is not optimized,- computes distances to all nodes in network. To speed up - when we find target node - we can only check these nodes which distance is less than target node distance. Also this code suits equally well for solving project euler problem no 81 - only change shortestpath() function parameter "directions". And finally - this solution shows shortest path visually,- prints it to screen. Program code:
___________________________________________________________

import time as t

def readmatrix(filename):
    f=open(filename, 'r')
    m=[map(lambda x: int(x),r.strip().split(",")) for r in f.readlines()]
    return m

def nextcell(dctresults):
    minw = 10000000
    mink = None
    for k in dctresults:
        if not dctresults[k][2]:
            if dctresults[k][1] < minw:
                minw = dctresults[k][1]
                mink = k
    return mink

def analyzecell(lstmat, tupcell, lstdirec, dctresults):
    i,j = tupcell
    for dr in lstdirec:
        di,dj = dr
        if (i+di) in range(len(lstmat)) and (j+dj) in range(len(lstmat)):
            if not dctresults.has_key((i+di,j+dj)):
                dctresults[(i+di,j+dj)] = [(i,j),dctresults[(i,j)][1]+lstmat[i+di][j+dj],False]
            else:
                oldw = dctresults[(i+di,j+dj)][1]
                neww = dctresults[(i,j)][1]+lstmat[i+di][j+dj]
                if neww < oldw:
                    dctresults[(i+di,j+dj)][0] = (i,j)
                    dctresults[(i+di,j+dj)][1] = neww
    dctresults[(i,j)][2] = True

def shortestpath(filename, directions):
    mat = readmatrix(filename)
    results = {}
    results[(0,0)] = [None,mat[0][0],False]
    cell = nextcell(results)
    while cell:
        analyzecell(mat,cell,directions,results)
        cell = nextcell(results)
    return results, results[(len(mat)-1,len(mat)-1)][1], mat

def showpath(res, inpmat):
    print ' '+'-'*len(inpmat)+' '
    pathcells = []
    curc = (len(inpmat)-1,len(inpmat)-1)
    while curc:
        pathcells += [curc]
        curc = res[curc][0]
    for i in range(len(inpmat)):
        s='|'
        for j in range(len(inpmat)):
            s+= '#' if (i,j) in pathcells else '.'
        print s+'|'
    print ' '+'-'*len(inpmat)+' '

t1 = t.time()*1000
resd, pathsum, mat = shortestpath('matrix.txt',[(-1,0),(+1,0),(0,-1),(0,+1)])
t2 = t.time()*1000
print 'Shortest path sum', pathsum, '; Calculation time is ', int(t2-t1), 'ms'
showpath(resd, mat)

___________________________________________________________
Have fun with shortest path algorithms !!!

Project Euler p92

2008-11-11T09:30:00.000-08:00

Problem description:
-------------------------------------------------------------------
A number chain is created by continuously adding the square of the digits in a number to form a new number until it has been seen before.

For example,

44 -> 32 -> 13 -> 10 -> 1 -> 1
85 -> 89 -> 145 -> 42 -> 20 -> 4 -> 16 -> 37 -> 58 -> 89

Therefore any chain that arrives at 1 or 89 will become stuck in an endless loop. What is most amazing is that EVERY starting number will eventually arrive at 1 or 89.

How many starting numbers below ten million will arrive at 89?
-------------------------------------------------------------------
Semi-bruteforce algorithm based on idea that we don`t need to construct full number chains from starting number until 1 or 89. Instead we will save in memory each starting number final result (1 or 89). And at the next starting number we will construct number chain only until previous calculated number.
Well solution is slow, but fits in "one minute rule" :-). Here it is-
_____________________________________________________________

import time as t

def nextnum(number):
 n = number
 s = 0
 while n!=0:
  m = n%10
  n = n/10
  s += m*m
 return s

def untilknown(number):
    start = number
    res = []
    while start >= number and start not in (1,89):
        start = nextnum(start)
        res += [start]
    return res

def endswith89(limit):
    seq = {1:1,89:89}
    ret = 0
    for x in range(1,limit):
        if not seq.has_key(x):
            k = untilknown(x)
            if k:
                kk = k[-1]
            else:
                kk = x
            seq[x] = seq[kk]
            if len(k) > 1:
                for i in range(len(k)-1):
                    seq[k[i]] = seq[kk]
        if seq[x] == 89:
            ret += 1
    return ret

t1 = t.time()
ans = endswith89(10000000)
t2 = t.time()

print ans, int(t2-t1),'s'

_____________________________________________________________

Have fun !!

Project Euler p99

2008-10-27T04:12:00.000-07:00

Project Euler problem 99 description:
------------------------------------------------------------------------
Comparing two numbers written in index form like 2^11 and 3^7 is not difficult, as any calculator would confirm that 2^11 = 2048 < 3^7 = 2187.

However, confirming that 632382^518061 > 519432^525806 would be much more difficult, as both numbers contain over three million digits.

Using base_exp.txt (right click and 'Save Link/Target As...'), a 22K text file containing one thousand lines with a base/exponent pair on each line, determine which line number has the greatest numerical value.

NOTE: The first two lines in the file represent the numbers in the example given above.
------------------------------------------------------------------------
The trick of solution is not to compare whole numbers (because raising to power will be too much CPU time and memory consuming operation), instead we can compare logarithms of numbers. So actually we will compare LOG(number^power). Ok, but we still have to do power operation ? Well, no more ! By applying known logarithm formula: LOG(number^power) = power*LOG(number). So we are comparing power*LOG(number). This is very fast operation.
Problem solution is found in a 3.7 ms. Solution Python code:
_____________________________________________________

import time as tm
import math as m

t1 = tm.time()*1000
f=open('base_exp.txt', 'r')
l=[i.strip().split(",")+[n+1] for (n,i) in enumerate(f.readlines())]
l = map(lambda i: (int(i[0]),int(i[1]), i[2]), l)
l = map(lambda (x,a,n): (a*m.log(x), n), l)
max = max(l)
t2 = tm.time()*1000

print max,t2-t1,'ms'

________________________________________________________
Have fun!

Calculation of Pi

2008-10-11T10:29:00.000-07:00

This time I tried to answer a question - what is the best method of calculating Pi to a required accuracy - for example,- 100 or 1000 or 10 000 decimal digits ?
So I chose 4 Pi formulas for testing:
1. Leibniz
2. Madhava transformed series
3. Machin like formula
4. Gauss algorithm
The idea for experiment - is to find out which algorithm gives more decimal digits of Pi in same number of iterations/terms. So measurement algorithm is this:
1. Calculate correct digits of Pi given from every Pi formula in each iteration.
2. Plot diagram - correct Pi decimal digits amount VS iteration number
(For this to work - you must have Matplotlib module).
Python code doing this stuff:
_____________________________________________________

import inspect, time, sys
from decimal import Decimal, getcontext
from pylab import plot,  show,  xlabel,  ylabel,  title,  grid,  legend

def correct_digits(calcpi):
 pi = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
 gd = 1000
 for i,s in enumerate(zip(pi,str(calcpi))):
  if s[0] != s[1]:
   gd = i - 2
   if gd < 0: gd = 0
   break
 return gd

def show_results(res):
 x=range(len(res[0][2]))
 plpar = ""
 legpar = ""
 for i in range(len(res)):
  if plpar: plpar +=","
  if legpar: legpar +=","
  plpar += "x, res["+str(i)+"][2]"
  legpar += "res["+str(i)+"][0][3:] +' ('+ str(int(res["+str(i)+"][1])) + 'ms)'"
 plpar = "plot("+plpar+", linewidth=2.0)"
 legpar = "legend(("+legpar+"), loc=0)"
 exec(plpar)
 exec(legpar)
 xlabel('Iterations', fontsize = 15)
 ylabel('Correct digits', fontsize = 15)
 title("Convergence of some $\pi$ formulas",fontsize=20)
 grid()
 show()

def test_pi_formulas(pi_func_list):
    res = []
    terms = 6
    print 'Calculating %d iterations...' % terms
    getcontext().prec = 1002
    for pif in pi_func_list:
        t1 = time.time()*1000
        name,r = pif(terms)
        t2 = time.time()*1000
        res.append((name,t2-t1,r))
    show_results(res)

def pi_Leibniz(terms):
 """
 http://en.wikipedia.org/wiki/Leibniz_formula_for_pi
 """
 st = inspect.stack()[0][3]
 ra = Decimal(0)
 ans = []
 
 for n in range(terms):
  ra += Decimal(4*(-1)**n)/Decimal(2*n+1)
  ans.append(correct_digits(ra))
 
 return st,ans

def pi_Madhava(terms):
 """
 http://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#The_value_of_.CF.80_.28pi.29
 """
 st = inspect.stack()[0][3]
 ra = Decimal(0)
 k = Decimal(12)**Decimal('0.5')
 ans = []
 
 for n in range(terms):
  ra += Decimal((-1)**n)/Decimal((2*n+1)*(3**n))
  ans.append(correct_digits(k*ra))

 return st, ans

def pi_Machin(terms):
 """
 http://en.wikipedia.org/wiki/Machin-like_formula
 """
 st = inspect.stack()[0][3]
 c1 = Decimal(1)/Decimal(5)
 c2 = Decimal(1)/Decimal(239)
 a = Decimal(4)
 atan1 = Decimal(0)
 atan2 = Decimal(0)
 ra = Decimal(0)
 ans = []
 
 for n in range(terms):
  atan1 += Decimal((-1)**n) * (c1**Decimal(2*n+1))/Decimal(2*n+1)
  atan2 += Decimal((-1)**n) * (c2**Decimal(2*n+1))/Decimal(2*n+1)
  ra = a * (a*atan1 - atan2)
  ans.append(correct_digits(ra))
 
 return st, ans

def pi_Gauss(terms):
 """
 http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
 """
 st = inspect.stack()[0][3]
 one = Decimal(1)
 ra = Decimal(0)
 val = (one,one/(Decimal(2)**Decimal('0.5')),one/Decimal(4),one)
 ans = []
 
 for n in range(terms):
  an, bn, tn, pn = val
  an1 = (an+bn)/Decimal(2)
  bn1 = (an*bn)**Decimal('0.5')
  tn1 = tn - pn * (an - an1)**Decimal(2)
  pn1 = Decimal(2) * pn
  val = (an1, bn1, tn1, pn1)
  ra = ((an+bn)**Decimal(2))/(Decimal(4)*tn)
  ans.append(correct_digits(ra))
 
 return st, ans

if __name__ == "__main__":
 mod = sys.modules[__name__]
 pifunc = []
 for key, value in mod.__dict__.items():
  if str(value)[:13] == '<function pi_':
   pifunc.append(value)
 test_pi_formulas(pifunc)

_____________________________________________________
BTW this code is very scalable - if there is a need to test additional Pi formula - you only need to add new Pi function code like this:
"def pi_new_required_function(terms):
[...]". (Function name MUST begin with "pi_"). The rest program will do it for you - measure execution times, measure decimal digits correctness, plot results...
This is thanks to Python reflection used in program -> program scans itself for ALL pi functions and dynamically invokes them.

So after execution of this test we get such graph:

Conclusion
From graph can be seen that fastest converges Gauss Pi algorithm (from the chosen set of Pi formulas). So if you need for some reason to calculate Pi in desired accuracy (desired decimal places amount) I suggest Gauss algorithm. Because you will get result in less iterations / CPU cycles compared with other algorithms from the set.

Have fun!

Detecting similar color zones in image

2008-09-21T12:33:00.000-07:00

After analyzing information entropy in image I started to think - what will be if we hash image not by it`s total entropy, but instead- by some other cumulative function ? After several tries, I found this hashing function interesting-

Where Ir,Ig,Ib - pixel value in red,green,blue channels respectively and N is number of pixels in image.
Python code doing this stuff of image zones comparision is very similar to the previous mentioned (information entropy) article code. It only differs on 2 aspects:
- Now we compare image blocks by above mentioned hash function
- There is no code of grouping image blocks.
Pseudo code of algorithm:
1. Split given image into 4x4 small image blocks.
2. Calculate color hash of each block.
3. Find blocks with similar hash value
4. Draw these blocks on top of original image as green color layer
Python code, which does that:
_______________________________________________________

from PIL import Image
from math import *

def colorhash(pixels):
 """
 Calculating color hash by summing
 intensity levels in all 3 channels
 separatelly and multiplying these
 accumulated intensities together
 also dividing by cube of pixel amount.
 """
 
 cr = sum([r for (r,g,b) in pixels])
 cg = sum([g for (r,g,b) in pixels])
 cb = sum([b for (r,g,b) in pixels])
 
 return ((cr+1) * (cg+1) * (cb+1)) / len(pixels)**3

def decompose(image, block_len):
 """
 Decomposing given image into some number of
 smaller images of size block_len*block_len
 """
 parts = []
 w, h = image.size
 
 for x in range(0, w, block_len):
  for y in range(0, h, block_len):
   locim = image.crop((x,y,x+block_len,y+block_len))
   acc = colorhash(list(locim.getdata()))
   parts.append((acc,x,y,locim))

 parts.sort()
 
 return parts

def similarparts(imagparts):
 """
 Detecting similar image blocks by comparing
 calculated color hashes of given images. 
 Two images considered being equal if 
 their color hash difference is not big.
 """
 dupl = []
 
 for i in range(len(imagparts)-1):
  acc1, x1, y1, im1 = imagparts[i]
   acc2, x2, y2, im2 = imagparts[i+1]
  
   if acc1 == acc2 == 0:
     diff = 0.0
   else:
     diff = 100.0 * (1.0 - float(acc1) / acc2)
  
   col = 0.3
 
  if diff > col:
    if imagparts[i] not in dupl:
     dupl.append(imagparts[i])
     if imagparts[i+1] not in dupl:
      dupl.append(imagparts[i+1])

 return dupl

def marksimilar(image, dparts):
 """
 Mark found similar image blocks on
 original image, by applying red layer
 on similar parts of image.
 """
 if dparts:
  colormask = Image.new('RGB', dparts[0][3].size,(0,255,0))
  for (acc,x,y,im) in dparts:
   im = Image.blend(im, colormask, 0.3)
   image.paste(im,(x,y))
 
 return image

if __name__ == '__main__':
 im = Image.open("1.jpg")
 ls = decompose(im, 4)
 dparts = similarparts(ls)
 im = marksimilar(im, dparts)
 im.show()

_________________________________________________________
Running this code on several images we get:

Conclusion
Seems in some cases by using such color hash function we can detect color similarities in image. Maybe even some symmetry in image (if it exists). Also this hash function should be improved, because somehow it "sticks" into finding similarities in darker zones of image.

Have fun!

Detecting similar entropy zones in image

2008-09-14T11:24:00.000-07:00

As I thought about information entropy one idea came to me - to write application which looks similar entropy zones in image. So after some time, I came with this algorithm (pseudo code):

1. Split image into 5x5 pixel image blocks
2. Calculate information entropy of these blocks (actually sum of entropy in 3 color channels)
3. Find similar entropy blocks.
4. [...] filter out small groups of blocks (seems like noise, huh ? ). Blah, blah...
5. Picture these similar entropy blocks on top of original image as red color layer.

Now real part (as you know already) - Python code which does the job [you need PIL module to run this]:
______________________________________________________________

from PIL import Image
from math import *

def entropysum(pixels):
 """
 Calculating information entropy for image
 region and returning entropy sum for
 all 3 color channels
 """
 cr = [r for (r,g,b) in pixels]
 cg = [g for (r,g,b) in pixels]
 cb = [b for (r,g,b) in pixels]
 
 er = 0.0
 eg = 0.0
 eb = 0.0
 
 for r in set(cr):
  p = float(cr.count(r))/len(cr)
  if p > 0.0: er += -p * log(p,2)

 for g in set(cg):
  p = float(cg.count(g))/len(cg)
  if p > 0.0: eg += -p * log(p,2)

 for b in set(cb):
  p = float(cb.count(b))/len(cb)
  if p > 0.0: eb += -p * log(p,2)

 return er + eg + eb

def decompose(image, block_len):
 """
 Decomposing given image into some number of
 smaller images of size block_len*block_len
 """
 parts = []
 w, h = image.size
 
 for x in range(0, w, block_len):
  for y in range(0, h, block_len):
   locim = image.crop((x,y,x+block_len,y+block_len))
   acc = entropysum(list(locim.getdata()))
   parts.append((acc,x,y,locim))
 
 parts.sort()
 
 return parts

def similarparts(imagparts):
 """
 Detecting similar image blocks by comparing
 entropy of given images. Two images considered
 being equal if entropy difference is not big.
 """
 dupl = []
 
 for i in range(len(imagparts)-1):
  acc1, x1, y1, im1 = imagparts[i]
  acc2, x2, y2, im2 = imagparts[i+1]
  
  if acc1 == acc2 == 0:
   gain = 0.0
  else:
   gain = 100.0 * (1.0 - acc1 / acc2)

  if 0.01 < gain < 0.1 :
   if imagparts[i] not in dupl:
    dupl.append(imagparts[i])
   if imagparts[i+1] not in dupl:
    dupl.append(imagparts[i+1])

 return dupl

def clusterparts(parts):
 """
 Grouping nearest images into groups.
 This is done, because we need to
 filter out very small groups. We
 want to know only big differences.
 """
 
 filtparts = []
 clust = {}
 belongs = {}
 w,h = parts[0][3].size
 
 # assign all parts to clusters
 for i in range(len(parts)):
  acc, x, y, im = parts[i]
  sides = []
  sides.append(str(x)+str(y)+str(x+w)+str(y))
  sides.append(str(x+w)+str(y)+str(x+w)+str(y+h))
  sides.append(str(x)+str(y+h)+str(x+w)+str(y+h))
  sides.append(str(x)+str(y)+str(x)+str(y+h))
  
  # detect side already in cluster
  fc = None
  for s in sides:
   if belongs.has_key(s):
    fc = belongs[s]
    break
  
  # if this is new cluster
  if fc == None:
   fc = len(clust) + 1
   clust[fc] = 1
  else:
   clust[fc] += 1
  
  # set cluster for rectangle sides
  for s in sides:
   if not belongs.has_key(s):
    belongs[s] = fc

 # filter out small clusters
 for i in range(len(parts)):
  acc, x, y, im = parts[i]
  side = str(x)+str(y)+str(x+w)+str(y)
  cl = belongs[side]
  if clust[cl] > 2:
   filtparts.append(parts[i])
 
 return filtparts

def marksimilar(image, dparts):
 """
 Mark found similar image blocks on
 original image, by applying red layer
 on similar parts of image.
 """
 if dparts:
  colormask = Image.new('RGB', dparts[0][3].size,(255,0,0))
  for (acc,x,y,im) in dparts:
   im = Image.blend(im, colormask, 0.4)
   image.paste(im,(x,y))
 
 return image

if __name__ == '__main__':
 im = Image.open("1.jpg")
 ls = decompose(im, 5)
 dparts = similarparts(ls)
 cparts = clusterparts(dparts)
 im = marksimilar(im, cparts)
 im.show()

________________________________________________________
So these are the results after running this script on several images:

Conclusion

So this algorithm is an interesting tool for exploration of information entropy in image. Maybe in some cases it could be a tool for analyzing very similar texture zones. BTW information entropy may be used for hashing image. Hashing image is useful, because it lets us to search similar images in database (for example) by its hash.

Have fun !

Number of runs in random sequence

2008-08-31T12:32:00.000-07:00

While reading this article about human ability on outputing random data-
http://faculty.rhodes.edu/wetzel/random/mainbody.html
I started to think - How we can simply detect that entered sequence is NOT random ( or entered by human) ? Of course to make some real-world measurements we need statistics a lot. But if we need just some straightfoward method for preliminary results ? Well, the idea was to check the number of runs in sequence and to compare the result with the expected value. (Number of runs is total number of subsequences in sequence).
So that if number_of_runs != expected_number_of_runs then we can conclude that data is not much random. OK. Number of runs we can measure simply. But what about expected_number_of_runs ? Somehow I suspected that this expected runs count should be dependent on the random variable values amount. For example sequences 'azazaz' and 'azxazx' are the same length, but one is composed of a,z (2 values) and second- a,z,x (3 values). So if we know what is going on with the number_of_runs when random sequences is composed from different values amount - we may calculate expected_number_of_runs given any sequence. Second point - for eliminating need to recalculate everything when sequence length changes - better calculate expected_number_of_runs / sequence_length.
So I made such experiment for detection of such dependancy.
Experiment Python code (to run it you need matplotlib python module):
_________________________________________________________________

from random import *
from pylab import plot,  show,  xlabel,  ylabel,  title,  grid,  legend

def runs(lst):
    return sum([int(not lst[i]==lst[i+1]) for i in range(len(lst)-1)])+1

def testruns(n):
    s=[randint(1, n) for x in range(10000)]
    return float(runs(s)) / len(s)

max = 20
xrange = range(2, max)
yruns = [testruns(x) for x in xrange] # Generating random data and counting runs
yreg = [1.0-1.0/x for x in xrange]  # Function for fitting data
plot(xrange, yruns ,  'ro',  xrange, yreg)
xlabel('Variable values amount')
ylabel('Number of runs per item')
title('Number of runs dependance on variable values amount')
legend(('Experimental' , 'Data fit:  ' r'$1-n^{-1}$'), loc=7)
grid()
show()

_________________________________________________________________

So after running this random data test- we get such plot:

So it tells us that expected_number_of_runs we need to calculate like that-
expected_number_of_runs = sequence_length * (1 - 1/n);
n- unique element count in sequence. So now we have some tool for detecting sequence being not much random.
First- we measure real number_of_runs in sequence.
Second - we calculate expected_number_of_runs.
Third - if these two values is not approximately equal - then we can conclude that data is not random. Of course we can`t be sure about randomness just from this test.
We need a lot of different measures - information entropy and as such... But
the goal of this article was to play with number_of_runs. Maybe later - I will write something about information entropy. It is very interesting also.

Have fun!!!

Codegolf - Vigenere Cipher

2008-08-27T12:34:00.000-07:00

Problem definition:
---------------------------------------------------
The Vigenere cipher is a simple form of a polyalphabetic substitution, using a series of different Caesar ciphers based on the letters of a keyword. It's over 450 years old now, but that shouldn't stop us having a bit of golfing fun with it now.

The following information is largely taken from Wikipedia's Vigenere cipher page. You might find some interesting information that isn't included here, so I recommend you check it out.

Your job is to write some code which takes a keyword and some plaintext, and encrypts it according to the Vigenere cipher.

In a Caesar cipher, each letter of the alphabet is shifted along some number of places; for example, in a Caesar cipher of shift 3, A would become D, B would become E and so on. The Vigenere cipher consists of several Caesar ciphers in sequence with different shift values.

To encipher, a table of alphabets can be used, termed a tabula recta, Vigenere square, or Vigenere table. It consists of the alphabet written out 26 times in different rows, each alphabet shifted cyclically to the left compared to the previous alphabet, corresponding to the 26 possible Caesar ciphers. At different points in the encryption process, the cipher uses a different alphabet from one of the rows. The alphabet used at each point depends on a repeating keyword.

Click here to see the Vigenere square.

See the examples section below for a worked example.
More Information

* Your code will be given both the keyword and the plaintext on stdin. It will be in the format

KEYWORD
PLAINTEXT

with newlines at the end of each line.

* You should print the encrypted plaintext on stdout.
* Both the keyword and the plaintext will contain only the uppercase letters A through Z.
* Your code will be run three times and will need to pass all three tests to be deemed successful.
---------------------------------------------------
Code golf challenge is interesting because it measures your ability to write code as short as possible. My try to give shortest solution in Python:

Run no 1 -> code size 111 bytes:

k=raw_input()*100;p=raw_input()
print ''.join(map(lambda x, y: chr(65+(ord(x)+ord(y)-130)%26), k[:len(p)], p))

Run no 2 -> code size 84 bytes:

i=raw_input;print ''.join(chr(65+(ord(x)+ord(y)-130)%26) for x,y in zip(i()*10,i()))

Have fun!
Code golf

Project Euler p22

2008-08-25T07:20:00.001-07:00

Problem description:
-----------------------------------------------------------
Using names.txt (right click and 'Save Link/Target As...'), a 46K text file containing over five-thousand first names, begin by sorting it into alphabetical order. Then working out the alphabetical value for each name, multiply this value by its alphabetical position in the list to obtain a name score.

For example, when the list is sorted into alphabetical order, COLIN, which is worth 3 + 15 + 12 + 9 + 14 = 53, is the 938th name in the list. So, COLIN would obtain a score of 938 * 53 = 49714.

What is the total of all the name scores in the file?
-----------------------------------------------------------
Solution in Python:

def ReadFile():
    """ Read data input file"""
    f = open('C:/ProjectEuler/names.txt','r')
    str = f.readlines()
    f.close()
    str[0] = str[0].replace('"','')
    ret = str[0].split(',')
    ret.sort()
    return ret

def AlphabeticalValue(string):
    alph = list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')
    values = map(lambda x: alph.index(x)+1, list(string))
    return sum(values)

def TotalScore():
    names = ReadFile()
    namsc = map(lambda (x,y): AlphabeticalValue(y)*(x+1), list(enumerate(names)))
    return sum(namsc)

if __name__ == '__main__':
    print TotalScore()

What is missing ? Of course - Have fun :-)

Google codejam 2008 contest practice test

2008-07-11T12:56:00.000-07:00

This time we will probe Google famous CodeJam contest. We will try to solve codejam 2008 one of practice tests. Test problem is defined as:
-----------------------------------------------------------------------
Problem

You're interested in writing a program to classify triangles. Triangles can be classified according to their internal angles. If one of the internal angles is exactly 90 degrees, then that triangle is known as a "right" triangle. If one of the internal angles is greater than 90 degrees, that triangle is known as an "obtuse" triangle. Otherwise, all the internal angles are less than 90 degrees and the triangle is known as an "acute" triangle.

Triangles can also be classified according to the relative lengths of their sides. In a "scalene" triangle, all three sides have different lengths. In an "isosceles" triangle, two of the sides are of equal length. (If all three sides have the same length, the triangle is known as an "equilateral" triangle, but you can ignore this case since there will be no equilateral triangles in the input data.)

Your program must determine, for each set of three points, whether or not those points form a triangle. If the three points are not distinct, or the three points are collinear, then those points do not form a valid triangle. (Another way is to calculate the area of the triangle; valid triangles must have non-zero area.) Otherwise, your program will classify the triangle as one of "acute", "obtuse", or "right", and one of "isosceles" or "scalene".

Input

The first line of input gives the number of cases, N. N test cases follow. Each case is a line formatted as

x1 y1 x2 y2 x3 y3

Output

For each test case, output one line containing "Case #x: " followed by one of these strings:

* isosceles acute triangle
* isosceles obtuse triangle
* isosceles right triangle
* scalene acute triangle
* scalene obtuse triangle
* scalene right triangle
* not a triangle

Limits

1 ≤ N ≤ 100,
x1, y1, x2, y2, x3, y3 will be integers.

Sample Input

0 0 0 4 1 2
1 1 1 4 3 2
2 2 2 4 4 3
3 3 3 4 5 3
4 4 4 5 5 6
5 5 5 6 6 5
6 6 6 7 6 8
7 7 7 7 7 7

Sample Output

Case #1: isosceles obtuse triangle
Case #2: scalene acute triangle
Case #3: isosceles acute triangle
Case #4: scalene right triangle
Case #5: scalene obtuse triangle
Case #6: isosceles right triangle
Case #7: not a triangle
Case #8: not a triangle

-----------------------------------------------------------------------
This time solution is written in Python language. To test this solution you need to have Python interpreter installed. You can install it from www.python.org. For trying this script- paste Sample Input lines into file named 'A-large.in' and place this file into the same folder where below mentioned Python script will be. Script is somewhat fast - 100 sample triangles is analyzed in just 2.2 ms !!!. Script code:

#!/usr/bin/env python

import sys, time

def ReadFile():
    """ Read data input file"""
    path = sys.path[0]
    f = open(path + '/A-large.in','r')
    str = f.readlines()
    f.close()
    return str

def IsTriangle(p):
    """Defines triangle by it`s area"""
    area = (p[0]*(p[3]-p[5])+p[2]*(p[5]-p[1])+p[4]*(p[1]-p[3]))
    if area == 0:
        return 'not a triangle'
    else:
        return 'triangle'

def SideType(p):
    """Defines triangle by sides lengths"""
    s1 = ((abs(p[0]-p[2])**2)+(abs(p[1]-p[3])**2))**0.5
    s2 = ((abs(p[0]-p[4])**2)+(abs(p[1]-p[5])**2))**0.5
    s3 = ((abs(p[2]-p[4])**2)+(abs(p[3]-p[5])**2))**0.5
    if s1!=s2 and s1!=s3 and s3!=s2:
        return 'scalene'
    if s1==s2 or s1==s3 or s3==s2:
        return 'isosceles'

def AngleType(p):
    """Defines triangle type by angle.
       Idea based on Pythagorean theorem and somewhat
       empirical hack about 'not perfect' right triangles ->
       for which (a^2+b^2)/c^2 is only approximately 1.
       And the more angle shifts from 90 degrees, the
       more (a^2+b^2)/c^2 relation shifts from 1"""
    s1 = ((abs(p[0]-p[2])**2)+(abs(p[1]-p[3])**2))**0.5
    s2 = ((abs(p[0]-p[4])**2)+(abs(p[1]-p[5])**2))**0.5
    s3 = ((abs(p[2]-p[4])**2)+(abs(p[3]-p[5])**2))**0.5
    k1 = round(float((s1**2)+(s2**2))/float(s3**2),3)
    k2 = round(float((s1**2)+(s3**2))/float(s2**2),3)
    k3 = round(float((s3**2)+(s2**2))/float(s1**2),3)
    
    if k1==1.0 or k2==1.0 or k3==1.0:
        return 'right'
    elif k1<1.0 or k2<1.0 or k3<1.0:
        return 'obtuse'
    else:
        return 'acute'

def ComputeResults():
    out = []
    for i,s in enumerate(ReadFile()):
        pt = map(int,s.rstrip().split())
        if pt.__len__() != 1:
            tr = IsTriangle(pt)
            st = ''
            at = ''
            if tr == 'triangle':
                st = SideType(pt) + ' '
                at = AngleType(pt) + ' '
            tot = 'Case #' + str(i) + ': ' + st + at + tr + '\n'
            out.append(tot)
    return out

def WriteFile(lst):
    """ Write results to data file"""
    path = sys.path[0]
    f = open(path + '/A-large.out','w')
    f.writelines(lst)
    f.close()

if __name__ == '__main__':
    t1 = time.time()*1000
    res = ComputeResults()
    t2 = time.time()*1000
    WriteFile(res)
    print 'Done.'
    print 'Computation time ',round(t2-t1,3),'ms'
    print 'Results are written in file A-large.out'

Have fun !

Project Euler p96 - solving sudoku

2008-06-10T09:15:00.000-07:00

This is interesting one. Problem description:
-----------------------------------------------
Su Doku (Japanese meaning number place) is the name given to a popular puzzle concept. Its origin is unclear, but credit must be attributed to Leonhard Euler who invented a similar, and much more difficult, puzzle idea called Latin Squares. The objective of Su Doku puzzles, however, is to replace the blanks (or zeros) in a 9 by 9 grid in such that each row, column, and 3 by 3 box contains each of the digits 1 to 9. Below is an example of a typical starting puzzle grid and its solution grid.


0 0 3  0 2 0  6 0 0
9 0 0  3 0 5  0 0 1
0 0 1  8 0 6  4 0 0

0 0 8  1 0 2  9 0 0
7 0 0  0 0 0  0 0 8
0 0 6  7 0 8  2 0 0

0 0 2  6 0 9  5 0 0
8 0 0  2 0 3  0 0 9
0 0 5  0 1 0  3 0 0
===================
4 8 3  9 2 1  6 5 7
9 6 7  3 4 5  8 2 1
2 5 1  8 7 6  4 9 3

5 4 8  1 3 2  9 7 6
7 2 9  5 6 4  1 3 8
1 3 6  7 9 8  2 4 5

3 7 2  6 8 9  5 1 4
8 1 4  2 5 3  7 6 9
6 9 5  4 1 7  3 8 2

A well constructed Su Doku puzzle has a unique solution and can be solved by logic, although it may be necessary to employ "guess and test" methods in order to eliminate options (there is much contested opinion over this). The complexity of the search determines the difficulty of the puzzle; the example above is considered easy because it can be solved by straight forward direct deduction.

The 6K text file, sudoku.txt (right click and 'Save Link/Target As...'), contains fifty different Su Doku puzzles ranging in difficulty, but all with unique solutions (the first puzzle in the file is the example above).

By solving all fifty puzzles find the sum of the 3-digit numbers found in the top left corner of each solution grid; for example, 483 is the 3-digit number found in the top left corner of the solution grid above.
-----------------------------------------------

Well, back to more understandable language like C# :-). Solution is done as semi brute-force, iterative method. Solution is quick enough - fifty puzzles are solved in just 122 ms. So this is it:


class Program
    {

        static int[][][] ReadSudoku(string file)
        {
            StreamReader sr = new StreamReader(file);
            int[][][] res = new int[50][][];
            int row = 0;
            int sud = -1;
            string line;

            while (!sr.EndOfStream)
            {
                line = sr.ReadLine();
                if (line.ToLower().IndexOf("grid") > -1)
                {
                    row = 0;
                    sud++;
                    res[sud] = new int[9][];
                }
                else
                {
                    res[sud][row] = Array.ConvertAll(line.ToCharArray(), c => 
                                    { return int.Parse(c.ToString()); });
                    row++;
                }
            }
            sr.Close();

            return res;
        }

        static void ShowSudoku(int[][] sud)
        {
            string o;
            o = sud.Aggregate("", (agg, next) =>
            {
                return agg +
                    ((agg == "") ? "":"\n") +
                    next.Aggregate("|", (agg2, next2) => { 
                        return agg2 + next2.ToString(); }) + "|";
            });
            o = "-----------\n" + o + "\n-----------";
            o = o.Replace("0", ".");
            Console.Write(o);
        }

        static bool ValueIsValid(int[][] grid, int val, int row, int col)
        {
            int regc, regr;
            int lc, lr;

            if (val == 0)
                return false;

            // checking row and column
            for (int i = 0; i < 9; i++)
                if ((grid[row][i] == val && i != col) || 
                    (grid[i][col] == val && i != row))
                    return false;
            
            // checking 3x3 region
            regr = 3*(row / 3);
            regc = 3*(col / 3);
            for (int i = 0; i < 3; i++)
            {
                lc = regc + i;
                for (int j = 0; j < 3; j++)
                {
                    lr = regr + j;
                    if (lc != col || lr != row)
                        if (grid[lr][lc] == val)
                            return false;
                }
            }

            return true;
        }

        static void CellWithMinValid(int[][] grid, 
                                    out int rowo, 
                                    out int colo)
        {
            int row = -1, col = -1;
            int minvalid = 9;
            int valcount;

            for (int i = 0; i < grid.Length; i++)
            {
                for (int j = 0; j < grid[i].Length; j++)
                {
                    if (grid[i][j] == 0)
                    {
                        valcount = 0;
                        for (int k = 1; k < 10; k++)
                            if (ValueIsValid(grid, k, i, j))
                                valcount++;
                        if (valcount < minvalid)
                        {
                            minvalid = valcount;
                            row = i;
                            col = j;
                        }
                    }
                    if (minvalid == 1) break;
                }
                if (minvalid == 1) break;
            }

            rowo = row;
            colo = col;
        }

        static void SolveSudoku(int[][] grid)
        {
            int row, col;
            int curval;
            int iter;
            bool force;
            int[,] steps = new int[81, 2];

            CellWithMinValid(grid, out row, out col);
            grid[row][col] = 1;
            steps[0, 0] = row;
            steps[0, 1] = col;
            force = false;
            iter = 0;

            while (row > -1)
            {
                curval = grid[steps[iter, 0]]
                             [steps[iter, 1]];
                if (ValueIsValid(grid, curval, 
                                steps[iter, 0], 
                                steps[iter, 1]) && !force)
                {
                    CellWithMinValid(grid, out row, out col);
                    if (row != -1)
                    {
                        force = false;
                        iter++;
                        grid[row][col] = 1;
                        steps[iter, 0] = row;
                        steps[iter, 1] = col;
                    }
                }
                else
                {
                    curval++;
                    if (curval < 10)
                    {
                        grid[steps[iter, 0]]
                            [steps[iter, 1]] = curval;
                        force = false;
                    }
                    else
                    {
                        grid[steps[iter, 0]]
                            [steps[iter, 1]] = 0;
                        force = true;
                        iter--;
                    }
                }
            }
        }

        static void SolveAll(string file)
        {
            int[][][] sud = ReadSudoku(file);
            int tot = 0;

            for (int i = 0; i < sud.Length; i++)
            {
                SolveSudoku(sud[i]);
                tot += (sud[i][0][0] * 100) + 
                    (sud[i][0][1] * 10) + 
                    (sud[i][0][2]);
            }

            Console.WriteLine("Answer to projectEuler {0}\nLast sudoku solution:\n",tot);
            ShowSudoku(sud[49]);
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch sw = 
                new System.Diagnostics.Stopwatch();
            sw.Start();
            SolveAll("C:\\ProjectEuler\\sudoku.txt");
            sw.Stop();
            Console.WriteLine("\n50-ty sudoku`s solved in {0} ms", 
                            sw.ElapsedMilliseconds);
            Console.ReadLine();
        }
    }

As always - have fun!.

J and Project Euler p35

2008-05-29T08:49:00.000-07:00

Problem definition:
------------------------------------------
The number, 197, is called a circular prime because all rotations of the digits: 197, 971, and 719, are themselves prime.

There are thirteen such primes below 100: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, and 97.

How many circular primes are there below one million?
------------------------------------------

This time again solution in J language:


   rot =: |: ". (i.@:#) ((({:,}:) ^: )) (@: ":)
   ffac =: 0{|: @: q: @: rot
   circ =: ((+/@: (ffac = rot)) = (#@:":))
   tot =: +/ @: (circ"0) @: (2&+)
   tot i.999998

Learning J language - euler problem36

2008-05-27T11:27:00.001-07:00

Project euler problem 36 definition:
--------------------------------------
The decimal number, 585 = 1001001001 (binary), is palindromic in both bases.

Find the sum of all numbers, less than one million, which are palindromic in base 10 and base 2.

(Please note that the palindromic number, in either base, may not include leading zeros.)
--------------------------------------

So, I has some interest in array programming language J. I try to learn it a bit. And as such, I tried to solve this euler problem with the help of J language. What I learned is that this J language is VERY powerful and also it`s programs is very short, but also it is very hard to learn this language and of course - hard to understand the code (if you don`t know the language). But as soon as you proceed - more things are showing up from the fog :-).

So this is the solution:

+/@:((0:`(0&+))@.((+/@:(=|.)=#)@:":*.(+/@:(=|.)=#)@:#:)"0)i.1000000

One shot - One Kill, oh sorry, I mean: One line - One Goal :-)

As you see this code is not much understandable, but I have a little bit more readable code, acctually code is the same- just with defined function names.
This is it:

equalbits=: +/@:(=|.)
decpalind =: (equalbits=#)@:":
binpalind =: (equalbits=#)@:#:
bothpalind =: (decpalind*.binpalind)"0
palindsum =: +/@:((0:`(0&+))@.bothpalind)
palindsum i.1000000

Also just 6 lines of code. And a bit more readable, but not more if you don`t know some J basics. This language is very powerfull - especially at data processing.
If you want to learn J some basics, I would recommend this online book (also comes for free as documentation in J installation pack). Note - J interpreter is also free. This book is very good introduction, without it I would not take a look at J.
So all in all J - powerfull,interesting,strange,hard language in the same time, IMHO.
Have fun !!

Project Euler p102 - finding origin

2008-05-21T05:11:00.000-07:00

Problem description:
-----------------------------------------------------
Three distinct points are plotted at random on a Cartesian plane, for which -1000 ≤ x, y ≤ 1000, such that a triangle is formed.

Consider the following two triangles:

A(-340,495), B(-153,-910), C(835,-947)

X(-175,41), Y(-421,-714), Z(574,-645)

It can be verified that triangle ABC contains the origin, whereas triangle XYZ does not.

Using triangles.txt (right click and 'Save Link/Target As...'), a 27K text file containing the co-ordinates of one thousand "random" triangles, find the number of triangles for which the interior contains the origin.

NOTE: The first two examples in the file represent the triangles in the example given above.
-----------------------------------------------------

I`ve exploited the idea that for origin to be in the triangle interior, triangle lines segments should cross x and y axis 2 times: x+ / x- and y+ / y-.
Basically it is enough to check does triangle cross y+ and y- axis parts, but I didn`t know that and checked both axis. This worked too :-)

Solution ( C# as always):


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace ConsoleApplication1
{
    class Program
    {
        static int[][] ReadTriangleCoordinates(string file)
        {
            int[][] res;
            StreamReader sr = new StreamReader(file);
            res = Array.ConvertAll(sr.ReadToEnd().TrimEnd().Split("\n".ToCharArray()), 
                  s => {return Array.ConvertAll(s.Split(",".ToCharArray()),
                                         n => { return int.Parse(n); }) ;});
            sr.Close();
            return res;
        }

        static void CrossesAxisAt(int x1, 
                                  int y1, 
                                  int x2, 
                                  int y2, 
                                  out double? xAt, 
                                  out double? yAt)
        {
            int dx = x2 - x1;
            int dy = y2 - y1;
            int xmin = Math.Min(x1, x2);
            int ymin = Math.Min(y1, y2);
            int xmax = Math.Max(x1, x2);
            int ymax = Math.Max(y1, y2);
            double? xcross;
            double? ycross;
            double a;
            double b;

            if (dx == 0)
            {
                xcross = x1;
                ycross = null;
            }
            else if (dy == 0)
            {
                xcross = null;
                ycross = y1;
            }
            else
            {
                a = (double)dy / (double)dx;
                b = y1 - (a * x1);
                xcross = -(b / a);
                ycross = b;
            }

            if (xcross != null)
            {
                if (xcross < xmin || xcross > xmax 
                                  || ymin > 0 || ymax < 0)
                    xcross = null;
            }

            if (ycross != null)
            {
                if (ycross < ymin || ycross > ymax 
                                  || xmin > 0 || xmax < 0)
                    ycross = null;
            }

            xAt = xcross;
            yAt = ycross;
        }

        static bool ContainsOrigin(int[] abc)
        {
            double? xc1, xc2, xc3;
            double? yc1, yc2, yc3;
            int xposc = 0, xnegc = 0, yposc = 0, ynegc = 0;
            bool res;

            CrossesAxisAt(abc[0], abc[1], abc[4], abc[5], out xc1, out yc1);
            CrossesAxisAt(abc[0], abc[1], abc[2], abc[3], out xc2, out yc2);
            CrossesAxisAt(abc[2], abc[3], abc[4], abc[5], out xc3, out yc3);

            xposc = ((xc1 > 0) ? 1 : 0) + ((xc2 > 0) ? 1 : 0) + ((xc3 > 0) ? 1 : 0);
            xnegc = ((xc1 < 0) ? 1 : 0) + ((xc2 < 0) ? 1 : 0) + ((xc3 < 0) ? 1 : 0);
            yposc = ((yc1 > 0) ? 1 : 0) + ((yc2 > 0) ? 1 : 0) + ((yc3 > 0) ? 1 : 0);
            ynegc = ((yc1 < 0) ? 1 : 0) + ((yc2 < 0) ? 1 : 0) + ((yc3 < 0) ? 1 : 0);

            if (xposc > 0 && xnegc > 0 && yposc > 0 && ynegc > 0)
                res = true;
            else
                res = false;

            return res;
        }

        static int TrianglesWithOrigin(string file)
        {
            int[][] tr = ReadTriangleCoordinates(file);
            int trcount = 0;

            for (int i = 0; i < tr.Length; i++)
            {
                if (ContainsOrigin(tr[i]))
                    trcount++;
            }

            return trcount;
        }

        static void Main(string[] args)
        {
            Console.WriteLine(TrianglesWithOrigin("C:\\ProjectEuler\\triangles.txt"));
            Console.ReadLine();
        }
    }
}

IBM Ponder This, May 2008 - cellular automaton

2008-05-18T12:21:00.000-07:00

From the IBM project "ponder this", may 2008:
-------------------------------------------------
This month's puzzle is about Conway's Game of life.
The game is a cellular automaton on a grid where in each generation, every cell determines its state (dead or alive) based on its eight neighbors: a dead cell will change its state to alive if and only if exactly three of its neighbors are alive; A live cell will stay alive if and only if two or three of its neighbors are alive.
All the changes are done simultaneously.
Find a possible state of the universe whose next generation is depicted below.


....................
....................
..###..###...#...#..
...#....#.#..##.##..
...#....##...#.#.#..
...#....#.#..#...#..
..###..###...#...#..
....................
....................

"." denotes dead cells and "#" denotes live cells.
-------------------------------------------------

So after a bit of thinking i chose Selfish Gene Algorithm, which I also presented on my blog article. Please note - I will present solution which will use 4 classes from the given my blog link above. So if you need to test this solution - download the code from that article also and import these classes:
- SGAallele.cs
- SGAgenotype.cs
- SGAlocus.cs
- SGAlocusGroup.cs
Now about solution. In general - solution with the help of selfish gene algorithm, tries to evolve cellular automaton configuration which next generation is similar to the one described in problem definition. Fitness is calculated as wrong cell amount - that is cell, which is not in place of target configuration, amount. Algorithm tries to minimize that wrong cell amount. So solution is this:


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SGA;

namespace ConsoleApplication18
{
    class Program
    {
        static int counter = 0;
        static bool[][] grid = InitialGrid();
        static SGAgenotype gen;

        static void ShowCells(bool[][] cells)
        {
            string outp;
            outp = cells.Aggregate("", (acc1, ba) => 
                                   acc1 + ba.Aggregate("", (acc2, b) => 
                                        acc2 + ((b) ? "#" : ".")) + "\n");
            Console.Write(outp);
        }

        static bool[][] InitialGrid()
        {
            bool[][] c = new bool[9][];
            c[0] = new bool[20];
            c[1] = new bool[20];
            c[2] = new bool[] { false, false, true, true, true, false, false, 
                                true, true, true, false, false, false, true, 
                                false, false, false, true, false, false };
            c[3] = new bool[] { false, false, false, true, false, false, false, 
                                false, true, false, true, false, false, true, 
                                true, false, true, true, false, false };
            c[4] = new bool[] { false, false, false, true, false, false, false, 
                                false, true, true, false, false, false, true, 
                                false, true, false, true, false, false };
            c[5] = new bool[] { false, false, false, true, false, false, false, 
                                false, true, false, true, false, false, true, 
                                false, false, false, true, false, false };
            c[6] = new bool[] { false, false, true, true, true, false, false, true, 
                                true, true, false, false, false, true, false, 
                                false, false, true, false, false };
            c[7] = new bool[20];
            c[8] = new bool[20];

            return c;
        }

        static int LiveCells(bool[][] grid, int row, int col)
        {
            int live = 0;

            for (int i = row - 1; i <= row + 1; i++)
            {
                for (int j = col - 1; j <= col + 1; j++)
                {
                    if (i != row || j != col)
                    {
                        if (grid[i][j])
                            live++;
                    }
                }
            }

            return live;
        }

        static bool WillLive(bool IsLive, int LiveNeighbours)
        {
            bool state;

            if (IsLive)
                state = (LiveNeighbours < 2 || LiveNeighbours > 3) ? false : true;
            else
                state = (LiveNeighbours == 3) ? true : false;
           
            return state;
        }

        static bool[][] ComputeStates(bool[][] grid)
        {
            bool[][] res = new bool[9][];
            res[0] = new bool[20];
            res[1] = new bool[20];
            res[2] = new bool[20];
            res[3] = new bool[20];
            res[4] = new bool[20];
            res[5] = new bool[20];
            res[6] = new bool[20];
            res[7] = new bool[20];
            res[8] = new bool[20];

            for (int i = 1; i < grid.Length - 1; i++)
            {
                for (int j = 1; j < grid[i].Length - 1; j++)
                {
                    res[i][j] = WillLive(grid[i][j], LiveCells(grid, i, j));
                }
            }

            return res;
        }

        static bool[][] CellGridFromGenotype(SGAgenotype gen, SGAgenotype.GenomeNo num)
        {
            int row, col;
            bool v;
            int gene = -1;
            bool[][] res = new bool[9][];
            res[0] = new bool[20];
            res[1] = new bool[20];
            res[2] = new bool[20];
            res[3] = new bool[20];
            res[4] = new bool[20];
            res[5] = new bool[20];
            res[6] = new bool[20];
            res[7] = new bool[20];
            res[8] = new bool[20];

            for (int i = 0; i < gen.locusGroups.Length; i++)
            {
                row = (i) / 20;
                col = (i) - (row * 20);

                if (num == SGAgenotype.GenomeNo.FirstGenome)
                    gene = gen.locusGroups[i].loci[0].FirstGeneIndex;
                else if (num == SGAgenotype.GenomeNo.SecondGenome)
                    gene = gen.locusGroups[i].loci[0].SecondGeneIndex;
                else if (num == SGAgenotype.GenomeNo.None)
                    gene = gen.locusGroups[i].loci[0].BestGeneIndex;

                v = Convert.ToBoolean(gen.locusGroups[i].loci[0].Alleles[gene].Value);

                res[row][col] = v;
            }

            return res;
        }

        static int WrongCellCount(bool[][] previousStates)
        {
            int res = 0;
            bool[][] nextStates = ComputeStates(previousStates);

            for (int i = 0; i < grid.Length; i++)
            {
                for (int j = 0; j < grid[i].Length; j++)
                {
                    if (nextStates[i][j] != grid[i][j])
                        res++;
                }
            }

            return res;
        }

        static SGAgenotype.CompareResults CompareCellAutomata()
        {
            SGAgenotype.CompareResults cr;
            bool[][] grid1 = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.FirstGenome);
            bool[][] grid2 = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.SecondGenome);
            bool[][] gridb = CellGridFromGenotype(gen, SGAgenotype.GenomeNo.None);

            int wcellsGene1 = WrongCellCount(grid1);
            int wcellsGene2 = WrongCellCount(grid2);
            int wcellsGeneb = WrongCellCount(gridb);

            cr.CompareEachOther = SGAgenotype.GenomeNo.None;
            cr.CompareWithBest = SGAgenotype.GenomeNo.None;

            if (wcellsGene1 < wcellsGene2)
                cr.CompareEachOther = SGAgenotype.GenomeNo.FirstGenome;
            else if (wcellsGene2 < wcellsGene1)
                cr.CompareEachOther = SGAgenotype.GenomeNo.SecondGenome;

            if (wcellsGene1 < wcellsGeneb)
                cr.CompareWithBest = SGAgenotype.GenomeNo.FirstGenome;
            else if (wcellsGene2 < wcellsGeneb)
                cr.CompareWithBest = SGAgenotype.GenomeNo.SecondGenome;

            counter++;

            if (counter%1000==0)
            {
                Console.Clear();
                ShowCells(gridb);
                Console.WriteLine("====================");
                ShowCells(ComputeStates(gridb));
                Console.WriteLine("====================\nPlease wait...");
            }

            return cr;
        }
        
        static void Main(string[] args)
        {
            gen = new SGAgenotype(new SGAlocus[] { new SGAlocus(0, 1, 1) }, 
                                       180, CompareCellAutomata);
            gen.AlleleAffectValue *= 0.1;
            gen.MutationProbability *= 1.5;
            gen.Evolve(400000);
            Console.WriteLine("DONE.");
            Console.ReadLine();
        }
    }
}

After execution we get this picture-


...............#....
.#....#.....#.....##
...##...##..........
...#....#....##.##.#
....#....##..#...#.#        ANSWER
..#.....#.....#..#..
...#....##..#.....#.
...#..#.....#...#...
....................
====================
....................
....................
..###..###...#...#..
...#....#.#..##.##..
...#....##...#.#.#..    NEXT GENERATION OF CELLS
...#....#.#..#...#..
..##...###...#...#..
....................
....................
====================

First picture is answer and second (below) is next generation of cells.
We see that this is approximate answer because next generation differs by 1 cell from the given situation in problem definition. Well, this is because I`ve used genetic algorithm which by definition is probabilistic and doesn`t guarantee to give exact solution. But very well calculates approximate solutions.
So all in all - it was fun to play with selfish gene algorithm and cellular automaton.

Traveler's dilemma and bonus factor

2008-05-13T00:05:00.000-07:00

Traveler's dilemma description is following (taken from wikipedia):
-------------------------------------------------------
An airline loses two suitcases belonging to two different travelers. Both suitcases happen to be identical and contain identical antiques. An airline manager tasked to settle the claims of both travelers explains that the airline is liable for a maximum of $100 per suitcase, and in order to determine an honest appraised value of the antiques the manager separates both travelers so they can't confer, and asks them to write down the amount of their value at no less than $2 and no larger than $100. He also tells them that if both write down the same number, he will treat that number as the true dollar value of both suitcases and reimburse both travelers that amount. However, if one writes down a smaller number than the other, this smaller number will be taken as the true dollar value, and both travelers will receive that amount along with a bonus/malus: $2 extra will be paid to the traveler who wrote down the lower value and a $2 deduction will be taken from the person who wrote down the higher amount. The challenge is: what strategy should both travelers follow to decide the value they should write down?
-------------------------------------------------------

Experiment was conducted as following:
- BONUS was chosen between 2 and 20
- Possible travelers request was between BONUS and 100.
- Travelers chooses cost randomly (for no-strategy imitation)
- Cost random selection was repeated 5000000 times
- After that optimal request was calculated as travelers choise which gives the greatest profit.
- And finally these optimal requests are plot in diagram as function of BONUS.

C# code which conducted experiments on random travelers choise is this:


    class Program
    {
        static void Payoff(int Aoffer, int Boffer, int extra, out int Apayoff, out int Bpayoff)
        {
            int Aextra, Bextra;
            int minoffer = Math.Min(Aoffer, Boffer);

            if (Aoffer == Boffer)
            {
                Aextra = Bextra = 0;
            }
            else
            {
                Aextra = (Aoffer == minoffer) ? extra : -extra;
                Bextra = (Boffer == minoffer) ? extra : -extra;
            }
            Apayoff = minoffer + Aextra;
            Bpayoff = minoffer + Bextra;
        }

        static long OptimalSelection(int costFrom, int costTo, int goods)
        {
            long[] payoffA = new long[costTo + 1];
            long[] payoffB = new long[costTo + 1];
            Random r = new Random();
            int Ac, Bc;
            int Pa, Pb;
            long Aopt = -1, Bopt = -1;
            long Amax = 0, Bmax = 0;

            for (int i = 1; i <= goods; i++)
            {
                Ac = r.Next(costFrom, costTo + 1);
                Bc = r.Next(costFrom, costTo + 1);
                Payoff(Ac, Bc, costFrom, out Pa, out Pb);
                payoffA[Ac] += Pa;
                payoffB[Bc] += Pb;
            }

            for (int j = costFrom; j <= costTo; j++)
            {
                if (payoffA[j] > Amax) { Amax = payoffA[j]; Aopt = j; }
                if (payoffB[j] > Bmax) { Bmax = payoffB[j]; Bopt = j; }
            }

            return (Aopt + Bopt) / 2;
        }

        static void MeasureBonusInfluence()
        {
            double[] bonus = Array.ConvertAll(Enumerable.Range(2, 20).ToArray(),
                             v => (double)(v));
            double[] optimal = Array.ConvertAll(bonus,
                             v => (double)OptimalSelection((int)v, 100, 5000000));

            GoogleChart.PlotLineChart(bonus,
                                      optimal,
                                      "Bonus ($)",
                                      "Optimal request ($)",
                                      "Traveler's dilemma - bonus factor",
                                      "400x200");
        }

        static void Main(string[] args)
        {
            MeasureBonusInfluence();
        }
    }

    class GoogleChart
    {
        public static void PlotLineChart(double[] Xvalues,
                                        double[] Yvalues,
                                        string Xlabel,
                                        string Ylabel,
                                        string Title,
                                        string imgsize)
        {
            string urltemplate = "http://chart.apis.google.com/chart?cht=lxy&chs=@imgsize@&chd=t:@xvalues@%7C@yvalues@&chds=@xmin@,@xmax@,@ymin@,@ymax@&chxr=0,@xmin@,@xmax@%7C1,@ymin@,@ymax@&chxt=x,y,x,y&chxl=2:%7C@xlabel@%7C3:%7C@ylabel@%7C&chxp=2,50%7C3,50&chtt=@title@";
            string xvalagg = Xvalues.Aggregate("", (old, next) =>
                             old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));
            string yvalagg = Yvalues.Aggregate("", (old, next) =>
                             old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));

            if (imgsize == null) imgsize = "300x200";

            urltemplate = urltemplate.Replace("@imgsize@", imgsize);
            urltemplate = urltemplate.Replace("@xvalues@", xvalagg);
            urltemplate = urltemplate.Replace("@yvalues@", yvalagg);
            urltemplate = urltemplate.Replace("@xmin@", Xvalues.Min().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@xmax@", Xvalues.Max().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@ymin@", Yvalues.Min().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@ymax@", Yvalues.Max().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@xlabel@", Xlabel);
            urltemplate = urltemplate.Replace("@ylabel@", Ylabel);
            urltemplate = urltemplate.Replace("@title@", Title);

            System.Diagnostics.Process.Start(urltemplate);
        }
    }

Finally after execution of this program, we get this picture of optimal request relation with bonus amount (if travelers chooses randomly):

Google chart API in .NET

2008-05-07T10:26:00.000-07:00

Yes I know - there are good .NET libraries for using google chart API, but I was interested to dig a little bit deeper into google chart API. And I have to say that google chart is worse than opensource program GNUPLOT (in the sense of functionality), but ok whatever... And second target was - to have as less code (including libraries) as possible for charting VERY simple XY graphs. So after studing a bit google chart documentation- I wrote very primitive class for plotting simple XY graph. For plotting graphs with this code you need:
- Internet connection
- .NET framework 3.5
- Internet browser (whatever you like)
So this is the code:


    class Program
    {
        static void Main(string[] args)
        {
            double[] xval;
            double[] ysin, yexp;

            xval = new double[] { 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5 };
            ysin = Array.ConvertAll(xval, v => { return Math.Sin(v); }); 
            yexp = Array.ConvertAll(xval, v => { return Math.Exp(v); }); 
            
            GoogleChart.PlotLineChart(xval, ysin, "x values", "sin(x)", "sin(x) function", null);
            GoogleChart.PlotLineChart(xval, yexp, "x values", "exp(x)", "exp(x) function", null);
        }
    }

    class GoogleChart
    {
        public static void PlotLineChart(double[] Xvalues, 
                                        double[] Yvalues, 
                                        string Xlabel, 
                                        string Ylabel, 
                                        string Title, 
                                        string imgsize)
        {
            string urltemplate = "http://chart.apis.google.com/chart?cht=lxy&chs=@imgsize@&chd=t:@xvalues@%7C@yvalues@&chds=@xmin@,@xmax@,@ymin@,@ymax@&chxr=0,@xmin@,@xmax@%7C1,@ymin@,@ymax@&chxt=x,y,x,y&chxl=2:%7C@xlabel@%7C3:%7C@ylabel@%7C&chxp=2,50%7C3,50&chtt=@title@";
            string xvalagg = Xvalues.Aggregate("", (old, next) => 
                             old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));
            string yvalagg = Yvalues.Aggregate("", (old, next) => 
                             old + ((old == "") ? "" : ",") + next.ToString(CultureInfo.InvariantCulture));

            if (imgsize == null) imgsize = "300x200";

            urltemplate = urltemplate.Replace("@imgsize@", imgsize);
            urltemplate = urltemplate.Replace("@xvalues@", xvalagg);
            urltemplate = urltemplate.Replace("@yvalues@", yvalagg);
            urltemplate = urltemplate.Replace("@xmin@", Xvalues.Min().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@xmax@", Xvalues.Max().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@ymin@", Yvalues.Min().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@ymax@", Yvalues.Max().ToString(CultureInfo.InvariantCulture));
            urltemplate = urltemplate.Replace("@xlabel@", Xlabel);
            urltemplate = urltemplate.Replace("@ylabel@", Ylabel);
            urltemplate = urltemplate.Replace("@title@", Title);

            System.Diagnostics.Process.Start(urltemplate);
        }
    }

With this sample code there are generated these nice graphs:

This is it. Have fun with charting !

Project Euler p59 - breaking encryption

2008-05-06T05:44:00.000-07:00

Problem definition:
------------------------------------------------
Each character on a computer is assigned a unique code and the preferred standard is ASCII (American Standard Code for Information Interchange). For example, uppercase A = 65, asterisk (*) = 42, and lowercase k = 107.

A modern encryption method is to take a text file, convert the bytes to ASCII, then XOR each byte with a given value, taken from a secret key. The advantage with the XOR function is that using the same encryption key on the cipher text, restores the plain text; for example, 65 XOR 42 = 107, then 107 XOR 42 = 65.

For unbreakable encryption, the key is the same length as the plain text message, and the key is made up of random bytes. The user would keep the encrypted message and the encryption key in different locations, and without both "halves", it is impossible to decrypt the message.

Unfortunately, this method is impractical for most users, so the modified method is to use a password as a key. If the password is shorter than the message, which is likely, the key is repeated cyclically throughout the message. The balance for this method is using a sufficiently long password key for security, but short enough to be memorable.

Your task has been made easy, as the encryption key consists of three lower case characters. Using cipher1.txt (right click and 'Save Link/Target As...'), a file containing the encrypted ASCII codes, and the knowledge that the plain text must contain common English words, decrypt the message and find the sum of the ASCII values in the original text.
------------------------------------------------

C# brute-force solution, running in 42 seconds:


        static byte[] ReadCipher(string file)
        {
            StreamReader sr = new StreamReader(file);
            string[] buf = sr.ReadToEnd().Trim().Split(",".ToCharArray());
            byte[] res = Array.ConvertAll(buf, s => { return byte.Parse(s); });
            sr.Close();

            return res;
        }

        static byte[] TryDecipher(byte[] cip, string key)
        {
            int ind = 0;
            byte[] dec = new byte[cip.Length];

            for (int i = 0; i < cip.Length; i++)
            {
                dec[i] = (byte)((int)cip[i] ^ (int)key[ind]);
                ind = (ind >= 2) ? 0 : ind + 1;
            }

            return dec;
        }

        static int KeywordsInText(string text, string[] keywords)
        {
            int kc = 0;

            for (int i = 0; i < keywords.Length; i++)
            {
                if (text.ToLower().IndexOf(" " + keywords[i] + " ") > -1) kc++;
            }

            return (kc*100) / keywords.Length;
        }

        static void Decipher(string file)
        {
            string[] keywords = new string[] {"the","be","is","are","was","to","of","and","a","in","into","that","have","has","had","I","it","for","not","on","with","he","as","you","do","did","done","at","this","but","his","by","from","they","we","say","her","she","or","an","will","shall","my","one","all","would","should","there","here","their","what","so","up","out","if","about","who","get","got","which","go","went","gone","me"};
            string passchars = "qwertyuioplkjhgfdsazxcvbnm";
            string curkey = "";
            byte[] cipher = ReadCipher(file);
            byte[] decoded = new byte[cipher.Length];
            string dectext;
            int keywordcount = 0;
            int maxkeywords = 0;
            string bestkey = "";
            byte[] bestbytes = new byte[cipher.Length];
            string besttext = "";
            long bytesum = 0;

            for (int i1 = 0; i1 < passchars.Length; i1++)
            {
                for (int i2 = 0; i2 < passchars.Length; i2++)
                {
                    for (int i3 = 0; i3 < passchars.Length; i3++)
                    {
                        curkey = passchars[i1].ToString() + passchars[i2].ToString() + passchars[i3].ToString();
                        decoded = TryDecipher(cipher, curkey);
                        dectext = ASCIIEncoding.ASCII.GetString(decoded);
                        keywordcount = KeywordsInText(dectext, keywords);
                        
                        if (keywordcount > maxkeywords)
                        {
                            maxkeywords = keywordcount;
                            bestkey = curkey;
                            bestbytes = decoded;
                            besttext = dectext;
                        }
                    }
                }
            }

            for (int i = 0; i < bestbytes.Length; i++)
            {
                bytesum += (long)bestbytes[i];
            }

            Console.WriteLine("Byte sum: {0}\nPassword: {1}\nDecoded message 50 chars: {2}\nMessage has {3}% common words",
                                      bytesum,bestkey,besttext.Substring(0,50),maxkeywords);
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
            sw.Start();
            Decipher("C:\\ProjectEuler\\cipher1.txt");
            sw.Stop();
            Console.WriteLine("Decryption time: {0} seconds",sw.Elapsed.Seconds);
            Console.ReadLine();
        }

Project Euler p79 - cracking passcode

2008-04-29T14:18:00.000-07:00

This one is very interesting. Problem description:
-------------------------------------------
A common security method used for online banking is to ask the user for three random characters from a passcode. For example, if the passcode was 531278, they may asked for the 2nd, 3rd, and 5th characters; the expected reply would be: 317.
The text file, keylog.txt, contains fifty successful login attempts.
Given that the three characters are always asked for in order, analyse the file so as to determine the shortest possible secret passcode of unknown length.
-------------------------------------------

So solution in C# is this, not optimized, but works (in a 16 seconds):


        static int[] ReadKeyLog(string file)
        {
            int[] dig;
            string buf;
            string[] barr;
            StreamReader sr = new StreamReader(file);
            buf = sr.ReadToEnd().TrimEnd("\n".ToCharArray());
            barr = buf.Split("\n".ToCharArray());
            dig = Array.ConvertAll(barr, s => { return Int32.Parse(s); });
            sr.Close();
            return dig;
        }

        static bool KeyPartExists(long testnumber, int keypart)
        {
            int n1 = keypart / 100;
            int n2 = (keypart / 10) % 10;
            int n3 = keypart % 10;
            bool n1e = false;
            bool n2e = false;
            bool n3e = false;

            for (long i = testnumber; i > 0; i /= 10)
            {
                if (!n3e)
                    n3e = (i % 10) == n3;
                else
                    if (!n2e)
                        n2e = (i % 10) == n2;
                    else
                        if (!n1e)
                            n1e = (i % 10) == n1;
                        else
                            break;
            }
            if (n1e && n2e && n3e)
                return true;
            else
                return false;
        }

        static bool HasAllKeyParts(long testnum, int[] keylog)
        {
            bool res = true;

            for (int i = 0; i < keylog.Length; i++)
            {
                if (!KeyPartExists(testnum, keylog[i]))
                {
                    res = false;
                    break;
                }
            }
            return res;
        }

        static long FindPassCode(string file)
        {
            int[] keylog = ReadKeyLog(file);
            long num = 10000;
            
            while (true)
            {
                if (HasAllKeyParts(num, keylog))
                    return num;
                num++;
            }
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
            sw.Start();
            Console.WriteLine("Bruteforced passcode is {0}",
                                FindPassCode("C:\\ProjectEuler\\keylog.txt")
            );
            sw.Stop();
            Console.WriteLine("Cracking time {0} seconds", sw.Elapsed.Seconds);
            Console.ReadLine();
        }

Project Euler Problem 17

2008-04-24T11:46:00.000-07:00

Definition:
---------------------------------------------
If the numbers 1 to 5 are written out in words: one, two, three, four, five, then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total.
If all the numbers from 1 to 1000 (one thousand) inclusive were written out in words, how many letters would be used?
NOTE: Do not count spaces or hyphens. For example, 342 (three hundred and forty-two) contains 23 letters and 115 (one hundred and fifteen) contains 20 letters. The use of "and" when writing out numbers is in compliance with British usage.

---------------------------------------------

C# code solution:

static long LettersOfNumber(long num, long[,] bases)
{
long count = 0;
long tempnum;
if (num > 1000 num < 1)
count = -1;
else if (num == 1000)
count = 11;
else
{
tempnum = num % 100;
if (tempnum < 20)
count += bases[tempnum, 0];
else
count += bases[num % 10, 0] + bases[(num / 10) % 10, 1];
if (num > 99)
count += bases[(num / 100), 0] + 7 + ((tempnum == 0) ? 0 : 3);
}
return count;
}
static long LettersUntilNumber(int num)
{
long letcount = 0;
long[,] bas = new long[20, 2] { { 0, 0 }, { 3, 0 }, { 3, 6 }, { 5, 6 }, { 4, 5 }, { 4, 5 }, { 3, 5 }, { 5, 7 }, { 5, 6 }, { 4, 6 }, { 3, 0 }, { 6, 0 }, { 6, 0 }, { 8, 0 }, { 8, 0 }, { 7, 0 }, { 7, 0 }, { 9, 0 }, { 8, 0 }, { 8, 0 } };
for (int i = 1; i <= num; i++)
letcount += LettersOfNumber(i, bas);
return letcount;
}
static void Main(string[] args)
{
Console.WriteLine(LettersUntilNumber(1000));
Console.ReadLine();
}