detailed_simulation_code_multi_critic.jl

########## Main simulation functions #############

# putting noise updates in a function (which must be called!)
#  rather than in the post() function, for debugging reasons
function update_noise()
  #global ksi = rand(Normal(0,output_noise), no_post_neurons);
  global ksi = rand(Normal(0,1), no_post_neurons) .* sqrt(output_noise_variance);
  global population_ksi = rand(Normal(0,1), no_post_neurons) .* sqrt(output_noise_variance);
end


function initialise_pre_population(tuning_type::linear_tc)
  # a and b only get set once for this paper
  #   they are receptive fields rather than 'noise'
  #global a = rand(Normal(input_baseline, input_baseline_variance), (no_pre_neurons_per_task, no_input_tasks));
  global a = rand(Normal(0,1), no_pre_neurons_per_task, no_input_tasks) .* sqrt(input_baseline_variance) + input_baseline;
  global b = zeros(no_pre_neurons_per_task, no_input_tasks);

  for i = 1:no_input_tasks
    #b[:,i] = rand(Normal(0, task_tuning_slope_variance[i]), no_pre_neurons_per_task);
    b[:,i] = rand(Normal(0,1), no_pre_neurons_per_task) .* sqrt(task_tuning_slope_variance[i]);
  end
end

function initialise_pre_population(tuning_type::gaussian_tc)
  # Multiple dispatch should use this function rather than the linear
  #   receptive fields function if a boolean is passed as a parameter
  global a;
  # Note: there is no b generated by this function!

  #CONSIDER: should probably extend mu beyond (-1,1) for complete coverage of interval
  #CONSIDER: should I reduce height by no_tuning_curves_per_input_neuron?
  a = Array(gaussian_tc_type, (no_pre_neurons_per_task, no_input_tasks) );
  for i = 1:no_input_tasks;
    #figure()
    for j=1:no_pre_neurons_per_task;
      tuning_mu = rand(Uniform(gaussian_tuning_mu_lower_bound,gaussian_tuning_mu_upper_bound), no_tuning_curves_per_input_neuron);
      if (fix_tuning_gaussian_width)
        tuning_sigma = ones(no_tuning_curves_per_input_neuron);
        tuning_sigma *= gaussian_tuning_sigma_mean;
      else
        #tuning_sigma = rand(Normal(gaussian_tuning_sigma, gaussian_tuning_sigma_width), no_tuning_curves_per_input_neuron);
        tuning_sigma = rand(Normal(0, 1), no_tuning_curves_per_input_neuron) .* sqrt(gaussian_tuning_sigma_variance) + gaussian_tuning_sigma_mean;
      end
      tuning_height = rand(Normal(0,1), no_tuning_curves_per_input_neuron) .* sqrt(gaussian_tuning_height_variance) + gaussian_tuning_height_mean;
      if ( normalise_height_of_multiple_gaussian_inputs )
        tuning_height = tuning_height ./ no_tuning_curves_per_input_neuron;
      end
      a[j,i] = gaussian_tc_type(no_tuning_curves_per_input_neuron, tuning_mu, tuning_sigma, tuning_height);

      #scatter(tuning_mu, tuning_height, c="r");
      #scatter(tuning_mu, tuning_sigma, c="b");
    end
  end
end


function plot_gaussian_tuning_single_input(neuron_id::Int=1, task_id::Int=1)
  #figure();
  x = linspace(-1,1,101);
  y = zeros(101);
  for i = 1:101
    #local_pre = pre(x[i], task_id, gaussian_tc() );
    #y[i] = sum(local_pre);
    for j = 1:a[neuron_id, task_id].no_curves
      f(x) = a[neuron_id, task_id].height[j] .* exp( -(x - a[neuron_id,task_id].mu[j]).^2 ./ (2 * ( a[neuron_id, task_id].sigma[j] .^2 ) ) );
      y[i] += f(x[i]);
    end
  end

  plot(x,y);
end

function plot_gaussian_tuning_multi_inputs(task_id::Int=1, begin_id::Int=1, end_id::Int=no_pre_neurons_per_task)
  figure()
  for i = begin_id:end_id
    plot_gaussian_tuning_single_input(i, task_id);
  end
  xlim([-1,1])
  ylim([0,3])
  title("Input layer tuning curves, direction $task_id")
  xlabel("Input range")
  ylabel("Neuronal firing rate")
end


# pre-synaptic firing rate upon presentation of pattern x
function pre(x::Float64, task_id::Int, tuning_type::linear_tc)
  local_pre = zeros(no_pre_neurons_per_task, no_input_tasks);
  local_pre[:,task_id] = collect(a[:,task_id] + b[:,task_id] .* x); # weird vcat error in Julia v0.4
  return local_pre;
end

# pre-synaptic firing rate upon presentation of pattern x
#   using gaussian based tuning curves
#   current version still maintains task specific populations
function pre(x::Float64, task_id::Int, tuning_type::gaussian_tc)
  local_pre = zeros(no_pre_neurons_per_task, no_input_tasks);

  for neuron_id = 1:no_pre_neurons_per_task
    for j = 1:a[neuron_id, task_id].no_curves
      f(x) = a[neuron_id, task_id].height[j] .* exp( -(x - a[neuron_id,task_id].mu[j]).^2 ./ (2 * ( a[neuron_id, task_id].sigma[j] .^2 ) ) );
      local_pre[neuron_id, task_id] += f(x);
    end
  end

  return local_pre;
end


function initialise_weight_matrix(tuning_type::gaussian_tc)
  # Remember: always call this after a and b have already been initialised!
  #set initial weights
  global w = rand(Uniform(0,1), (no_pre_neurons_per_task, no_post_neurons, no_input_tasks));
  #w = ones(no_pre_neurons_per_task, no_post_neurons, no_input_tasks);

  centre_of_mass = zeros(no_pre_neurons_per_task, no_input_tasks);
  #CONSIDER: adding height to centre of mass calculation (sum(mu.*h)/no_curves)
  for i = 1:no_input_tasks
    for j = 1:no_pre_neurons_per_task
      centre_of_mass[j, i] = sum(a[j,i].mu) ./ a[j,i].no_curves ;
    end
    w[:,1,i] += -gaussian_weight_bias.*centre_of_mass[:,i];
    w[:,2,i] += gaussian_weight_bias.*centre_of_mass[:,i];
  end

  # hard bound weights at +/- 10
  w[w .> weights_upper_bound] = weights_upper_bound;
  w[w .< weights_lower_bound] = weights_lower_bound;
end

function initialise_weight_matrix(tuning_type::linear_tc)
  # Remember: always call this after a and b have already been initialised!
  #set initial weights
  global w; #Array{Float64,3,w};# :: Array{Float64,3};

  if (!use_defined_performance_setup)
    w = rand(Uniform(0,1), no_pre_neurons_per_task, no_post_neurons, no_input_tasks);
    for i = 1:no_input_tasks
      w[:,1,i] += -initial_weight_bias.*b[:,i];
      w[:,2,i] += initial_weight_bias.*b[:,i];
    end
  else
    w = ones(no_pre_neurons_per_task, no_post_neurons, no_input_tasks);
    w *= 0.5;
    D_rev(p) = 2.0 * sqrt(output_noise_variance) * erfinv(2 * p - 1.0);
    positive_difference_in_outputs_1 = D_rev(defined_performance_task_1);
    positive_difference_in_outputs_2 = D_rev(defined_performance_task_2);
    local_bias = zeros(no_pre_neurons_per_task, no_post_neurons, no_input_tasks);
    for i = 1:no_input_tasks
      # currently assuming only two input tasks (+/-1)
      local_pre_1 = pre(-1.0, i, linear_tc());
      local_pre_2 = pre(1.0, i, linear_tc());

      #TODO: dimension mismatch, need to realign these matrices
      local_bias[:,1,i] = positive_difference_in_outputs_1 ./ ( 2 * local_pre_1[:,i] .* b[:,i] );
      local_bias[:,2,i] = positive_difference_in_outputs_1 ./ ( 2 * local_pre_2[:,i] .* b[:,i] ) ;

      w[:,1,i] += local_bias[:,1,i] .* b[:,i];
      w[:,2,i] += local_bias[:,2,i] .* b[:,i];
    end
  end

  # calculate square root of average squared weight value, want to keep this constant during weight normalisation process
  #   Henning and I left out sqrt() when we were writing this together
  global subject_initial_weight_scale = sqrt( mean(w[:,:,:].^2) ) :: Float64; #sqrt(sum(w.^2))
  # as a test, we look at normalisation separately per task
  #   these normalisation factors must still combine the two outputs as the real learning
  #   is to discriminate between left and right
  global subject_task_1_initial_weight_scale = sqrt( mean(w[:,:,1].^2) ) :: Float64;
  global subject_task_2_initial_weight_scale = sqrt( mean(w[:,:,2].^2) ) :: Float64;
end


function initialise()
  srand(random_seed);

  if(use_gaussian_tuning_function)
    # use gaussian basis functions
    tuning_type = gaussian_tc();
  elseif(use_linear_tuning_function)
    # use linear tuning functions
    tuning_type = linear_tc();
  else
    print("ERROR: you need to define a tuning function\n");
    error(1);
  end

  initialise_pre_population(tuning_type);
  update_noise();
  initialise_weight_matrix(tuning_type);

  if !use_hard_coded_critic
    initialise_critic_parameters()
  end

  global decision_criterion_monitor = zeros(no_decision_monitors,1);

  global average_delta_reward = 0.0;
  global average_choice = 0.0;
  #global n = 0 :: Int; # use this to monitor trial ID per block (very important: this is a block level counter!)
  global n_within_block = 0 :: Int; # use this to monitor trial ID per block (very important: this is a block level counter!)
  global n_task_within_block = zeros(Int, no_input_tasks) :: Array{Int,1};
  # changing to multi-critic model
  #   critic can be per block or over entire learning history
  global n_critic = zeros(Int, no_task_critics, no_choices_per_task_critics); # use this to monitor trial ID per critic
  global average_reward = zeros(no_task_critics, no_choices_per_task_critics); # running average, stored values represent end of a block value
  global average_block_reward = 0.0;
  global instance_correct = 0;
  global instance_incorrect = 0;

  # intrinsic plasticity requires running averages of post-synaptic firing rates
  global average_post = zeros(no_post_neurons);
  global n_post = 0 :: Int;

  global proportion_1_correct = 0.0;
  global proportion_2_correct = 0.0;

  global exp_results = Array{RovingExperiment}(0);

  global enable_weight_updates = true;
  typeassert(enable_weight_updates, Bool);

  global use_fixed_external_bias = false;
  typeassert(use_fixed_external_bias, Bool);

  print("RND seeded $random_seed\n")
end


__init__ = initialise();


# winner takes all
function wta(left::Float64, right::Float64, debug_on::Bool = false)
  # debugging code to highlight negative post firing rates
  if(debug_on)
    if(verbosity > 0)
      if (left < 0)
        print("Flag left -ve\n")
      end
      if (right < 0)
        print("Flag right -ve\n")
      end
    end
  end

	if (left > right)
    if(verbosity > 0)
      if(debug_on)
        print("Left!\n")
      end
    end
    right = 0.
	else
    if(verbosity > 0)
      if(debug_on)
        print("Right!\n")
      end
    end
    left = 0.
	end
	return [left right]
end

# population pooling of post-synaptic decision neurons leads to a feedback of WTA from decision layer
#   to weight learning neurons
function wta(local_post::Array{Float64,2}, pop_rate::Array{Float64,2})
  local_pop_rate = wta(pop_rate[1], pop_rate[2]);
  if (local_pop_rate[1] > local_pop_rate[2])
    return [local_post[1] 0.0];
  else
    return [0.0 local_post[2]];
  end
end


function noise_free_post(x::Float64, task_id::Int, tuning_type::TuningSelector)
  local_pre = pre(x, task_id, tuning_type)

  noise_free_left = sum(local_pre[:,task_id] .* w[:,1,task_id]);
  noise_free_right = sum(local_pre[:,task_id] .* w[:,2,task_id]);

  # intrinsic plasticity subtracts a running average of post
  #   we're probably ignoring this for now, but I'll leave the code hear and add a note
  #   to be careful enabling the parameter
  if( use_intrinsic_plasticity )
    # we use intrinsic_baseline to offset post synaptic firing from zero
    left = left - average_post[1] + intrinsic_baseline[1];
    right = right - average_post[2] + intrinsic_baseline[2];
  end

  # alternative way of biasing decisions 50:50 Left:Right
  if( use_decision_criterion_learner )
    local_monitor_task_id = min(no_decision_monitors, task_id);
    noise_free_left = noise_free_left + decision_criterion_monitor[local_monitor_task_id];
    noise_free_right = noise_free_right - decision_criterion_monitor[local_monitor_task_id];
  end

  return (noise_free_left, noise_free_right);
end


# update average_post variable
#   note that intrinsic_baseline is not included in the running average
function update_intrinsic_excitability(x::Float64, task_id::Int, tuning_type::TuningSelector)
  global average_post;
  typeassert(average_post, Array{Float64,1});
  global n_post;
  typeassert(n_post, Int);

  (noise_free_left, noise_free_right) = noise_free_post(x, task_id, tuning_type);

  if(use_intrinsic_plasticity_with_noise)
    left = noise_free_left + ksi[1]
  	right = noise_free_right+ ksi[2]
  end

  # we use intrinsic_baseline to offset post synaptic firing from zero, so don't add it here
  #   or we will get zero firing rate on average
  left = left - average_post[1];
  right = right - average_post[2];

  # hack: putting a lower bound on post synaptic firing
  if (left < floor_on_post)
    left = floor_on_post;
  end
  if (right < floor_on_post)
    right = floor_on_post;
  end

  # we need to decide whether we really want the average_post to use wta() form of output or not
  if( (!disable_winner_takes_all) && (use_intrinsic_plasticity_with_wta_form) )
    (left,right) = wta(left,right);
  end

  # update average_post
  n_post += 1;
  tau = min(intrinsic_plasticity_window_length, n_post);

  #print("$average_post \n")

  average_post[1] = ( (tau - 1) * average_post[1] + left ) / tau;
  average_post[2] = ( (tau - 1) * average_post[2] + right ) / tau;
end


# criterion learner is capable of accounting for biases in presentation frequency between left and right
function update_decision_criterion_learner(local_post, task_id::Int)
  global decision_criterion_monitor;

  # following the logic that a critic which cannot distinguish the two tasks
  #     will only have a single criterion learner... (bad idea?)
  local_monitor_task_id = min(task_id, no_decision_monitors);

  # associate 1 with output 2 and 0 with output 1
  update_choice = (local_post[2] > local_post[1] ? 1. : 0.) :: Float64;
  # I previously forgot to include the following term, which maintains the 50:50 ratio
  update_choice = update_choice - criterion_learner_expectation;

  decision_criterion_monitor[local_monitor_task_id] = decision_criterion_monitor[local_monitor_task_id] + (decision_criterion_timescale * update_choice);

  # the multiplicative decision criterion monitor was clamped between 0 and 1, this is currently not of use
#=  if (decision_criterion_monitor > 1)
    decision_criterion_monitor = 1;
  elseif (decision_criterion_monitor) < 0
    decision_criterion_monitor = 0;
  end=#
end


# post-synaptic firing rate upon presentation of pattern x
#  no longer generating a new noise value (ksi) on each call,
#  this must be done externally to allow for repeatibility during debug
# Note: local_post returns a tuple where one value is 0. All comparisons to find the non zero value should use absolute comparison.
function post(x::Float64, task_id::Int, tuning_type::TuningSelector, debug_on::Bool=false)
	#local_pre = pre(x, task_id, tuning_type)
  #noise_free_left = sum(local_pre[:,task_id] .* w[:,1,task_id]);
  #noise_free_right = sum(local_pre[:,task_id] .* w[:,2,task_id]);
  (noise_free_left, noise_free_right) = noise_free_post(x, task_id, tuning_type);

  # noise is scaled via the size of an imaginary pool of post neurons, per output decision category
  #   this is single neuron code, so I think the noise should only be scaled by output_noise_variance (in ksi)
  #   and not by the number of scaling neurons in the post population
  left = noise_free_left + ksi[1] * sqrt(no_pop_scaling_post_neurons)
	right = noise_free_right+ ksi[2] * sqrt(no_pop_scaling_post_neurons)

  # moved intrinsic plasticity application to noise_free_post()
  # moved decision criterion monitor application to noise_free_post()

  # calculated probability of getting this result given de-noised results and error size
  #   TODO: finish this code
  trial_probability_left = 0.5 + erf( (noise_free_left - noise_free_right) / (sqrt(output_noise_variance) * 2.0) ) * 0.5;

  # hack: putting a lower bound on post synaptic firing
  if (left < floor_on_post)
    left = floor_on_post;
  end
  if (right < floor_on_post)
    right = floor_on_post;
  end

  if(debug_on)
    if(verbosity > 0)
      print("n_within_block: $n_within_block, x: $x, left: $left, right: $right,\n noise_free_left: $noise_free_left, noise_free_right: $noise_free_right, trial_probability_left: $trial_probability_left ")
    end
  end
  if(disable_winner_takes_all)
    return [left right];
  else
    if ( no_pop_scaling_post_neurons > 1)
      # population pooling of post-synaptic decision neurons leads to a feedback of WTA from decision layer
      #   to weight learning neurons
      # there is a clear performance punishment by re-calculating the pooled_population_post_rate here
      #     (post() gets called multiple times!)
      #   but it seems cleaner than most alternative implementations without re-writing substantial
      #   pieces of code. It also maintains invariance of most functions and a logical implementation
      #   of the WTA function.
      local_post = [left right];
      pop_rate = pooled_population_post_rate(x, task_id, tuning_type, local_post);
      return wta(local_post, pop_rate);
    else
      return wta(left,right, debug_on);
    end
  end
end


# calculate the noise-free difference in the outputs for a given input
#   return value is positively associated with correct classification
function noise_free_output_positive_difference(x::Float64, task_id::Int, tuning_type::TuningSelector)
  #local_pre = pre(x, task_id, tuning_type)
  #noise_free_left = sum(local_pre[:,task_id] .* w[:,1,task_id]);
  #noise_free_right = sum(local_pre[:,task_id] .* w[:,2,task_id]);

  (noise_free_left, noise_free_right) = noise_free_post(x, task_id, tuning_type);

  output_sign_selector = 1;
  if (x > 0)
    # then output right should be larger
    output_sign_selector = -1;
  end

  return ( output_sign_selector * (noise_free_left - noise_free_right) );
end


# probability output x is chosen given that input x is given
function probability_correct(x::Float64, task_id::Int, tuning_type::TuningSelector)
  local_noise_free_output_positive_difference = noise_free_output_positive_difference(x, task_id, tuning_type);

  return (0.5 + 0.5 * erf( ( local_noise_free_output_positive_difference ) / (sqrt(output_noise_variance) * 2) ) );
end


function post_hoc_calculate_thresholds(tuning_type::TuningSelector, subjects::Array{Subject,2}, split_output::Bool=false)
  # globals required for correct processing of pre()
  global a,b,w;
  global no_pre_neurons_per_task;
  global no_input_tasks;

  # The second dimension of subjects is basically a separate experimental
  #   protocol. Call this variable no_experimental_tasks_dimension
  #   to make it stand out
  (local_no_subjects, local_no_experimental_tasks_dimension) = size(subjects);
  no_points = 30;
  #x = linspace(0,1,no_points); # not sure if this is memory safe (scope rules)

  for j = 1:local_no_experimental_tasks_dimension
    #print("j: $j")
    for i = 1:local_no_subjects
      #print(", i: $i")
      a = deepcopy(subjects[i,j].a);
      if( isa(tuning_type, linear_tc) )
        b = deepcopy(subjects[i,j].b);
      end

      x = linspace(0,1,no_points); #check scope rules before switching to out of loop declaration of x
      # Calculate pre() for an entire linspace of inputs for this subject
      #   This is the heavy part of the processing which I wanted to reduce
      (no_pre_neurons_per_task, no_input_tasks) = size(a);
      local_pre_pos = zeros(no_pre_neurons_per_task, no_input_tasks, no_points);
      local_pre_neg = zeros(no_pre_neurons_per_task, no_input_tasks, no_points);
      #print(", generating pre(xi)...")
      for task_id = 1:no_input_tasks
        for m = 1:no_points
          # Assuming that pre returns zero entries for all non task specific entries
          local_pre_pos[:,task_id,m] = pre(x[m], task_id, tuning_type)[:,task_id];
          local_pre_neg[:,task_id,m] = pre(-x[m], task_id, tuning_type)[:,task_id];
        end
      end
      #print("done\n");

      # Calculate threshold for this subject using his trial by trial
      #   weight matrix and the associated task_id
      local_no_blocks_per_experiment = length(subjects[i,j].blocks);
      local_no_trials_per_block = length(subjects[i,j].blocks[1].trial);
      for k = 1:local_no_blocks_per_experiment
        #print("block: $k, ")
        # block level statistical variables
        local_average_threshold = 0.0;
        local_average_task_threshold = zeros(no_input_tasks);
        local_n_task_within_block = zeros(no_input_tasks);

        for l = 1:local_no_trials_per_block
          #print("trial: $l \n")
          # finally we get to processing a single threshold calculation
          task_id = subjects[i,j].blocks[k].trial[l].task_type;
          w = deepcopy(subjects[i,j].blocks[k].trial[l].w);

          x = linspace(0,1,no_points); # re-initialise due to sorting in loop below
          error_rate = zeros(no_points);
          split_error = zeros(no_points,2)

          # Loop over the xi in x, current implementation not really amenable
          #   to non-loop slicing
          for m = 1:no_points
            # calculate noise free post for xi
            local_noise_free_post_pos_left = sum(local_pre_pos[:,task_id,m].*w[:,1,task_id]);
            local_noise_free_post_pos_right = sum(local_pre_pos[:,task_id,m].*w[:,2,task_id]);

            # calculate noise free post for -xi
            local_noise_free_post_neg_left = sum(local_pre_neg[:,task_id,m].*w[:,1,task_id]);
            local_noise_free_post_neg_right = sum(local_pre_neg[:,task_id,m].*w[:,2,task_id]);
            if(verbosity > 2)
              print("DEBUG: $local_noise_free_post_pos_left, $local_noise_free_post_pos_right, ")
              print("$local_noise_free_post_neg_left, $local_noise_free_post_neg_right, ")
            end

            # probability, for a positive input (i) that we choose left
            p_pos_left = 0.5 + 0.5 * erf( (local_noise_free_post_pos_left - local_noise_free_post_pos_right) / (sqrt(output_noise_variance) * 2.0) );
            p_pos_right = (1. - p_pos_left);

            if(verbosity > 2)
              print("p_pos_left: $p_pos_left, p_pos_right: $p_pos_right ")
            end

            # probability, for a negative input (-i) that we choose left
            p_neg_left = 0.5 + 0.5 * erf( (local_noise_free_post_neg_left - local_noise_free_post_neg_right) / (sqrt(output_noise_variance) * 2.0) );
            p_neg_right = (1. - p_neg_left);
            if(verbosity > 2)
              print("p_neg_left: $p_neg_left, p_neg_right: $p_neg_right,")
            end

            # probability of choosing the wrong side
            error_rate[m] = ( p_pos_left + p_neg_right ) / 2.0;
            split_error[m,1] = p_pos_left;
            split_error[m,2] = p_neg_right; # prob for a negative input we falsely choose right of zero

            if(verbosity > 1)
              print(" m: $m, error_rate: ", error_rate[m],", error_left: ", split_error[m,1], ", error_right: ", split_error[m,2], "\n")
            end
          end # loop over m:no_points

          # do a sort to enforce increasing error rates
          A = [error_rate x];
          A = sortrows(A, by=x->(x[1],-x[2]), rev=false)

          # linear interpolator from target error rate back to value on input space producing this error rate
          z = InterpIrregular(A[:,1], A[:,2], 1, InterpLinear); # could also use BCnan as non-interp value

          if(!split_output)
            if(verbosity > 1)
              print("n_within_block: $n_within_block, z: ", z[detection_threshold], "\n");
            end
            subjects[i,j].blocks[k].trial[l].error_threshold = z[detection_threshold];
            #return z[detection_threshold];
          else # not currently in use
            Al = [split_error[:,1] x];
            Al = sortrows(Al, by=x->(x[1],-x[2]), rev=false); #sort by ascending error, and descending distance from 0
            zl = InterpIrregular(Al[:,1],Al[:,2], 1, InterpLinear);

            Ar = [split_error[end:-1:1,2] x[end:-1:1]];
            Ar = sortrows(Ar, by=x->(x[1],-x[2]), rev=false);
            zr = InterpIrregular(Ar[:,1],Ar[:,2], 1, InterpLinear);

            print("Al: ",Al,"\n")
            print("Ar: ",Ar,"\n")
            if(verbosity > 1)
              print("n_within_block: $n_within_block, z: ", z[detection_threshold], ", zl: ", zl[detection_threshold], ", zr: ", zr[detection_threshold],"\n");
            end

            #return [zl[detection_threshold] zr[detection_threshold]];
          end

          # update statistics
          local_average_threshold += z[detection_threshold];
          local_average_task_threshold[task_id] += z[detection_threshold];
          local_n_task_within_block[task_id] += 1;
        end # loop over trials per block

        #print("av_threshold: $local_average_threshold, ")
        # save statistics at block level
        local_average_threshold = local_average_threshold / local_no_trials_per_block;
        #print("local_no_trials_per_block: $local_no_trials_per_block, av_threshold: $local_average_threshold\n")
        local_average_task_threshold = local_average_task_threshold ./ local_n_task_within_block;
        subjects[i,j].blocks[k].average_threshold = local_average_threshold;
        subjects[i,j].blocks[k].average_task_threshold = local_average_task_threshold;
      end # loop over blocks per experiment

    end # loop over subjects in an experiment class
  end # loop over experiments classes for a group of subjects
  print("Post-hoc calculation of thresholds completed.\n\nNote: subject currently in memory has changed\n");
end

function detect_threshold(tuning_type::TuningSelector, task_id::Int=1, split_output::Bool=false)
  # find the detection threshold with current weight matrix and current subject
  no_points = 30;
  error_rate = zeros(no_points);
  split_error = zeros(no_points,2)
  x = linspace(0,1,no_points); # linspace of x values for detection threshold
  i = 1;
  for xi in x
    # calculate pre for +/- xi
    local_pre_pos = pre(xi, task_id, tuning_type);
    local_pre_neg = pre(-xi, task_id, tuning_type);

    #print("DEBUG: $local_pre_pos, $local_pre_neg ")

    # calculate noise free post for xi
    #TODO: there is now a function for this, switch to using it
    local_noise_free_post_pos_left = sum(local_pre_pos[:,task_id].*w[:,1,task_id]);
    local_noise_free_post_pos_right = sum(local_pre_pos[:,task_id].*w[:,2,task_id]);

    # calculate noise free post for -xi
    local_noise_free_post_neg_left = sum(local_pre_neg[:,task_id].*w[:,1,task_id]);
    local_noise_free_post_neg_right = sum(local_pre_neg[:,task_id].*w[:,2,task_id]);
    if(verbosity > 2)
      print("DEBUG: $local_noise_free_post_pos_left, $local_noise_free_post_pos_right, ")
      print("$local_noise_free_post_neg_left, $local_noise_free_post_neg_right, ")
    end

    # probability, for a positive input (i) that we choose left
    p_pos_left = 0.5 + 0.5 * erf( (local_noise_free_post_pos_left - local_noise_free_post_pos_right) / (sqrt(output_noise_variance) * 2.0) );
    p_pos_right = (1. - p_pos_left);

    if(verbosity > 2)
      print("p_pos_left: $p_pos_left, p_pos_right: $p_pos_right ")
    end

    # probability, for a negative input (-i) that we choose left
    p_neg_left = 0.5 + 0.5 * erf( (local_noise_free_post_neg_left - local_noise_free_post_neg_right) / (sqrt(output_noise_variance) * 2.0) );
    p_neg_right = (1. - p_neg_left);
    if(verbosity > 2)
      print("p_neg_left: $p_neg_left, p_neg_right: $p_neg_right,")
    end

    # probability of choosing the wrong side
    error_rate[i] = ( p_pos_left + p_neg_right ) / 2.0;
    split_error[i,1] = p_pos_left;
    split_error[i,2] = p_neg_right; # prob for a negative input we falsely choose right of zero

    if(verbosity > 1)
      print(" i: $i, xi: $xi, error_rate: ", error_rate[i],", error_left: ", split_error[i,1], ", error_right: ", split_error[i,2], "\n")
    end
    i += 1;
  end

  # do a sort to enforce increasing error rates
  A = [error_rate x];
  A = sortrows(A, by=x->(x[1],-x[2]), rev=false)

  # linear interpolator from target error rate back to value on input space producing this error rate
  z = InterpIrregular(A[:,1], A[:,2], 1, InterpLinear); # could also use BCnan as non-interp value

  if(!split_output)
    if(verbosity > 1)
      print("n_within_block: $n_within_block, z: ", z[detection_threshold], "\n");
    end
    return z[detection_threshold];
  else
    Al = [split_error[:,1] x];
    Al = sortrows(Al, by=x->(x[1],-x[2]), rev=false); #sort by ascending error, and descending distance from 0
    zl = InterpIrregular(Al[:,1],Al[:,2], 1, InterpLinear);

    Ar = [split_error[end:-1:1,2] x[end:-1:1]];
    Ar = sortrows(Ar, by=x->(x[1],-x[2]), rev=false);
    zr = InterpIrregular(Ar[:,1],Ar[:,2], 1, InterpLinear);

    print("Al: ",Al,"\n")
    print("Ar: ",Ar,"\n")
    if(verbosity > 1)
      print("n_within_block: $n_within_block, z: ", z[detection_threshold], ", zl: ", zl[detection_threshold], ", zr: ", zr[detection_threshold],"\n");
    end

    return [zl[detection_threshold] zr[detection_threshold]];
  end
end


# a single post-synaptic neuron only contributes (1/N) to the decision process
#   the purpose is to decouple the decision making somewhat from an individual
#   cell firing rate
function pooled_population_post_rate(x::Float64, task_id::Int, tuning_type::TuningSelector, local_post::Array{Float64,2})
  (left,right) = noise_free_post(x, task_id, tuning_type);
  pop_noise_free_post = [left right];
  # pop_rate = local_post + (no_pop_scaling_post_neurons - 1) .* pop_noise_free_post + ( sqrt(no_pop_scaling_post_neurons * (no_pop_scaling_post_neurons - 1)) .* transpose(population_ksi) );
  # I think the following is the correct scaling of population noise
  pop_rate = local_post + ( (no_pop_scaling_post_neurons - 1) .* pop_noise_free_post ) + ( sqrt( no_pop_scaling_post_neurons * (no_pop_scaling_post_neurons - 1) ) .* transpose(population_ksi) );

  return pop_rate;
end


# this is the only function which actually knows if things went right or wrong
# instance_correct = 0;
# instance_incorrect = 0;
function reward(x::Float64, task_id::Int, tuning_type::TuningSelector)
	local_post = post(x, task_id, tuning_type, true);

  # pooling of decision across a population of, per decision class, post-synaptic neurons
  if(use_pooled_scaling_of_post_population_for_decisions)
    local_post = pooled_population_post_rate(x, task_id, tuning_type, local_post);
    # local_post is getting overwritten her by pooled post rate variable for decision making
    #   this will not influence the post-synaptic firing rate in the weight update equation
  end

  if(disable_winner_takes_all)
    # Disabling winner-takes-all requires a change in logic: we just
    #   take the output neuron with the higher firing rate (even if
    #   it's negative).
    if( (x > 0) && (local_post[2] > local_post[1]) )
      return (1);
    elseif( (x <= 0) && (local_post[1] >= local_post[2]) )
      return (1);
    else
      return (-1);
    end
  else # maintaining winner-takes-all logic
    # I've had some trouble with the logic here due to wta() accepting negative inputs
    #   originally used abs() > 0 here for wta logic [which worked]
    #   changed to abs(a) > abs(b) when developing non-wta logic
	 if ( (x > 0) && (abs(local_post[2]) > abs(local_post[1]) ) ) #(local_post[2] > local_post[1]) )#right
    if(verbosity > 1)
      global instance_correct += 1;
      print("Greater than zero (x: $x)\n")
    end
		return (1);
	 elseif ( (x <= 0) && (abs(local_post[1]) >= abs(local_post[2]) ) ) #(local_post[1] > local_post[2]) )#left
    if(verbosity > 1)
      instance_correct += 1;
      print("Less than zero (x: $x)\n")
    end
		return (1);
	 else
    if(verbosity > 1)
      global instance_incorrect += 1;
      print("WRONG\n")
    end
	 	return (-1);
	 end
  end
end


# individual critics for running rewards
# no_task_critics = 2
# no_choices_per_task_critics = 2
# initialise
#  n = 0
#  average_reward = 0
function multi_critic_running_av_reward(R, task_critic_id::Int, choice_critic_id::Int)
  # note regarding julia deprecation of :: Notation
  # my use of the notation should come back into use in
  # future versions of the language, so don't bother to remove
  global n_critic;
  typeassert(n_critic, Array{Int, 2});
  global average_reward;
  typeassert(average_reward, Array{Float64,2});

  tau_r = running_av_window_length;

  if (n_critic[task_critic_id, choice_critic_id] < tau_r)
    n_critic[task_critic_id, choice_critic_id] += 1;
  end

  tau = min(tau_r, n_critic[task_critic_id, choice_critic_id]);

  Rn = ( (tau - 1) * average_reward[task_critic_id, choice_critic_id] + R ) / tau;

  # update average_reward monitor
  average_reward[task_critic_id, choice_critic_id] = Rn;

  return Rn :: Float64;
end


function update_weights(x::Float64, task_id::Int, tuning_type::TuningSelector, trial_dat::Trial)
  if(verbosity > 3)
    global instance_reward;
    global instance_average_reward;
  end
  global n_within_block += 1;
  global n_task_within_block;
  n_task_within_block[task_id] += 1;

  # don't forget to update noise externally to this function on separate iterations
  local_pre = pre(x, task_id, tuning_type);
  # Note: local_post returns a tuple where one value is 0 [in wta mode]. All comparisons to find the non zero value should use absolute comparison.
  local_post = post(x, task_id, tuning_type) :: Array{Float64,2};
  # reward() is now invariant again
  local_reward = reward(x, task_id, tuning_type) :: Int; # it is important that noise is not updated between calls to post() and reward()

  if(use_pooled_scaling_of_post_population_for_decisions)
    pop_rate = pooled_population_post_rate(x, task_id, tuning_type, local_post);
    # we're keeping local_post and pop_rate separate to allow for different decision processes below
  else
    pop_rate = local_post;
  end

  if(perform_detection_threshold)
    local_threshold = detect_threshold(tuning_type, task_id);
    trial_dat.error_threshold = local_threshold;
  end
  if(verbosity > 3)
    instance_reward[n] = local_reward;
  end

  # Save some data for later examination
  trial_dat.task_type = task_id; #(is_problem_1 ? 1 : 2);
  trial_dat.correct_answer = x #(x > 0 ? 1 : -1);
  if(disable_winner_takes_all)
    #=trial_dat.chosen_answer = ( (local_post[2] > local_post[1]) ? 1 : -1)
    local_choice = ( (local_post[2] > local_post[1]) ? 2 : 1);=#
    trial_dat.chosen_answer = ( (pop_rate[2] > pop_rate[1]) ? 1 : -1)
    local_choice = ( (pop_rate[2] > pop_rate[1]) ? 2 : 1);
    if(debug_print_now)
      print("chosen_answer: ", trial_dat.chosen_answer, "\n");
    end
  else
    #print("Error\n")
    #=trial_dat.chosen_answer = ( (abs(local_post[2]) > abs(local_post[1])) ? 1 : -1)
    local_choice = ( (abs(local_post[2]) > abs(local_post[1])) ? 2 : 1);=#
    trial_dat.chosen_answer = ( (abs(pop_rate[2]) > abs(pop_rate[1])) ? 1 : -1)
    local_choice = ( (abs(pop_rate[2]) > abs(pop_rate[1])) ? 2 : 1);
  end
  trial_dat.got_it_right = ((local_reward > 0) ? true : false);

  # monitor average choice per block here
  #   using n independent of critic, for now
  #local_choice = (abs(local_post[1]) > 0 ? 1 : 2);
  global average_choice = ( (n_within_block-1) * average_choice + local_choice ) / (n_within_block);
  global average_block_reward = ( ( n_within_block - 1) * average_block_reward + local_reward ) / (n_within_block);
  global average_task_reward;
  average_task_reward[task_id] = ( (n_task_within_block[task_id] - 1) * average_task_reward[task_id] + local_reward) / (n_task_within_block[task_id]);


  ## The Critic: Reward Predictor
  #running_av_reward(local_reward); # Nicolas is doing this before dw update, so first timestep is less unstable...
  # TODO: need to improve critic response axis and task type bin logic here
  # Binning along input-output/selection choice axis
  local_critic_response_bin = 1::Int;
  if(no_choices_per_task_critics > 1) #current logic has max 2 critics on this axis of choice
    local_correct_answer = (x > 0 ? 2 : 1);
    local_critic_response_bin = local_correct_answer;
  end
  # Binning along task type axis
  local_critic_task_bin = 1::Int;
  if(no_task_critics > 1)
    local_critic_task_bin = trial_dat.task_type;
  end
  if use_hard_coded_critic
    # as in Fremaux paper
    multi_critic_running_av_reward(local_reward, local_critic_task_bin, local_critic_response_bin);
  else
    #TODO: link to my new critic representation learner
    # needs task_id, to form representation
    # and reward received, to form association
    update_critic_representation(task_id, local_reward, true);
  end

  # decide which form of average reward we're using here:
  if use_hard_coded_critic
    if( use_multi_critic )
      #   multi critic:
      #     currently expanding logic to use number of critics declared in params
      local_average_reward = average_reward[local_critic_task_bin, local_critic_response_bin];
    elseif (use_single_global_critic)
      #   single global critic:
      local_sum_reward = 0.;
      local_sum_critics = 0;
      for i=1:no_task_critics
        for j = 1:no_choices_per_task_critics
          local_sum_reward += average_reward[i, j] * n_critic[i,j];
          local_sum_critics += n_critic[i,j];
        end
      end
      local_average_reward = ( local_sum_reward / local_sum_critics );
    else
      # Rmax (no running average):
      #TODO: actually this is not Rmax, that would also require (post-\bar{post}) in the dw formula
      local_average_reward = 0.;
    end
  else # use critic representation learner version
    local_average_reward = get_reward_prediction(task_id,true)[1]; # note this returns a single element array
  end

  # the weight update matrix
  dw = zeros(no_pre_neurons_per_task, no_post_neurons, no_input_tasks);
  if(binary_outputs_mode) # chosen output is 1, other output is 0
    binary_post = wta(local_post[1], local_post[2]);
    if (maximum(binary_post) != 0)
      binary_post /= maximum(binary_post);
    else
      #print("DEBUG: NaN avoided\n")
    end
    dw[:,1,task_id] = learning_rate * local_pre[:,task_id] * binary_post[1] * (local_reward - local_average_reward);
    dw[:,2,task_id] = learning_rate * local_pre[:,task_id] * binary_post[2] * (local_reward - local_average_reward);
  elseif(rescaled_outputs_mode)
    big_post = maximum(local_post);
    if( !(big_post == 0) ) # avoid divide by zero
      dw[:,1,task_id] = learning_rate * local_pre[:,task_id] * (local_post[1] / big_post) * (local_reward - local_average_reward);
      dw[:,2,task_id] = learning_rate * local_pre[:,task_id] * (local_post[2] / big_post) * (local_reward - local_average_reward);
    else
      # CONSIDER: what if big_post == 0 ?
      #   could use epsillon
      #   could just assume both outputs are 0 and do no weight change
      # print("big_post == 0 occurred\n");
    end
  else # original unbounded post-synaptic firing mode
    #this one stays
    dw[:,1,task_id] = learning_rate * local_pre[:,task_id] * local_post[1] * (local_reward - local_average_reward);
    dw[:,2,task_id] = learning_rate * local_pre[:,task_id] * local_post[2] * (local_reward - local_average_reward);
  end
  # Save some data for later examination
  trial_dat.reward_received = (local_reward - local_average_reward);
  trial_dat.w = deepcopy(w);
  trial_dat.dw = deepcopy(dw);


  if (verbosity > 3)
    instance_average_reward[n_within_block] = local_average_reward;
    reward_signal = (local_reward - local_average_reward)
    print("local_reward-average_reward: $local_reward - $average_reward = $reward_signal\n")
    #TODO: can I find a better measure of dw here?
    l1_norm_dw = sum(abs(dw));
    el_sum_dw = sum(dw);
    print("l1_norm_dw: $l1_norm_dw, element sum dw: $el_sum_dw\n")
    l1_norm_w = sum(abs(w))
    el_sum_w = sum(w)
    print("l1 norm w: $l1_norm_w, element sum w: $el_sum_w\n")
    left_sum_w = sum(w[:,1]);
    right_sum_w = sum(w[:,2]);
    print("before weight change, sum w left: $left_sum_w, sum w right: $right_sum_w\n")
  end

  trial_dat.mag_dw = sum(abs.(dw));

  # the weight update
  if(enable_weight_updates)
    if (!use_soft_bounded_weight_rule)
      global w += dw;
    else
      update_array = deepcopy(dw);
      update_array /= learning_rate;
      # be careful modifying the following logic
      update_array[update_array.>0.0] .*= (weights_upper_bound - w[update_array.>0.0]);
      update_array[update_array.<0.0] .*= (weights_lower_bound - w[update_array.<0.0]);
      global w += (update_array * learning_rate) / weights_upper_bound;
    end

    if (verbosity > 3)
      #for now at least these sums are across all tasks, could make them task specific
      left_sum_w = sum(w[:,1,:]);
      right_sum_w = sum(w[:,2,:]);
      print("after weight change, sum w left: $left_sum_w, sum w right: $right_sum_w\n")
    end


    # weight normalisation according to quadratic norm, multiplicative rule
    if(use_multiplicative_weight_normalisation)
      output_1_weights_norm = sqrt( mean(w[:,1,:].^2) ) #sqrt( sum(w[:,1,1].^2 ) + sum( w[:,1,2].^2) )
      output_2_weights_norm = sqrt( mean(w[:,2,:].^2) ) #sqrt( sum(w[:,2,1].^2 ) + sum( w[:,2,2].^2) )
      w[:, 1, :] = ( w[:, 1, :] ./ output_1_weights_norm ) * subject_initial_weight_scale;
      w[:, 2, :] = ( w[:, 2, :] ./ output_2_weights_norm ) * subject_initial_weight_scale;
    elseif(use_per_task_multiplicative_weight_normalisation)
      output_1_task_1_weights_norm = sqrt( mean(w[:,1,1].^2) ) #sqrt( sum(w[:,1,1].^2 ) + sum( w[:,1,2].^2) )
      output_2_task_1_weights_norm = sqrt( mean(w[:,2,1].^2) ) #sqrt( sum(w[:,2,1].^2 ) + sum( w[:,2,2].^2) )
      output_1_task_2_weights_norm = sqrt( mean(w[:,1,2].^2) ) #sqrt( sum(w[:,1,1].^2 ) + sum( w[:,1,2].^2) )
      output_2_task_2_weights_norm = sqrt( mean(w[:,2,2].^2) ) #sqrt( sum(w[:,2,1].^2 ) + sum( w[:,2,2].^2) )
      w[:, 1, 1] = ( w[:, 1, 1] ./ output_1_task_1_weights_norm ) * subject_task_1_initial_weight_scale;#  subject_initial_weight_scale;
      w[:, 2, 1] = ( w[:, 2, 1] ./ output_2_task_1_weights_norm ) * subject_task_1_initial_weight_scale;#  subject_initial_weight_scale;
      w[:, 1, 2] = ( w[:, 1, 2] ./ output_1_task_2_weights_norm ) * subject_task_1_initial_weight_scale;#  subject_initial_weight_scale;
      w[:, 2, 2] = ( w[:, 2, 2] ./ output_2_task_2_weights_norm ) * subject_task_1_initial_weight_scale;#  subject_initial_weight_scale;
    elseif(use_subtractive_weight_normalisation)
      # normalise via sigma w
      output_1_weights_norm = mean(dw[:,1,:]);
      output_2_weights_norm = mean(dw[:,2,:]);
      w[:,1,:] = w[:,1,:] - output_1_weights_norm;
      w[:,2,:] = w[:,2,:] - output_2_weights_norm;
    end


    # hard bound weights at +/- 10
    w[w .> weights_upper_bound] = weights_upper_bound;
    w[w .< weights_lower_bound] = weights_lower_bound;
  end # enable_weight_updates


  # update the running average of post-synaptic firing rates (only once per trial and either before or after all calls to post() )
  update_intrinsic_excitability(x, task_id, tuning_type);
  # decision_criterion_learner tries to ensure a 50:50 left:right ratio
  #update_decision_criterion_learner(local_post);
  update_decision_criterion_learner(pop_rate, task_id);
  trial_dat.decision_criterion_monitor = deepcopy(decision_criterion_monitor);

  return (local_reward+1); # make it 0 or 2, rather than +/-1
end