Hi,

I have this monster program. But, mostly I need help in Main and in Derivs. I know some is written in a juvenille form, but I'm not a C++ student!

I have defined all these arrays in the Main function and I need to write an "if" statement in the Derivs function to define a value for ngen (genotype), for sex (both, male, female), and risk level (1, 0.5, 0.25). I then need to write a small function to calculate a transition intensity using each of the given values for ngen, sex, and risk inside the if statement. I need to have a loop in the main function where I can change the values of ngen, sex, and male so I can run the program for only one set of values, then easily change the values and run it again.

Does any of this make any sense? Any help, even how to write an if statement that has three conditions would be appreciated.

I've tried to do an if statement in Derivs, but don't know if it's correct.

Thanks!!

Here's my program: (remember, Main and Derivs only. Don't be frightened off)

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <fstream.h>


//GLOBAL VARIABLES

const int nvar = 4; //No. of states
const int nmom = 3; //No. of moments
const int ngen = 5; //No. of ApoE genotypes
const int sex = 3; //both, female, and male
const int risk = 3; //No. of levels of relative risk
int x1; //Starting value for t
int x2; //Ending value for t
double h; //Step size

double delta;

double A[nvar+1][nvar+1];
double B[nvar+1][nvar+1];
double C[nvar+1][nvar+1];

double D[nvar+1];
double E[ngen+1][sex+1];
double F[ngen+1][sex+1];
double G[ngen+1][sex+1];
double H[ngen+1][sex+1];
double I[nvar+1][nvar+1];
double J[nvar+1];
double K[nvar+1][nvar+1];
double k1[ngen+1][sex+1];
double k2[ngen+1][sex+1];
double m[risk+1];
double r_m[sex+1][risk+1];

double dydx[nvar+1][nmom+1];

double mu[nvar+1][nvar+1];
double f[ngen+1][sex+1];

double vstart[nvar+1][nmom+1], cmom[nvar+1][nmom+1], mom[nvar+1][nmom+1] ;
double v[nvar+1][nmom+1], vout[nvar+1][nmom+1], dv[nvar+1][nmom+1];
double dym[nvar+1][nmom+1], dyt[nvar+1][nmom+1], yt[nvar+1][nmom+1];


//PROGRAM HEADERS

void dumpMat1(double a[nvar+1][nvar+1]);

double fact (int); //factorial program
double comb (int,int); //combinatorial program
double power (double i,double j); //to the power program

void derivs (double t, double yt[nvar+1][nmom+1], double dyt[nvar+1][nmom+1]);

void rk4(double y[nvar+1][nmom+1], double dydx[nvar+1][nmom+1], double x,
double yout[nvar+1][nmom+1]);

void rkdumb (double vstart[nvar+1][nmom+1], double cmom[nvar+1][nmom+1]);

void rkmom(double cmoments[nvar+1][nmom+1], double moments[nvar+1][nmom+1]);

void dumpMat(double a[nvar+1][nmom+1]);


//SPECIFY OUTPUT FILE (NEEDS TO BE CHANGED)

ofstream outfile ("h:/test.txt", ios::out);


//////////////////////////////////////////////////////////////////////////////////////
// //
// RUNGE KUTTA ALGORITHM ADAPTED FROM 'NUMERICAL RECIPES IN C', PRESS ET AL. (1989) //
// CAMBRIDGE UNIVERSITY PRESS. PROGRAM TO SOLVE NORBERGS EQUATIONS (1995) FOR THE //
// MOMENTS OF THE PRESENT VALUE OF BENEFITS IN A GENERAL MARKOV MODEL. //
// //
// THE FUNCTION 'DERIVS' NEEDS TO BE COMPLETED IN THIS VERSION - IT SHOULD BE //
// PROGRAMMED TO RETURN THE DERIVATIVE VALUES GIVEN BY NORBERG'S EQUATIONS. //
// //
// DELME PRITCHARD 28/7/2004 Version 1.1 //
// //
//////////////////////////////////////////////////////////////////////////////////////


// MAIN PROGRAM

main()
{
int i, j, k, t;
delta = 0.05;

A[1][2] = 0; //describes the transition intensities
A[1][3] = 0;
A[1][4] = 0;
A[2][1] = 0;
A[2][3] = 0.18896;
A[2][4] = 0;
A[3][1] = 0;
A[3][2] = 0;
A[3][4] = 0.173;


B[1][2] = 0.000000131275; //describes the transition intensitites
B[1][3] = 0;
B[1][4] = 0.65*0.000025934;
B[2][1] = 0;
B[2][3] = 0;
B[2][4] = 0.33502*0.65*0.000025934;
B[3][1] = 0;
B[3][2] = 0;
B[3][4] = 0.65*0.000025934;

C[1][2] = 0.145961; //describes the transition intensitities
C[1][3] = 0;
C[1][4] = 0.093605;
C[2][1] = 0;
C[2][3] = 0;
C[2][4] = 0.093605;
C[3][1] = 0;
C[3][2] = 0;
C[3][4] = 0.093605;

A[1][1] = 0; //describes the transition intensitities
B[1][1] = 0;
C[1][1] = 0;
A[2][2] = 0;
B[2][2] = 0;
C[2][2] = 0;
A[3][3] = 0;
B[3][3] = 0;
C[3][3] = 0;
A[4][4] = 0;
B[4][4] = 0;
C[4][4] = 0;

D[1] = 0; //describes the benefits
D[2] = 0;
D[3] = 0.05;
D[4] = 0;

I[1][2] = 0; //describes the benefits
I[1][3] = 0;
I[1][4] = 0;
I[2][1] = 0;
I[2][3] = 0;
I[2][4] = 0;
I[3][1] = 0;
I[3][2] = 0;
I[3][4] = 0;

J[1] = 0; //describes the benefits
J[2] = 0;
J[3] = 1;
J[4] = 0;

K[1][2] = 0; //describes the benefits
K[1][3] = 0;
K[1][4] = 0;
K[2][1] = 0;
K[2][3] = 0;
K[2][4] = 0;
K[3][1] = 0;
K[3][2] = 0;
K[3][4] = 0;
K[1][1] = 0;
K[2][2] = 0;
K[3][3] = 0;
K[4][4] = 0;

E[1][1] = 0.754;
E[1][2] = 0.675;
E[1][3] = 0.434;
E[3][1] = 2.87;
E[3][2] = 4.21;
E[3][3] = 1.42;
E[4][1] = 2.98;
E[4][2] = 3.68;
E[4][3] = 1.92;
E[5][1] = 13.5;
E[5][2] = 10.4;
E[5][3] = 8.94;

F[1][1] = 0;
F[1][2] = 0;
F[1][3] = 0;
F[3][1] = 0.00938;
F[3][2] = 0.0102;
F[3][3] = 0.00506;
F[4][1] = 0.00312;
F[4][2] = 0.00319;
F[4][3] = 0.00103;
F[5][1] = 0.00529;
F[5][2] = 0.00504;
F[5][3] = 0.00656;

G[1][1] = 0.00859;
G[1][2] = 0.00692;
G[1][3] = 0.0160;
G[3][1] = 0;
G[3][2] = 0;
G[3][3] = 0;
G[4][1] = 0;
G[4][2] = 0;
G[4][3] = 0;
G[5][1] = 0;
G[5][2] = 0;
G[5][3] = 0;
H[1][1] = 0;
H[1][2] = 0;
H[1][3] = 0;
H[3][1] = 1;
H[3][2] = 1;
H[3][3] = 0;
H[4][1] = 1;
H[4][2] = 1;
H[4][3] = 0;
H[5][1] = 1;
H[5][2] = 1;
H[5][3] = 1;

k1[1][1] = 0;
k1[1][2] = 0;
k1[1][3] = 0;
k1[3][1] = 68;
k1[3][2] = 68;
k1[3][3] = 67;
k1[4][1] = 62;
k1[4][2] = 62;
k1[4][3] = 51;
k1[5][1] = 60;
k1[5][2] = 60;
k1[5][3] = 60;

k2[1][1] = 60;
k2[1][2] = 60;
k2[1][3] = 60;
k2[3][1] = 0;
k2[3][2] = 0;
k2[3][3] = 0;
k2[4][1] = 0;
k2[4][2] = 0;
k2[4][3] = 0;
k2[5][1] = 0;
k2[5][2] = 0;
k2[5][3] = 0;

m[1] = 1;
m[2] = 0.5;
m[3] = 0.25;

r_m[1][1] = 0.93;
r_m[1][2] = 0.96;
r_m[1][3] = 0.97;
r_m[2][1] = 0.88;
r_m[2][2] = 0.94;
r_m[2][3] = 0.97;
r_m[3][1] = 1.27;
r_m[3][2] = 1.11;
r_m[3][3] = 1.05;

E[2][1] = 0;
E[2][2] = 0;
E[2][3] = 0;
F[2][1] = 0;
F[2][2] = 0;
F[2][3] = 0;
G[2][1] = 0;
G[2][2] = 0;
G[2][3] = 0;
H[2][1] = 0;
H[2][2] = 0;
H[2][3] = 0;
k1[2][1] = 0;
k1[2][2] = 0;
k1[2][3] = 0;
k2[2][1] = 0;
k2[2][2] = 0;
k2[2][3] = 0;


h = -0.0005;

x1 = 120, x2 = 60;

for (int i = 0; i <= nvar; i++)
{
for(int j = 0; j <= nmom; j++)
{
vstart[i][j] = 0;
vstart[i][0] = 1;
cmom[i][j] = 0;
}
}




rkdumb (vstart, cmom);

rkmom (cmom, mom);

dumpMat(mom);

return 0;
}



// START OF RKDUMB

void rkdumb (double vstart[nvar+1][nmom+1], double cmom[nvar+1][nmom+1])
{

int i, j, q;
double x, nstep;
long int k;

for (q = 0; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
v[i][q] = vstart[i][q];
cmom[i][q] = vstart[i][q];
}
}

x = x1;
nstep =(x2 - x1)/h;

for (k = 1; k <= nstep; k++)
{
derivs(x, v, dv);
rk4(v, dv, x, vout);

if (x + h == x) outfile << "ERROR"<<endl;

x += h;

for (j = 1; j <= nmom; j++)
{
for (i = 1; i <= nvar; i ++)
{
v [i][j] = vout[i][j];

if (k > nstep-1)
{
cmom[i][j] = vout[i][j];
cmom[i][0] = 1;
}
}
}
}
}


//START OF RK4

void rk4(double y[nvar+1][nmom+1], double dydx [nvar+1][nmom+1], double x,
double yout [nvar+1][nmom+1])
{

int i, q;

double hh = h*0.5, h6=h/6.0, xh=x+hh;


for (q = 1; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
dyt[i][q] = 0;
}
}

for (q = 1; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
yt[i][q] = y[i][q] + hh * dydx[i][q];
}
}


derivs (xh, yt, dyt);


for (q = 1; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
yt[i][q] = y[i][q] + hh * dyt[i][q];
}
}


derivs (xh, yt, dym);


for (q = 1; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
yt[i][q] = y[i][q] + h * dym[i][q];
dym[i][q] += dyt[i][q];
}
}

derivs (x+h, yt, dyt);

for (q = 1; q <= nmom; q++)
{
for (i = 1; i <= nvar; i++)
{
yout[i][q] = y[i][q] + h6 * (dydx[i][q] + dyt[i][q] + 2.0 * dym[i][q]);
}
}

}


//START OF DERIVS

void derivs (double t, double yt[nvar+1][nmom+1], double dyt[nvar+1][nmom+1])
{
int i, j, r, q, k;

double summation[nvar +1][nmom +1];
double mu_start[nvar+1];
double sum_part[nvar+1][nvar+1][nmom+1];
double mu[nvar+1][nvar+1];
double beni[nvar+1];
double benij[nvar+1][nvar+1];










for(i=1; i <= nvar; i++)
{
yt[i][0] = 1; // set 0th conditional moments equal to zero.
for(j=1; j <= nvar; j++)
{
mu[i][j] = A[i][j] + B[i][j] * exp(C[i][j]*t);
//cout<<"mu"<<i<<j<<" = "<<mu[i][j]<<endl;
// char w; cin>>w;
}

}



if(i = 1)
{
if(j = 1)
{
if(k = 1)
{
f[i][j] = E[i][j] * (exp(-F[i][j]*(power((t-k1[i][j]), 2) - G[i][j]*(t - k2[i][j])))) + H[i][j];
mu[1][2] = r_m[j][k] * ((f[i][j] - 1) * m[k] + 1) * mu[1][2];
}
}
}




for(i=1; i <= nvar; i++)
{
beni[i] = J[i] * exp(D[i]*(t-x2));

for(j=1; j <= nvar; j++)
{
benij[i][j] = K[i][j] * exp(I[i][j]*(t-x2));
//cout<<"ben"<<i<<" = "<<beni[i]<<endl;
//cout<<"ben"<<i<<j<<" = "<<benij[i][j]<<endl;
//char w; cin>>w;

}
}



for(i=1; i <= nvar; i++)
{
mu_start[i] = 0;

for(j=1; j <= nvar; j++)
{
mu_start[i] += mu[i][j];
// cout<<"mu "<<i<<" = "<<mu_start[i]<<endl;
// char w; cin>>w;
}
}

for(i=1; i <= nvar; i++) //sums up the second part of the summation in Norberg's equation
{
for(q=1; q <= nmom; q++)
{
for(k=1; k <= nvar; k++)
{
sum_part[i][k][q] = 0;

for(r=0; r <= q; r++)
{
sum_part[i][k][q] += comb(q, r) * power(benij[i][k], r) * yt[k][q-r];
}
}
}
}

for(i = 1; i <= nvar; i++) //multiplies and stores the total summation in Norberg's equation
{
for(q = 1; q <= nmom; q++)
{
summation[i][q] = 0;

for(k = 1; k <= nvar; k++)
{
summation[i][q] += mu[i][k] * sum_part[i][k][q];
}
}
}

for(i=1; i <= nvar; i++)
{
for(q=1; q <= nmom; q++)
{
dyt[i][q] = (q * delta + mu_start[i]) * yt[i][q] - q * beni[i] * yt[i][q-1] - summation[i][q];
}

}
}


// OTHER PROGRAMS

void rkmom(double cmoments[nvar+1][nmom+1], double moments[nvar+1][nmom+1])
{
for (int i = 0; i <= nvar; i++)
{
for(int j = 0; j <= nmom; j++)
{
moments[i][j] = 0;
}
}

for (int q = 0; q <= nmom; q++)
{
for (int i = 1; i <= nvar; i++)
{
for (int p = 0; p <= q; p++)
{
if (q == 0 || q == 1) moments[i][q] = cmoments[i][q];
else
moments[i][q] += comb(q,p) * power(-1, q-p) * cmoments[i][p] *
power(cmoments[i][1], q-p);
}
}
}
}

void dumpMat1(double a[nvar+1][nvar+1])
{
for (int i = 1; i <= nvar; i++)
{
for (int j = 1; j<= nvar; j++)
{
outfile << a[i][j] << " \t";
}
outfile<<endl;
}
outfile << endl;
}


void dumpMat(double a[nvar+1][nmom+1])
{
for (int i = 1; i <= nvar; i++)
{
for (int j = 1; j<= nmom; j++)
{
outfile << a[i][j] << " \t";
}
outfile<<endl;
}
outfile << endl;
}


double fact(int n)
{
if (n < 0) return 0;
double f = 1;
while (n > 1) f *= n--;

return f;
}


double comb(int i,int j)
{
double c = 1;
if (j==0) return c;
else c = (fact (i)) / (fact (j) * fact (i-j));

return c;
}


double power(double i,double j)
{
double p=1;
if (j == 0) return 1;
else if (i==0) return 0;
else p = pow(i,j);

return p;
}