Added to the synch process an intial computation to estimate the optimal threshold (which a asych process would use) in order to comput the expected accuracy such process would have, and finally use the expected accuracy to set the interrogation time.

Andreagiovanni Reina · Andreagiovanni Reina · commit 251aa4302bf2 · 2019-08-14T10:44:03.000+01:00
diff --git a/conf/DecNet.config b/conf/DecNet.config
@@ -1,7 +1,7 @@
 [experiment]
 randomSeed=2289
 # bad case for space with size 1 and range 0.3 seed: 1670054
-numberOfExperiments=1
+numberOfExperiments=3
 repetitionsPerDDM=1
 outputTxtFile=../../data/out.txt
 outputPdfFile=../../data/output
@@ -25,15 +25,15 @@ truncatePoor=true
 
 [DDM]
 # possible driftDistributions are [ 'uniform', 'normal', 'from_accuracy' ]
-driftDistribution=from_accuracy
+driftDistribution=normal
 baseDrift=0.1
 driftStdDev=0.5
 randomDriftRangePlusMinus=1
 noiseStdDev=0.5
 dt=0.01
 threshold=1
 interrogationTimeFromAccuracy=true
-interrogationTime=0.6
+interrogationTime=10
 prior=0.5
 # costMatrix accepts two values [ timeCost, errorCost ]
 costMatrix=[ 1 , 20 ]
diff --git a/conf/DecNet.template.config b/conf/DecNet.template.config
@@ -31,7 +31,10 @@ randomDriftRangePlusMinus=DDM_DRIFT_RANGE_PLUS_MINUS
 noiseStdDev=DDM_NOISE_STD_DEV
 dt=0.01
 threshold=1
+interrogationTimeFromAccuracy=INTERR_TIME_FROM_ACC
 interrogationTime=INTERR_TIME
+prior=PRIOR_DIST
+costMatrix=[ COST_MATRIX_T , COST_MATRIX_E ]
 
 [Network]
 number_of_nodes=NUM_NODES
diff --git a/scripts/staticSPERAB.sh b/scripts/staticSPERAB.sh
@@ -7,64 +7,75 @@ mkdir -p ${CONF_DIR}
 TEMPLATE_SETTINGS="${PROJECT_HOME}/conf/DecNet.template.config"
 PYPATH="${PROJECT_HOME}/src/"
 EXEC_FILE="${PROJECT_HOME}/src/DecNet/DecisionProcess.py"
-OUTPUT_DATA_DIR="${PROJECT_HOME}/static_lin6/"
+OUTPUT_DATA_DIR="${PROJECT_HOME}/synch_1/"
 
 SEED=221189
 NUM_EXP=100
 NUM_RUN=1
 MAX_TIME=1000
 
-AGENT_TYPE='simple'
-DEC_MODEL_LIST=( 'conf-perfect' )
+AGENT_TYPE='DDM'
+DEC_MODEL_LIST=( 'log-odds-distr' )
 #UPDATE_CONF_LIST=( 'optim-up' 'belief-up')
-UPDATE_CONF_LIST=( 'finite-time')
+UPDATE_CONF_LIST=( 'no-up' 'optim-up' )
 #BELIEF_EPSILON_LIST=($(seq 0.07 0.001 0.1))
 #BELIEF_EPSILON_LIST=($(seq 0.01 0.005 0.15))
 #BELIEF_EPSILON_LIST=(0.01 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00)
 #BELIEF_EPSILON_LIST=($(seq 0.01 0.01 0.10))
 #BELIEF_EPSILON_LIST=($(seq 0.001 0.003 0.016))
-BELIEF_EPSILON_LIST=( 0 0.05 )
+BELIEF_EPSILON=0
 FINITE_TIME_EXPONENT=0.5
 
 ### Unnecessary for this analysis 
 DDM_DRIFT_DIST='normal'
-DDM_BASE_DRIFT=0.1
+#DDM_BASE_DRIFT=0.1
+DDM_BASE_DRIFT_LIST=( 0.05 0.1 0.2 )
 DDM_DRIFT_RANGE_PLUS_MINUS=1
-DDM_DRIFT_STD_DEV=0.3
+#DDM_DRIFT_STD_DEV=0.3
+DDM_DRIFT_STD_DEV_LIST=($(seq 0.10 0.10 0.50))
 DDM_NOISE_STD_DEV=0.5
-#THRESHOLD=1
+THRESHOLD=1
+INTERR_TIME_FROM_ACC='true'
 INTERR_TIME=10
+PRIOR_DIST=0.5
+COST_MATRIX_T=1
+COST_MATRIX_E=20
+
 
 #NUM_NODES=20
 NUM_NODES_LIST=($(seq 20 10 100))
 #NUM_NODES_LIST=($(seq 30 20 90))
 #NUM_NODES_LIST=($(seq 50 50 200))
 #NUM_NODES_LIST=($(seq 11 12 47))
 #NUM_NODES_LIST=( 40 )
-NET_TYPE_LIST=('erdos-renyi' 'barabasi-albert' 'space')
+NET_TYPE_LIST=('erdos-renyi' 'barabasi-albert' 'rgg-fixed-degree')
 #NET_TYPE_LIST=('space')
 #NET_TYPE_LIST=('erdos-renyi')
 #LINK_PROBABILITY_LIST=($(seq 0.2 0.1 0.8))
-LINK_PROBABILITY_LIST=(0.2 0.4 0.6 0.8)
+#LINK_PROBABILITY_LIST=(0.2 0.4 0.6 0.8)
+LINK_PROBABILITY_LIST=(0.2 0.6)
 #LINK_PROBABILITY_LIST=(0.2)
 #NUM_EDGES=3
-NUM_EDGES_LIST=($(seq 3 3 12))
+#NUM_EDGES_LIST=($(seq 3 3 12))
 #NUM_EDGES_LIST=($(seq 5 5 20))
+NUM_EDGES_LIST=(3 9)
 #INTERACTION_RADIUS_LIST=($(seq 0.05 0.05 0.3))
-INTERACTION_RADIUS_LIST=($(seq 0.20 0.05 0.35))
+#INTERACTION_RADIUS_LIST=($(seq 0.20 0.05 0.35))
+INTERACTION_RADIUS_LIST=( 10 )
 #INTERACTION_RADIUS_LIST=( 0.05 )
-ENV_SIZE_LIST=( 1 )
+ENV_SIZE=1
+#ENV_SIZE_LIST=( 1 )
 #ENV_SIZE_LIST=(0.5 2 3 4 5 6 7 8 9 15 20)
 #ENV_SIZE_LIST=(0.01 0.1 0.5 1 2 3 4 5 10 15 20)
 
 DYNAMICNET='false'
 #AGENTSPEED_LIST=($(seq 0.05 0.05 0.2))
-AGENTSPEED_LIST=( 0.05 )
+AGENTSPEED=0.05
 MOVE_STDDEV=45
 PERIODICBOUND='false'
 
-ACC_MEAN_LIST=( 0.6 )
-ACC_STD_DEV_LIST=( 0.12 )
+ACC_MEAN=0.6
+ACC_STD_DEV=0.12
 ACC_TRUNCATED='true'
 
 COUNT=0
@@ -81,106 +92,101 @@ do
 				do
 					for NUM_EDGES in ${NUM_EDGES_LIST[*]}
 					do
-						for ACC_MEAN in ${ACC_MEAN_LIST[*]}
+						for DDM_BASE_DRIFT in ${DDM_BASE_DRIFT_LIST[*]}
 						do
-							for ACC_STD_DEV in ${ACC_STD_DEV_LIST[*]}
+							for DDM_DRIFT_STD_DEV in ${DDM_DRIFT_STD_DEV_LIST[*]}
 							do
 								for UPDATE_CONF in ${UPDATE_CONF_LIST[*]}
 								do
-									for BELIEF_EPSILON in ${BELIEF_EPSILON_LIST[*]}
-									do
-										for AGENTSPEED in ${AGENTSPEED_LIST[*]}
-										do
-											for ENV_SIZE in ${ENV_SIZE_LIST[*]}
-											do
-												
-												### Skipping unnecessary iterations
-												if [ "$NET_TYPE" != "erdos-renyi" ] && [ $LINK_PROBABILITY != ${LINK_PROBABILITY_LIST[0]} ];
-												then
-													continue
-												fi
-												if [ "$NET_TYPE" != "barabasi-albert" ] && [ $NUM_EDGES != ${NUM_EDGES_LIST[0]} ];
-												then
-													continue
-												fi
-												if [ "$NET_TYPE" != "space" ] && [ $INTERACTION_RADIUS != ${INTERACTION_RADIUS_LIST[0]} ];
-												then
-													continue
-												fi
-												#if [ "$DEC_MODEL" == "belief-init" ] && [ $UPDATE_CONF != "belief-up" ];
-												#then
-												#	continue
-												#fi
-												if [ "$UPDATE_CONF" == "optim-up" ] && [ $BELIEF_EPSILON != ${BELIEF_EPSILON_LIST[0]} ];
-												then
-													continue
-												fi
+											
+									### Skipping unnecessary iterations
+									if [ "$NET_TYPE" != "erdos-renyi" ] && [ $LINK_PROBABILITY != ${LINK_PROBABILITY_LIST[0]} ];
+									then
+										continue
+									fi
+									if [ "$NET_TYPE" != "barabasi-albert" ] && [ $NUM_EDGES != ${NUM_EDGES_LIST[0]} ];
+									then
+										continue
+									fi
+									if [ "$NET_TYPE" != "space" ] && [ $INTERACTION_RADIUS != ${INTERACTION_RADIUS_LIST[0]} ];
+									then
+										continue
+									fi
+									#if [ "$DEC_MODEL" == "belief-init" ] && [ $UPDATE_CONF != "belief-up" ];
+									#then
+									#	continue
+									#fi
+									if [ "$UPDATE_CONF" == "optim-up" ] && [ $BELIEF_EPSILON != ${BELIEF_EPSILON_LIST[0]} ];
+									then
+										continue
+									fi
+									
+									NET_PAR=0
+									if [ "$NET_TYPE" == "erdos-renyi" ];
+									then
+										NET_PAR=${LINK_PROBABILITY}
+									fi
+									if [ "$NET_TYPE" == "barabasi-albert" ];
+									then
+										NET_PAR=${NUM_EDGES}
+									fi
+									if [ "$NET_TYPE" == "space" ];
+									then
+										NET_PAR=${INTERACTION_RADIUS}
+									fi
+									
+									JOB_PARAM="net-${NET_TYPE}_nodes-${NUM_NODES}_link-${NET_PAR}_model-${DEC_MODEL}_up-${UPDATE_CONF}_driftbase-${DDM_BASE_DRIFT}_driftrange-${DDM_DRIFT_STD_DEV}"
+									OUT_TXT="${OUTPUT_DATA_DIR}out_${JOB_PARAM}.txt"
+										
+									CONF_FILE="${CONF_DIR}/decnet_${JOB_PARAM}.config"
+									
+									sed -e "s|SEED|${SEED}|" \
+										-e "s|NUM_EXP|${NUM_EXP}|" \
+										-e "s|NUM_RUN|${NUM_RUN}|" \
+										-e "s|MAX_TIME|${MAX_TIME}|" \
+										-e "s|AGENT_TYPE|${AGENT_TYPE}|" \
+										-e "s|ACC_MEAN|${ACC_MEAN}|" \
+										-e "s|ACC_STD_DEV|${ACC_STD_DEV}|" \
+										-e "s|ACC_TRUNCATED|${ACC_TRUNCATED}|" \
+										-e "s|DEC_MODEL|${DEC_MODEL}|" \
+										-e "s|UPDATE_CONF|${UPDATE_CONF}|" \
+										-e "s|BELIEF_EPSILON|${BELIEF_EPSILON}|" \
+										-e "s|FINITE_TIME_EXPONENT|${FINITE_TIME_EXPONENT}|" \
+										-e "s|DDM_DRIFT_DIST|${DDM_DRIFT_DIST}|" \
+										-e "s|DDM_BASE_DRIFT|${DDM_BASE_DRIFT}|" \
+										-e "s|DDM_DRIFT_RANGE_PLUS_MINUS|${DDM_DRIFT_RANGE_PLUS_MINUS}|" \
+										-e "s|DDM_DRIFT_STD_DEV|${DDM_DRIFT_STD_DEV}|" \
+										-e "s|DDM_NOISE_STD_DEV|${DDM_NOISE_STD_DEV}|" \
+										-e "s|THRESHOLD|${THRESHOLD}|" \
+										-e "s|INTERR_TIME|${INTERR_TIME}|" \
+										-e "s|INTERR_TIME_FROM_ACC|${INTERR_TIME_FROM_ACC}|" \
+										-e "s|PRIOR_DIST|${PRIOR_DIST}|" \
+										-e "s|COST_MATRIX_T|${COST_MATRIX_T}|" \
+										-e "s|COST_MATRIX_E|${COST_MATRIX_E}|" \
+										-e "s|NET_TYPE|${NET_TYPE}|" \
+										-e "s|NUM_NODES|${NUM_NODES}|" \
+										-e "s|NUM_EDGES|${NUM_EDGES}|" \
+										-e "s|LINK_PROBABILITY|${LINK_PROBABILITY}|" \
+										-e "s|INTERACTION_RADIUS|${INTERACTION_RADIUS}|" \
+										-e "s|DYNAMICNET|${DYNAMICNET}|" \
+										-e "s|AGENTSPEED|${AGENTSPEED}|" \
+										-e "s|MOVE_STDDEV|${MOVE_STDDEV}|" \
+										-e "s|ENV_SIZE|${ENV_SIZE}|" \
+										-e "s|PERIODICBOUND|${PERIODICBOUND}|" \
+										-e "s|OUT_TXT|${OUT_TXT}|" \
+											${TEMPLATE_SETTINGS} > ${CONF_FILE}
 												
-												NET_PAR=0
-												if [ "$NET_TYPE" == "erdos-renyi" ];
-												then
-													NET_PAR=${LINK_PROBABILITY}
-												fi
-												if [ "$NET_TYPE" == "barabasi-albert" ];
-												then
-													NET_PAR=${NUM_EDGES}
-												fi
-												if [ "$NET_TYPE" == "space" ];
-												then
-													NET_PAR=${INTERACTION_RADIUS}
-												fi
-												
-												JOB_PARAM="net-${NET_TYPE}_nodes-${NUM_NODES}_link-${NET_PAR}_bound-${PERIODICBOUND}_env-${ENV_SIZE}_model-${DEC_MODEL}_up-${UPDATE_CONF}_acc-${ACC_MEAN}_acstdv-${ACC_STD_DEV}_eps-${BELIEF_EPSILON}"
-												OUT_TXT="${OUTPUT_DATA_DIR}out_${JOB_PARAM}.txt"
-													
-												CONF_FILE="${CONF_DIR}/decnet_${JOB_PARAM}.config"
-												
-												sed -e "s|SEED|${SEED}|" \
-													-e "s|NUM_EXP|${NUM_EXP}|" \
-													-e "s|NUM_RUN|${NUM_RUN}|" \
-													-e "s|MAX_TIME|${MAX_TIME}|" \
-													-e "s|AGENT_TYPE|${AGENT_TYPE}|" \
-													-e "s|ACC_MEAN|${ACC_MEAN}|" \
-													-e "s|ACC_STD_DEV|${ACC_STD_DEV}|" \
-													-e "s|ACC_TRUNCATED|${ACC_TRUNCATED}|" \
-													-e "s|DEC_MODEL|${DEC_MODEL}|" \
-													-e "s|UPDATE_CONF|${UPDATE_CONF}|" \
-													-e "s|BELIEF_EPSILON|${BELIEF_EPSILON}|" \
-													-e "s|FINITE_TIME_EXPONENT|${FINITE_TIME_EXPONENT}|" \
-													-e "s|DDM_DRIFT_DIST|${DDM_DRIFT_DIST}|" \
-													-e "s|DDM_BASE_DRIFT|${DDM_BASE_DRIFT}|" \
-													-e "s|DDM_DRIFT_RANGE_PLUS_MINUS|${DDM_DRIFT_RANGE_PLUS_MINUS}|" \
-													-e "s|DDM_DRIFT_STD_DEV|${DDM_DRIFT_STD_DEV}|" \
-													-e "s|DDM_NOISE_STD_DEV|${DDM_NOISE_STD_DEV}|" \
-													-e "s|THRESHOLD|${THRESHOLD}|" \
-													-e "s|INTERR_TIME|${INTERR_TIME}|" \
-													-e "s|NET_TYPE|${NET_TYPE}|" \
-													-e "s|NUM_NODES|${NUM_NODES}|" \
-													-e "s|NUM_EDGES|${NUM_EDGES}|" \
-													-e "s|LINK_PROBABILITY|${LINK_PROBABILITY}|" \
-													-e "s|INTERACTION_RADIUS|${INTERACTION_RADIUS}|" \
-													-e "s|DYNAMICNET|${DYNAMICNET}|" \
-													-e "s|AGENTSPEED|${AGENTSPEED}|" \
-													-e "s|MOVE_STDDEV|${MOVE_STDDEV}|" \
-													-e "s|ENV_SIZE|${ENV_SIZE}|" \
-													-e "s|PERIODICBOUND|${PERIODICBOUND}|" \
-													-e "s|OUT_TXT|${OUT_TXT}|" \
-														${TEMPLATE_SETTINGS} > ${CONF_FILE}
-															
-												export PYTHONPATH=${PYPATH}
-												#COMMAND="python3 ${EXEC_FILE} ${CONF_FILE}"
-												#COMMAND="./run_job.sh ${EXEC_FILE} ${CONF_FILE}"
-												COMMAND="qsub run_job.sh ${EXEC_FILE} ${CONF_FILE}"
-												#${COMMAND}
-												while ! ${COMMAND}
-												do
-													sleep 2
-												done
-													
-												COUNT=$((COUNT + 1))
-											done
-										done
+									export PYTHONPATH=${PYPATH}
+									#COMMAND="python3 ${EXEC_FILE} ${CONF_FILE}"
+									#COMMAND="./run_job.sh ${EXEC_FILE} ${CONF_FILE}"
+									COMMAND="qsub run_job.sh ${EXEC_FILE} ${CONF_FILE}"
+									#${COMMAND}
+									while ! ${COMMAND}
+									do
+										sleep 2
 									done
+										
+									COUNT=$((COUNT + 1))
 								done
 							done
 						done
diff --git a/src/AsynchKicks/DecisionNet.py b/src/AsynchKicks/DecisionNet.py
@@ -80,14 +80,17 @@ def setDDMAgent(self, driftDistribution, baseDrift, noiseStdDev, threshold, prio
             self.accuracyStdDev = args[1]
         if useBayesRisk:
             costMatrix = args[-1]
-            def compute_BR_opt_thresh(thresh): # from Eq. (5.6) of Bogacz et al. Psy.Rev. 2006
-                return costMatrix[1]/costMatrix[0] * 2 * (meanDrift**2) / (self.noiseStdDev**2) - 4 * meanDrift * thresh / (self.noiseStdDev**2) + np.exp( - 2 * meanDrift * thresh / (self.noiseStdDev**2) ) -  np.exp( 2 * meanDrift * thresh / (self.noiseStdDev**2) )            
-            self.threshold = scipy.optimize.fsolve( compute_BR_opt_thresh, 0.25*meanDrift*costMatrix[1]/costMatrix[0], maxfev=1000 )[0] # second parameter is the starting point which is set to the approximation of Eq. (5.7) of Bogacz et al. Psy.Rev. 2006 
-            if (self.DEBUG): print("Computed threshold with BR cost matrix is: " + str(self.threshold) )
-            #if (self.DEBUG): print("CHECK: " + str( compute_BR_opt_thresh(self.threshold)) )
-            integrResult = integrate.quad( self.integrationFuncDdmThreshNormal, -np.inf, np.inf, args=( costMatrix, self.noiseStdDev ) )
-            self.threshold = integrResult[0]
-            if (self.DEBUG): print("Computed (with integral with error: " + str(integrResult[1]) + ") threshold with BR cost matrix is: " + str(self.threshold) )
+            integrateOverAllValues = True
+            if integrateOverAllValues:
+                integrResult = integrate.quad( self.integrationFuncDdmThreshNormal, -np.inf, np.inf, args=( costMatrix, self.noiseStdDev ) )
+                self.threshold = integrResult[0]
+                if (self.DEBUG): print("Computed (with integral with error: " + str(integrResult[1]) + ") threshold with BR cost matrix is: " + str(self.threshold) )
+            else:
+                def compute_BR_opt_thresh(thresh): # from Eq. (5.6) of Bogacz et al. Psy.Rev. 2006
+                    return costMatrix[1]/costMatrix[0] * 2 * (meanDrift**2) / (self.noiseStdDev**2) - 4 * meanDrift * thresh / (self.noiseStdDev**2) + np.exp( - 2 * meanDrift * thresh / (self.noiseStdDev**2) ) -  np.exp( 2 * meanDrift * thresh / (self.noiseStdDev**2) )            
+                self.threshold = scipy.optimize.fsolve( compute_BR_opt_thresh, 0.25*meanDrift*costMatrix[1]/costMatrix[0], maxfev=1000 )[0] # second parameter is the starting point which is set to the approximation of Eq. (5.7) of Bogacz et al. Psy.Rev. 2006 
+                if (self.DEBUG): print("Computed threshold with BR cost matrix is: " + str(self.threshold) )
+                #if (self.DEBUG): print("CHECK: " + str( compute_BR_opt_thresh(self.threshold)) )
         else:
             self.threshold = threshold
         
diff --git a/src/DecNet/DecisionProcess.py b/src/DecNet/DecisionProcess.py
@@ -9,7 +9,11 @@
 import os
 import configparser
 import json
+import numpy as np
 import numpy.random as rand
+import scipy.optimize
+import scipy.integrate as integrate
+from scipy.stats import norm
 import math
 # import matplotlib.pyplot as plt
 # from statsmodels.graphics import plottools
@@ -310,6 +314,17 @@
                 args.append( accuracyStdDev )
             if interrogationTimeFromAccuracy:
                 args.append(costMatrix)
+                def integrationFuncDdmThreshNormal(drift, costMatrix, noiseStdDev, baseDrift, driftStdDev):
+                    prob = norm.pdf(drift, baseDrift, driftStdDev)
+                    drift = abs(drift)
+                    def compute_BR_opt_thresh(thresh): # from Eq. (5.6) of Bogacz et al. Psy.Rev. 2016
+                        return costMatrix[1]/costMatrix[0] * 2 * (drift**2) / (noiseStdDev**2) - 4 * drift * thresh / (noiseStdDev**2) + np.exp( - 2 * drift * thresh / (noiseStdDev**2) ) -  np.exp( 2 * drift * thresh / (noiseStdDev**2) )             
+                    pdf = scipy.optimize.fsolve( compute_BR_opt_thresh, 0.25*drift*costMatrix[1]/costMatrix[0], maxfev=1000 )[0] # second parameter is the starting point which is set to the approximation of Eq. (5.7) of Bogacz et al. Psy.Rev. 2016
+                    return pdf*prob
+                integrResult = integrate.quad( integrationFuncDdmThreshNormal, -np.inf, np.inf, args=( costMatrix, noiseStdDev, baseDrift, driftStdDev ) )
+                thresh = integrResult[0]
+                expectedAccuracy = 1 - ( 1.0/ ( 1 + np.exp(2*baseDrift*thresh/(noiseStdDev**2)) ))
+                interrogationTime = expectedAccuracy
             decNet.setDDMAgent(driftDistribution, baseDrift, noiseStdDev, interrogationTimeFromAccuracy, interrogationTime, prior, args)
         if (DEBUG): print("Initialising agents")
         agentParams = []
@@ -367,3 +382,4 @@
     outFile.close()   
     if DEBUG: print("Process Ended")
     
+    
