Hi all, I'm currently using pmap to distribute to a bunch of cores some work to do. I was just wondering how you guys are structuring this.
Basically I have a function written in C that I need to pass a bunch of parameters to so that it can do its work, so I first create a tuple with all the parameters needed for one worker to find results, for every combination of parameters: (I'm solving an infinite dimensional programming problem using bellman equations and brute force to maximize, the C code searches through the control space and finds the best choices to maximize the objective) input = [(w,z,mu,model,currentValue[:,:,nextPromise[mu]],squeeze(squeeze(capitals[w,z,:,:,:,:],2),1)) for w in 1:nW, z in 1:nZ, mu in 1:nMu]; And then I simply pass it with pmap to the wrapper of my C function: results = pmap(findPolicyC,input) And use linear indexing to recuperate the answers: anArray = results[w + nW*((z-1) + nZ*(mu-1))] Is this how you guys would structure this parallel computation? I've attached the code for any that are interested. You need GSL to run the C part: Eg. compile with, where the directories point to GSL: gcc -I/home/gcam/china/gsl -L/home/gcam/china/gsl/.libs -shared -fPIC findPolicySOE.c -lgsl -lm -o findPSoe.so Thanks for the input!
aux.jl
Description: Binary data
//
//
//
// C functions for heavy duty
//
// gcc -shared -fPIC findPolicySOE.c -lgsl -lm -lgslcblas -o findPSoe.so
// gcc -I/home/gcam/china/gsl -L/home/gcam/china/gsl/.libs -shared -fPIC findPolicySOE.c -lgsl -lm -o findPSoe.so
//
//
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_roots.h>
struct rootFinderParams { double theta; double delta; double q1Value; double q2Value; double q3Value; double wealth;
double betaGov; double p1; double p2; double p3;};
int linearIndex2d(int x, int nX, int y){ return x + nX*y; }
int linearIndex3d(int x, int nX, int y, int nY, int z){ return x + nX*(y + nY*z); }
int linearIndex4d(int x, int nX, int y, int nY, int z, int nZ, int last){ return x + nX*(y + nY*(z + nZ*last)); }
double* findPolicySOEInC(int w, int z, double mu, double alpha, double theta, double delta, double beta, double* vW, int nW, double* vK,
int nK, double* vQ, int nQ,double* vZ, double* currentValueZ1, double* currentValueZ2, double* currentValueZ3, double* transition,double* capital, double betaGov){
double highestVal = -1000000000009;
int policyK = -1;
int policyQ1 = -1;
int policyQ2 = -1;
int policyQ3 = -1;
//Generate interpolation functions
//gsl_set_error_handler_off();
gsl_interp_accel *accZ1 = gsl_interp_accel_alloc ();
gsl_interp_accel *accZ2 = gsl_interp_accel_alloc ();
gsl_interp_accel *accZ3 = gsl_interp_accel_alloc ();
// One for each state
gsl_spline *splineZ1 = gsl_spline_alloc (gsl_interp_cspline, nW);
gsl_spline *splineZ2 = gsl_spline_alloc (gsl_interp_cspline, nW);
gsl_spline *splineZ3 = gsl_spline_alloc (gsl_interp_cspline, nW);
gsl_spline_init (splineZ1, vW, currentValueZ1, nW);
gsl_spline_init (splineZ2, vW, currentValueZ2, nW);
gsl_spline_init (splineZ3, vW, currentValueZ3, nW);
int k,q1,q2,q3;
for(k = 0; k < nK; k++){
for(q1 = 0; q1 < nQ; q1++){
for(q2 = 0; q2 < nQ; q2++){
for(q3 = 0; q3 < nQ; q3++){
double newVal = 1.0;
double capitalFromPrecompute = capital[linearIndex4d(k,nK,q1,nQ,q2,nQ,q3)];
double kValue = capitalFromPrecompute <= 0.01 ? vK[k] : capitalFromPrecompute;
double q1Value = theta*(1-delta)*kValue - vQ[q1] < 0.0 ? theta*(1-delta)*kValue : vQ[q1];
double q2Value = theta*(1-delta)*kValue - vQ[q2] < 0.0 ? theta*(1-delta)*kValue : vQ[q2];
double q3Value = theta*(1-delta)*kValue - vQ[q3] < 0.0 ? theta*(1-delta)*kValue : vQ[q3];
double w1 = vZ[0]*pow(kValue,alpha) + (1-delta)*kValue - q1Value;
double w2 = vZ[1]*pow(kValue,alpha) + (1-delta)*kValue - q2Value;
double w3 = vZ[2]*pow(kValue,alpha) + (1-delta)*kValue - q3Value;
double expectedDebt = q1Value*transition[linearIndex2d(z,3,0)] + q2Value*transition[linearIndex2d(z,3,1)] +
q3Value*transition[linearIndex2d(z,3,2)];
double d = vW[w] - kValue + betaGov*expectedDebt;
if(w1 < vW[0] || w2 < vW[0] || w3 < vW[0] || capitalFromPrecompute < -1){ // last check is for states where no k' makes dividends non-negative
newVal = -10000000000;
} else{
double valueTomZ1 = 0.0;
double valueTomZ2 = 0.0;
double valueTomZ3 = 0.0;
// To compute tommorrow's value, interpolate the value function
// Take care of off-grid values with linear extrapolation
if(w1 < vW[0]){
valueTomZ1 = currentValueZ1[0] + ((w1 - vW[0])/(vW[1]-vW[0]))*(currentValueZ1[1] - currentValueZ1[0]);
} else if(w1 > vW[nW-1]){
valueTomZ1 = currentValueZ1[nW-1] + ((w1 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ1[nW-1] - currentValueZ1[nW-2]);
} else{
valueTomZ1 = gsl_spline_eval(splineZ1, w1, accZ1);
}
if(w2 < vW[0]){
valueTomZ2 = currentValueZ2[0] + ((w2 - vW[0])/(vW[1]-vW[0]))*(currentValueZ2[1] - currentValueZ2[0]);
} else if(w2 > vW[nW-1]){
valueTomZ2 = currentValueZ2[nW-1] + ((w2 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ2[nW-1] - currentValueZ2[nW-2]);
}
else{
valueTomZ2 = gsl_spline_eval(splineZ2, w2, accZ2);
}
if(w3 < vW[0]){
valueTomZ3 = currentValueZ3[0] + ((w3 - vW[0])/(vW[1]-vW[0]))*(currentValueZ3[1] - currentValueZ3[0]);
} else if(w3 > vW[nW-1]){
valueTomZ3 = currentValueZ3[nW-1] + ((w3 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ3[nW-1] - currentValueZ3[nW-2]);
} else{
valueTomZ3 = gsl_spline_eval(splineZ3, w3, accZ3);
}
//printf("cap %f intV %f for w1 %f \n",capitalFromPrecompute,valueTomZ1,w1);
double continuation = valueTomZ1*transition[linearIndex2d(z,3,0)] + valueTomZ2*transition[linearIndex2d(z,3,1)] +
valueTomZ3*transition[linearIndex2d(z,3,2)];
newVal = d*(1+mu) + beta*continuation;
}
if(newVal > highestVal){
highestVal = newVal;
policyK = k+1; //Julia 1 based indexing
policyQ1 = q1 +1;
policyQ2 = q2 +1;
policyQ3 = q3 +1;
}
}
}
}
}
// Now release all splines
gsl_spline_free (splineZ1);
gsl_interp_accel_free (accZ1);
gsl_spline_free (splineZ2);
gsl_interp_accel_free (accZ2);
gsl_spline_free (splineZ3);
gsl_interp_accel_free (accZ3);
double *results = (double *)malloc(sizeof(double)*5);
results[0] = policyK;
results[1] = policyQ1;
results[2] = policyQ2;
results[3] = policyQ3;
results[4] = highestVal;
return results;
}
int c1(double k, void* p){
struct rootFinderParams *params = (struct rootFinderParams *)p;
int answer = params->theta *(1-params->delta)*k - params->q1Value <= 0.0 ? 1 : 0;
return answer;
}
int c2(double k, void* p){
struct rootFinderParams *params = (struct rootFinderParams *)p;
int answer = params->theta *(1-params->delta)*k - params->q2Value <= 0.0 ? 1 : 0;
return answer;
}
int c3(double k, void* p){
struct rootFinderParams *params = (struct rootFinderParams *)p;
int answer = params->theta *(1-params->delta)*k - params->q3Value <= 0.0 ? 1 : 0;
return answer;
}
double findRootInC(double k,void* p){
struct rootFinderParams *params = (struct rootFinderParams *)p;
double answer = k - (params->wealth + params->betaGov*((1-c1(k,p))*params->q1Value*params->p1 +
(1-c2(k,p))*params->q2Value*params->p2 +(1-c3(k,p))*params->q3Value*params->p3 ))/(
1 - params->betaGov*params->theta*(1-params->delta)*
(c1(k,p)*params->p1 +
c2(k,p)*params->p2 +
c3(k,p)*params->p3 ));
return answer;
}
/*
This function precomputes all possible state-policy combinations
to check whether or not the non-negative dividend constraint is satisfied.
Basically for every state (w,z), go through each policy (k,q1,q2,q3) and check
1. Given this k, do I violate the collateral constraints for any qi?
If yes, make them binding, if not ignore them
2. Are the dividends that I give in this period non-negative?
If d >= 0 then assign a 0 to the capitals matrix for this slot
If d<0 find the choice of k' that make dividends exactly zero and assign
to the capitals matrix
*/
double* findZeroInC(double wealth, int z,int k, double theta, double delta, double* vK,
int nK, double* vQ, int nQ, double* transition, double betaGov){
int q1,q2,q3;
double *capitalChoices = (double *)malloc(sizeof(double)*nQ*nQ*nQ);
for(q1 = 0; q1 < nQ; q1++){
for(q2 = 0; q2 < nQ; q2++){
for(q3 = 0; q3 < nQ; q3++){
double kValue = vK[k];
double q1Value = vQ[q1];
double q2Value = vQ[q2];
double q3Value = vQ[q3];
double collateral1 = theta*(1-delta)*kValue -q1Value;
double collateral2 = theta*(1-delta)*kValue - q2Value;
double collateral3 = theta*(1-delta)*kValue -q3Value;
if(collateral1 < 0.0){
q1Value = theta*(1-delta)*kValue;
}
if(collateral2 < 0.0){
q2Value = theta*(1-delta)*kValue;
}
if(collateral3 < 0.0){
q3Value = theta*(1-delta)*kValue;
}
double expectedDebt = q1Value*transition[linearIndex2d(z,3,0)] + q2Value*transition[linearIndex2d(z,3,1)] +
q3Value*transition[linearIndex2d(z,3,2)];
double d = wealth - kValue + betaGov*expectedDebt;
if(d < 0){
struct rootFinderParams params = {theta, delta, q1Value, q2Value, q3Value, wealth, betaGov, transition[linearIndex2d(z,3,0)],
transition[linearIndex2d(z,3,1)],transition[linearIndex2d(z,3,2)]};
int signLow = signbit(findRootInC(0.0,¶ms)); // 0 if positive, nonzero if negative
int signHigh = signbit(findRootInC(vK[nK-1],¶ms));
if((signLow == 0 && signHigh == 0) || (signLow != 0 && signHigh != 0)) { // This is the case when wealth - debt is negative, so no K>=0 gives zero divs
capitalChoices[linearIndex3d(q1,nQ,q2,nQ,q3)] = -5;
} else{
// Must solve the non-linear equation for k'
int status;
int iter = 0, max_iter = 100;
double root = 0;
double diff = 0;
gsl_function F;
F.function = &findRootInC;
F.params = ¶ms;
//printf("ELSE: At w %f beta*expD %f, fval0 %f, fvalTop %f, q1 %f, q2 %f, q3 %f k %f\n",wealth,betaGov*expectedDebt,findRootInC(0.0,¶ms),findRootInC(vK[nK-1],¶ms),q1Value,q2Value,q3Value,kValue);
const gsl_root_fsolver_type * T = gsl_root_fsolver_bisection;
gsl_root_fsolver * s = gsl_root_fsolver_alloc (T);
double kLow = 0.0;
double kHigh = vK[nK-1];
gsl_root_fsolver_set (s,&F,kLow,kHigh);
do{
iter++;
status = gsl_root_fsolver_iterate (s);
root = gsl_root_fsolver_root (s);
diff = findRootInC(root,¶ms);
status = gsl_root_test_residual(diff,0.000001);
}
while (status == GSL_CONTINUE && iter < max_iter);
//printf("the root is %f \t diff, %f \t iterations %d \t q2Value in params %f coll %f actualq2 %f \n",root,findRootinC(root,¶ms),iter,params.q2Value,theta*(1-delta)*kValue - vQ[q2],vQ[q2] );
capitalChoices[linearIndex3d(q1,nQ,q2,nQ,q3)] = root;
gsl_root_fsolver_free (s);
}
} else{
capitalChoices[linearIndex3d(q1,nQ,q2,nQ,q3)] = 0;
}
}
}
}
return capitalChoices;
}
double findPolicySocialPlannerInC(int w, int z, double alpha, double theta, double delta, double beta, double* vW, int nW, double* vK,
int nK, double* vQ, int nQ,double* vZ, double* currentValueZ1, double* currentValueZ2, double* currentValueZ3, double* transition, double betaGov){
double highestVal = -1000000000009;
//Generate interpolation functions
//gsl_set_error_handler_off();
//Some tests to make sure it's all OK
gsl_interp_accel *accZ1 = gsl_interp_accel_alloc ();
gsl_interp_accel *accZ2 = gsl_interp_accel_alloc ();
gsl_interp_accel *accZ3 = gsl_interp_accel_alloc ();
// One for each state
gsl_spline *splineZ1 = gsl_spline_alloc(gsl_interp_cspline, nW); // gsl_spline_alloc(gsl_interp_akima,nW);
gsl_spline *splineZ2 = gsl_spline_alloc(gsl_interp_cspline, nW); //(gsl_interp_akima,nW);
gsl_spline *splineZ3 = gsl_spline_alloc(gsl_interp_cspline, nW); //(gsl_interp_akima,nW);
gsl_spline_init (splineZ1, vW, currentValueZ1, nW);
gsl_spline_init (splineZ2, vW, currentValueZ2, nW);
gsl_spline_init (splineZ3, vW, currentValueZ3, nW);
double valueTomZ1 = 0.0;
double valueTomZ2 = 0.0;
double valueTomZ3 = 0.0;
int k,q1,q2,q3;
for(k = 0; k < nK; k++){
for(q1 = 0; q1 < nQ; q1++){
for(q2 = 0; q2 < nQ; q2++){
for(q3 = 0; q3 < nQ; q3++){
double newVal = 1.0;
double kValue = vK[k];
double q1Value = vQ[q1];
double q2Value = vQ[q2];
double q3Value = vQ[q3];
double w1 = vZ[0]*pow(kValue,alpha) + (1-delta)*kValue - q1Value;
double w2 = vZ[1]*pow(kValue,alpha) + (1-delta)*kValue - q2Value;
double w3 = vZ[2]*pow(kValue,alpha) + (1-delta)*kValue - q3Value;
double expectedDebt = q1Value*transition[linearIndex2d(z,3,0)] + q2Value*transition[linearIndex2d(z,3,1)] +
q3Value*transition[linearIndex2d(z,3,2)];
double d = vW[w] - kValue + betaGov*expectedDebt;
if(d < 0){
kValue = vW[w] + betaGov*expectedDebt;
w1 = vZ[0]*pow(kValue,alpha) + (1-delta)*kValue - q1Value;
w2 = vZ[1]*pow(kValue,alpha) + (1-delta)*kValue - q2Value;
w3 = vZ[2]*pow(kValue,alpha) + (1-delta)*kValue - q3Value;
d = 0;
}
if(w1 < vW[0] || w2 < vW[0] || w3 < vW[0]){//(w1 < 0 || w2 < 0 || w3 < 0){
newVal = -10000000000;
} else{
// To compute tommorrow's value, interpolate the value function
// Take care of off-grid values with linear extrapolation
if(w1 < vW[0]){
valueTomZ1 = currentValueZ1[0] + ((w1 - vW[0])/(vW[1]-vW[0]))*(currentValueZ1[1] - currentValueZ1[0]);
} else if(w1 > vW[nW-1]){
valueTomZ1 = currentValueZ1[nW-1] + ((w1 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ1[nW-1] - currentValueZ1[nW-2]);
} else{
valueTomZ1 = gsl_spline_eval(splineZ1, w1, accZ1);
}
if(w2 < vW[0]){
valueTomZ2 = currentValueZ2[0] + ((w2 - vW[0])/(vW[1]-vW[0]))*(currentValueZ2[1] - currentValueZ2[0]);
} else if(w2 > vW[nW-1]){
valueTomZ2 = currentValueZ2[nW-1] + ((w2 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ2[nW-1] - currentValueZ2[nW-2]);
}
else{
valueTomZ2 = gsl_spline_eval(splineZ2, w2, accZ2);
}
if(w3 < vW[0]){
valueTomZ3 = currentValueZ3[0] + ((w3 - vW[0])/(vW[1]-vW[0]))*(currentValueZ3[1] - currentValueZ3[0]);
} else if(w3 > vW[nW-1]){
valueTomZ3 = currentValueZ3[nW-1] + ((w3 - vW[nW-1])/(vW[nW-1]-vW[nW-2]))*(currentValueZ3[nW-1] - currentValueZ3[nW-2]);
} else{
valueTomZ3 = gsl_spline_eval(splineZ3, w3, accZ3);
}
//printf("cap %f intV %f for w1 %f \n",capitalFromPrecompute,valueTomZ1,w1);
double continuation = valueTomZ1*transition[linearIndex2d(z,3,0)] + valueTomZ2*transition[linearIndex2d(z,3,1)] +valueTomZ3*transition[linearIndex2d(z,3,2)];
//printf("k %d q1 %d q2 %d q3 %d c1 %f c2 %f c3 %f w1 %f w2 %f w3 %f \n",k+1,q1+1,q2+1,q3+1,valueTomZ1,valueTomZ2,valueTomZ3,w1,w2,w3);
newVal = d + beta*continuation;
}
if(newVal > highestVal){
highestVal = newVal;
}
}
}
}
}
// Now release all splines
gsl_spline_free (splineZ1);
gsl_interp_accel_free (accZ1);
gsl_spline_free (splineZ2);
gsl_interp_accel_free (accZ2);
gsl_spline_free (splineZ3);
gsl_interp_accel_free (accZ3);
return highestVal;
}
soe.jl
Description: Binary data
