21/10/2014

Graph library design choices

I am currently writing the documentation for my library. The effort of structuring and explaining the underlying logic of it is definitely a great help in making apparent what the library lacks and what could be done to make it better.

Unfortunately, time is limited, and I will never have the time to make the perfect Java graph library. Here are nevertheless some tips and facts to keep in mind in order to do so.


Graphs are about:
- relation between objects

And possibly:
- information about objects
- information about relation between objects

yes, but:

What?

i.e. which information?

* about the relations:
- directed or undirected
- normal graph or multigraph
- normal graph or hypergraph
- both directed and directed edges?

* ad-hoc information:
- weighted? (e.g. for min paths)
- weighted with several weights? (e.g. for max flows)
- with spatial information?
- with any other information?

Now, where, when and how to store these information?

Where?

* within the objects: in an information field of the objects

n.name = "toto";
n.name;

+: do not use hashmap
-: increases the size of a node

* outside of the graph: in a map object-> information that stores information separately from the objects and the graph

names.put(n,"toto");
names.get(n);

+: the algorithms that modify information sets leaves the graph untouched (the graph does not contains state information)
-: uses hashtables (heavy, centralised)

* within the graph but outside the objects: the same map, but stored at the level of the graph

g.names.put(n,"toto");
g.names.get(n);

+: the info about the graph is stored within the graph (easy to access)
-: the graph contains state information
-: uses hashtable

When?

* If the info associated with a node is frequently accessed it is better to store this info in a dedicated field.

+: fast
-: need to implement ad-hoc classes

* If several algorithms are going to be used on a given graph, it is therefore not possible or desirable to modify the graph, i.e. to change its state. As such it is better to store the information relevant only to an algorithm outside of the graph.

+: respect functional decomposition
+: easy to implement
-: state info is sometimes interesting or necessary, solution: return it with the standard return in a tupple

How?

* use ad-hoc node object to contain the information (information directly into the node)

+: all the info there, plus the neighbours in one place
-: need to implement specific nodes for each new application

e.g.

class CityNode {
   String name;
   int population;
   List<CityNode> neighbours;
}


CityNode n;


* use generic node objects to contains the ad-hoc information (information pointed from the node)

+: generic, does not need to develop new nodes
-: the ad-hoc information and neighbour information are separated, from the former it is not possible to access the later

e.g.

class GenericNode<K> {
   K data;
   List<GenericNode<K>> neighbours;
}

class CityNodeData {
   String name;
   int population;
}

GenericNode<CityNodeData> n;

* use no node objects only identifiers and separated maps to contain the information (information independent from the nodes)

+: totally respect function decomposition
-: easier to store data in one place sometimes
-: edge information need two levels of indirection with hashmas, can be cumbersome

e.g.

class Node {
   List neighbours;
}

Map<Node, String>  mapNodeCityName;
Map<Node, Integer> mapNodeCityPopulation;

Node n;

BUT

adjacency information is node information... and the same questions applies: shall adjacency:
- be encoded within the nodes with an adjacency list field: requiring the creation of additional objects
- be encoded in map node->adjacency list : node object are not even needed in this case

And shall lists be encoded with linked lists or arrays lists?
And shall maps be encoded with hashmaps or carefully managed within arrays?

Usage of arrays is extremely compact, lightning fast for iteration thanks to caching but slower for modification.
Usage of linked list and hashmap is not compact (require to save pointer information) and is slower for iteration, but faster for modification.

two distinct use cases appear:
- high performance graph library: like grph which will take the pain of optimising in the minute details to the point of coding in C in order to save time. In order to go further, data should be arranged to match the cache access pattern within each implemented algorithm.
- user and developer friendly graph library: like any other Java library, provide slower than C highly reusable Java code based on Java concepts. That's the point of using Java. (and it is strongly typed and faster than Python: that's the point of not using Python)

It is clear that graphs that are static once loaded could benefit from a speedup over dynamic one if this is known in advance. However as these encompass two different style of coding algorithms it is not as transparent as using two different implementation for on abstract data type.

* Ease of development

The fact is that it is much easier to implement only a bare relation structure and additional information encoded outside of it in maps than implementing the ad-hoc graphical structure that precisely fit one's purpose.

e.g. for a graph of cities with two edge information, three node information including spatial ones, it is easier to specify:

Map mapNames = new Hashmap()
mapNames.put(name1)
etc...

Graph g = new Graph(g, {mapNames, mapPopulation, mapCoordinates}, {mapDistance, mapCapacity})  // not correct Java code btw, but let your imagination fly to understand

Than:

AdHocGraph g = new AdHocGraph()
g.addNode(name1, population1, coordinates1)
etc...

because this second solution require the developement of an AdHocGraph with AdHocNode able to hold precisely the information of a name, a population and a coordinate

The solution, in order to make development easier, is to be more generic and provide a class of node able to all contain an information of the same parametrized type, e.g. City, where a City object contains the precise three information we are interested in.

So the question is:
Information Within Nodes : for static information of the graph
Information Within Maps : for ad-hoc algorithms

* Sequential or parallel

But you know, we are entering the age of distributed and parallel algorithms... What we just discussed is fine for centralized and sequential processing... What about the future?

Distributed algorithms are piece of information spread over a network (...a graph...) which interact by exchanging information. If there is not central point of control, the algorithm is more than distributed, it is decentralised.

To be perfectly clear about the meanings:
- sequential: fixed workflow + 1 thread
- parallel: fixed workflow + N threads (same computer)
- distributed: fixed workflow + N threads (different computer) + communication delays + faults
- dencentralised: distributed with several independent workflows

Not all algorithms can be distributed, but also not all algorithms even need to be distributed, what we would be more interested into is harnessing the power of parallel computation than coping with the full set of difficulties of distributed computing (server faults, communication faults) and methods (data replication, process replication, consensus algorithms).

Sometimes an algorithm can not be fully distributed, but can be partially parallelized, with some part of the problem being computed independently, and possibly combined. For instance any matrix summation can be fully parallelized, with worker jobs working on independent parts of the problem and combining their solution by writing them in the result matrix. Matrix multiplication can be only partially parallelised, using Strassen algorithm for instance, because the computation over independent parts of the problem require additional computation to be combined, therefore introducing sequential dependencies over the computation that needs to be performed, and that can't therefore be all executed in parallel.

When using a decentralised algorithm there is only one choice for the location of data: inside the nodes of the graph. However for distribution and parallelisation, such is not the case. As in the sequential case, data can be stored inside the graph at the level of nodes or outside of the graph in a separate map. In distribution, storing the data makes it simpler to transmit info to processors: they just process the nodes, while storing the data in maps requires to first extract the relevant parts and transmit them along the graph to the processors. In parallelisation, all the processors can access the same memory making the split irrelevant.

The catch here is that what parts of a computation on a graph can be parallelised are determined by the graph structure itself, but all the code that parallelise, compute and aggregate piece of solution has nothing to do with the design of a graph library, still we need to take the distribution at  the core of the algorithms.

Let us consider different kinds of distribution in the problem of computing shortest paths with dijkstra:
* - node level: each node execute a local dijkstra algorithms that updates its information given the message he receives and updates other nodes. The heap, which is a global data must be managed throught a consensus algorithm. Such a solution is slow due to the communication overhead and should be used when strictly necessary
* - group of nodes level: each processor is responsible for the computation of a subparts of the graph: this can provide significant speedup if communication costs or memory contention are low. Note that this woudln't be the exact dijkstra algorithms as separating the graph and assembling parts of the solution require additional work to be done in a principled way.
* - graph level: no distribution, no speedup.






24/09/2014

Functional programming or Efficient programming?


The expressivity of functional programming is awesome: in a few lines, you can write complex functions. Plus it is much more secure than imperative programming. However this comes at a price, a performance price: everything has to be the result of a function call.

This means that if you want to compute both a shortest path between two node AND the cost of this path, you should theoretically call two different functions, that would do almost the exact same job. This is rudely inefficient and inelegant.

To overcome this issue, functional programming allows you to return tupples. The only catch with this approach is that the name of the function is less adequate: to be correct in the shortest path example it would be compute-shortest-path-and-associated-cost.

Imperative programming however does not elegantly solve this problem: to return multiple argument in C you have to use the arguments of the function - which is as simple and efficient that it is awful. On a side note in a logical programming language, i.e. Prolog, this is the only way possible to return arguments - no wonder I fucking hated it back at the uni.


In order to make my Java code better understandable, elegant and secure, I inserted in my graphical framework several usual features of functional programming language, absent in java. Notably: tupples and lambda functions. As a consequence I tend to use iterators more.

use of tupples

Tupples are used to return several argument from one function call, for instance again, a shortest path request. The alternatives would be to:
- return in the arguments (unthinkable: against code readability and clarity);
- store the relevant data (e.g. node parent lists and node costs) of the first function call in a global data structure which can be access by another function call (even more unthinkable: would be a total mess);
- returning a dedicated result structure for each function (unthinkable but a bit less: will bloat the framework with barely used classes);
- returning a hashtable containing the different information accessed by different keys (why not?).

A templated tupple structure is therefore a good compromise solution between the last two so-so solutions, it is a data structure that can be used by several functions, and the type information make it dedicated to the return data of this function.
However, because template have a precise length, it is necessary to create Tupple classes for each number of argument. We will have therefore Tupple2<K1, K2> for tupple of two objects, Tupple3<K1, K2, K3>, for tupples of three object, etc. In my implementation I stopped at three.
The only issue with it is that templated code is very verbose, but this is a problem of Java, not of this approach to tackling the problem.

use of lambda functions

The purpose of the support of lambda function is not to be fun to play around with when you don't know them but to provide code bricks that can be easily combined to perform complex and secure (and slow) operations.
In the framework lambda function have been used extensively to build the math.algebra section, whose most used member is the Semiring class.

This class is actually the very reason why lambda function where introduced: because the operation perform to optimize constraints, compute shortest path, make probabilistic inference, compute reachability all share the exact same mathematical structure, the one of a semiring, they just differ in the valuation and combination they use.
Therefore by providing a generic Semiring object we save a lot of coding by factoring it all in a few key functions and instantiations, this is exactly what I did in the framework.

The drawback of all this theoretical power is pragmatical slowness. I benchmarked simple usage of lambda functions (map/reduce + list/iterator + ad hoc lambda function)  against the same functionality written with a proper loop (for + list/iterator + ad hoc loop code): lambda function execute about 3 times slower than the standard loop. This is why when I latter read that Scala was 3 times slower than Java I was everything but surprised: I knew exactly why.

the choice: lambda or not

I first knew about lambda calculus after I finished my studies. I was unofficially following a class about "mathematical foundations of logic" at Paris 7 university. This class blew my mind: we had an awesome teacher (Gille Dowek) for an awesome topic. (On a side note again, a very good profane read about this very topic is the Logicomix book).
Since then I wanted to study it further, and this meant implement it and play with it. This is what I did, and by thinking intensely in a functional way it gave me an even deeper understanding of what was going on in code, making structure apparent.

While the abstraction is powerful to handle an combine different pieces of code, it is performance-wise quite catastrophic. This is why I removed almost all usage of map/reduce loops in my code and replace them by straight loops. It was a bit of a hearthbreaker at the beginning, but I knew I was doing the right thing.
Now it is used in the code of the framework only in two place: one as a souvenir where it won't hinder performances, and one where the power of the abstraction is more interesting than speedup of regular loops as it positively strongly impact on code reusability, readability and quality.

In the end of the day, it is always setting the trade-off among what matters to you, but you need to have it clear first.








03/08/2014

Aventures grammaticales : le mystère de l'érgativité

Rappel: l'érgativité c'est une propriété d'une langue de marquer de la meme manière le sujet d'un verbe intransitif et l'objet d'un verbe transitif. C'est simple à dire, très complexe à comprendre, c'est pourquoi c'est fascinant.

Pour mieux fixer les idées, prenons un exemple simpliste en marquant A le sujet d'un verbe transitif et B le reste:

le chat-B mange
le chat-A mange la souris-B

Il existe différente manière de réaliser

Il existe plusieurs langues naturelles qui possèdent cette propriété à différents degrés. Qui réalisent ce marquage de différentes facon, le plus souvent pas l'usage de cas: changement de la terminaison d'un mot, dont un cas peut ou non etre marqué ou les deux. Mais cela peut également etre marqué par la position par rapport au verbe.

Bref, la question que j'aborde dans ce post est: qu'est ce que l'érgativité dans les languages informatique peut nous apporter sur le compréhension de cette notion.

Et plus particulièrement: quels seraient les equivalents en langage de programmation de la notion d'érgativité dans les langues naturelles? Et qu'est ce que cela peut nous apprendre sur celles-ci?

Revoyons d'abord un peut plus en détails ce qui se passe dans les langues naturelles:

* quand le verbe se passe par défaut d'acteur, le complément précise celui qui controlle le processus

le plat cuit
le chef cuit le plat

ou il peut s'exprimer par tournure reflexive:

le plat se cuit

ou il peut s'exprimer de manière factitive

le chef fait cuire le plat

* quand le verbe se passe par défaut d'objet, le complément précise l'objet précise qui subit le processus

le chef cuisine
le chef cuisine le plat

la tournure reflexive ne marche pas:

le chef se cuisine (le sens est différent: le chef ne cuisine pas un plat, mais lui meme)

la tournure factitive ne marche pas:

le chef fait cuisiner le plat (le sens est different: ce n'est pas le chef qui cuisine, mais qqun d'autre)

Donc le "truc" de l'érgativité n'as pas avoir avec les marqueurs  des cas, mais avec la nature des verbes.


Voyons à présent quels pourraient en etre les équivalent dans les langages informatiques:

= le verbe cuire

* en version predicative (logique), impérative (programmation) :

le numero ou le nom de la position sert à marquer le cas

cuire(quoi="plat") : predicat
cuire(quoi="plat", qui="chef") : le plus naturel, l'argument optionel est en deuxieme position
cuire(qui="chef", quoi="plat")

* en version objet:

. et () servent ici à marquer les cas

lequel est correct est une question d'interpretation:

- correct si avant le . est indiqué l'acteur:

chef.cuire() : le chef fait cuire qqchose  
chef.cuire(plat)

- correct si avant le . est ce qui subit le processus de cuison

plat.cuire() : le plat cuit 
plat.cuire(chef)

Une vrai transposition de la grammaire langue nat // marqueurs serait:

Langue nomitaive-accusative:

chef.cuit()
chef.cuit(plat)

[ ou, de manière extremement inhabituelle mais pas de manière inconcevable:

(plat)cuire
(plat)cuire.chef ]

langue érgative-absolutive:

plat.cuit()
chef.cuit(plat)

ou

(plat)cuire
(chef)cuire.plat

autre représentation:

cuire(x="plat", y="chef") : predicatif
y="chef".cuire(x="plat")  : objet nominatif-accusatif
y="plat".cuire(x="chef")  : objet ergatif-absolutif

= le verbe cuisiner

cuisine(chef) : prédicat
cuisine(chef, plat) : le plus naturel, l'argument optionel est en deuxième position
cuisine(plat, chef) : moins naturel, les arguments se rajoutent en file
chef.cuisine(plat) : objet

lequel est correct est une question d'interpretation:

chef.cuisine() : le chef cuisine
plat.cuisine() : le plat est cuisiné

"chef".cuisine(?) : le chef cuisine qqchose
?.cuisine("plat") : qqun cuisine le plat

Donc les acteurs et objets du verbe, comment on les marque en langue naturelle et comment on les exprime en langue logique ou informatique, dependent totallement du sens fondamental de ce verbe, pas que de sa transitivité, mais plutot du fait que celui-ci indique un acteur ou subisseur naturel, un changement d'état ou non.

C'est comme si certain verbe exprimaient natrellement un acteur ou un objet, et que donc naturellement le complement le plus direct du sens à ce verbe, sera d'exprimer respectivement l'objet ou l'acteur manquant par défaut.

x cuit => (Acteur=?, Objet=x)
x cuisine => (Acteur=x, Objet=?)

x cuit y => (Acteur=x, Objet=y)
x cuisine y => (Acteur=x, Objet=y)

de ce point de vue la on peut dire qu'un verbe ergatif dans une langue marquant le cas par position par rapport au verbe, inverse la position entre les cas intransitif/transitif, tandis qu'un verbe non ergatif (noté ici "nominatif") ne l'inverse pas.

noté de manière prédicative cela serait équivalent à:

predicat_ergatif(x)
predicat_ergatif(y, x)

predicat_nominatif(x)
predicat_nominatif(x, y)


* en mathématiques

"-" est un opérateur unaire et binaire, notons le "minus"

minus(x) : -x
minus(x,y) : x - y

ici le role de l'argument x est inversé dans la version binaire par rapport à la version unaire: en unaire, x est retiré à un 0 qui n'est pas représenté, tandis qu'en binaire on retire de x. Cela ressemble à un construction ergative du predicat "minus".






25/07/2014

The Difficulties of Designing a Graph Library


I am currently building a graph library for the language Java, and as such I am faced with several tantalizing design questions, to which I'll try to give an overview here and the response I came up with.

Since my Master thesis till my PhD I have been working with my very own libraries to support my inference algorithms, which I kept expending, improving and redesigning, to a point I believe is ripe for diffusion and use.

I here review over some of the difficulties and choices one is faced in order to answer this fundamental and intricate question:

How to represent a graph?

while being:
- computationally efficient
- easily reusable by the end-users
- easily reusable by the developers

Here I am not talking about choosing a data structure to represent a graph ADT between an adjacency list graph, and incidence edge graph or an edge list graph.

Well, not only.

I mean how to design any implementation such that the three design constraints hold.

First and foremost: what is a graph?

Short and easy answer: It is a set of relations between objects.

Second: how to represent these relations?

We can have:
- for each object associated a set of all the related objects, this is an adjacency-based representation
- for each relation associated the set of objects within this relation, this is an incidence-based representation

Third: how to implement these representations?

adjacency-based representation:
- a list of objects for each objects: adjacency list
- a square matrix indexed by the objects: adjacency matrix

incidence-based representation:
- a list of objects for each relation: edge list
- a rectangular matrix indexed by the objects and the relations: incidence matrix

Fourth: where to store these information?

This is the actually the most fundamental question, because the whole design of the library will revolve around that. Where to store the information depends on how the information is used.

We could:
- store the list at the level of the objects: one list for each objects, nothing in the graph structure
- store the list at the level of the graph: nothing in the objects, one map object->list<object> in the graph structure

Fifth: how do we choose among all these options?

Implementation-wise, regarding representations and implementation, we don't: we should provide each representation, and let the users choose which one he wants to use according to its needs.

Usage-wise it's up to the performance for a particular usage, or the way a particular algorithm is expressed.

For instance Kruskal's algorithm for maximal spanning tree is based on edge list, Dijkstra's algorithm for shortest paths is based on adjacency list and the Floyd-Warshall algorithm for shortest paths is based on adjacency matrix.

05/07/2014

Fan control on ASUS UX32VD under linux

Out of the box, the fan is either full speed or off. This is very annoying, but there exists several practical solutions. I here report my experience and how you can safely and easily fix the issue.

My configuration is an ASUS ux32vd with first ubuntu 13.10 and now 14.04.

In a nutshell

the problem disapear after running the provided program for a few days.

1) Introduction

The most comprehensive thread on the net about the issue is a very interesting read for technical details, but that you may want to skip in order to just get a working solution. In a nutshell it is possible to control the fan speed by writing in the memory of the embedded controller, however not doing it properly through the ACPI interface is hazardous.

Some integrated work around have been developed: among the most complete I found on the internet are one solution based on a kernel module, and here is another solution based on a shell script and acpi_call. They unfortunately lack documentation which can be an obstacle for the layman. I tried, used and modified the last one by Emil Lind, which is based on a shell script and the acpi_call module. I here report more information about how to properly use it, as well as my own modification.

2) What you need to know

First it has to be noted that contrary to what has been written elsewhere, what is controlled by writing to the embedded controller memory through ACPI calls is not the actual speed of the fan, but the maximal speed of the fan.

Secondly, the danger of running these programs is of not exiting them in a safe state, i.e. not to restore the automatic control of the fan and a maximal fan speed when they exit, which could theoretically bring the laptop to overheat.

If you have this in mind an ascertain that each time you use it it leaves properly you can safely use this program.

Thirdly, it has to be noted that after experimenting with this program, the problem disappeared: the fan where able to react proportionally to heat without any ad-hoc fan control loop in the background. So you can expect experimenting the same after a while.

3) How you can fix the fan

I here provide a heavily modified version of the shell program asusfanctrld of Emil Lind: it has a control loop which automatically switches between a set of presets RPM temperature thresholds. The thresholds have been determined empirically in order to minimise the switching between states and provide and overall unnoticeable and cool experience ;) The real fix is the following: for an unknown reason after a few days of use of the following program, the problem disapear by itself. It may reappear if you boot windows. In which case, running again the provided program will fix the issue.

3.1) In more details:
there are four predefined temperature thresholds and four predefined fan RPM values. When a threshold is hit from the upper of lower values, the fan speed is set to the corresponding RPM until another threshold, lower or higher, is reached. ACPI commands are written through the special device /proc/acpi/call provided by the acpi_call module kernel module that must be installed.

3.2) Install:
- step 1 : the acpi_call module has to be downloaded, compiled and loaded (from a root shell: insmod acpi_call.ko) prior to running the fan control loop, which is the most complicated task to be done to get a silent zenbook.
- step 2 : the following shell script as to be copied/pasted to a file, let us call it asusfanctrlnico.sh

3.3) Usage:
The script has to be run as root, leave it running in a windows to check everything runs smoothly until you either get confident enough in it or the fan problems goes away. In order to allow the hard drive to sleep you need to comment out the file logging commands in the main function.

It is not meant to be used by startup scripts. It however can with the appropriate kill commands added to the init.d or start-stop-daemon services (which I didn't as the problem went away after using the program for some days).

3.4) The code of asusfanctrlnico.sh :


ACPI_CALL=/proc/acpi/call

# average temperature without control:
# sleep ~ 27 - 30
# on, idle ~ 45 - 50
# internet browsing, text editing ~ 50 - 58
# temperature threshold (Celsius)

HIGH=70
MID=60
LOW=55
LOWEST=50
MARGIN=1

# corresponding fan RPM
LOWESTRPM=0x35
LOWRPM=0x50
MIDRPM=0x70
HIGHRPM=0xFF

DEBUGACPI=0
DEBUGBEHV=1

ME="$(basename $0)"

LOGFILE=asusfanlog
LOGFILETEMP=fantemplog

fatal() {
    logger -id -t $ME -s "FATAL: $@"
}
info() {
    logger -id -t $ME "INFO: $@  `date +"%H %M %S %s"`"
    echo "INFO $@  `date +"%H %M %S %s"`" >> $LOGFILE
    if [ "$DEBUGBEHV" -gt 0 ]; then
        echo ""
    fi
    echo -en "INFO $@  `date +"%H %M %S %s"`\n"
}
debug_acpi() {
    if [ $DEBUGACPI -gt 0 ]; then
        echo "DEBUG: $@" >> $LOGFILE
    fi
}
debug_behaviour() {
    if [ $DEBUGBEHV -gt 0 ]; then
        echo "DEBUG: $@" >> $LOGFILE
    fi
}


if ! [ -e $ACPI_CALL ]; then
    modprobe acpi_call
fi
if ! [ -e $ACPI_CALL ]; then
    fatal "You must have acpi_call kernel module loaded."
    fatal modprobe acpi_call
    exit 1
fi

if [ "$(id -u)" != "0" ]; then
    echo "Must be run as root."
    exit 1
fi

set_fanspeed() {
    RPM=$1
    command='\_SB.PCI0.LPCB.EC0.ST98 '$1
    echo "$command" > "${ACPI_CALL}"
    debug_acpi $(cat "${ACPI_CALL}")
}

set_fanstate() {
    current=$STATE
    asked=$1
    if [ "$current" != "$asked" ]; then
        info "reset fan state: $current => $asked"
        if [ "$asked" = "high" ]; then
            set_fanspeed $HIGHRPM
        elif [ "$asked" = "mid" ]; then
            set_fanspeed $MIDRPM
        elif [ "$asked" = "low" ]; then
            set_fanspeed $LOWRPM
        elif [ "$asked" = "lowest" ]; then
            set_fanspeed $LOWESTRPM
 else
      fail "error unknown state: $asked"
 fi
 STATE=$asked
    fi
}

set_auto() {
    command='\_SB.ATKD.QMOD 0x02'
    echo "$command" > "${ACPI_CALL}"
    debug_acpi $(cat "${ACPI_CALL}")
    command='\_SB.PCI0.LPCB.EC0.SFNV 0 0'
    echo "$command" > "${ACPI_CALL}"
    debug_acpi $(cat "${ACPI_CALL}")
    CONFIG="auto"
    STATE="auto"
}

set_manual() {
 TEMP=$1
 if [ $TEMP -ge $HIGH -o \( "$STATE" = "high" -a $TEMP -ge $[$MID+$MARGIN] \) ]; then
  set_fanstate high
 elif [ $TEMP -ge $MID -o \( "$STATE" = "mid" -a $TEMP -ge $[$LOW+$MARGIN] \) ]; then
  set_fanstate mid
 elif [ $TEMP -ge $LOW -o \( "$STATE" = "low" -a $TEMP -ge $[$LOWEST+$MARGIN] \)  ]; then
  set_fanstate low
 else
  set_fanstate lowest
 fi
}

set_fanmode() {
 asked=$1
 if [ "$asked" = "auto" ]; then
  set_auto
  CONFIG="auto"
  STATE="auto"
  info "set fan mode: auto"
 elif [ "$asked" = "manual" ]; then
  CONFIG="manual"
  STATE="auto"
  info "set fan mode: manual"
 else
      fail "error unknown mode: $asked"
 fi
}
fail(){
 info "$@"
 set_fanmode auto
 exit 1
}

fatal_callback() {
 fail "stop, fatal signal"
}
sigusr1_callback() {
 set_fanmode manual
 logger -id -s -t $ME "got SIGUSR1, invoking manual control of fans"
 info "siguser1 mode : auto => manual"
}
sigusr2_callback() {
 set_fanmode auto
 logger -id -s -t $ME "got SIGUSR2, setting mode to auto"
 info "siguser1 mode : manual => auto"
}

main() {
    trap "fatal_callback" INT TERM EXIT KILL
    trap "sigusr1_callback" SIGUSR1
    trap "sigusr2_callback" SIGUSR2

    info start
    set_fanmode manual
    if [ $DEBUGBEHV -gt 0 ]; then
        echo ""
    fi
    echo "`date +"%s"` CONFIG $LOWEST $LOW $MID $HIGH  $LOWESTRPM $LOWRPM $MIDRPM $HIGHRPM" >> $LOGFILETEMP
    LASTTEMPSTATE=""
    while :; do
        TEMP0=$[$(</sys/devices/virtual/thermal/thermal_zone0/temp)/1000]
        TEMP1=$[$(</sys/devices/virtual/thermal/thermal_zone1/temp)/1000]
 TEMP=$TEMP1
 if [ $DEBUGBEHV -gt 0 ]; then
            echo -en "\r${TEMP0} ${TEMP1} mode=$CONFIG state=$STATE"
 fi
        if [ "$CONFIG" = "manual" ]; then
  set_manual $TEMP
 fi
 TEMPSTATE="$TEMP0 $TEMP1 $CONFIG $STATE $(cat /proc/loadavg | cut -d\  -f 1,2,3)"
 if [ "$TEMPSTATE" != "$LASTTEMPSTATE" ]; then
  echo "`date +"%s"` $TEMPSTATE" >> $LOGFILETEMP
  LASTTEMPSTATE=$TEMPSTATE
 fi
        sleep 1
    done
}

set -e
main 








04/07/2014

Goodness of fit test for log-normal and exponential distributions by bootstrapping


Stackexchange users be welcome to this page :), if you find the program useful or have suggestions, please feel free to leave a comment.


The poweRlaw R library provides the bootstrap_p function which allows to test the goodness of fit of a power law to the data using bootstrapping.  Unfortunately the library does not provides such methods for other distributions. We here propose such functions for log-normal and exponential models.
The code is based on the principle of the paper by Clauset et al. 2009 to test goodness of fit by bootstrapping for the log-normal and exponential distributions and reuse some of the powerLaw functions.

I wrote it after realizing that the code provided on the page http://users.dimi.uniud.it/~massimo.franceschet/R/fit.html meant for the same purpose was incorrect: the way the goodness of fit is computed is not correct as it does not strictly follow the procedure described in the paper, and implemented in the poweRlaw package.
Specifically: the synthetic data generation procedure is half implemented, it does not search for the best xmin as provided by the extimate_xmin function of the poweRlaw package, and finally the ks.test discards all the ties, while powerLaw doesn't with its built-in KS test.


Usage: a poweRlaw discrete distribution object has to be given as the first parameter. Prints the p-value on the output. Note that following the paper, about 2500 simulations should be generated to get a 2 points precision.

Interpretation of the p-value: a p-value < 0.1 rejects the hypothesis that the tested law generated the data, while a value > 0.1 does not reject it, but does not disprove the model either.

Note that model comparison based on the likelihood ratio (like the Vuong test, which is implemented in the compare_distribution function of poweRlaw) is better able to tell which of two models better fits the data and if there are enough data to confidently claim so.

Code for log-normal distributions:

bootstrap_ln = function(m,nbsims=5) {
 xmin = m$getXmin()
 ntot = length(m$dat)
 n1 = sum(m$dat<xmin)
 n2 = ntot - n1
 print(n1)
 print(n2)
 print(ntot)

 start_time = Sys.time()

 m_est = estimate_xmin(m)

 end_time = Sys.time()
 total_time = difftime(end_time, start_time, units="secs")

 mean = m_est$pars[1]
 sd = m_est$pars[2]

 print(m_est$KS)
 print(unlist(m_est))

 print(total_time)

 count=0
 for(i in 1:nbsims) {
  # generate synthetic data set
  if(xmin>1) {
   syn1 = sample(m$dat[m$dat<xmin],n1,replace=TRUE)
  } else {
   syn1 = c()
  }
  # sample from distribution
  syn2 = c()
  while(length(syn2)<n2) {
   n2bis = n2-length(syn2)
   syn2bis = round(rlnorm(n2bis,meanlog=mean,sdlog=sd))
   syn2bis = syn2bis[syn2bis>0]
   syn2 = c(syn2,syn2bis)
  }  
  syn = c(syn1,syn2)
  print(length(syn))
  # estimate gof of syn data  
  msyn = dislnorm$new(syn)
  syn_est = estimate_xmin(msyn)
  print(syn_est$KS)
  print(unlist(syn_est))
  if(syn_est$KS >= m_est$KS) {
   count = count + 1
  }
 }
 print(count/nbsims)
}

example:

dsyn = round(rlnorm(1000,2.4,0.98))
dsyn = dsyn[dsyn>0]
length(dsyn)
e_lnsyn = dislnorm(dsyn)
e_lnsyn$setXmin(estimate_xmin(e_lnsyn))
bootstrap_ln(e_lnsyn)


Code for exponential distributions:

bootstrap_exp = function(m,nbsims=5) {
 xmin = m$getXmin()
 ntot = length(m$dat)
xmin
 n1 = sum(m$dat<xmin)
 n2 = ntot - n1
 print(n1)
 print(n2)
 print(ntot)

 start_time = Sys.time()

 m_est = estimate_xmin(m)

 end_time = Sys.time()
 total_time = difftime(end_time, start_time, units="secs")

 rate = m_est$pars[1]

 print(m_est$KS)
 print(unlist(m_est))

 print(total_time)

 count=0
 for(i in 1:nbsims) {
  # generate synthetic data set
  if(xmin>1) {
   syn1 = sample(m$dat[m$dat < xmin],n1,replace=TRUE)
  } else {
   syn1 = c()
  }
  # sample from distribution
  syn2 = c()
  while(length(syn2)<n2) {
   n2bis = n2-length(syn2)
   syn2bis = round(rexp(n2bis,rate=rate))
   syn2bis = syn2bis[syn2bis>0]
   syn2 = c(syn2,syn2bis)
  }  
  syn = c(syn1,syn2)
  # estimate gof of syn data  
  msyn = disexp$new(syn)
  syn_est = estimate_xmin(msyn)
  print(syn_est$KS)
  print(unlist(syn_est))
  if(syn_est$KS >= m_est$KS) {
   count = count + 1
  }
 }
 print(count/nbsims)
}