Rhesus Monkey Young and Aged L3 PFC Pyramidal Neurons (Rumbell et al. 2016)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:184497
A stereotypical pyramidal neuron morphology with ion channel parameter combinations that reproduce firing patterns of one young and one aged rhesus monkey L3 PFC pyramidal neurons. Parameters were found through an automated optimization method.
Reference:
1 . Rumbell TH, Draguljic D, Yadav A, Hof PR, Luebke JI, Weaver CM (2016) Automated evolutionary optimization of ion channel conductances and kinetics in models of young and aged rhesus monkey pyramidal neurons. J Comput Neurosci 41:65-90 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neocortex L2/3 pyramidal GLU cell;
Channel(s): I Na,p; I Na,t; I A; I K; I M; I h; I K,Ca; I Sodium; I Calcium; I Potassium; I_AHP; I Cl, leak;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Ion Channel Kinetics; Parameter Fitting; Detailed Neuronal Models; Aging/Alzheimer`s;
Implementer(s):
Search NeuronDB for information about:  Neocortex L2/3 pyramidal GLU cell; I Na,p; I Na,t; I A; I K; I M; I h; I K,Ca; I Sodium; I Calcium; I Potassium; I_AHP; I Cl, leak;
//params and set up for DE

pop_size=256				// size of the population.	
Cross_prob=0.9				// probability to crossover
Num_Total_Gen = 1			// Total number of generations to run
Best=0						// curent best individual
Elapsed=0					// Total time of run
F=0.9						// scale factor for differential
rand_max_pop=16 // maximum number of parameter sets to run in random mode
max_time=144000		// 40 hours // maximum time to run for, in seconds
longest_gen=0				// hold longest generation so far, in s

DEBUG=0	//OUTPUT: 0 - no debugging msgs; 1 - all debugging msgs

objref pop,res,mut_pop,mut_res,pop_vec,mut_pop_vec,ffs,mut_ffs,pool,pool_ffs,errVec
objref pool_ind,transvec,paravec,changed
objref p_file,pc_file,param_file,g_file,fv_file,sa_file,randomData,scratch_file
objref iv_file,fi_file,ivout_file,fiout_file,nsde_ff_combs_file

strdef cmd,tstr,parstr

p_file = new File()
pc_file = new File()
param_file = new File()
g_file = new File()
fv_file = new File()
sa_file = new File()
randomData = new File()
scratch_file = new File()
iv_file = new File()
fi_file = new File()
ivout_file = new File()
fiout_file = new File()
nsde_ff_combs_file = new File()

//open main DE log file
g_file.wopen("../output/de_log.dat")

//open fitnessVal log file
fv_file.wopen("../output/fvevo.dat")

//open sensitivity analysis log file
sa_file.wopen("../output/SA.dat")

//File containing random seeds to use
randomData.ropen("../data/rands.dat")

pop = new Matrix( pop_size , NP )
mut_pop = new Matrix( pop_size , NP )
res = new Vector( pop_size , 1e50 )
mut_res = new Vector( pop_size , 1e50 )
ffs = new Matrix( pop_size , totalFFs )
mut_ffs = new Matrix( pop_size , totalFFs )
pool = new Matrix( ( 2 * pop_size ) , NP )
pool_ffs = new Matrix( ( 2 * pop_size ) , totalFFs )
pool_ind = new Vector()
transvec = new Vector( NP )
changed = new Vector( pop_size , 1 )
errVec = new Vector( totalFFs )

// random numbers and initial boundaries for the population
objref  r[NP],ru[pc.nhost()]

// initialise a random number generator for each parameter
for j=0,(NP-1){
	r[j] = new Random(randomData.scanvar())
	r[j].uniform(minvec.x[j],maxvec.x[j])
}

// initialise a random number generator for each parallel host
for j=0,(pc.nhost()-1){
	ru[j] = new Random(randomData.scanvar())
	ru[j].uniform(0,1)
}

// close the random number file
randomData.close()

LOAD_GENERATION=0	// if zero generates random initial population, if one load generation from file

num_gen=0

//=======================MAIN LOOP FOR DE================================
proc DE() { local i,j,k

	init_pop()
	
	// make file of DE parameters for producing output figures
	param_file.wopen( "output/metaparams.dat" )
	param_file.printf( "%s: %d\n" , "population" , pop_size ) 
	param_file.printf( "%s: %d\n" , "parameters" , NP )
	param_file.printf( "%s: %d\n" , "generators" , NGens )
	param_file.printf( "%s: %d\n" , "totalFFs" , totalFFs )
	param_file.printf( "%s: %d\n" , "totalT" , (Tend.max()/voltvecdt) )
	param_file.printf( "%s: %d\n" , "MRF" , MRFflag )
	param_file.printf( "%s: %g\n" , "dt" , voltvecdt)

	begintime = startsw()
	
	Evaluate_pop(0)
	print_pop()
	num_gen += 1

	//while(num_gen<=Num_Total_Gen){	// termination criterion besed on fixed # Generations
    //while(Best>10){				  // termination criterion besed on target score
	while ( ( Elapsed + ( 2 * longest_gen ) ) < max_time ) {
		startgentime = startsw()
		//Evaluate_pop(0)
       	Mutants()
        Evaluate_pop(1)
       	Merge_pops()
       	gen_time = startsw() - startgentime
       	if ( gen_time > longest_gen ) { longest_gen = gen_time }
       	find_best()
		print_things()
		print_pop()
		num_gen += 1
        Elapsed = startsw() - begintime
	}
	Num_Total_Gen = num_gen-1
	
	param_file.printf( "%s: %d\n" , "generations" , Num_Total_Gen )
	param_file.close()
 	print "-------------------------------end of Differential Evolution----------------------------------"
	for (j=0;j<=NP-1;j+=1) g_file.printf("%10.10f ",pop.x[Best][j])
	g_file.printf("%10.10f ",res.x[Best])
    g_file.printf("\n")
	g_file.close()
	fv_file.close()
}

//=======================MAIN LOOP FOR NSDE================================
proc NSDE() { local i,j,k
	
	init_pop()
	
	// make file of DE parameters for producing output figures
	param_file.wopen( "output/metaparams.dat" )
	param_file.printf( "%s: %d\n" , "population" , pop_size ) 
	param_file.printf( "%s: %d\n" , "parameters" , NP )
	param_file.printf( "%s: %d\n" , "generators" , NGens )
	param_file.printf( "%s: %d\n" , "totalFFs" , totalFFs )
	param_file.printf( "%s: %d\n" , "totalT" , (Tend.max()/voltvecdt) )
	param_file.printf( "%s: %d\n" , "MRF" , MRFflag )
	param_file.printf( "%s: %g\n" , "dt" , voltvecdt)

	begintime=startsw()
	
	Evaluate_pop(0)
	print_pop()
	num_gen += 1

	//while(num_gen<=Num_Total_Gen){	// termination criterion besed on fixed # Generations
    //while(Best>10){				  // termination criterion besed on target score
	while ( ( Elapsed + ( 2 * longest_gen ) ) < max_time ) {
		startgentime = startsw()
		//Evaluate_pop(0)
       	Mutants()
       	Evaluate_pop(1)
       	non_dominated_sorting()
       	Merge_pops_NSDE()
       	gen_time = startsw() - startgentime
       	if ( gen_time > longest_gen ) { longest_gen = gen_time }
       	find_best()
		print_things()
		print_pop()
		num_gen += 1
        Elapsed = startsw() - begintime
	}
	Num_Total_Gen = num_gen-1
	
	param_file.printf( "%s: %d\n" , "generations" , Num_Total_Gen )
	param_file.close()
 	print "-------------------------------end of Differential Evolution----------------------------------"
	for (j=0;j<=NP-1;j+=1) g_file.printf("%10.10f ",pop.x[Best][j])
	g_file.printf("%10.10f ",res.x[Best])
    g_file.printf("\n")
	g_file.close()
	fv_file.close()
}

/***************************************
 * NSDE processes and functions
 **************************************/
// Assign fitness values to individuals based on their level of non-domination
// within the pool of both population and mutants
proc non_dominated_sorting() { local i , j , pool_ind_size	localobj curr_rank_ind
	DEBUG=0
	starttime = startsw()
	
	printf("Non-dominated sorting: start\n")
	
	// set up curr_rank_ind as an empty Vector
	curr_rank_ind = new Vector()
	
	// fill pool and pool_ffs with the values from pop, mut_pop, ffs, mut_ffs
	for ( i = 0 ; i < ( 2 * pop_size ) ; i+=1 ) {
		if ( i < pop_size ) {
			pool.setrow( i , pop.getrow( i ) )
			pool_ffs.setrow( i , ffs.getrow( i ) )
		} else {
			pool.setrow( i , mut_pop.getrow( i - pop_size ) )
			pool_ffs.setrow( i , mut_ffs.getrow( i - pop_size ) )
		}
	}
	if (DEBUG) pool.printf
	if (DEBUG) pool_ffs.printf
	
	// set pool_ind vector to every index in pool
	pool_ind.indgen( 0 , ( ( pop_size * 2 ) - 1 ) , 1 )
	if (DEBUG)pool_ind.printf

	// initialise rank variable
	rank = 1
	
	// start while loop for having individuals left in the pool
	while ( pool_ind.size > 0 ) {
		
		// set the pool size for this iteration
		pool_ind_size = pool_ind.size
		// reset curr_rank_ind
		curr_rank_ind.resize(0)
		if (DEBUG) printf( "NS rank = %d; pool size = %d\n" , rank , pool_ind_size )
		
		// non-domination testing:
		// for each member remaining in the pool:
		for ( i = 0 ; i < pool_ind_size ; i+=1 ) {
			// is this individual dominated by any others in the pool?
			if ( is_dominated( pool_ind.x[ i ] ) == 0 ) {
				curr_rank_ind.append( pool_ind.x[ i ] )
			}
		}
		
		// check if there are any currently non-dominated individuals...
		if ( curr_rank_ind.size == 0 ) {
			// if there are no dominant individuals, set all to current rank
			curr_rank_ind = pool_ind.c
		}
		
		if (DEBUG) curr_rank_ind.printf
		
		// for each currently non-dominated individual, set fitness to rank
		for ( i = 0 ; i < ( pop_size * 2 ) ; i+=1 ) {
			if ( curr_rank_ind.contains( i ) ) {
				if ( i < pop_size ) {
					res.x[ i ] = rank
				} else {
					mut_res.x[ i - pop_size ] = rank
				}
			}
		}
		
		// adjust ranks using a specified method for curr_rank_ind mini-pool
		// Methods:	1 = crowdedness in fitness space
		//			2 = total combined fitness (normalised per FF)
		//adjust_res_for_mini_pool( curr_rank_ind , 1 )
		
		// remove members of mini-pool from pool_ind
		for ( i = 0 ; i < curr_rank_ind.size ; i+=1 ) {
			if ( pool_ind.contains( curr_rank_ind.x[ i ] ) ) {
				pool_ind.remove( pool_ind.indwhere( "==" , curr_rank_ind.x[ i ] ) )
			}
		}
		
		// increment the rank variable
		rank += 1
	}
	
	// actually we want to calculate rankings for e.g. crowdedness using the 
	// entire pool, not just the mini-pools. So here we submit the whole pool 
	// to the process:
	pool_ind.indgen( 0 , ( ( pop_size * 2 ) - 1 ) , 1 )
	adjust_res_for_mini_pool( pool_ind , 1)
	
	if (DEBUG) res.printf
	if (DEBUG) mut_res.printf
	
	DEBUG=0
	printf("Non-dominated sorting: time taken: %f\n", startsw() - starttime)
}

// Check is an individual is dominated by any other members of the pool - takes 
// index i into 'pool' as input, and returns 1 if true and 0 if false
func is_dominated() { local i , j , output_is_dominated
	i = $1
	output_is_dominated = 0
	if (DEBUG) printf( "Domination test %d:\n" , i )
	for ( j = 0 ; j < pool_ind.size ; j+=1 ) {
		if ( dominates( pool_ind.x[ j ] , i ) ) {
			output_is_dominated = 1
		}
	}
	if (DEBUG) printf( "%d is_dominated = %d\n" , i , output_is_dominated )
	return output_is_dominated
}

// Check if an individual dominates another individual - takes indices to two 
// individuals indexed into 'pool' as two input integers, and returns 1 if $1 
// dominates $2 and 0 if not
// NSDE_COMB_FFS tells us whether we want to look up the file 'nsde_ff_combs.dat'
// to get sets of FFs that should be summed prior to determining dominance.
func dominates() { local i , j , better , worse , output_dominates , num_groups	, \
				   groupsum1 , groupsum2	\
				   localobj groups_vec
	better = 0
	worse = 0
	DEBUG=0
	// better is if $1 < $2 - this direction is caused by use of 'error' func.
	// i.e. error should be minimised
	if ( NSDE_COMB_FFS == 0 ) {
		for ( i = 0 ; i < totalFFs ; i+=1 ) {
			if ( pool_ffs.x[ $1 ][ i ] < pool_ffs.x[ $2 ][ i ] ) {
				better+=1
			} else if ( pool_ffs.x[ $1 ][ i ] > pool_ffs.x[ $2][ i ] ) {
				worse+=1
			}
		}
	} else if ( NSDE_COMB_FFS == 1 ) {
		if (DEBUG) printf( "Combining FFs into groups...\n" )
		nsde_ff_combs_file.ropen("setup/nsde_ff_combs.dat")
		groups_vec = new Vector()
		groups_vec.scanf( nsde_ff_combs_file )
		num_groups = groups_vec.max
		if (DEBUG) printf( "num_groups = %d\n" , num_groups )
		for ( j = 1 ; j <= num_groups ; j+=1 ) {
			if (DEBUG) printf( "group %d:\n" , j )
			groupsum1 = 0
			groupsum2 = 0
			for ( i = 0 ; i < totalFFs ; i+=1 ) {
				if ( groups_vec.x[ i ] == j ) {
					if (DEBUG) printf( "FF %d in group %d: adding %f to group 1 and %f to group 2\n" , \
										i , j , pool_ffs.x[ $1 ][ i ] , pool_ffs.x[ $2 ][ i ] )
					groupsum1 += pool_ffs.x[ $1 ][ i ]
					groupsum2 += pool_ffs.x[ $2 ][ i ]
				}
			}
			if (DEBUG) printf( "groupsum1: %f : groupsum2 : %f\n" , groupsum1 , groupsum2 )
			if ( groupsum1 < groupsum2 ) {
				better += 1
			} else if ( groupsum1 > groupsum2 ) {
				worse += 1
			}
		}
	}
	
	if ( ( worse == 0 ) && ( better > 0 ) ) {
		output_dominates = 1
	} else {
		output_dominates = 0
	}
	if (DEBUG) printf( "Dominates : %d vs %d : better = %d ; worse = %d : dominates = %d\n" , \
								$1 , $2 , better , worse , output_dominates )
	DEBUG=0
	return output_dominates
}

// Calculate the crowdedness of each individual in a pool_ffs, with the indices 
// into the 'pool_ffs' object provided as argument. Returns vector of crowdedness.
obfunc simple_crowdedness() { 	local i , j , curr_val, above , above_ind , below , \
						below_ind , distance	localobj crowdedness_vec , \
						ind_pool_ffs , ind_dest , pool_ffs_col , ind_pool_ffs_col , \
						ind_pool_ffs_col_sorted_index
	ind_pool_ffs = $o1
	if (DEBUG) ind_pool_ffs.printf
	ind_dest = new Vector( ind_pool_ffs.size )
	if ( ind_dest.size == 1 ) {
		ind_dest.x[ 0 ] = 0
	} else {
		ind_dest.indgen( 0 , ( ind_pool_ffs.size - 1 ) , 1 )
	}
	if (DEBUG) ind_dest.printf
	crowdedness_vec = new Vector( ind_pool_ffs.size )
	pool_ffs_col = new Vector( pop_size * 2 )
	ind_pool_ffs_col = new Vector( ind_pool_ffs.size )
	ind_pool_ffs_col_sorted_index = new Vector( ind_pool_ffs.size )
	
	// simple crowdedness metric is calculated according to:
	// 	for an individual, for an objective function:
	//		find the adjacent objective function values in the rest of the pool_ffs
	//		subtract the value one above from the value one below
	//		add this to the total crowdedness score
	//		lower is worse
	//		if there is no bigger or smaller value in the pool_ffs, 
	//			crowdedness score gets an arbitrarily large value added to it
	
	for ( i = 0 ; i < ind_pool_ffs.size ; i+=1 ) {
		crowdedness_vec.x[ i ] = 0
		for ( j = 0 ; j < totalFFs ; j+=1 ) {
			// get a column from pool_ffs of the current ff values (j), 
			// according to the indices in ind
			pool_ffs_col = pool_ffs.getcol( j )
			if (DEBUG) pool_ffs_col.printf
			// copy members of the ff col that are part of the source pool into
			// the target pool at indices given by ind_dest,0:size of the vector-1
			ind_pool_ffs_col.copy( pool_ffs_col , ind_pool_ffs , ind_dest )
			if (DEBUG) ind_pool_ffs_col.printf
			// scale so that all distances are normalised - so each ff is the same
			ind_pool_ffs_col.scale( 0 , 1 )
			if (DEBUG) ind_pool_ffs_col.printf
			// current value for finding crowdedness
			//curr_val = ind_pool_ffs_col.x[ i ]
			// sort the pool_ffs so we can find one above and one below
			// get the sorted indices so that we can find our original value
			ind_pool_ffs_col_sorted_index = ind_pool_ffs_col.sortindex()
			ind_of_orig = ind_pool_ffs_col_sorted_index.x[ i ]
			if (DEBUG) ind_pool_ffs_col.printf
			// check if we're at the bottom of the pile:
			if ( ind_of_orig == 0 ) {
				below_ind = -1
			} else {
				below_ind = ind_pool_ffs_col_sorted_index.x[ i - 1 ]
			}
			//printf("i=%d j=%d\n",i,j)
			// also check if we're at the top of the pile:
			if ( ind_of_orig == ( ind_pool_ffs.size - 1 ) ) {
				above_ind = -1
			} else {
				//printf("made it \n")
				above_ind = ind_pool_ffs_col_sorted_index.x[ i - 1 ]
			}
			// if ind_of_orig is at the top or bottom, add a big number
			if ( ( above_ind == -1 ) || ( below_ind == -1 ) ) {
				crowdedness_vec.x[ i ] += totalFFs
			} else {
				// else find how crowded it is
				above = ind_pool_ffs_col.x[ above_ind ]
				below = ind_pool_ffs_col.x[ below_ind ]
				crowdedness_vec.x[ i ] += ( above - below )
			}
		}
	}
	if (DEBUG) crowdedness_vec.printf
	// scale the crowdedness values so that they don't change rank by too much
	if ( crowdedness_vec.min == crowdedness_vec.max ) {
		crowdedness_vec.fill( 1 )
	} else {
		crowdedness_vec.scale( 	( 1 - ( 1 / ( ( 2 * pop_size ) + 1 ) ) ) , \
								( 1 + ( 1 / ( ( 2 * pop_size ) + 1 ) ) ) )
	}
	if (DEBUG) crowdedness_vec.printf
	return crowdedness_vec
}

// Calculate the total fitness of each individual in a pool_ffs, with the 
// indices into the 'pool_ffs_ffs' object provided as argument. Returns vector 
// of total fitnesses
obfunc total_fitness() { 	local i , j , curr_val, above , above_ind , below , \
						below_ind , distance	localobj total_fitness_vec , \
						ind_pool_ffs , ind_dest , pool_ffs_col , ind_pool_ffs_col
	ind_pool_ffs = $o1
	if (DEBUG) ind_pool_ffs.printf
	ind_dest = new Vector( ind_pool_ffs.size )
	if ( ind_dest.size == 1 ) {
		ind_dest.x[ 0 ] = 0
	} else {
		ind_dest.indgen( 0 , ( ind_pool_ffs.size - 1 ) , 1 )
	}
	if (DEBUG) ind_dest.printf
	total_fitness_vec = new Vector( ind_pool_ffs.size )
	total_fitness_vec.fill( 0 )
	pool_ffs_col = new Vector( pop_size * 2 )
	ind_pool_ffs_col = new Vector( ind_pool_ffs.size )
	
	// total fitness metric is calculated according to:
	// 	for an FF:
	//		normalise the FF range to 0:1
	//		for an individual:
	//			add the current FF score to his total fitness score
	
	for ( j = 0 ; j < totalFFs ; j+=1 ) {
		// get a column from pool_ffs of the current FF values (j), 
		// according to the indices in ind
		pool_ffs_col = pool_ffs.getcol( j )
		if (DEBUG) pool_ffs_col.printf
		ind_pool_ffs_col.copy( pool_ffs_col , ind_pool_ffs , ind_dest )
		if (DEBUG) ind_pool_ffs_col.printf
		// scale so that the FF is normalised, with better values as bigger numbers
		ind_pool_ffs_col.scale( 1 , 0 )
		// and set to 0 if the scale didn't work
		if ( ind_pool_ffs_col.min == ind_pool_ffs_col.max ) {
			ind_pool_ffs_col.fill( 0 )
		}
		if (DEBUG) ind_pool_ffs_col.printf
		// for each individual:
		for ( i = 0 ; i < ind_pool_ffs.size ; i+=1 ) {
			total_fitness_vec.x[ i ] += ind_pool_ffs_col.x[ i ]
		}
	}
	
	if (DEBUG) total_fitness_vec.printf
	// check if all fitness values are identical, in which case, fine, set to 1
	// scale fitness values in a range appropriate to the amount of potential
	// non-domination ranks that might be present, according to population size
	if ( total_fitness_vec.min == total_fitness_vec.max ) {
		total_fitness_vec.fill( 1 )
	} else {
		total_fitness_vec.scale( 	( 1 - ( 1 / ( ( 2 * pop_size ) + 1 ) ) ) , \
								( 1 + ( 1 / ( ( 2 * pop_size ) + 1 ) ) ) )
	}
	if (DEBUG) total_fitness_vec.printf
	return total_fitness_vec
}

// Adjust res and mut_res values based on the crowdedness of a mini-pool with 
// the indices into 'pool' for the mini-pool supplied as an argument
proc adjust_res_for_mini_pool() { local i , j , method	localobj ind , modifier
	
	ind = $o1
	method = $2
	
	// call appropriate method here - modifier should be high for strong indivs
	if ( method == 1 ) {	
		modifier = simple_crowdedness( ind )
	} else if ( method == 2 ) {
		modifier = total_fitness( ind )
	}
	
	j = 0
	
	for ( i = 0 ; i < ( pop_size * 2 ) ; i+=1 ) {
		if ( ind.contains( i ) ) {
			// use divide by modifier because modifier of greater than 1 
			// indicates better result, and res.x[] values should be low
			if ( i < pop_size ) {
				res.x[ i ] = ( res.x[ i ] / modifier.x[ j ] )
			} else {
				mut_res.x[ i - pop_size ] = ( mut_res.x[ i - pop_size ] / \
														modifier.x[ j ] )
			}
			j+=1
		}
	}
}

/************************************
 * DE processes and functions
 ***********************************/

// Main DE function to generate a new mutant population using the 
// standard DE algorithm
proc Mutants() {local i,j,inside,outside,starttime,x,y,z,num_muts_gend

	starttime = startsw()
	num_muts_gend=0
	printf("Mutants: start\n")
	//parallel version isn't faster unless > 1000 cpus are available...
	//even then it is barely faster...
	// serial version:

	if (pc.nhost<1000) {
		if (DEBUG) printf("Mutants: made it, serial\n")
		for (i=0;i<pop_size;i+=1) {
			if (DEBUG) printf("Population member: %d\n",i)
			inside = 0 
			while ( inside == 0 ) {
				mutate( i )
				num_muts_gend += 1
				if (DEBUG) printf("    Generated a mutant\n")
				if (DEBUG) printf("    Generated %d mutants so far\n",num_muts_gend)
				outside = 0
				// loop through each parameter and check that none are outside boundaries
				for ( j = 0; j < NP; j += 1 ) {
					if (DEBUG) printf("param: %f; max: %f; min: %f\n",mut_pop.x[i][j],maxvec.x[j],minvec.x[j])
					if ( mut_pop.x[i][j] < minvec.x[j] ) {
						outside = 1
						if (DEBUG) printf("    Parameter %d outside lower boundary\n",j)
					} else if ( mut_pop.x[i][j] > maxvec.x[j] ) {
						outside = 1
						if (DEBUG) printf("    Parameter %d outside upper boundary\n",j)
					}
				}
				if ( outside == 0 ) {
					inside = 1
					if (DEBUG) printf("All parameters inside boundaries!\n")
				}
			}
		} // end of pop_size for loop (i)
	} else { // bulletin board version:

		send_pop_from_master()

		for (i=0;i<pop_size;i+=1) {
			pc.submit("mutPar", i)
		}
		while ((id = pc.working()) != 0) {
			pc.look_take(id)
			ind = pc.upkscalar()
			transvec = pc.upkvec()
			for (j=0;j<NP;j+=1) {
				mut_pop.x[ind][j] = transvec.x[j]
			}
		}
	}

	printf("Mutants: mutants generated: %d\n", num_muts_gend)
	printf("Mutants: time taken: %f\n", startsw() - starttime)
}

proc mutate() { local i, j, r0, r1, r2, jrand
	r0 = $1	
	r1 = $1	
	r2 = $1	

	while ( r0 == $1 ) {
		r0 = int( ru[0].repick() * pop_size )
	}
	while ( r1 == $1 ) {
		r1 = int( ru[0].repick() * pop_size )
	}
	while ( r2 == $1 ) {
		r2 = int( ru[0].repick() * pop_size )
	}

	// pick parameter to make sure at least one is mutated
	jrand = int( ru[0].repick()*NP )

	// loop through each parameter and generate the mutant vector
	for ( j = 0; j < NP; j += 1 ) {
		if ( ru[0].repick() <= Cross_prob || j == jrand ) {
			mut_pop.x[$1][j] = ( pop.x[r0][j] + ( F * ( pop.x[r1][j] - pop.x[r2][j] ) ) )
		} else {
			mut_pop.x[$1][j] = pop.x[$1][j]
		}
	}
}

// pair of functions to transfer the pop matrix  
// to the workers by encoding as a vector
func get_pop_from_master() { local i, j
  for i=0, pop.nrow()-1 {
    for j=0, pop.ncol()-1 {
      pop.x[i][j] = $o1.x[i*pop.ncol() + j]
    }
  }
  return 0
}

func send_pop_from_master() { local i, j
  pop_vec = new Vector(pop.nrow()*pop.ncol())
  for i=0, pop.nrow()-1 {
    for j=0, pop.ncol()-1 {
      pop_vec.x[i*pop.ncol() + j] = pop.x[i][j]
    }
  }
  pc.context("get_pop_from_master", pop_vec)
  return 0
}

// pair of functions to transfer the pop matrix  
// to the workers by encoding as a vector
func get_mut_pop_from_master() { local i, j
  for i=0, mut_pop.nrow()-1 {
    for j=0, mut_pop.ncol()-1 {
      mut_pop.x[i][j] = $o1.x[i*mut_pop.ncol() + j]
    }
  }
  return 0
}

func send_mut_pop_from_master() { local i, j
  mut_pop_vec = new Vector(mut_pop.nrow()*mut_pop.ncol())
  for i=0, mut_pop.nrow()-1 {
    for j=0, mut_pop.ncol()-1 {
      mut_pop_vec.x[i*mut_pop.ncol() + j] = mut_pop.x[i][j]
    }
  }
  pc.context("get_pop_from_master", mut_pop_vec)
  return 0
}

func mutPar(){local ind,key,r0,r1,r2,i,j	localobj mutvec
	key = hoc_ac_
    ind = $1

	mutvec = new Vector(NP)

	//pick base indices for the mutant to make the while loops work:
	r0=ind
	r1=ind
	r2=ind
	while (r0==ind) {
		r0 = int(ru[pc.id].repick()*pop_size)
	}
	while (r1==ind || r1==r0) {
		r1 = int(ru[pc.id].repick()*pop_size)
	}
	while (r2==ind || r2==r0 || r2==r1) {
		r2 = int(ru[pc.id].repick()*pop_size)
	}
	//pick parameter to make sure at least one comes from mutant 
	jrand=int(ru[pc.id].repick()*NP)

	// loop through each parameter and generate the mutant vector
	for (j=0;j<NP;j+=1)	{
		if (ru[pc.id].repick()<=Cross_prob || j==jrand) {
			mutvec.x[j]=(pop.x[r0][j]+(F*(pop.x[r1][j]-pop.x[r2][j])))
			if (mutvec.x[j]<=minvec.x[j]) {//if mut param is -ve...
				mutvec.x[j]=minvec.x[j] //...make it +ve again
			} else if (mutvec.x[j]>maxvec.x[j]) {//if mut param is huge
				mutvec.x[j]=maxvec.x[j] //reign it in
			}
		} else {
			mutvec.x[j]=pop.x[ind][j]
		}
	}//end of NP for loop (j)

	pc.pack(ind)
	pc.pack(mutvec)
    pc.post(key)
    return key
}

// For every index in the population, checks which is fitter out of 
// the parent and the mutant - saves the fitter one in the pop matrix
proc Merge_pops() {local i
	for ( i = 0 ; i < pop_size ; i += 1 ) {
		//printf("Merging : %d : pop: %f vs %f mut\n",i,res.x[i],mut_res.x[i])
		if ( mut_res.x[ i ] < res.x[ i ] ) {
			res.x[ i ] = mut_res.x[ i ]
			ffs.setrow( i , mut_ffs.getrow( i ) )
			for ( j = 0 ; j <= NP - 1 ; j += 1 ) {
				pop.x[ i ][ j ] = mut_pop.x[ i ][ j ]
			}
			changed.x[ i ] = 1
		} else {
			changed.x[ i ] = 0
		}
	}
}

// Combine the two pools, pop and mut, and then take top pop_size members of 
// this combined pool to be the new pop
proc Merge_pops_NSDE() { local i , ind	localobj combPops_res , combPops_sort_ind , tempPop , tempPop_res , tempPop_ffs
	// temporary vectors for holding things while sorting
	combPops_res = new Vector()
	combPops_sort_ind = new Vector()
	tempPop = new Matrix( pop_size , NP )
	tempPop_res = new Vector()
	tempPop_ffs = new Matrix( pop_size , totalFFs )
	
	// combine the results lists into one big one and then sort them all
	combPops_res.append( res )
	combPops_res.append( mut_res )
	combPops_res.sortindex( combPops_sort_ind )
	
	// now we're taking the top pop_size members of the population, and putting 
	// them into one new population
	for ( i = 0 ; i < pop_size ; i += 1 ) {
		ind = combPops_sort_ind.x[ i ]
		if ( ind < pop_size ) {
			tempPop_res.append( res.x[ ind ] )
			tempPop_ffs.setrow( i , ffs.getrow( ind ) )
			tempPop.setrow( i , pop.getrow( ind ) )
		} else {
			tempPop_res.append( mut_res.x[ ( ind - pop_size ) ] )
			tempPop_ffs.setrow( i , mut_ffs.getrow( ( ind - pop_size ) ) )
			tempPop.setrow( i , mut_pop.getrow( ( ind - pop_size ) ) )
		}
	}
	pop = tempPop
	res = tempPop_res
	ffs = tempPop_ffs
}

//---------------------------------initialpopulation, random or from file -----------------	
proc init_pop(){local i,j

	if(LOAD_GENERATION==0){
		for (i=0;i<=pop_size-1;i+=1){
			for (j=0;j<NP;j+=1) {
				pop.x[i][j]=r[j].repick()
			}
		}
	} else {
		p_file.ropen("output/curr_population")
		pop.scanf(p_file)
		res.scanf(p_file)
		p_file.close()
	}
}

proc print_pop(){local i , j
	for ( i = 0 ; i <= pop_size-1 ; i += 1 ) {
		fv_file.printf( "%d %d " , num_gen , i )
		for ( j = 0 ; j <= NP-1 ; j += 1 ) {
			fv_file.printf( "%1.16f " , pop.x[ i ][ j ] )
		}
		errVec = ffs.getrow( i )
		for ( j = 0 ; j <= totalFFs-1 ; j += 1 ) {
			fv_file.printf( "%g " , ffs.x[ i ][ j ] )
		}
		if (MULOBJ) {
			fv_file.printf( "%g\n" , errVec.sum() )
		} else {
			fv_file.printf( "%g\n" , res.x[ i] )
		}
		fv_file.flush()
	}
}

// population being evaluated is either the parent population or the 
// population of mutants - input variable 0 is parents, 1 is mutants
proc Evaluate_pop(){local i,j,k,ind,fitVal,which_pop,starttime

	starttime = startsw()

	which_pop = $1
	
	printf("Evaluating pop!: %d %d\n", which_pop, num_gen)
	
	if (pc.nhost == 1) {    // use the serial form
		for (i=0;i<=pop_size-1;i+=1){
			if (which_pop == 0) {
				for(j=0;j<=NP-1;j+=1) {
					transvec.x[j]=pop.x[i][j]
				}
				res.x[i]=tfunk()
			} else {
				for(j=0;j<=NP-1;j+=1) {
					transvec.x[j]=mut_pop.x[i][j]
				}
				mut_res.x[i]=tfunk()
			}
		}
	} else {                  // use the bulletin board form
		for (i=0;i<=pop_size-1;i+=1){
			if (which_pop == 0) {
				for(j=0;j<=NP-1;j+=1) {
					transvec.x[j]=pop.x[i][j]
				}
			} else {
				for(j=0;j<=NP-1;j+=1) {
					transvec.x[j]=mut_pop.x[i][j]
				}
			}
			pc.submit( "tfunkpar" , i , transvec )
		}
		while ((id = pc.working()) != 0) {
			pc.look_take(id)
			if (DEBUG) printf("made it into while\n")
			ind = pc.upkscalar()
			if (DEBUG) printf("ind=%d\n",ind)
			fitVal = pc.upkscalar()
			if (DEBUG) printf("made it past upkscalars\n")
			errVec = pc.upkvec()
			// normalise FFs if desired
			if ( NORMFFS ) {
				xopen( "setup/normffs.hoc" )
				fitVal = errVec.sum
			}
			if ( DEBUG ) printf( "collected %d error values\n" , errVec.size() )
			if (which_pop == 0) {
				res.x[ind]=fitVal
				ffs.setrow( ind , errVec )
			} else {
				mut_res.x[ind]=fitVal
				mut_ffs.setrow( ind , errVec )
			}
			if (DEBUG) printf("made it past set fitVal\n")
		}
	}
	if (DEBUG) printf("made it out of evaluate pop...\n")
	printf("Pop evaluated!: %d %d\n", which_pop, num_gen)
	printf("EvaluatePop: time taken = %f\n", startsw() - starttime)
}

func tfunk(){local k,fitnessVal
	fitnessVal=0
	
	//set parameters
	sprint( xstr , "set_n_params(" )
	for ( k = 0 ; k < NP ; k += 1 ) {
		if ( k == 0 ) {
			sprint( xstr , "%stransvec.x[%d]" , xstr , k )
		} else {
			sprint( xstr , "%s,transvec.x[%d]" , xstr , k )
		}
	}
	sprint( xstr , "%s)" , xstr )
	execute(xstr)

	//adjust conductances
	set_conds()
	//adjust kinetics
	set_kins()
	
	MRF.p.run()
	
	fitnessVal=MRF.p.pf.errval
	errVec = get_error_values()
	if ( NORMFFS ) {
		xopen( "setup/normffs.hoc" )
		fitnessVal = errVec.sum
	}
	
	fv_file.printf("%d %d ",1,1)
	for ( j = 0 ; j <= NP-1 ; j += 1 ) {
		fv_file.printf( "%1.16f " , transvec.x[ j ] )
	}
	for ( j = 0 ; j <= errVec.size()-1 ; j += 1 ) {
		fv_file.printf( "%g " , errVec.x[ j ] )
	}
	fv_file.printf("%10.2f\n",fitnessVal)
	fv_file.flush()
	
	return fitnessVal
}

func tfunkpar(){local i, k, fitnessVal, ind, key	localobj locerrVec

	key = hoc_ac_
    ind = $1
    transvec = $o2
    
    if (DEBUG) printf("%d: ind=%d\n",pc.id,ind)

	fitnessVal=0

	//set parameters
	sprint( parstr , "set_n_params(" )
	for ( k = 0 ; k < NP ; k += 1 ) {
		if ( k == 0 ) {
			sprint( parstr , "%stransvec.x[%d]" , parstr , k )
		} else {
			sprint( parstr , "%s,transvec.x[%d]" , parstr , k )
		}
	}
	sprint( parstr , "%s)" , parstr )
	execute( parstr  )

	//adjust conductances
	set_conds()
	//adjust kinetics
	set_kins()
	
	if (DEBUG) printf("made it...1\n")
	MRF.p.run()
	if (DEBUG) printf("made it...2\n")
	fitnessVal=MRF.p.pf.errval
	locerrVec = get_error_values()
	
	pc.pack(ind)
    pc.pack(fitnessVal)
    pc.pack(locerrVec)
    pc.pack(transvec)
    pc.post(key)

    return key
}

// output function
proc print_things() { local j

	// Add best result of this generation to the list

	//printf("[Generation # %g]  ",num_gen)
	g_file.printf( "%g\t  " , num_gen )
	for ( j = 0 ; j <= NP - 1 ; j+=1 ) {
		//printf("%10.5f \t",best_indiv.x[j]) 
		g_file.printf( "%1.16f \t" , pop.x[ Best ][ j ] )

        //prax_par[j]=best_indiv.x[j]
	}
	//printf("%15.7e\n",Best)
	//g_file.printf("%10.2f\n",res.x[Best])
	g_file.printf( "%10.2f\n" , ffs.getrow( Best ).sum )
	g_file.flush()
	//pop.printf("%5.6f ")
	
	// Overwrite the whole of the current population for resuming
	// or for making vtraces.dat
	p_file.wopen( "output/curr_population" )
	pop.fprint( p_file , "%g\t" , "\n" )
	res.printf( p_file , "%g\n" )
	p_file.close()
}

proc find_best(){local i,j
	Best = 0
	for ( i = 1 ; i < pop_size ; i+=1 ) {
		if ( res.x[ i ] < res.x[ Best ] ) {
			Best = i
		}
	}
}

/************************************
 * Output processes and functions
 ***********************************/

proc save_voltages() { local i, j, k , ind	localobj paramcombs

	// First entry in vtrace file is Best
	find_best()
	printf("Saving vTrace %d\n",i)
	starttimevolt=startsw()	
	for( j = 0 ; j <= ( NP - 1 ) ; j += 1 ) transvec.x[j] = pop.x[Best][j]
	sprint( xstr , "outputVtrace(" )
	for ( k = 0 ; k < NP ; k += 1 ) {
		if ( k == 0 ) {
			sprint( xstr , "%stransvec.x[%d]" , xstr , k )
		} else {
			sprint( xstr , "%s,transvec.x[%d]" , xstr , k )
		}
	}
	sprint( xstr , "%s)" , xstr )
	printf(xstr)
	execute(xstr)
	printf("Save took %f\n",startsw()-starttimevolt)

	// Then open up parameter combination file and set paramcombs Vector
	// to contain the parameter sets selected for saving vtraces of
	sprint( xstr, "%soutput/paramcombs.dat", PARENTDIR )
	pc_file.ropen( xstr )
	paramcombs = new Vector()
	//paramcombs.append(1,29,58,86,114,143,171,199,228,256)
	paramcombs.scanf(pc_file)
	pc_file.close()
	
	for ( i = 0 ; i < ( paramcombs.size ) ; i += 1 ) {
		ind = paramcombs.x[i]
		printf("Saving vTrace %d\n",i)
		starttimevolt=startsw()
		for( j = 0 ; j <= ( NP - 1 ) ; j += 1 ) transvec.x[j] = pop.x[ind-1][j]
		sprint( xstr , "outputVtrace(" )
		for ( k = 0 ; k < NP ; k += 1 ) {
			if ( k == 0 ) {
				sprint( xstr , "%stransvec.x[%d]" , xstr , k )
			} else {
				sprint( xstr , "%s,transvec.x[%d]" , xstr , k )
			}
		}
		sprint( xstr , "%s)" , xstr )
		execute(xstr)
		printf("Save took %f\n",startsw()-starttimevolt)
	}
	vt_file.close()
}

proc save_iv() { local i , j , k , l , stimlvl , numtowrite , ind	localobj targetV , modelV , paramcombs

	// load up the target voltage vector
	sprint( xstr , "data/%s/iv.dat" , cellname )
	iv_file.ropen( xstr )
	targetV = new Vector()
	targetV.scanf( iv_file )
	
	// open up the output dat file
	ivout_file.wopen( "output/iv.dat" )
	
	// Use the best parameter set for testing IV relationship
	find_best()
	printf( "Saving IV curve data \n" )

	// save the best one (index 0... hopefully...) and then save every one in 
	// paramcombs.dat, then write the number of spikes to the ilfe...
	
	// Open parameter combination file and set paramcombs vector to contain
	// the parameter sets selected for saving data
	sprint( xstr, "%soutput/paramcombs.dat", PARENTDIR )
	pc_file.ropen( xstr )
	paramcombs = new Vector()
	paramcombs.scanf(pc_file)
	pc_file.close()
	
	// modelV matrix will be width of paramcombs.size + 1 and length 9...
	modelV = new Matrix( 9 , ( paramcombs.size + 1 ) )
	
	for ( i = -1 ; i < ( paramcombs.size ) ; i += 1 ) {
		if( i==-1 ) {
			ind = 1
		} else {
			ind = paramcombs.x[i]
		}
		// we save the best fit first, and assumption here is a sorted 
		// curr_population file with best fit as the first entry
		printf("Saving IV %d\n",ind)
		for( j = 0 ; j <= ( NP - 1 ) ; j += 1 ) transvec.x[j] = pop.x[ind][j]
		sprint( xstr , "set_n_params(" )
		for ( k = 0 ; k < NP ; k += 1 ) {
			if ( k == 0 ) {
				sprint( xstr, "%stransvec.x[%d]", xstr, k )
			} else {
				sprint( xstr, "%s,transvec.x[%d]", xstr, k )
			}
		}
		sprint( xstr, "%s)", xstr )
		execute( xstr  )
		//adjust conductances
		set_conds()
		//adjust kinetics
		set_kins()
		
		// 9 levels of current injection tested:
		// -160, -120, -80, -40, 0, 40, 80, 120, 160
		for( l = 0 ; l < 9 ; l += 1 ) {
			stimlvl = ( ( l * 0.04 ) - 0.16 )
			sprint( xstr , "%s=%g" , stimname , stimlvl )
			execute( xstr  )
			
			// init and run, take the voltage at 150ms...
			init()
			tstop=150
			run()
			
			// take the voltage
			modelV.x[l][i+1] = v( 0.5 )
		}
	}

	// write to file: stimlvl , target , each model
	// write all 9 if using cell Dec15, but only 6 if using Jun24
	if ( CELL == 2 ) {
		numtowrite = 9
	} else if ( CELL == 3 ) {
		numtowrite = 6
	} else {
		numtowrite = 6
	}
	
	for( l = 0 ; l < numtowrite ; l += 1 ) {
		stimlvl = ( ( l * 0.04 ) - 0.16 ) * 1000
		sprint( xstr , "%g\t%g\t" , stimlvl , targetV.x[ l ] )
		ivout_file.printf( xstr )
		for( i = 0 ; i < (paramcombs.size+1) ; i += 1 ) {
			sprint( xstr , "%g\t" , modelV.x[l][i] )
			ivout_file.printf( xstr )
		}
		sprint( xstr , "\n" )
		ivout_file.printf( xstr )
	}
	
	ivout_file.close()
}

proc save_fi() { local i , j , k , l , stimlvl  ,ind	localobj targetFR , voltVecSpikes , spikeHist , modelFR , paramcombs
	// load vectors
	voltVecSpikes = new Vector()
	spikeHist = new Vector()
	
	// load up the target voltage vector
	sprint( xstr , "data/%s/fi.dat" , cellname )
	fi_file.ropen( xstr )
	targetFR = new Vector()
	targetFR.scanf( fi_file )
	
	// open up the output dat file
	fiout_file.wopen( "output/fi.dat" )
	
	printf( "Saving FI curve data \n" )
	
	// save the best one (index 0... hopefully...) and then save every one in 
	// paramcombs.dat, then write the number of spikes to the ilfe...
	
	// Open parameter combination file and set paramcombs vector to contain
	// the parameter sets selected for saving data
	sprint( xstr, "%soutput/paramcombs.dat", PARENTDIR )
	pc_file.ropen( xstr )
	paramcombs = new Vector()
	paramcombs.scanf(pc_file)
	pc_file.close()
	
	// modelFR matrix will be width of paramcombs.size + 1 and length 8...
	modelFR = new Matrix( 8 , ( paramcombs.size + 1 ) )
	
	for ( i = -1 ; i < ( paramcombs.size ) ; i += 1 ) {
		if( i==-1 ) {
			ind = 0
		} else {
			ind = paramcombs.x[i]
		}
		// we save the best fit first, and assumption here is a sorted 
		// curr_population file with best fit as the first entry
		printf("Saving FI %d\n",ind)
		for( j = 0 ; j <= ( NP - 1 ) ; j += 1 ) transvec.x[j] = pop.x[ind-1][j]
		sprint( xstr , "set_n_params(" )
		for ( k = 0 ; k < NP ; k += 1 ) {
			if ( k == 0 ) {
				sprint( xstr, "%stransvec.x[%d]", xstr, k )
			} else {
				sprint( xstr, "%s,transvec.x[%d]", xstr, k )
			}
		}
		sprint( xstr, "%s)", xstr )
		printf( "i = %d ; ind = %d ; %f\n" , i , ind , transvec.x[1] )
		execute( xstr  )
		//adjust conductances
		set_conds()
		//adjust kinetics
		set_kins()
		
		// 8 levels of current injection tested:
		// 30, 80, 130, 180, 230, 280, 330, 380
		for( l = 0 ; l < 8 ; l += 1 ) {
			stimlvl = ( ( l * 0.05 ) + 0.03 )
			sprint( xstr , "%s=%g" , stimname , stimlvl )
			execute( xstr  )
			
			// init and run, take the voltage at 200ms...
			init()
			run()
			
			// count the spikes
			voltVecSpikes.spikebin( voltVec , -20 )
			spikeHist = voltVecSpikes.histogram( 0 , 1 , 1 )
			modelFR.x[l][i+1] = spikeHist.x[ 2 ]
		}
	}

	// write to file: stimlvl , target , each model
	for( l = 0 ; l < 8 ; l += 1 ) {
		stimlvl = ( ( l * 0.05 ) + 0.03 )
		sprint( xstr , "%g\t%g\t" , stimlvl , targetFR.x[ l ] )
		fiout_file.printf( xstr )
		for( i = 0 ; i < (paramcombs.size+1) ; i += 1 ) {
			sprint( xstr , "%g\t" , modelFR.x[l][i] )
			fiout_file.printf( xstr )
		}
		sprint( xstr , "\n" )
		fiout_file.printf( xstr )
	}
	
	fiout_file.close()
}

/* process sensitivity_analysis()
 * Generates fitness values for all fitness functions in all currently loaded 
 * generators for a range of parameters centred around the curr_population.
 * This means the process takes the current best fitting parameter values for 
 * the model, tests them with +- a % to each parameter in turn, and saves the 
 * associated parameter values and fitness value from every function in an 
 * output file.
 * 
 * Runs in parallel using the same code the Evaluate_Pop() uses to distribute 
 * jobs.
 * 
 * Currently unsure how exactly to generate the populations. Can use:
 * 		N_processes == number of permutations of parameter values, then get 
 * 		all of the permutations for one parameter in the population at a time, 
 * 	OR
 * 		N_processes == population size, and do one permutation of each 
 * 		parameter at a time, generating that permutation for every member of 
 * 		the population.
 * 
 */

proc sensitivity_analysis(){local i,j,k,numPerms,fitnessVal,ind	localobj permMat
	/* Procesing order:
	 * 1) Generate populations to be submitted to workers one at a time, 
	 * 	  sending them off to run, waiting for the results, gathering results 
	 * 	  as they come back, and writing to the output file.
	 */
	 
	 printf( "Sensitivity Analysis commencing in 3...2...1... GO!\n" )
	 
	//	1) going to use first methoed where n_processes = number of permutations
	// initialise vector holding permutations, and work out number of 
	// permutations based on number of parameters (NP)
	
	numPerms = ( ( NP * 4 ) + 1 )
	permMat = new Matrix( numPerms , NP )
	
	// set up permutations for members of population one at a time
	// Can set pop_size here to not run through lots of very similar parameters
	pop_size = 1
	
	for ( i = 0 ; i <= pop_size-1 ; i += 1 ) {
		printf( "Population member: %d\n" , i )
		printf( "Generating permutations...\n" )
		// fill the permMat with identical vectors, the current pop member:
		for ( j = 0 ; j <= numPerms-1 ; j += 1 ) {
			for ( k = 0 ; k <= NP-1 ; k += 1 ) {
				permMat.x[ j ][ k ] = pop.x[ i ][ k ]
			}
		}
		// now modulate one parameter by -10, -5, +5 and +10 % in 4 adjacent 
		// rows, for each paramter, starting at row index 1 (base is 0)...
		row = 1
		for ( j = 0 ; j <= NP-1 ; j += 1 ) {
			permMat.x[ row ][ j ] = ( pop.x[ i ][ j ] * 0.9 )
			row += 1
			permMat.x[ row ][ j ] = ( pop.x[ i ][ j ] * 0.95 )
			row += 1
			permMat.x[ row ][ j ] = ( pop.x[ i ][ j ] * 1.05 )
			row += 1
			permMat.x[ row ][ j ] = ( pop.x[ i ][ j ] * 1.1 )
			row += 1
		}
		permMat.printf()
		// now run simulations to get fitness values for each permutation, and 
		// save output to the sa_file()
		if ( pc.nhost == 1 ) {    // use the serial form
			// can't just use tfunk(), as this is always set to save to fv_file, 
			// so the guts of it are copied here:
			//for each permutation...
			for ( j = 0 ; j <= numPerms-1 ; j += 1 ) {
				printf( "...Simulating permutation %d...\n" , j )
				//...retrieve parameters from permMat...
				for ( k = 0 ; k <= NP-1 ; k += 1 ) {
					transvec.x[ k ] = permMat.x[ j ][ k ]
				}
				transvec.printf()
				//...set parameters...
				sprint( xstr , "set_n_params(" )
				for ( k = 0 ; k < NP ; k += 1 ) {
					if ( k == 0 ) {
						sprint( xstr , "%stransvec.x[%d]" , xstr , k )
					} else {
						sprint( xstr , "%s,transvec.x[%d]" , xstr , k )
					}
				}
				sprint( xstr , "%s)" , xstr )
				execute(xstr)
				//...adjust conductances...
				set_conds()
				//...adjust kinetics...
				set_kins()
				//...run simulation...
				MRF.p.run()
				//...retrieve total error...
				fitnessVal=MRF.p.pf.errval
				//...retrieve individual error values...
				errVec = get_error_values()
				//...write this stuff to the output file...
				sa_file.printf("%d %d ",1,1)
				for ( k = 0 ; k <= NP-1 ; k += 1 ) {
					sa_file.printf( "%1.16f " , transvec.x[ k ] )
				}
				for ( k = 0 ; k <= errVec.size()-1 ; k += 1 ) {
					sa_file.printf( "%g " , errVec.x[ k ] )
				}
				sa_file.printf("%10.2f\n",fitnessVal)
				sa_file.flush()
				//.. finally loop to next permutation
			}
		} else {	// use the bulletin board form
			// for each permutation...
			for ( j = 0 ; j <= numPerms-1 ; j += 1 ) {
				//...retrieve parameters from permMat...
				for ( k = 0 ; k <= NP-1 ; k += 1 ) {
					transvec.x[ k ] = permMat.x[ j ][ k ]
				}
				if ( DEBUG ) printf( "j=%d\n" , j )
				//...submit a job for this permutation to the worker nodes...
				pc.submit( "tfunkpar" , j , transvec )
				//...loop back and set up a job for the next permutation...
			}
			// While	 results aren't in from every job yet...
			while ((id = pc.working()) != 0) {
				//...when something gets returned, retrieve the output values...
				pc.look_take(id)
				if (DEBUG) printf("made it into while\n")
				ind = pc.upkscalar()
				if (DEBUG) printf("ind=%d\n",ind)
				fitVal = pc.upkscalar()
				if (DEBUG) printf("made it past upkscalars\n")
				errVec = pc.upkvec()
				if ( DEBUG ) printf( "collected %d error values\n" , errVec.size() )
				//...write the results to the output file...
				//fv_file.printf("%d %d %1.16f %1.16f %1.16f %10.2f\n",id,ind,pop.x[ind][0],pop.x[ind][1],pop.x[ind][2],fitVal)
				sa_file.printf("%d %d ",id,ind)
				for ( j = 0 ; j <= NP-1 ; j += 1 ) {
					sa_file.printf( "%1.16f " , pop.x[ ind ][ j ] )
				}
				for ( j = 0 ; j <= errVec.size()-1 ; j += 1 ) {
					sa_file.printf( "%g " , errVec.x[ j ] )
				}
				sa_file.printf("%10.2f\n",fitVal)
				sa_file.flush()
			}
		}
	}
	sa_file.close()
	if ( DEBUG ) printf( "made it out of sensitivity analysis...\n" )
	printf( "Sensitivity analysed!\n" )
}

proc random_testing(){local i , j , k , fitVal	localobj paramVec
	
	// make file of randomTesting parameters for producing output figures
	param_file.wopen( "output/metaparams.dat" )
	param_file.printf( "%s: %d\n" , "population" , rand_max_pop ) 
	param_file.printf( "%s: %d\n" , "parameters" , NP )
	param_file.printf( "%s: %d\n" , "generators" , NGens )
	param_file.printf( "%s: %d\n" , "totalFFs" , totalFFs )
	param_file.printf( "%s: %d\n" , "totalT" , (Tend.max()/voltvecdt) )
	param_file.printf( "%s: %d\n" , "MRF" , MRFflag )
	param_file.printf( "%s: %g\n" , "dt" , voltvecdt)
	param_file.printf( "%s: %d\n" , "generations" , 1 )
	param_file.close()

	begintime=startsw()
	
	// submit a load of jobs for the parallel workers to attack...
	for ( i = 0 ; i < rand_max_pop ; i += 1 ) {
		for ( j = 0 ; j < NP ; j += 1 ) {
			transvec.x[ j ] = r[ j ].repick()
		}
		pc.submit( "tfunkpar" , i , transvec )
	}
	//... wait for results back from the workers, and print them to a file
	while ((id = pc.working()) != 0) {
		pc.look_take(id)
		ind = pc.upkscalar()
		fitVal = pc.upkscalar()
		errVec = pc.upkvec()
		paramVec = pc.upkvec()
		fv_file.printf("%d %d ",id,ind)
		for ( i = 0 ; i < paramVec.size() ; i += 1 ) {
			fv_file.printf( "%1.16f " , paramVec.x[ i ] )
		}
		for ( i = 0 ; i < errVec.size() ; i += 1 ) {
			fv_file.printf( "%g " , errVec.x[ i ] )
		}
		fv_file.printf("%10.2f\n",fitVal)
		fv_file.flush()
	}
	Elapsed = startsw() - begintime
 	print "-------------------------------end of Random Testing----------------------------------"
	printf( "\nTime Taken = %d\n" , Elapsed )
	g_file.close()
	fv_file.close()
}