package proxel;
import java.util.LinkedList;
import java.util.List;
class Proxel {
int id; /* unique proxel id for searching */
int s; /* discrete state of SPN */
int tau1k; /* first supplementary variable */
int tau2k; /* second supplementary variable */
double val; /* proxel probability */
Proxel left, right; /* reference to child proxels in tree */
}
public class MTProxel {
static double[][] y; /* vectors for storing solution */
static double tmax ; /* maximum simulation time */
static int TAUMAX; /*the largest possible aging time */
static int totcnt; /* counts total proxels processed */
static int maxccp; /* counts max # concurrent proxels */
static int ccpcnt ; /* counts concurrent proxels */
static Proxel[] root=new Proxel[2]; /* trees for organizing proxels */
static List<Proxel> firstfree=new LinkedList<Proxel>(); ; /* linked list of free proxels */
static double eerror = 0; /* accumulated error */
static int sw = 0; /* switch for old and new time steps */
static int len;
static double dt; /* DELTA*/
final static double ENDTIME=10;
final static double DELTA=0.005;
final static double MINPROB=1.0e-12;
final static int HPM= 0;
final static int LPM= 2;
//outcome of the machine production
static double defective;
static double working;
/********************************************************/
/* distribution functions */
/* instantaneous rate functions */
/********************************************************/
/* returns weibull IRF */
static double weibullhrf(double x, double alpha, double beta, double x0) {
double y;
y = beta/alpha*Math.pow((x-x0)/alpha,beta-1);
return y;
}
/* returns deterministic IRF */
static double dethrf(double x, double d) {
double y;
if (Math.abs(x - d) < dt/2)
y = 1.0/dt;
else
y = 0.0;
return y;
}
/* returns uniform IRF */
public static double unihrf(double x, double a, double b) {
double y;
if ((x >= a) && (x < b))
y = 1.0/(b-x);
else
y = 0.0;
return y;
}
/********************************************************/
/* output functions */
/********************************************************/
/* print all proxels in tree */
static void printtree(Proxel p) {
if (p == null)
return;
System.out.printf("s %d t1 %d t2 %d val %f \n",p.s,p.tau1k,p.tau2k,p.val);
printtree(p.left);
printtree(p.right);
}
/* print out complete solution */
static void plotsolution(int kmax) {
System.out.printf("\n\n");
for(int k=0; k<=kmax; k++)
{
System.out.printf("%7.5f \t %7.5e \t %7.5e \n", k*dt, y[0][k], y[2][k]);
}
}
/* print out a proxel */
static void printproxel(Proxel c) {
System.out.printf("processing %2d %2d %7.5e \n", c.s, c.tau1k, c.val);
}
/********************************************************/
/* proxel manipulation functions */
/********************************************************/
/* compute unique id from proxel state */
static int state2id(int s, int t1k, int t2k) {
return(TAUMAX*(TAUMAX*s+t1k)+t2k);
}
/* compute size of tree */
static int size(Proxel p) {
int sl, sr;
if (p == null)
return(0);
sl = size(p.left);
sr = size(p.right);
return(sl+sr+1);
}
/* get a fresh proxel and copy data into it */
static Proxel insertproxel(int s, int tau1k, int tau2k, double val) {
Proxel temp;
/* create new proxel or grab one from free list */
if (firstfree.isEmpty()== true)
temp = new Proxel();
else {
/*temp = firstfree;
firstfree = firstfree.right;*/
temp=((LinkedList<Proxel>) firstfree).removeFirst();
}
/* copy values */
temp.id = state2id(s, tau1k, tau2k);
temp.s = s;
temp.tau1k = tau1k;
temp.tau2k = tau2k;
temp.val = val;
ccpcnt += 1;
//to count the max number of the concurrent Proxel
if (maxccp < ccpcnt) {
maxccp = ccpcnt;
//printf("\n ccpcnt=%d",ccpcnt);
}
return temp;
}
/*returns a proxel from the tree*/
static Proxel getproxel()
{
Proxel temp;
Proxel old;
final int LEFT = 0, RIGHT = 1;
int dir=1; //to record it is right or left child of its parent
int cont = 1; //a mark to label if it has been already the leaf
if (root[1-sw] == null)
return null;
temp = root[1-sw];
old = temp;
/*move down the tree to a leaf */
while (cont == 1)
{
/* go right */
if ((temp.right != null) && (temp.left == null))
{
old = temp;
temp = temp.right;
dir = RIGHT;
}
/* go left */
else if ((temp.right == null) && (temp.left != null))
{
old = temp;
temp = temp.left;
dir = LEFT;
}
/* choose right/left at random */
else if ((temp.right != null) && (temp.left != null))
{
if (Math.random() >0.5) // 0.5 is a random bound here,RAND_MAX/2
{
old = temp;
temp = temp.left;
dir = LEFT;
}
else
{
old = temp;
temp = temp.right;
dir = RIGHT;
}
}
else
cont = 0;
} //end while
//remove the proxel which has been returned from the tree
if (temp == root[1-sw])
root[1-sw] = null;
else
{
if (dir == RIGHT)
old.right = null;
else
old.left = null;
}
/*c: old = firstfree;
firstfree = temp;
temp.right = old;
ccpcnt -= 1;*/
((LinkedList<Proxel>)firstfree).addFirst(temp);
return temp ;
}
/* adds a new proxel to the tree */
static void addproxel(int s, int tau1k, int tau2k, double val) {
Proxel temp, temp2;
int cont = 1,id;
/* Alarm! TAUMAX overstepped! */
if (tau1k >= TAUMAX) {
// printf(">>> %3d %3d %3d %7.5le \n", s, tau1k, val, TAUMAX);
tau1k = TAUMAX - 1;
}
/* New tree, add root */
if (root[sw] == null) {
root[sw] = insertproxel(s,tau1k, tau2k, val);
root[sw].left = null;
root[sw].right = null;
return;
}
/* compute id of new proxel */
id = state2id(s,tau1k, tau2k);
/* Locate insertion point in tree */
temp = root[sw];
while (cont == 1) {
if ((temp.left != null) && (id < temp.id))
temp = temp.left;
else
if ((temp.right != null) && (id > temp.id))
temp = temp.right;
else
cont = 0;
}
/* Insert left leaf into tree */
if ((temp.left == null) && (id < temp.id)) {
temp2 = insertproxel(s, tau1k,tau2k, val);
temp.left = temp2;
temp2.left = null;
temp2.right = null;
return;
}
/* Insert right leaf into tree */
if ((temp.right == null) && (id > temp.id)) {
temp2 = insertproxel(s, tau1k,tau2k, val);
temp.right = temp2;
temp2.left = null;
temp2.right = null;
return;
}
/* Proxels have the same id, just add their vals */
if (id == temp.id) {
temp.val += val;
return;
}
System.out.printf("\n\n\n!!!!!! addproxel failed !!!!!\n\n\n");
}
/********************************************************/
/* model specific distribtuions */
/********************************************************/
/* INSTANTANEOUS RATE FUNCTION 1 */
static double hpm2lpm(double age) {
//return unihrf(age, 0.25, .5);
return weibullhrf(age,55,4,0);
}
/* INSTANTANEOUS RATE FUNCTION 2 */
static double lpm2hpm(double age) {
//return exphrf(age, 2);
return unihrf(age, 9,11);
}
static double manufacture(double age){
double temp;//to ensure the manufacture process happening every 1 minute
if(Math.round(age)>age)
{
temp=age-(Math.round(age)-1);
}
else
temp=age-Math.round(age);
return dethrf(temp,1);
}
/********************************************************/
/* main processing loop */
/********************************************************/
public static void main(String[] args){
int kmax; //necessary sim level in proxel tree including the initial level
Proxel currproxel;
double val, z,m;
int s, tau1k; //, tau2k;
/* initialize the simulation */
root[0] = null;
root[1] = null;
eerror=0.0;
totcnt = 0;
maxccp = 0;
tmax = ENDTIME;
dt = DELTA;
kmax=(int) (tmax/dt+1); //add one to include the simulation at time 0
y=new double[3][kmax+2];
for (int k = 0; k < 3; k++) {
//..........................
for (int j = 0; j < kmax+2; j++)
y[k][j] = 0.0;
}
TAUMAX =(int)(tmax/dt+1); //TAUMAX=kmax?
defective=0.0;
working=0.0;
/* set initial proxel */
addproxel(HPM, 0, 0, 1.0);
/* first loop: iteration over all kmax time steps*/
for (int k = 1; k <= kmax+1; k++) {
//if(k==79 || k==78) printtree(root[sw]);
//printf("\nSTEP %d\n",k);
/* current model time is k*dt */
/* print progress information */
if (k%100==0) {
System.out.printf("\nSTEP %d\n",k);
System.out.printf("Size of tree %d\n",size(root[sw]));
}
sw = 1 - sw;
/* second loop: iterating over all proxels of a time step */
while (root[1-sw] != null)
{
totcnt++;
currproxel = getproxel();
while ((currproxel.val < MINPROB) && (root[1-sw] != null)) {
val=currproxel.val;
eerror += val;
currproxel = getproxel();
}
val = currproxel.val;
tau1k = currproxel.tau1k;
/*tau2k = currproxel->tau2k;*/
s = currproxel.s;
y[s][k-1] += val;
/* create child proxels */
switch (s) {
case HPM:
z = dt*hpm2lpm(tau1k*dt);
m=dt*manufacture(tau1k*dt);
if(m==1){
defective+=0.05;//expected value of the products
working+=0.95;
}
if (z < 1.0) {
addproxel(LPM, 0, 0, val*z);
addproxel(HPM, tau1k+1, 0, val*(1-z));
} else
addproxel(LPM, 0, 0, val);
break;
case LPM :
z = dt * lpm2hpm(tau1k*dt);
m=dt*manufacture(tau1k*dt);
if(m==1){
defective+=0.2;
working+=0.8;
}
if (z < 1.0) {
addproxel(HPM, 0, 0, val*z);
addproxel(LPM, tau1k+1, 0, val*(1-z));
} else
addproxel(HPM, 0, 0, val);
}
}
}
System.out.printf("error = %7.5e\n", eerror);
System.out.printf("ccpx = %d\n", maxccp);
System.out.printf("count = %d\n", totcnt);
System.out.printf("count = %d\n", defective);
plotsolution(kmax);
//return 0;
}
}