Monday, December 28, 2009

Corrupted Git

Well, I have a directory with all my undergraduate project and I use git to ease my work, but recently I was not able to create a backup in a external disk because a error:


remote: Counting objects: 313, done.
remote: Compressing objects: 100% (116/116), done.
fatal: pack has bad object at offset 11893474: inflate returned 1
fatal: index-pack failed


I have files of almost 500MB on this repo, I suppose that the problem. I try to repack but it was useless, I also can check the hash of the object (with fsck), but I can't find where is that file, or what commit it belongs to.

This code will help in this case:


git checkout master
git branch new_master
git checkout new_master
for i in $!$(seq 1 100); do git checkout new_master~$!$i || break; done
for j in $!$(seq $!$i 100); do git checkout new_master~$!$j && break; done
# The error is between new_master~$!$j and new_master~($!$i-2)
git checkout master
git branch -D new_master
git branch new_master
git checkout new_master
git rebase --onto new_master~$!$j new_master~$!$((i-2)) new_master


It works (al least for my problem!!!), and now i want to rebase master:

git checkout master
git rebase new_master
git branch -d new_master


That's all.

Thursday, December 03, 2009

Compute the normal vector of a hyperplane in $R^n$

Well, basically any vector $N$ that holds $N \dot (p_i - p_j) = 0$ for for any vector $p_i, p_j \in H$ $H$ is the hyperplane.

Only $n$ different points are needed, call $P$ the set of $n$ points in $H$, and define $P^* = \{ p^*_i = p_i - p_{i+1} : p_i,p_{i+1} \in H , i \in [1,n-1]\}$.

$P^*$ is the set of vectors on $H$.

Any normal vector $N$ that holds that $N \dot p^*_i , i \in [1,n-1]$ is a normal vector.

$N$ has $n$ components, so we got here $n$ unknowns. And with the $n_1$ $p^*_i$ vectors we got $n-1$ equations.

An extra equation could be defined if we want a normalized vector, but this add a quadratic term in the systems, as the vector with all zero components can't be a normal vector of any hyperplane at least one of its components has a non-zero value, suppose that the $j$ component of $N$ has value $1$. and we have then a Ax = B, with n-1 equations and unknowns. we just have to try which j makes the systm non singular and solve it.


The system of equations:
$\begin{pmatrix}
p^*_{1,1} & p^*_{1,2} & \cdots & p^*_{1,j-1} & p^*_{1,j+1} & \cdots & p^*_{1,n} \\ p^*_{2,1} & p^*_{2,2} & \cdots & p^*_{2,j-1} & p^*_{2,j+1} & \cdots & p^*_{2,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
p^*_{n,1} & p^*_{n,2} & \cdots & p^*_{n,j-1} & p^*_{n,j+1} & \cdots & p^*_{n,n} \\
\end{pmatrix} \begin{pmatrix} N_1 \\ N_2 \\ \vdots \\ N_{j-1} \\ N_{j+1} \\ \vdots \\ N_n \end{pmatrix} = - \begin{pmatrix} p^*_{1,j} \\ p^*_{2,j} \\ \vdots \\ p^*_{(j-1),j} \\ p^*_{(j-1),j} \\ \vdots \\ p^*_{n,j} \end{pmatrix}
$

Start checking for $j \in [1,n]$ and you will get at much $n$ and at least 1 normal vector.

Python:


from numpy import array, dot
from numpy.linalg import solve, norm


def normal(P):
n = len(P)

if n < 1:
raise Exception, "O-size list."

m = len(P[0])

if array([len(p) != m for p in P]).any():
raise Exception, "Inconsistent set of points."

if n != m:
raise Exception

P_ = [p-q for p,q in zip(P[:-1],P[1:])]

N = []

for i in range(n):
try:
x = solve(
array([list(p)[:i] + list(p)[i+1:] for p in P_]),
array([-p[i] for p in P_])
)

x = array(list(x)[:i] + [1] + list(x)[i:])
N.append(x)

except:
pass

n = N[array([(dot(n,n)-1)**2 for n in N]).argmin()]

return n / norm(n)



if __name__ == '__main__':
print normal([array([1,0,0]), array([0,1,0]), array([1,1,0])])
print normal([array([1,123,0123]), array([0345,1,034]), array([132,134,034])])
print normal([array([1,0,0]), array([0,1,0]), array([1,123,1123])])
I'm gonna create a image with Debian and Django to upload it to AWS.



dd if=/dev/zero of=debian-ami count=1000 bs=1m
sudo mkfs.ext3 -F debian-ami
mkdir /tmp/chroot
sudo mount -o loop debian-ami /tmp/chroot
sudo debootstrap --arch i386 lenny /tmp/chroot/ http://ftp.debian.org
sudo chroot /tmp/chroot/
# Inside the chroot
mount -t proc none /proc
cd /dev
MAKEDEV console
MAKEDEV std
echo -e 'auto lo\niface lo inet loopback\nauto eth0\niface eth0 inet dhcp' >> /etc/network/interfaces
echo -e 'proc /proc proc defaults 0 0\n/dev/sda1 / reiserfs defaults 0 1\n/dev/sda2 swap swap defaults 0 0' > /etc/fstab
aptitude update
aptitude install locales-all
aptitude install ssh
exit

# In your machine
sudo umount -l /tmp/chroot
export EC2_PRIVATE_KEY=xxxxxxx.pem
export EC2_CERT=xxxxxxxxxxxxxxxxxx.pem
export EC2_ACCNO=xxxxxxxxxx
export ACCESS_KEY=xxxxxxxxxxxxxxxx
export SECRET_KEY=xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
ec2-bundle-image -i debian-ami --cert EC2_CERT --privatekey EC2_PRIVATE_KEY -u EC2_ACCNO
ec2-upload-bundle -b 20091203-linux-debian-lenny -m /tmp/debian-ami.manifest.xml -a ACCESS_KEY -s SECRET_KEY
ec2-register -n 20091203-linux-debian-lenny 20091203-linux-debian-lenny/debian-ami.manifest.xml

Sunday, August 23, 2009

Algebra

I'm studying algebra with some friends, and i will try to write here what subjects are we studying right now.

We start with Abstract Algebra using Hungerford book, but i want to follow the MITOpenCourseWare about Algebra I (http://ocw.mit.edu/OcwWeb/Mathematics/18-701Fall-2007/Readings/index.htm).

Composition Law:
A composition law $f$ select a pair of elements on a set $S$ and return an element of $S$, ie:
$f : S \times S \rightarrow S$, I will write $f(a,b)$ as $a*b$ for seek of simplicity.

If $*$ holds that $(a*b)*c = a*(b*c) \quad \forall a,b,c \in S$ then we say that is an associative law.

If $*$ holds that $a*b=b*a \quad \forall a,b \in S$ we say that is a commutative law.

Semigroups

Let $S$ be a set and $*$ be a composition law on $S$, if $*$ is associative then $(S,*)$ is a semigroup.

Monoids

Let $(S,*)$ be a semigroup if $S$ has an identity element $e$ such that:

$\forall x \in S$ $x*e = x = e*x$,

then $(S,*)$ is a Monoid.

Groups

Let $(S,*)$ be a Monoid if:

$\forall x \in S$ $\exists x^{-1} \in S$ s.t. $x*x^{-1} = e = x^{-1}*x$,

then $(S,*)$ is a Group.

if $*$ is commutative then $(S,*)$ is a commutative or Abelian group.

Subgroups

Let $(G,*)$ be a group, if a proper subset $H$ of $G$ is also a group with $*$ then $H$ is a subgroup of $G$

Theorem

If $(G,*)$ is a group and $H$ is a nonempty subset of $G$:

$(H,*)$ is a subgroup of $(G,*)$ if and only if $a*b^{-1} \in H \quad \forall a,b \in H$

Relation

Let $A$ and $B$ be sets and $R$ a map from $A$ to $B$ s.t.:

$R : A \rightarrow B$
with $a\in A$ and $b=R(a) \in B$ and
$(a,b) = (a,R(a)) \in A \times B$
for seek of simplicity we write $(a,R(a)) \in A \times B$ as $aRb$.

Equivalence Relations

Let $S$ be a set and $R$ a relation in $S \times S$, if:

i) $aRa \quad \forall a \in S$ (Reflexive)

ii) $aRb \implies bRa$ (Symmetric)

iii) $aRb \wedge bRc \implies aRc$ (Transitive)

then $R$ is an Equivalence Relation on $S$. (When we speak of an equivalence relation we'll use the symbol $=$)

Partial Order Relations

A partial order relation is a relation $R$ in a set $S$ with a equivalence relation $=$ if $R$ is reflexive and transitive and:

$aRb \wedge bRa \implies a=b$ (antisymmetric)

We write a partial order relation as $\leq$. The above oreder is a non-strict relation.

A strict partial relation order $R$ in $S$ is:

i) $a\not R a \quad \forall a \in S$ (irreflexive)

ii) $aRb \wedge bRa \implies a=b$ (antisymmetric)

iii) $aRb \wedge bRc \implies aRc$ (Transitive)

We write a strict partial order relation as $<$.

Complete Order Relations

The above partial relations don't hold $\forall a,b \in S$, if they do then the order is a complete order.

Congruence Relation

Let $(G,*)$ be a Monoid, if a equivalence relation $=$ on $G$ also holds that:

$a_1=a_2 \wedge b_1 = b_2 \implies a_1*b_1 = a_2*b_2$,

then $=$ is also congruence relation.

Equivalence Classes

Let $S$ be a set with an equivalence relations $=$, the equivalence class $\bar{a}$ (or $[a]$) of an element $a \in S$ is:

$\bar{a} = \{ x \in S | x = a \}$

The class of all the equivalence classes on $S$ is denoted as $(S/=)$, and is called the quotient class of $S$ by $=$.

The union of all the equivalence classes is $S$, i.e.

$\bigcup_{a \in S} \bar{a} = A = \bigcup_{\bar{a} \in (S/=)} \bar{a}$ or $\bar{a} = \bar{b}$

If $a,b \in S$ either $\bar{a} \cap \bar{b} = \emptyset$

Theorem:
Let $(G,*)$ be a Monoid and $=$ a congruence relation on $G$, then:

$((G/=), *)$ is a Monoid with

$* : (G/=) \times (G/=) \rightarrow (G/=)$ and $[a]*[b] = [a*b] \quad [a],[b] \in G/=$ (sorry for the change on notation but the bar don't expand enough in the a*b)

If $(G,*)$ is Abelian then $(G/=,*)$ is Abelian too.

Group Homomorphism

Let $(G,\circ)$ and $(H,\diamond)$ be two groups and a function $f : G \rightarrow H$ such that:

$f(x \circ y ) = f(x) \diamond f(y) \quad \forall x,y \in G$,

then $f$ is a Group homomorphism.

If $f$ is injective then $f$ is a Group Monomorphism.

If $f$ is surjective then $f$ is a Group Epimorphism.

If $f$ is bijective then $f$ is a Group Isomorphism.

If a Group Homomorphism is $f: G \rightarrow G$ then $f$ is a Group Endomorphism of $G$.

If a Group Isomorphism is $f: G \rightarrow G$ then $f$ is a Group Automorphism of $G$.

Kernel

Let $(G,\circ)$ and $(H,\diamond)$ be two groups and a homomorphism $f : G \rightarrow H$. $e_G$ and $e_H$ are the ideintity element of each group, then, the kernel of $f$ is:

$Ker f = \{a \in G | f(a) = e_H\}$,

The kernel is a subset of $G$ that maps to the identity element on $H$ with the homomorphism $f$.

Theorem:
If $f : G \rightarrow H$ is a monomorphism $\Longleftrightarrow$ $Ker f = \{e_G\}$

Theorem:
If $f : G \rightarrow H$ is a isomorphism $\Longleftrightarrow$ $f^{-1} : H \rightarrow G$ is a homomorphism and $f(f^{-1}(a)) = a \quad \forall a \in H$ and $f^{-1}(f(a)) = a \quad \forall a \in G$

Cyclic groups

Let $(G,*)$ be a group and $X$ a subset of $G$.

Let $\{ H_i | i \in I \}$ be the family of all the subgroups ($H_i$) of $G$ that contains $X$ i.e. $X \subset H_i \forall i \in I$

$\bigcap_{i \in I} H_i$ is a subgroup of $G$ generated by $X$ and denoted $\langle X\rangle$

The elements of $X$ are the generators of $\langle X\rangle$.

Several different subsets on $G$ can generate the same $\langle X\rangle$, so, in general $\langle X\rangle = \langle Y \rangle$ with $X \neq Y$.

If $X$ is a finite set such that $X = \{a_1, a_2, \cdots , a_n \}$, we write $\langle a_1, a_2, \cdots, a_n \rangle$ instead of $\langle X\rangle$.

If $G = \langle a_1,a_2, \cdots, a_n \rangle$ with $a_i \in G$, then $G$ is said to be finitely generated.

If $a \in G$, then the subgroup $\langle a\rangle$ is the Cyclic Group or Cyclic Subgroup generaterd by a.

Theorem

Let $(G,*)$ be a group and $X$ a nonempty subset of $G$, then the subgroup $\langle X \rangle$ generated by $X$ consists of all finite products $a_1^{n_1}a_2^{n_2}\cdots a_1t^{n_t} \quad a_i \in X; \quad n_i \in \mathbb{Z}$.

And $\forall a \in G, \langle a \rangle = \{ a^n | n \in \mathbb{Z}\}$

Theorem

Let $H$ be a cycli subgroup of $(\mathbb{Z},+)$, either $H=\langle 0 \rangle$ or $H=\langle m \rangle$, with $m$ the least positive integer in $H$, then if $H \neq \langle 0 \rangle$, then $H$ is infinite.

Theorem

Every infinite cyclic group is isomorphic to $(\mathbb{Z},+)$ and every finite cyclic group of order $m$ is isomorphic to $(\mathbb{Z}_m,+)$ (¿¿¿$\langle m \rangle$???).

Cosets

First we define that $a \equiv b (mod m) \Longleftrightarrow m | a - b \Longleftrightarrow a-b \in \langle m \rangle$.

Let $(G,*)$ be a group, and $(H,*)$ a subgroup of $(G,*)$ and $a,b \in G$.

a is right congruent to $b$ modulo $H$, denoted as $a \equiv_r b (mod H)$ if $a*b^{-1} \in H$.

a is left congruent to $b$ modulo $H$, denoted as $a \equiv_l b (mod H)$ if $a^{-1}*b \in H$.

If $G$ is abelian then if one of the congruencies hold then the other does:

$a*b^{-1} \in H \Longleftrightarrow (a*b^{-1})^{-1} \in H$

$(a*b^{-1})^{-1} = b*a^{-1} = a^{-1}*b$

and

$a^{-1}*b \in H \Longleftrightarrow (a^{-1}*b)^{-1} \in H$

$(a^{-1}*b)^{-1} = b^{-1}*a = a*b^{-1}$

Theorem

Let $(H,*)$ be a subgroup of $(G,*)$, then

i) Right(Left) congruence modulo $H$ is an equivalence relation on $G$.

ii) The equivalence class of $a \in G$ with the right (left) congruence is the set $Ha = \{h*a | h \in H\} ($aH = \{a*h | h \in H\})

iii) $|Ha| = |H| = |aH| \quad \forall a \in G$

$Ha$ is called a right coset of $H$ in $G$.
$aH$ is called a left coset of $H$ in $G$.


Corollary

i) $G$ is the union of the right(left) cosets of $H$ on $G$.

ii) Two right(left) cosets of $H$ on $G$ are either disjoint or equal.

iii) $\forall a,b \in G, \quad Ha = Hb \Longleftrightarrow a*b^{-1} \in H$ and $aH = bH \Longleftrightarrow a^{-1}*b \in H$

iv) If $R$ is the set of distinct right(left) cosets of $H$ in $G$, and $L$ is the set of distinct left cosets of $H$ in $G$, then

Definition

Let $(G,*)$ be a group, and $(H,*)$ a subgroup of $(G,*)$ the index of $H$ in $G$, denoted as $[G : H]$ is the cardinal number of the set of distinct right(left) cosets. of $H$ in $G$.

Saturday, August 22, 2009

Runge Kutta in Haskell

Well, Runge-Kutta is a better solver that Euler so:


-- Runge-Kutta method (t_o, t_f,h, x_o, function(t,x))
runge t_o t_f h x_o f = case t_o < t_f of
True ->
do
let k_1 = f t_o x_o
let k_2 = f (t_o + h/2) (x_o + h * k_1 / 2)
let k_3 = f (t_o + h/2) (x_o + h * k_2 / 2)
let k_4 = f (t_o + h) (x_o + h * k_3)
[(t_o,x_o)] ++ ( runge (t_o + h) t_f h (x_o + h * (k_1 + 2*k_2 + 2*k_3 + k_4) / 6) f)
False ->
[]


-- Function
f :: Double -> Double -> Double
f t x = - x

-- Print in column
show_col :: [(Double, Double)] -> String
show_col x = case x of
(a:b) ->
do
let (t, v) = a
(show t) ++ " , " ++ (show v) ++ "\n" ++ (show_col b)
[] ->
""

main = putStrLn(show_col(runge 0 5 0.01 10 f))

Friday, August 21, 2009

Euler Method in Haskell

I decide today to learn haskell, but i can't get out of my head how could it be used in simulation, so, i will try with several basic solvers, starting with Euler Method to solve: y'(t) = f(t, y(t))


--Euler method (t_o, t_f,h, x_o, function(t,x))
euler1 :: Double -> Double -> Double -> Double -> (Double -> Double -> Double) -> [(Double, Double)]
euler1 t_o t_f h x_o f = case t_o < t_f of
True ->
[(t_o,x_o)] ++ (euler1 (t_o+h) t_f h (x_o+(h * f t_o x_o)) f)
False ->
[]



--Euler method (t_o, t_f,steps, x_o, function(t,x))
euler2 :: Double -> Double -> Integer -> Double -> (Double -> Double -> Double) -> [(Double, Double)]
euler2 t_o t_f steps x_o f = case steps > 0 of
True ->
do
let h = (t_f - t_o) / (fromIntegral steps)
[(t_o, x_o)] ++ euler2 (t_o + h) t_f (steps - 1) (x_o + (h * f t_o x_o)) f
False ->
[]



-- Function
f :: Double -> Double -> Double
f t x = - x

-- Print in column
show_col :: [(Double, Double)] -> String
show_col x = case x of
(a:b) ->
do
let (t, v) = a
(show t) ++ " , " ++ (show v) ++ "\n" ++ (show_col b)
[] ->
""

--main = putStrLn(show_col(euler1 0 1 0.00001 10 f))
main = putStrLn(show_col(euler2 0 1 100000 10 f))

Thursday, August 13, 2009

Qucs

Mmmm, Creo que por fin encontre un buen programa para simular circuitos, trabaja con SPICE y ademas es libre.

http://qucs.sourceforge.net/
http://qucs.sourceforge.net/docs.html
http://qucs.sourceforge.net/download.html

Wednesday, August 12, 2009

Symbolic Integration (Manuel Bronstein)

Well, now that i have a new book, Symbolic Integration I of Manuel Bronstein, i will implement some of the code examples on the initial chapter: polydivide, euclidean algorithm, square free factorization, all of them using GiNaC library.

PolyDivide

This one find the unique q and r in a field K[z] that holds a = b*q + r


lst PolyDivide(ex A, ex B, symbol x)
{
ex Q = 0;
ex R = A;

int d = R.degree(x) - B.degree(x);

while ((not R.is_zero()) and (d >=0))
{
ex T = R.lcoeff(x) / B.lcoeff(x) * pow(x,d);
Q = expand(Q + T);
R = expand(R - B*T);
d = R.degree(x) - B.degree(x);
}
return lst(Q,R);
}


PolyPseudoDivide
Same that last, but in a integral domain, like Z.


lst PolyPseudoDivide(ex A, ex B, symbol x)
{
ex b = B.lcoeff(x);
int N = A.degree(x) - B.degree(x) + 1;
ex Q = 0;
ex R = A;

int d = R.degree(x) - B.degree(x);

while ((not R.is_zero()) and (d >=0))
{
ex T = R.lcoeff(x) * pow(x,d);
N -= 1;
Q = expand(b * Q + T);
R = expand(b * R - T * B);
d = R.degree(x) - B.degree(x);
}

b = pow(b, N);
return lst( b * Q, b * R);
}


Well, testing both algorithm against the quo and rem that do the same computation in GiNaC, we have this:



The test is a division between two polynomials of the same n degree randomly generated, in x-axis is n and in y-axis is ellapsed time of the compute code in sec.


int main()
{
symbol x("x");
int i = 0;
ex a=0, b=0;

clock_t s, e;
double t1,t2, t3, t4, t5;

for(int n = 1; n < 100000 ; n += 1000)
{
srand ( time(NULL) );

for(; i < n; i++)
{
a += (rand() % 200 - 100) * pow(x,i);
b += (rand() % 200 - 100) * pow(x,i);
}

s = clock();
PolyDivide(a,b,x);
e = clock();
t1 = (e-s);

s = clock();
PolyPseudoDivide(a,b,x);
e = clock();
t2 = (e-s);

s = clock();
quo(a,b,x); rem(a,b,x);
e = clock();
t3 = (e-s);

ex q;
s = clock();
q = quo(a,b,x);a - b*q;
e = clock();
t4 = (e-s);

ex r;
s = clock();
r = rem(a,b,x);(a-r)/b;
e = clock();
t5 = (e-s);


t1 /= CLOCKS_PER_SEC;
t2 /= CLOCKS_PER_SEC;
t3 /= CLOCKS_PER_SEC;
t4 /= CLOCKS_PER_SEC;
t5 /= CLOCKS_PER_SEC;

cout << n << " " << t1 << " " << t2 << " " << t3 << " " << t4 << " " << t5 << endl;
}

return 0;
}

Wednesday, July 22, 2009

Create a new toolbox for Scilab (II)

Now i gonna follow the instructions at: http://www.scilab.org/contrib/index_contrib.php?page=toolbox_guide

1) First i need to create all the files needed to the toolbox:


toolbox$ cd scilab/scilab/contrib/
toolbox/scilab/scilab/contrib$ mkdir symbolic
toolbox/scilab/scilab/contrib$ cd symbolic/
toolbox/scilab/scilab/contrib/symbolic$ mkdir macros src sci_gateway help etc unit_tests demos includes
toolbox/scilab/scilab/contrib/symbolic$ touch readme.txt builder.sce loader.sce license.txt


2) Macros, well, basically i only have one macro at the symbolic module, but i gonna test the macro at the instructions, so:


toolbox/scilab/scilab/contrib/symbolic$ cd macros/
toolbox/scilab/scilab/contrib/symbolic/macros$ cat << EOF > foo1.sci
function [X]=foo1(A)
// This function returns the positive components of the A diagonal

// Check the type and the size of A
if type(A)<>1 then
error("type of input argument must be a double");
end
if size(A,1)<>size(A,2) then
error("input argument must be a square matrix");
end
//Extraction of the positive components
X=[];
for i=1:size(A,1)
if A(i,i)>0 then
X($+1)=A(i,i);
end
end
endfunction

EOF
toolbox/scilab/scilab/contrib/symbolic/macros$ cat << EOF > buildmacros.sce
mode(-1)
toolboxname='symbolic'
pathB=get_absolute_file_path('buildmacros.sce')
disp('Building macros in ' +pathB)
genlib(toolboxname+'lib',pathB,%t)
clear pathB genlib toolboxname

EOF
toolbox/scilab/scilab/contrib/symbolic/macros$ cat << EOF > loadmacros.sce
mode(-1)
pathL=get_absolute_file_path('loadmacros.sce')
disp('Loading macros in ' +pathL)
load(pathL+'/lib')
clear pathL

EOF


3) We're going to create some primitives now, basically is like the module, so, there is a sci_gateway path and a src path, in the first we define the gateways to our functions, and in the second we write the main code of the primitives.


toolbox/scilab/scilab/contrib/symbolic$ cd src/
toolbox/scilab/scilab/contrib/symbolic/src$ cat << EOF > vectsum.c
void vectsum(int n, double * a, double * b, double * y)
{
int k;
for (k = 0; k < n; ++k)
y[k] = a[k] + b[k];
}

EOF
jcardona@terminus:~/stuff/personal/symbolic/main/toolbox/scilab/scilab/contrib/symbolic/src$ cd ..
toolbox/scilab/scilab/contrib/symbolic$ cd sci_gateway/
toolbox/scilab/scilab/contrib/symbolic/sci_gateway$ cat << EOF > sci_sumab.c
#include "stack-c.h"
extern int vectsum(int n, double * a, double * b, double * y);

void sci_sumab(char *fname){
int l1, m1, n1, l2, m2, n2, l3, n;

/* 1 - Check the number of inputs/outputs arguments */
int minlhs=1, maxlhs=1, minrhs=2, maxrhs=2;
CheckRhs(minrhs,maxrhs) ;
CheckLhs(minlhs,maxlhs) ;

/* 2 - Check inputs arguments type, and get the size
and the address in the Scilab stack of the inputs
arguments
*/
GetRhsVar(1, "d", &m1, &n1, &l1);
GetRhsVar(2, "d", &m2, &n2, &l2);

/* 3 - Check that the inputs arguments have the same size */
/* it's possible to use the chekdims and getscalar
functions to make these checks
*/
n=m2*n2;
if( n1!=n2 || m1!=m2)
{
cerro("inputs arguments must have the same size");
return 0;
}
if(n1!=0 && m1!=0)
if(n1!=1 && m1!=1)
{
cerro("inputs arguments must be vectors");
return(0);
}


/* 4 - Create a new variable corresponding to the output argument */
CreateVar(3,"d",&m2,&n2,&l3);

/* 5 -call vectsum routine: returns in stk(l3) the sum of a and b*/
vectsum(n,stk(l1),stk(l2),stk(l3));

/* 6 - Specif ouput argument */
LhsVar(1) = 3;
return 0;
}

EOF



The builder for the primitives, we need a buildsrc.sce and a buildsci_gateway.sce, in the respective path. The first one creates a shared library that is linked to the second shared library (the one of the gateway).


toolbox/scilab/scilab/contrib/symbolic/src$ cat << EOF > buildsrc.sce
//ilib_for_link('symbolicsrc',['fun1.o','fun2.o','vectsum.o'],[],"c")
ilib_for_link('symbolicsrc',['vectsum.o'],[],"c")

EOF
toolbox/scilab/scilab/contrib/symbolic/sci_gateway$ cat << EOF >buildsci_gateway.sce
// must be run from this directory
ilib_name = 'libsymbolic' // interface library name
//files = ['sci_fun.o','sci_sumab.o']; // objects files
files = ['sci_sumab.o']; // objects files
libs = ["../src/libsymbolicsrc"] // other libs needed for linking
//table = [ /*'fun', 'sci_fun';*/
table = ['sumab','sci_sumab']; // table of (scilab_name,interface-name)
// do not modify below
ilib_build(ilib_name,table,files,libs)

EOF



4) The help files are located at the help path with a xml format, and a dtd defined by scilab, we only add a template for the files:


toolbox/scilab/scilab/contrib/symbolic/help$ cat << EOF > sumab.xml
<?xml version="1.0" encoding="UTF-8"?>
<refentry version="5.0-subset Scilab" xml:id="sumab" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<info>
<pubdate>$LastChangedDate: 2008-03-26 09:50:39 +0100 (mer., 26 mars 2008)$</pubdate>
</info>

<refnamediv>
<refname>sumab</refname>

<refpurpose>Purpose</refpurpose>
</refnamediv>

<refsynopsisdiv>
<title>Calling Sequence</title>

<synopsis>sequence</synopsis>
</refsynopsisdiv>

<refsection>
<title>Description</title>

<para>Do something</para>

<para>Add here a paragraph of the function description </para>
</refsection>

<refsection>
<title>Examples</title>

<programlisting role="example">exmaple</programlisting>
</refsection>

<refsection>
<title>Authors</title>

<simplelist type="vert">
<member>YOUR NAME</member>
</simplelist>
</refsection>
</refentry>

EOF
toolbox/scilab/scilab/contrib/symbolic/help$ cat sumab.xml | sed -e 's?sumab?foo1?' > foo1.xml
toolbox/scilab/scilab/contrib/symbolic/help$ cat << EOF > buildhelp.sce
mode(-1) //force silent execution
path=get_absolute_file_path('buildhelp.sce');//get the absolute path of this file
add_help_chapter("Symbolic",path);//add help chapter
xmltohtml(path,"Symbolic")
//clear the variable stack
clear path add_help_chapter get_absolute_file_path

EOF
toolbox/scilab/scilab/contrib/symbolic/help$ cat << EOF > loadhelp.sce
mode(-1) //force silent execution
path=get_absolute_file_path('loadhelp.sce');//get the absolute path of this file
add_help_chapter("Symbolic",path);//add help chapter
clear path add_help_chapter get_absolute_file_

EOF


5) Finally we set the global builder and loader:


toolbox/scilab/scilab/contrib/symbolic$ cat << EOF > builder.sce
mode(-1);
mainpathB=get_absolute_file_path('builder.sce');
chdir(mainpathB);
if isdir('src') then
chdir('src');
exec('buildsrc.sce');
chdir('..');
end
if isdir('sci_gateway') then
chdir('sci_gateway');
exec('buildsci_gateway.sce');
chdir('..');
end
if isdir('macros') then
chdir('macros');
exec('buildmacros.sce');
chdir('..');
end
if isdir('help') then
chdir('help');
exec('buildhelp.sce');
chdir('..');
end
clear mainpathB

EOF
toolbox/scilab/scilab/contrib/symbolic$ cat << EOF > loader.sce
mode(-1);
mainpathL=get_absolute_file_path('loader.sce');
chdir(mainpathL);
if isdir('sci_gateway') then
chdir('sci_gateway');
exec('loader.sce');
chdir('..');
end
if isdir('macros') then
chdir('macros');
exec('loadmacros.sce');
chdir('..');
end
if isdir('help') then
chdir('help');
exec('loadhelp.sce');
chdir('..');
end
clear mainpathL

EOF


6) Build it at scilab:


toolbox/scilab/scilab$ bin/scilab

-->exec("contrib/symbolic/builder.sce")

-->mode(-1);
Generate a loader file
Generate a Makefile
ilib_gen_Make: Copy compilation files (Makefile*, libtool...) to TMPDIR
ilib_gen_Make: Copy vectsum.c to TMPDIR
ilib_gen_Make: Modification of the Makefile in TMPDIR.
Running the Makefile
Generate a cleaner file
ans =

libsymbolicsrc.so
ilib_name =

libsymbolic
libs =

../src/libsymbolicsrc
Generate a gateway file
Generate a loader file
Generate a Makefile
ilib_gen_Make: Copy compilation files (Makefile*, libtool...) to TMPDIR
ilib_gen_Make: Copy sci_sumab.c to TMPDIR
ilib_gen_Make: Copy libsymbolic.c to TMPDIR
ilib_gen_Make: Modification of the Makefile in TMPDIR.
Running the makefile
Generate a cleaner file

Building macros in /home/jcardona/stuff/personal/symbolic/main/toolbox/scilab/scilab/contrib/symbolic/macros/

Building the master document:
SCI/contrib/symbolic/help

Building the manual file [html] in SCI/contrib/symbolic/help. (Please wait building ... this can take a while)
Warning : redefining function: get_absolute_file_path . Use funcprot(0) to avoid this message

-->exec("contrib/symbolic/loader.sce")

-->mode(-1);
Shared archive loaded.
Link done.
Shared archive loaded.
Link done.

Loading macros in /home/jcardona/stuff/personal/symbolic/main/toolbox/scilab/scilab/contrib/symbolic/macros/

-->sumab(1,2)
ans =

3.

-->sumab(2,2)
ans =

4.


We have then a brand new toolbox enjoy summing.

Tuesday, July 21, 2009

Create a new toolbox for Scilab (I)

As i have to change my symbolic module to a symbolic toolbox, i want to write this to use as a future reference.

1) Get Scilab:


$ mkdir toolbox
$ cd toolbox/
toolbox$ git clone git://git.scilab.org/scilab


2) First compilation:



toolbox$ cd scilab/scilab

toolbox/scilab/scilab$ ./configure
...
checking hdf5.h usability... no
checking hdf5.h presence... no
checking for hdf5.h... no
configure: error: Cannot find headers (hdf5.h) of the library HDF5. Please install the dev package



3) HDF5 Headers:

sudo aptitude install libhdf5-serial-dev



4) Installing JHDF5 :


wget http://www.hdfgroup.org/ftp/HDF5/hdf-java/src/hdf-java-2.5-src.tar
tar -xvf hdf-java-2.5-src.tar
cd hdf-java/
./configure --with-jdk=/usr/lib/jvm/java-6-sun-1.6.0.14/include/,/usr/lib/jvm/java-6-sun-1.6.0.14/lib/
make
sudo mkdir -p /usr/share/java/jar
sudo cp lib/jhdf5.jar /usr/share/java/jar


5) Configuring again:


toolbox/scilab/scilab$ ./configure
toolbox/scilab/scilab$ make
toolbox/scilab/scilab$ bin/scilab


Well, we now have our scilab base.
We have to set all the toolboxes files, but that will be tonight.

Wednesday, May 13, 2009

Lock free programming

http://en.wikipedia.org/wiki/Lock-free_and_wait-free_algorithms
http://en.wikipedia.org/wiki/Compare-and-swap

http://www.google.com/search?hl=en&client=iceweasel-a&rls=org.debian%3Aen-US%3Aunofficial&q=lock+free+&btnG=Search

Jack en gumstix.

MMM, necesito un framework de simulacion por bloques que este ligado la kernel, se me vino a la cabeza jack y su forma de pasar datos entre plugins, y simplemente vi esto en unos slides:

Although not a requirement, Jack should support any streaming data type, not just audio.

Podria ser una buena opcion meter jack en el gumstix y crear los bloques de control como si fueran plugins y el paso de datos serian rapidito, el llamaod a los codigos de los bloques. Suena interesante la cosa. Veremos que sale.

Tuesday, April 28, 2009

Drivers en linux (II)

Ahora crearemos un dispositivo de caracteres en /dev/algo, para esto debemos definir cuatro nuevas funciones: para abrir y cerrar el dispositivo, y para leer y escribirlo. Hay mas operaciones, pero con esas es suficiente por ahora (http://www.linuxtopia.org/online_books/linux_kernel/linux_kernel_module_programming_2.6/x569.html).

El archivo test.h queda asi:

#ifndef _TEST_SPI_H_
#define _TEST_SPI_H_

#define DEVICE_NAME "test-spi"

MODULE_AUTHOR("Jorge Eduardo Cardona <jorgeecardona@gmail.com>");
MODULE_DESCRIPTION("SPI test driver.");
MODULE_LICENSE("GPL");

// Character device functions.
static int this_open(struct inode *inode, struct file *file);
static int this_release(struct inode *inode, struct file *file);
static ssize_t this_read(struct file *file, char __user *buffer, size_t len, loff_t *offset);
static ssize_t this_write(struct file *file, const char __user *buffer, size_t len, loff_t *offset);

// File Operations structure.
struct file_operations this_fops = {
.open = this_open,
.release = this_release,
.read = this_read,
.write = this_write
};

// Linux module functions.
static int __init this_init(void);
static void __exit this_exit(void);
module_init(this_init);
module_exit(this_exit);

#endif



Para crear el dispositivo de caracteres debemos llamar a la funcion: register_chrdev (http://kerneltrap.org/man/linux/man9/register_chrdev.9)(http://www.silicontao.com/ProgrammingGuide/GNU_function_list/register_chrdev.html)

#include <linux/module.h>
#include <linux/fs.h>
#include <linux/err.h>


#include "test.h"

// Global data
static int this_major;


// Function that is called every time we try to open a device.
static int this_open(struct inode *inode, struct file *file)
{
printk(DEVICE_NAME": Open device.\n");
return 0;
}

// Function that is called when we close a device.
static int this_release(struct inode *inode, struct file *file)
{
printk(DEVICE_NAME": Close device.\n");
return 0;
}

// Function that is called when we read from a device.
static ssize_t this_read(struct file *file, char __user *buffer, size_t len, loff_t *offset)
{
printk(DEVICE_NAME": Read device.\n");
return 0;
}

// Function that is called when we write to a device.
static ssize_t this_write(struct file *file, const char __user *buffer, size_t len, loff_t *offset)
{
printk(DEVICE_NAME": Write device.\n");
return 0;
}

// Function that is called when we load the module.
static int __init this_init(void)
{
int ret;

printk(DEVICE_NAME ": Init.\n");

// Register char device: 0 for dynamically allocation of the major number.
ret = register_chrdev(0, DEVICE_NAME, &this_fops);
if (~IS_ERR_VALUE(ret))
{
printk(DEVICE_NAME ": Character device registered, with major=%d.\n",ret);
this_major = ret;

return 0;
}

unregister_chrdev(this_major, DEVICE_NAME);
printk(DEVICE_NAME ": Character device unregistered.\n");

return ret;
}

// Function that is called when we remove the module.
static void __exit this_exit(void)
{
unregister_chrdev(this_major, DEVICE_NAME);
printk(DEVICE_NAME ": Character device unregistered.\n");
printk(DEVICE_NAME ": Exit.\n");
}

(http://tomoyo.sourceforge.jp/cgi-bin/lxr/ident?i=IS_ERR_VALUE)

Bueno, igual que antes, configuramos, compilamos y tenemos:

root@gumstix-custom-connex:~$ modprobe test
root@gumstix-custom-connex:~$ dmesg | grep "test-spi"
<4>test-spi: Init.
<4>test-spi: Character device registered, with major=254.
root@gumstix-custom-connex:~$ mknod /dev/test c 254 0
root@gumstix-custom-connex:~$ cat /dev/test
root@gumstix-custom-connex:~$ dd bs=1 count=0 if=/dev/zero of=/dev/test
0+0 records in
0+0 records out
root@gumstix-custom-connex:~$ echo "ja" > /dev/test

root@gumstix-custom-connex:~$ rmmod test
root@gumstix-custom-connex:~$ dmesg | grep "test-spi"
<4>test-spi: Init.
<4>test-spi: Character device registered, with major=254.
<4>test-spi: Open device.
<4>test-spi: Read device.
<4>test-spi: Close device.
<4>test-spi: Open device.
<4>test-spi: Close device.
<4>test-spi: Open device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Write device.
<4>test-spi: Close device.
<4>test-spi: Character device unregistered.
<4>test-spi: Exit.
root@gumstix-custom-connex:~$


Los varios mensajes de Write se deben a que retornamos 0 esto quiere decir que no escribimos ningun caracter y la funcion es llamada hasta que vacie el buffer, pero aqui toca matarla con un Ctrl-C.

Bueno, un paso mas.

Drivers en linux (I)

Bueno iba a hacer esto en un solo post pero por partes es mejor, primero lo básico, cargar y bajar el modulo. En mi caso necesito un modulo para SPI, entonces lo ubico en drivers/spi/test.c

Comenzando con lo basico, entrar y bajar el modulo.

test.c

#include <linux/module.h>
#include "test.h"

static int __init this_init(void)
{
// Initialize module.

printk(DEVICE_NAME ": Init.\n");
}

static void __exit this_exit(void)
{
// Exit module.

printk(DEVICE_NAME ": Exit.\n");
}


test.h

#ifndef _TEST_SPI_H_
#define _TEST_SPI_H_

#define DEVICE_NAME "test-spi"

MODULE_AUTHOR("Jorge Eduardo Cardona <jorgeecardona@gmail.com>");
MODULE_DESCRIPTION("SPI test driver.");
MODULE_LICENSE("GPL");

// Linux module functions.
static int __init this_init(void);
static void __exit this_exit(void);
module_init(this_init);
module_exit(this_exit);

#endif


Tenemos que manipular dos archivos mas: Kconfig y Makefile, ambos en la misma carpeta:

Kconfig

.
.
.
config SPI_TEST
tristate "Test SPI driver."
depends on SPI_MASTER && EXPERIMENTAL
help
Test SPI driver.
.
.
.


Makefile

.
.
.
obj-$(CONFIG_SPI_TEST) += test.o
.
.
.


Al configurar el kernel nos saldra la opcion, la seleccionamos y podremos hacer esto:

root@gumstix-custom-connex:~$ modprobe test
root@gumstix-custom-connex:~$ dmesg | tail -n1
<4>; test-spi: Init.
root@gumstix-custom-connex:~$ rmmod test
root@gumstix-custom-connex:~$ dmesg | tail -n1
<4> test-spi: Exit.
root@gumstix-custom-connex:~$


Primer paso listo.

Wednesday, April 15, 2009

Contando racionales.

Un numero racional es todo aquel que puede ser escrito de la forma:

racional = a/b : siendo a y b numeros enteros.

Usando el metodo como Georg Ferdinand Ludwig Philipp Cantor demostro que los racionales son contables podemos crear un pequeñito programa para sacar el i-esimo racional de la tabla de cantor:


from math import sqrt

def rational(i):
b = int((sqrt(1+8*i)-1)/2)
n = i - (b**2 + b)/2 + 1
d = b + 2 - n

return (n,d)


for i in xrange(100):
(n,d) = rational(i)

print "(%d) %d/%d \t= %f"%(i,n,d,float(n)/float(d))



Bueno, La tabla de Cantor es sencilla, el numerador aumenta columna a columna desde 1 hasta aleph_0 y el denumenador fila a fila de 1 a aleph_0, de esta forma todos los nuemros racionales son cubiertos por la tabla, y por tanto son contables, se les puede asociar uno a uno un elemento de los enteros.

La tabla:

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


Y el conteo:

0 2 5 9 14 20 27 35 44 ...
1 4 8 13 19 26 34 43 53 ...
3 7 12 18 25 33 42 52 ...
6 11 17 24 32 41 51 ...
10 16 23 31 40 50 ...
15 22 30 39 49 ...
21 29 38 48 ...
28 37 47 ...
36 46 ...
45 ...
...


Resultado de la ejecucion del codigo:

(0) 1/1 = 1.000000
(1) 1/2 = 0.500000
(2) 2/1 = 2.000000
(3) 1/3 = 0.333333
(4) 2/2 = 1.000000
(5) 3/1 = 3.000000
(6) 1/4 = 0.250000
(7) 2/3 = 0.666667
(8) 3/2 = 1.500000
(9) 4/1 = 4.000000
(10) 1/5 = 0.200000
(11) 2/4 = 0.500000
(12) 3/3 = 1.000000
(13) 4/2 = 2.000000
(14) 5/1 = 5.000000
(15) 1/6 = 0.166667
(16) 2/5 = 0.400000
(17) 3/4 = 0.750000
(18) 4/3 = 1.333333
(19) 5/2 = 2.500000
(20) 6/1 = 6.000000
(21) 1/7 = 0.142857
(22) 2/6 = 0.333333
(23) 3/5 = 0.600000
(24) 4/4 = 1.000000
(25) 5/3 = 1.666667
(26) 6/2 = 3.000000
(27) 7/1 = 7.000000
(28) 1/8 = 0.125000
(29) 2/7 = 0.285714
(30) 3/6 = 0.500000
(31) 4/5 = 0.800000
(32) 5/4 = 1.250000
(33) 6/3 = 2.000000
(34) 7/2 = 3.500000
(35) 8/1 = 8.000000
(36) 1/9 = 0.111111
(37) 2/8 = 0.250000
(38) 3/7 = 0.428571
(39) 4/6 = 0.666667
(40) 5/5 = 1.000000
(41) 6/4 = 1.500000
(42) 7/3 = 2.333333
(43) 8/2 = 4.000000
(44) 9/1 = 9.000000
(45) 1/10 = 0.100000
(46) 2/9 = 0.222222
(47) 3/8 = 0.375000
(48) 4/7 = 0.571429
(49) 5/6 = 0.833333
(50) 6/5 = 1.200000
(51) 7/4 = 1.750000
(52) 8/3 = 2.666667
(53) 9/2 = 4.500000
(54) 10/1 = 10.000000
(55) 1/11 = 0.090909
(56) 2/10 = 0.200000
(57) 3/9 = 0.333333
(58) 4/8 = 0.500000
(59) 5/7 = 0.714286
(60) 6/6 = 1.000000
(61) 7/5 = 1.400000
(62) 8/4 = 2.000000
(63) 9/3 = 3.000000
(64) 10/2 = 5.000000
(65) 11/1 = 11.000000
(66) 1/12 = 0.083333
(67) 2/11 = 0.181818
(68) 3/10 = 0.300000
(69) 4/9 = 0.444444
(70) 5/8 = 0.625000
(71) 6/7 = 0.857143
(72) 7/6 = 1.166667
(73) 8/5 = 1.600000
(74) 9/4 = 2.250000
(75) 10/3 = 3.333333
(76) 11/2 = 5.500000
(77) 12/1 = 12.000000
(78) 1/13 = 0.076923
(79) 2/12 = 0.166667
(80) 3/11 = 0.272727
(81) 4/10 = 0.400000
(82) 5/9 = 0.555556
(83) 6/8 = 0.750000
(84) 7/7 = 1.000000
(85) 8/6 = 1.333333
(86) 9/5 = 1.800000
(87) 10/4 = 2.500000
(88) 11/3 = 3.666667
(89) 12/2 = 6.000000
(90) 13/1 = 13.000000
(91) 1/14 = 0.071429
(92) 2/13 = 0.153846
(93) 3/12 = 0.250000
(94) 4/11 = 0.363636
(95) 5/10 = 0.500000
(96) 6/9 = 0.666667
(97) 7/8 = 0.875000
(98) 8/7 = 1.142857
(99) 9/6 = 1.500000