text:
using DataFrames
using CSV
import Base.getindex
export Chromosome, print_chromosome, getindex, random_chromosome, mutate_chromosome!, mutate_all, PredType
export num_mutate_locations, set_active_to_false, fraction_active, check_recursive, node_values
export output_values, number_active, number_active_gates, hamming_distance, ihamming_distance, hamming, deactivate_chromosome!
export copy_chromosome!, mutational_robustness, fault_tolerance_fitness, number_active_old
export build_chromosome, Input_node, Int_node, Output_node, print_build_chromosome, circuit_code, circuit_int
export circuit, print_circuit
export circuit_distance, remove_inactive, count_circuits
export code_to_circuit   # outdate 9/13/21
export insert_gate!, delete_gate!
export enumerate_circuits, chromosome_to_int, gate_int, gate_int_list, int_to_gate, int_to_chromosome

PredType = Int64
#= Commented out.  The included definition is in aliases.jl.
mutable struct Chromosome
    params::Parameters
    inputs::Vector{InputNode}
    interiors::Vector{InteriorNode}
    outputs::Vector{OutputNode}
    fitness::Float64
    robustness::Union{Float64,PredType}
end
=#

CPopulation = Vector{Chromosome}

function Chromosome(p::Parameters; ident::PredType=PredType(0))
    inputs = Array{InputNode}( undef, p.numinputs)
    interiors = Array{InteriorNode}( undef, p.numinteriors )
    outputs = Array{OutputNode}( undef, p.numoutputs)
    fitness = 0.0
    if ident == PredType(0)
      robustness = 0.0
    else
      robustness = ident
    end
    return Chromosome(p, inputs, interiors, outputs, fitness, robustness )
end

function random_chromosome(p::Parameters; ident::PredType=PredType(0))
  funcs = default_funcs(p.numinputs)
  random_chromosome(p,funcs, ident=ident)
end

function random_chromosome(p::Parameters, funcs::Vector{Func}; ident::PredType=PredType(0))
    #println("Creating random Chromosome")
    c = Chromosome(p)

    for index = 1:p.numinputs
        c.inputs[index] = InputNode(index)
        c.inputs[index].active = false
        c.inputs[index].cache = 0
    end
    for index = (p.numinputs +1):(p.numinputs+p.numinteriors)
        maxindex = index - 1
        minindex = max(1,index-p.numlevelsback)
        #println("index: ",index,"  minindex: ",minindex,"  maxindex: ",maxindex)
        func = funcs[rand(1:end)]
        inputs = Array{Int64}(undef, func.arity)
        for i = 1:func.arity
            inputs[i] = rand(minindex:maxindex)
            #println("i: ",i,"  inputs[i]: ",inputs[i])
        end
        c.interiors[index - p.numinputs] = InteriorNode(func, inputs)
    end
    minindex = max(1,p.numinputs + p.numinteriors - p.numlevelsback + 1)
    maxindex = p.numinputs + p.numinteriors
    #println("output (minindex, maxindex): ",(minindex,maxindex))
    for i = 1:length(c.outputs)
        #index = rand(minindex:maxindex)
        index = p.numinputs + p.numinteriors + i - p.numoutputs   # use the last numoutputs interiors
        #println("output ",i," index: ",index)
        c.outputs[i] = OutputNode(index)
        c[index].active = true  #  Output nodes are always active
    end
    set_active_to_false(c)
    if ident != PredType(0)
      c.robustness = ident
    end
    return c
end

function set_active_to_false( c::Chromosome )
  for index = 1:(c.params.numinputs+c.params.numinteriors)
    c[index].active = false
    c[index].cache = MyInt(0)
  end
end

# mutates chromosome c by changing a random node or function
# if length(funcs) == 1, then only connections are mutated
# Returns new chromosome c and Bool active which is true if the node mutated was active
# BUG!! as of 3/3/21:  Does not modify c, instead returns the mutated c.  Probably OK:  3/8/21
function mutate_chromosome!( c::Chromosome, funcs::Vector{Func}, mutate_location::Int64=0;
    insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  #println("mutate_chromosome! insert_gate_prob: ",insert_gate_prob)
  if rand() < insert_gate_prob && c.params.numinteriors < MyIntBits( MyInt )
    return(insert_gate!(c),false)
  elseif rand() < delete_gate_prob && c.params.numinteriors > c.params.numinputs
    # Don't reduce number of gates less than number of interiors
    return (delete_gate!(c),true)  # active set to true, but has no meaning here
  end
  num_mutate_locs = num_mutate_locations( c, funcs )
  if mutate_location == 0
    mutate_location = rand(1:num_mutate_locs)
  end
  #println("num_mutate_locations: ",num_mutate_locs,"  mutate_location: ",mutate_location)
  interiors_inputs_list = [ c.interiors[i].inputs for i = 1:c.params.numinteriors ]
  num_inputs_list = map(length, interiors_inputs_list)  # A list of the number of inputs of the interior node functions
  #println("num_inputs_list: ",num_inputs_list)
  num_funcs_to_mutate = length(funcs) > 1 ? c.params.numinteriors : 0
  if mutate_location <= num_funcs_to_mutate   # mutate a func
    #println("mutate a func ml:",mutate_location)
    active = c[mutate_location].active
    new_func_index = rand(1:length(funcs))
    #println("new_funct_index: ",new_func_index,"  funcs[new_func_index].func: ",funcs[new_func_index].func)
    #println("  c.interiors[mutate_location].func: ",c.interiors[mutate_location].func.func)
    #println("condition: ", funcs[new_func_index].func == c.interiors[mutate_location].func.func )
    while( funcs[new_func_index].func == c.interiors[mutate_location].func.func )
      new_func_index = rand(1:length(funcs))
      #println("while new_funct_index: ",new_func_index)
    end
    #println("new_func_index: ",new_func_index)
    old_func_arity = c.interiors[mutate_location].func.arity
    old_inputs = c.interiors[mutate_location].inputs
    new_func = funcs[new_func_index]
    new_inputs = Array{Int64}(undef, new_func.arity)
    for i = 1:new_func.arity
      if i <= length(old_inputs)
        new_inputs[i] = old_inputs[i]
      else
        maxindex = mutate_location + c.params.numinputs - 1
        minindex = max(1,mutate_location + c.params.numinputs - c.params.numlevelsback)
        #println("mut_location: ",mutate_location,"  minindex: ",minindex,"  maxindex: ",maxindex)
        new_inputs[i] = rand(minindex:maxindex)
      end
    end
    c.interiors[mutate_location] = InteriorNode(new_func,new_inputs)
    #println("new interior node: ",c.interiors[mutate_location])
    #println("new interior node: ",c[mutate_location+c.params.numinputs])
    #println("func ",rand(1:100))
  elseif mutate_location <= num_funcs_to_mutate + sum(num_inputs_list)  # mutate an interior
    ml = mutate_location
    #println("mutate a interior node ml:",mutate_location)
    interior_mutate_location = mutate_location - num_funcs_to_mutate
    #println("interior mutate_location: ",interior_mutate_location)
    i = 1  # i will be the index of the interior node one of whose inputs will be mutated
    num_inputs_sum = num_inputs_list[1]
    #println("i: ",i,"  num_inputs_sum: ",num_inputs_sum)
    while(num_inputs_sum < interior_mutate_location)
      i += 1
      num_inputs_sum += num_inputs_list[i]
      #println("i: ",i,"  num_inputs_sum: ",num_inputs_sum)
    end
    #println("i: ",i,"  c.interiors[i]: ",c.interiors[i])
    active = c.interiors[i].active
    j  = num_inputs_sum - interior_mutate_location + 1  # j is the index of the input to mutate
    #println("(i,j): ",(i,j),"  num_inputs_sum: ",num_inputs_sum)
    maxindex = i + c.params.numinputs -1
    minindex = max(1, i + c.params.numinputs - c.params.numlevelsback)
    intindex = rand(minindex:maxindex)
    if maxindex > minindex
      while( intindex == c.interiors[i].inputs[j] )
        intindex = rand(minindex:maxindex)
      end
    end
    #println("minindex: ",minindex,"  maxindex: ",maxindex,"  intindex: ",intindex)
    c.interiors[i].inputs[j] = intindex
    #println("(i,j): ",(i,j),"  c.interiors[i].inputs: ",c.interiors[i].inputs)
    #print("ml: ",mutate_location,"  ")
    #print_circuit(c)
    #println("inputs: ",rand(1:100))
  elseif mutate_location > num_funcs_to_mutate + sum(num_inputs_list)   # mutate an output
    # Mutating an output node should never happen because random chromosomes always use that last interior nodes as input,
    #   and num_mutate_locations() returns a value that does not include mutation of output nodes.
    error("mutate_chromosome!() is mutating an output node which should never happen")
    #println("mutate the output: ",mutate_location - num_funcs_to_mutate - sum(num_inputs_list ))
    active = true   # outputs are always active
    minindex = max(1, num_funcs_to_mutate + c.params.numinputs - c.params.numlevelsback + 1 )
    maxindex = num_funcs_to_mutate + c.params.numinputs
    if maxindex == mutate_location
      print_chromosome(c)
      exit()
    end
    #println("minindex: ",minindex,"  maxindex: ",maxindex)
    if maxindex > minindex   # no change if there is only one alternative for the output node index
      outindex = rand(minindex:maxindex)
      #println("outindex: ",outindex)
      #println("num_funcs_to_mutate: ", num_funcs_to_mutate )
      #println("sum(num_inputs_list): ", sum(num_inputs_list))
      while outindex == c.outputs[mutate_location - num_funcs_to_mutate - sum(num_inputs_list)].input
        outindex = rand(minindex:maxindex)
      end
      c.outputs[mutate_location - num_funcs_to_mutate - sum(num_inputs_list)].input = outindex
    end
  end
  #println("active: ",active)
  set_active_to_false(c)
  (c,active)
end

# Mutates c in all possible ways and returns a list of the outputs of the mutated chromosomes
#    or the pair (avg_robustness, evolvability)
# Deterministic if all functions in default_funcs() have the same arity
# If robustness_only==true returns the pair: (avg_robustness, evolvability)
# Otherwise returns vector of either outputs, or chromosomes, or (outputs, chromosomes) pairs
function mutate_all( c::Chromosome, funcs::Vector{Func};
      robustness_only::Bool=false, output_outputs::Bool=true, output_chromosomes::Bool=false )
  #println("mutate_all: numlevelsback: ", c.params.numlevelsback )
  #sav_c = deepcopy(c)
  if robustness_only
    output_outputs = output_chromosomes = false
    robustness_sum = 0.0; robustness_count = 0
    result = Vector{MyInt}[]   # Save outputs to give a measure of evolvability
  elseif output_outputs && !output_chromosomes
    result = Vector{MyInt}[]
  elseif !output_outputs && output_chromosomes
    result = Chromosome[]
  elseif output_outputs && output_chromosomes
    result = Tuple{Vector{MyInt},Chromosome}[]
  end
  context = construct_context(c.params.numinputs)
  orig_output = output_values(c)
  num_mutate_locs = num_mutate_locations( c, funcs )
  interiors_inputs_list = [ c.interiors[i].inputs for i = 1:c.params.numinteriors ]
  num_inputs_list = map(length, interiors_inputs_list)
  num_funcs_to_mutate = length(funcs) > 1 ? c.params.numinteriors : 0
  for mutate_location = 1:num_mutate_locs
    new_c = output_chromosomes ? deepcopy(c) : c
    if mutate_location <= num_funcs_to_mutate   # mutate a func
      #println("mutate a func ml:",mutate_location,"  node: ",c.interiors[mutate_location])
      active = c[mutate_location].active
      for new_func_index = 1:length(funcs)
        # Does not generate all alternative mutations in the very unusual case when the new function has greater arity
        #    than the old function
        #println("new_func: ",funcs[new_func_index].func,"  old_func: ", c.interiors[mutate_location].func.func )
        if funcs[new_func_index].func == c.interiors[mutate_location].func.func
          continue   # skip
        end
        #println("new_funct_index: ",new_func_index,"  funcs[new_func_index].func: ",funcs[new_func_index].func)
        #println("  c.interiors[mutate_location].func: ",c.interiors[mutate_location].func.func)
        #println("condition: ", funcs[new_func_index].func == c.interiors[mutate_location].func.func )
        #println("new_func_index: ",new_func_index)
        old_func_arity = c.interiors[mutate_location].func.arity
        old_inputs = c.interiors[mutate_location].inputs
        new_func = funcs[new_func_index]
        new_inputs = Array{Int64}(undef, new_func.arity)
        for i = 1:new_func.arity
          if i <= length(old_inputs)
            new_inputs[i] = old_inputs[i]
          else
            maxindex = mutate_location + c.params.numinputs - 1
            minindex = max(1,mutate_location + c.params.numinputs - c.params.numlevelsback)
            #println("mut_location: ",mutate_location,"  minindex: ",minindex,"  maxindex: ",maxindex)
            new_inputs[i] = rand(minindex:maxindex)
          end
        end
        new_c = deepcopy(c)  # This is necessary:  test done on 9/24/21
        #new_c = output_chromosomes ? deepcopy(c) : c
        new_c.interiors[mutate_location] = InteriorNode(new_func,new_inputs)
        deactivate_chromosome!(new_c)
        new_output = execute_chromosome(new_c,context)
        if robustness_only
          robustness_sum = (new_output == orig_output) ? robustness_sum + 1 : robustness_sum
          robustness_count+=1
          push!(result,new_output)
        elseif output_outputs && !output_chromosomes
          push!(result,new_output)
        elseif !output_outputs && output_chromosomes
          push!(result,new_c)
        elseif output_outputs && output_chromosomes
          push!(result,(new_output,new_c))
        end
      end
      #println("new interior node: ",new_c.interiors[mutate_location])
    elseif mutate_location <= num_funcs_to_mutate + sum(num_inputs_list)  # mutate an interior
      #println("mutate an interior node ml:",mutate_location, "  node: ",c.interiors[mutate_location-num_funcs_to_mutate])
      interior_mutate_location = mutate_location - num_funcs_to_mutate
      #println("interior mutate_location: ",interior_mutate_location)
      i = 1  # i will be the index of the interior node one of whose inputs will be mutated
      num_inputs_sum = num_inputs_list[1]
      #println("i: ",i,"  num_inputs_sum: ",num_inputs_sum)
      while(num_inputs_sum < interior_mutate_location)
        i += 1
        num_inputs_sum += num_inputs_list[i]
        #println("i: ",i,"  num_inputs_sum: ",num_inputs_sum)
      end
      #println("i: ",i,"  c.interiors[i]: ", c.interiors[i])
      j  = num_inputs_sum - interior_mutate_location + 1  # j is the index of the input to mutate
      maxindex = i + c.params.numinputs -1
      minindex = max(1, i + c.params.numinputs - c.params.numlevelsback)
      #println("(i,j): ",(i,j),"  num_inputs_sum: ",num_inputs_sum,"  minidex: ",minindex,"  maxindex: ",maxindex)
      for intindex = minindex:maxindex
        #println("intindex: ",intindex,"  c.interiors[i].inputs[j]: ",c.interiors[i].inputs[j])
        if intindex == c.interiors[i].inputs[j]
          #println("continue")
          continue
        end
        new_c = deepcopy(c)   # Is this necessary???  Maybe only when the chromosome is returned
        new_c.interiors[i].inputs[j] = intindex
        #println("intindex: ",intindex,"  (i,j): ",(i,j),"  new_c.interiors[i].inputs: ",new_c.interiors[i].inputs)
        deactivate_chromosome!(new_c)
        new_output = execute_chromosome(new_c,context)
        if robustness_only
          robustness_sum = (new_output == orig_output) ? robustness_sum + 1 : robustness_sum
          robustness_count+=1
          push!(result,new_output)
        elseif output_outputs && !output_chromosomes
          push!(result,new_output)
        elseif !output_outputs && output_chromosomes
          push!(result,new_c)
        elseif output_outputs && output_chromosomes
          push!(result,(new_output,new_c))
        end
      end
    elseif mutate_location > num_funcs_to_mutate + sum(num_inputs_list)   # mutate an output
      # Mutating an output node should never happen because random chromosomes always use that last interior nodes as input,
      #   and num_mutate_locations() returns a value that does not include mutation of output nodes.
      error("mutate_all() is mutating an output node which should never happen")
      #println("mutate the output: ",mutate_location - num_funcs_to_mutate - sum(num_inputs_list ))
      #active = true   # outputs are always active
      minindex = max(1, num_funcs_to_mutate + c.params.numinputs - c.params.numlevelsback + 1 )
      maxindex = num_funcs_to_mutate + c.params.numinputs
      if maxindex == mutate_location
        print_chromosome(c)
        exit()
      end
      #println("minindex: ",minindex,"  maxindex: ",maxindex)
      if maxindex > minindex   # no change if there is only one alternative for the output node index
        for outindex = minindex:maxindex
          #outindex = rand(minindex:maxindex)
          if outindex == c.outputs[mutate_location - num_funcs_to_mutate - sum(num_inputs_list)].input
            continue
          end
        end
        println("outindex: ",outindex)
        new_c.outputs[mutate_location - num_funcs_to_mutate - sum(num_inputs_list)].input = outindex
        new_output = execute_chromosome(new_c,context)
        push!(result,new_c)
      end
    end
  end
  if robustness_only
    (robustness_sum/robustness_count, length(unique(result))/length(result)) # pair of avg robustenss and evolvability
  else
    result
  end
end

function num_mutate_locations( c::Chromosome, funcs::Vector{Func} )
  num_funcs_to_mutate = length(funcs) > 1 ? c.params.numinteriors : 0
  interiors_inputs_list = [ c.interiors[i].inputs for i = 1:c.params.numinteriors ]
  outputs_input_list = [ c.outputs[i].input for i = 1:c.params.numoutputs ]
  #println("int_list: ", interiors_inputs_list, "   out_list: ",outputs_input_list)
  num_inputs_list = map(length, interiors_inputs_list)
  #num_inputs = sum(num_inputs_list) + c.params.numoutputs
  num_inputs = sum(num_inputs_list)
  #println("num_inputs_list: ",num_inputs_list,"  num_inputs: ",num_inputs)
  if length(funcs) == 1
    num_mutate_locations = num_inputs
  else
    num_mutate_locations = num_funcs_to_mutate + num_inputs
  end
  num_mutate_locations
end

function print_chromosome(f::IO, c::Chromosome )
  for i = 1:(c.params.numinputs + c.params.numinteriors + c.params.numoutputs)
    println(f,"i: ",i,"  ",c[i])
  end
  check_recursive(c)
end

function print_chromosome( c::Chromosome )
  for i = 1:(c.params.numinputs + c.params.numinteriors + c.params.numoutputs)
    println("i: ",i,"  ",c[i])
  end
  check_recursive(c)
end

function check_recursive( c::Chromosome )
  for i = (c.params.numinputs + 1) : c.params.numinteriors
    for j = 1:length(c[i].inputs)
      #println("i: ",i,"  inputs[",j,"] = ",c[i].inputs[j])
      if i == j
        println("recursive ")
        exit()
      end
    end
  end
end

function number_active_old( c::Chromosome )
  num_active = 0
  for i = 1:c.params.numinputs
    if c.inputs[i].active
      num_active += 1
    end
  end
  for i = 1:c.params.numinteriors
    if c.interiors[i].active
      num_active += 1
    end
  end
  num_active
end

function number_active( c::Chromosome )
  if !c[c.outputs[1].input].active   # if chromosome has not been executed
    context = construct_context(c.params.numinputs)
    execute_chromosome(c,context)
  end
  num_act_in = reduce(+,[c.inputs[i].active for i = 1:length(c.inputs)])
  num_act_int =  reduce(+,[c.interiors[i].active for i = 1:length(c.interiors)])
  num_act_in + num_act_int
end

function number_active_gates( c::Chromosome )
  if !c[c.outputs[1].input].active   # if chromosome has not been executed
    context = construct_context(c.params.numinputs)
    execute_chromosome(c,context)
  end
  num_act_int =  reduce(+,[c.interiors[i].active for i = 1:length(c.interiors)])
end

# Requires that cache values are set in the various versions of evaluat_node()
function node_values( c::Chromosome )
  output_values = [c[c.outputs[i].input].cache for i = 1:length(c.outputs)]
  input_values = [c.inputs[i].cache for i = 1:length(c.inputs)]
  interior_values = [c.interiors[i].cache for i = 1:(length(c.interiors))]
  #output_values = [c.interiors[i].cache for i = (length(c.interiors)-length(c.outputs)+1):length(c.interiors)]
  (input_values,interior_values,output_values)
end

function output_values( c::Chromosome )
  if !c[c.outputs[1].input].active   # if chromosome has not been executed
    context = construct_context(c.params.numinputs)
    return execute_chromosome(c,context)
  else
    return [c[c.outputs[i].input].cache for i = 1:length(c.outputs)]
  end
end

function deactivate_chromosome!( c::Chromosome )
  for in in c.inputs
    in.active = false
  end
  for int in c.interiors
    int.active = false
  end
  c
end

function fraction_active( c::Chromosome )
  number_active(c)/(c.params.numinputs+c.params.numinteriors)
end

# Returns a chromosome with the same output values as c and with all
#  inactive gate nodes removed.
function remove_inactive( c::Chromosome )
  if !c[c.outputs[1].input].active   # if chromosome has not been executed
    context = construct_context(c.params.numinputs)
    execute_chromosome(c,context)
  end
  new_ints = InteriorNode[]
  numinputs = length(c.inputs)
  numints = length(c.interiors)
  int_inds = collect(1:numints)
  # if i refers to a node in c, then map_indices[i] references
  #   the corresponding node in the new chromosome
  map_indices = zeros(Int64, numinputs+numints )
  active_gate_indices = int_inds[ [c.interiors[i].active for i = 1:numints]]
  #map_indices = ttt(numinputs,numints,active_gate_indices)
  map_indices = collect(1:numinputs)
  j = 1
  for i = 1:numints
    if active_gate_indices[j] > i
      push!(map_indices,0)
    else
      push!(map_indices,j+numinputs)
      j+= 1
    end
  end
  #println("map_indices: ",map_indices)
  for j in active_gate_indices
    new_inputs = map(x->map_indices[x],c.interiors[j].inputs)
    #println("j: ",j,"  new_inputs: ",new_inputs)
    new_int = InteriorNode( c.interiors[j].func, new_inputs )
    push!(new_ints,new_int)
  end
  new_outs = [ OutputNode( map_indices[c.outputs[i].input])  for i = 1:length(c.outputs) ]
  p = c.params
  new_p = Parameters( p.numinputs, p.numoutputs, length(new_ints), p.numlevelsback )
  nc = Chromosome( new_p, c.inputs, new_ints, new_outs, c.fitness, c.robustness )
  @assert output_values(nc) == output_values(c)
  nc
end

function ttt(numinputs,numints,active_gate_indices)
  map_indices = collect(1:numinputs)
  j = 1
  for i = 1:numints
    if active_gate_indices[j] > i
      push!(map_indices,0)
    else
      push!(map_indices,j+numinputs)
      j+= 1
    end
  end
  map_indices
end

function getindex(c::Chromosome, index::Integer)
    if index <= c.params.numinputs   # input node
        return c.inputs[index]
    end

    if index > c.params.numinteriors + c.params.numinputs  # output node
        return c.outputs[index - c.params.numinteriors - c.params.numinputs]
    end

    return c.interiors[index-c.params.numinputs]  # interior node
end

# Number of bits where x and y differ.  x and y are interpreted as bit strings
function hamming( x::MyInt, y::MyInt )
  xr = xor( x, y )
  result = 0
  my_one = convert(MyInt,1)
  my_zero = convert(MyInt,0)
  while xr != my_zero
    result +=  xr & my_one
    xr >>= 1
  end
  result
end

# Hamming distance between MyInts which is between 0 and 1.
# If hamming(x,y) >= 2^numinputs/2, then hamming_distance(x,y,numinputs) = 1.0
function hamming_distance( x::MyInt, y::MyInt, numinputs::Int64 )
  result = hamming(x,y)/2^numinputs
end

# Hamming distance between Goals which is between 0 and 1.
function hamming_distance( x::Vector{MyInt}, y::Vector{MyInt}, numinputs::Int64 )
  @assert length( x ) == length( y )
  sum( hamming_distance( x[i], y[i], numinputs ) for i = 1:length(x))/length(x)
end

# the minimum of the hamming distance and 1+hamming_distance( x, bitwise_complement(y))
function ihamming_distance( x::MyInt, y::MyInt, numinputs::Int64 )
  Ones = Main.CGP.construct_ones(numinputs)[numinputs]
  complement_y = xor(y,Ones)
  #@printf("0x0%0x\n",complement_y)
  #println("hd:  ",hamming_distance(x,y,numinputs))
  #println("chd: ",hamming_distance(x,complement_y,numinputs))
  #println("xhd: ",1.0/numinputs+hamming_distance(x,complement_y,numinputs))
  return min(hamming_distance(x,y,numinputs),1.0/2^numinputs+hamming_distance(x,complement_y,numinputs))
end

function ihamming_distance( x::Vector{MyInt}, y::Vector{MyInt}, numinputs::Int64 )
  sum( ihamming_distance( x[i], y[i], numinputs ) for i = 1:length(x))/length(x)
end

# Change the fields of chromosome c to be the fields of chromosom c_to_copy
function copy_chromosome!( c::Chromosome, c_to_copy::Chromosome )
  c.params = c_to_copy.params
  c.inputs = c_to_copy.inputs
  c.interiors = c_to_copy.interiors
  c.outputs = c_to_copy.outputs
  c.interiors = c_to_copy.interiors
end

# Computes the fraction of mutations that do not change the output
# Not really accurate because it counts one mutation per num_mutate_location(), while in fact
#    there may be more than one random choice per num_mutate_location()
# Superceded by calling   mutate_all( c, funcs, robustness_only=true ) which returns the pair (robustness,evolvability)
# If active_only==true, then only mutations that change active nodes are considered
# If active_only==false, then all mutations are considered
# Note that results are partially random since the result of mutations at a specific location can be different.
function mutational_robustness( c::Chromosome, funcs::Vector{Func}; active_only::Bool=false )
  context = construct_context( c.params.numinputs )
  prev_out = execute_chromosome( c, context )
  count_no_change = 0
  count_no_change = 0.0
  num_mut_locs = 0
  for j = 1:num_mutate_locations( c, funcs )
    cc = deepcopy(c)
    (cc,active) = mutate_chromosome!( cc, funcs, j )
    out_c = execute_chromosome( cc, context )
    #println("j: ",j,"  active: ",active)
    #print_chromosome( cc )
    if !active_only || active # consider all mutations if !active_only, only active mutations if active_only
      #out_c = execute_chromosome( cc, context )
      count_no_change += ((out_c == prev_out) ? 1 : 0)
      if active_only && active
        num_mut_locs += 1
      end
    end
    #println("j: ",j,"  active: ",active,"  no_change: ",(out_c == prev_out),"  count_no_change: ",count_no_change)
  end
  #println("num_mut_locs: ",num_mut_locs,"  num_mutate_locations(): ",num_mutate_locations(c,funcs),"  count_no_change: ",count_no_change)
  #@assert num_mut_locs == num_mutate_locations(c,funcs)
  if !active_only
    num_mut_locs = num_mutate_locations( c, funcs )
  end
  count_no_change/num_mut_locs
end

# Implements equation 3.3 of Macia and Sole (2009)
#  The average Hamming deviation of the output from the unperturbed output under perturbation of the output of each node
function fault_tolerance_fitness( c::Chromosome )
  numinputs = c.params.numinputs
  numoutputs = c.params.numoutputs
  numints = c.params.numinteriors
  #println("(numinputs,numoutputs,numints):",(numinputs,numoutputs,numints))
  context = construct_context( c.params.numinputs )
  outputs = output_values( c )
  #println("outputs: ",outputs)
  distances = fill(0.0,numints)
  for i = (1+numinputs):(numinputs+numints)
    out_ft = execute_chromosome_ft(c,context,i)
    #println("i: ",i,"  out_ft: ",out_ft)
    #println("hamming dis: ",hamming_distance(execute_chromosome_ft(c,context,i)[1], outputs[1], numinputs))
    distances[i-numinputs] =
      sum( hamming_distance(execute_chromosome_ft(c,context,i)[k], outputs[k], numinputs) for k = 1:numoutputs)
  end
  #println("distances: ",distances)
  1.0-sum(distances)/numints/numoutputs
end

#Assumes nodearity==2
mutable struct Gate_node
  node_index::Int64
  func::Func
  input1::Int64
  input2::Int64
end

# Creates a Chromosome (circuit) from the concise format.
# Example:
# julia>circuit((1,2,3), ((4,OR,1,2), (5,AND,2,3), (6,XOR,4,5)))
#    creates a circuit with 3 input nodes and 3 gate nodes and 1 implicit output node.
# gate nodes have 4 fields:
#    node_index:  should be the index of the node.  Not used in construcing the circuit, but checked for correctness
#    node_function:  See Func.jl for options
#    node_input1:  an integer in the range from 1 to the index of the node minus 1
#    node_input2:  an integer in the range from 1 to the index of the node minus 1
# Assumes nodearity==2 and numoutputs == 1
function circuit( inputs::Tuple, gates::Tuple;
    levsback::Int64=length(inputs)+length(gates))
  numinputs = length(inputs)
  errors = false
  for i = 1:length(gates)
    if gates[i][1] != i+numinputs
      println("index of gate ",i," not correct: it should be: ",i+numinputs)
      errors = true
    end
  end
  for i = 1:length(gates)
    if !((1 <= gates[i][3])  && (gates[i][3] < numinputs+i))
      error("illegal first gate input for gate: ",i)
    end
  end
  for i = 1:length(gates)
    if !((1 <= gates[i][4])  && (gates[i][4] < numinputs+i))
      error("illegal second gate input for gate: ",i)
    end
  end
  p = Parameters( numinputs=length(inputs), numoutputs=1, numinteriors=length(gates), numlevelsback=levsback )
  in_nodes = [InputNode(i) for i in inputs]
  gate_nodes = [InteriorNode(g[2], [g[3],g[4]] ) for g in gates ]
  out_nodes = [OutputNode( length(inputs)+length(gates ) ) ]
  c = Chromosome( p, in_nodes, gate_nodes, out_nodes, 0.0, 0.0 )
  if errors
    println("correct input to this function should have been: ")
    print_circuit(c)
  end
  c
end

# Outputs the concise format of Chromosome (circuit) c to IO stream f
# print_node_tuple() and gate_tuple() are defined in Node.jl
function print_circuit( f::IO, c::Chromosome; include_fitness::Bool=false, include_robustness::Bool=false,
    include_pheno::Bool=false )
  @assert c.params.nodearity==2
  @assert c.params.numinputs>1
  print(f,"circuit(")
  print_node_tuple(f, c.inputs )
  print(f,",")
  gate_tuple(f, c.interiors, c.params.numinputs )
  if c.params.numoutputs>1
    print_node_tuple(f, c.outputs )
  end
  if include_fitness
    @printf(f,", %5.3f",c.fitness)
  end
  if include_robustness
    if typeof(c.robustness) == PredType
      print(f,", ",c.robustness)
    elseif typeof(c.robustness) == Float64
      @printf(f,", %5.3f",c.robustness)
    end
  end
  if include_pheno
    if c.params.numoutputs == 1
      @printf(f,", 0x%x",output_values(c)[1])
    else
      print(f,", [")
      for i = 1:(c.params.numoutputs-1)
        @printf(f,"0x%x, ",output_values(c)[i])
      end
      @printf(f,"0x%x]",output_values(c)[c.params.numoutputs])
    end
  end
  println(f,")")
end

function print_circuit( c::Chromosome; include_fitness::Bool=false, include_robustness::Bool=false, include_pheno::Bool=false )
  print_circuit( Base.stdout, c, include_fitness=include_fitness, include_robustness=include_robustness,
    include_pheno=include_pheno )
end

# The following struct defintions and function definitions build_chromosome() and print_build_chromosome()
#   have mostly been replaced by the above definitions of circuit() and print_circuit()
mutable struct Int_node
  func::Func
  inputs::Vector{Int64}
end

mutable struct Input_node
  index::Int64
end

mutable struct Output_node
  input::Integer
end

# Example call: build_chromosome( [Input_node(1),Input_node(2)], [Int_node(OR,[1,2]),Int_node(AND,[2,3])],[Output_node(4)])
# Depreciated 12/18/20
function build_chromosome( input_nodes::Vector{Input_node}, interior_nodes::Vector{Int_node}, output_nodes::Vector{Output_node};
    levsback::Int64=length(input_nodes)+length(interior_nodes))
  num_in = length(input_nodes)
  num_ints = length(interior_nodes)
  num_outs = length(output_nodes)
  p = Parameters( numinputs=num_in, numoutputs=num_outs, numinteriors=num_ints, numlevelsback=levsback )
  in_nodes = [InputNode(in_node.index) for in_node in input_nodes]
  int_nodes = [InteriorNode(int_node.func, int_node.inputs) for int_node in interior_nodes]
  out_nodes = [OutputNode(outnode.input) for outnode in output_nodes]
  Chromosome( p, in_nodes, int_nodes, out_nodes, 0.0, 0.0 )
end

# Example calls:  build_chromosome((1,2), ((OR,[1,2]),(AND,[2,3])),(4,))
# Example calls:  build_chromosome((1,2), ((OR,[1,2]),(AND,[2,3])),(4,),0.0)
# Depreciated 12/18/20
function build_chromosome( inputs::Tuple, ints::Tuple, outs::Tuple, fitness::Float64=0.0;
    levsback::Int64=length(inputs)+length(ints))
  num_in = length(inputs)
  num_ints = length(ints)
  num_outs = length(outs)
  p = Parameters( numinputs=num_in, numoutputs=num_outs, numinteriors=num_ints, numlevelsback=levsback )
  in_nodes = [InputNode(in_index) for in_index in inputs]
  int_nodes = [InteriorNode(int_pair[1], int_pair[2]) for int_pair in ints]
  out_nodes = [OutputNode(out_index) for out_index in outs]
  Chromosome( p, in_nodes, int_nodes, out_nodes, fitness, 0.0 )
end

function print_build_chromosome( f::IO, c::Chromosome; include_fitness::Bool=false )
  print(f, "build_chromosome(")
  print_node_tuple(f, c.inputs )
  print(f,",")
  print_node_tuple(f, c.interiors )
  print(f,",")
  print_node_tuple(f, c.outputs )
  if include_fitness
    println(f,", ",c.fitness)
  end
  println(f,")")
  if typeof(f) == IOStream  # Don't close Base.stdout since this kills julia
    close(f)
  end
end

function print_build_chromosome( c::Chromosome )
  print_build_chromosome( Base.stdout, c)
end

# Generates a Int128 integer corresponding to chromosome ch.  The integer is unique for parameters p.
# Does not take into account the levelsback parameter
# Assumes that outputs come from the last numoutputs interior nodes and thus don't affect the number of chromosomes.
# See test/testChromosome.jl for a test----but only where all gates have arity max_arity.
function chromosome_to_int( ch::Chromosome, funcs::Vector{Func}=default_funcs(ch.params.numinputs); maxarity::Int64=2 )
  result = Int128(0)
  for i = 1:length(ch.interiors)
    result += gate_int( i, ch, maxarity, funcs )
    multiplier = i < length(ch.interiors) ? length(funcs)*(ch.params.numinputs+i+1)^maxarity : 1
    result *= multiplier
    #result = i < length(ch.interiors) ? result*length(funcs)*(ch.params.numinputs+i+1)^maxarity : result
    #println("i: ",i,"  gate_int: ",gate_int( i, ch, maxarity, funcs ),"  multiplier: ",multiplier, "  result: ",result)
    #println("ci i: ",i,"  gate_int: ",gate_int( i, ch, maxarity, funcs ),"  result: ",result)
  end
  #println("result: ",result)
  result
end

function gate_int( i::Int64, ch::Chromosome, maxarity::Int64, funcs::Vector{Func} )
  numinputs = ch.params.numinputs
  funcs_int = findfirst(x->x==ch.interiors[i].func,funcs)[1]-1
  gate_inputs = ch.interiors[i].inputs
  il = inputsList( maxarity, numinputs+i )
  ni = length(il)   # multiplier which should be (numinputs+i)^maxarity
  #println("ni: ",ni,"  maxarity: ",maxarity,"  numinputs+i: ",numinputs+i)
  @assert ni == (numinputs+i)^maxarity
  inputs_int = findfirst(x->x==gate_inputs,il)
  funcs_int*ni + inputs_int
end

# Generates a list of all Int64 vectors of length numinputs and values from 1:maxval
function inputsList( numinputs::Int64, maxval::Int64 )
  if numinputs == 1
    return [ [i] for i = 1:maxval ]
  end
  result = inputsList( numinputs-1, maxval )
  new_result = Vector{Int64}[]
  for r in result
    for i = 1:maxval
      dcr = deepcopy(r)
      push!(dcr,i)
      push!(new_result,dcr)
      #println("i: ",i,"  dcr: ",dcr)
    end
  end
  #println("numimnputs: ",numinputs,"  new_result: ",new_result)
  new_result
end

# Produces a list of all circuits corresponding the paramters p --- except that it ignores the levelsback parameter
function enumerate_circuits( p::Parameters, funcs::Vector{Func}=default_funcs(p.numinputs); maxarity::Int64=2 )
  enumerate_circuits( p, p.numinteriors, funcs, maxarity=maxarity )
end

# Recursive helper function
function enumerate_circuits( p::Parameters, numints::Int64, funcs::Vector{Func}; maxarity::Int64=2 )
  if numints == 0
    result = [Chromosome( p, map(i->InputNode(i),collect(1:p.numinputs)), InteriorNode[], [OutputNode(p.numinputs)], 0.0, 0.0 )]
    return result
  end
  prev_result = enumerate_circuits( p, numints-1, funcs )
  result = Chromosome[]
  for prev_ch in prev_result
    for func in funcs
      for inputs in inputsList( maxarity, p.numinputs+numints-1 )
        new_interiors = vcat(prev_ch.interiors,[InteriorNode(func,inputs)])
        new_ch = Chromosome( p, prev_ch.inputs,new_interiors,[OutputNode(p.numinputs+numints)],0.0,0.0)
        push!(result,new_ch)
      end
    end
  end
  result
end

function int_to_chromosome( ch_int::Int128, p::Parameters, funcs::Vector{Func}=default_funcs(p.numinputs); maxarity::Int64=2 )
  inputs = [ InputNode(i) for i = 1:p.numinputs ]
  g_ints = gate_int_list( ch_int, p, maxarity, funcs )
  interiors = InteriorNode[]
  for i = 1:p.numinteriors
    push!(interiors, int_to_gate( g_ints[i], i, p, maxarity, funcs ) )
  end
  outputs = [ OutputNode(i) for i = (p.numinputs + p.numinteriors -p.numoutputs + 1):p.numinputs + p.numinteriors]
  Chromosome( p, inputs, interiors, outputs, 0.0, 0.0 )
end

# list of the gate_ints used by chromsosome_int
function gate_int_list( ch_int::Integer, p::Parameters, maxarity::Int64, funcs::Vector{Func}=default_funcs(p.numinputs) )
  result = zeros(Int64,p.numinteriors)
  multipliers = [ length(funcs)*(p.numinputs+i+1)^maxarity for i = 1:p.numinteriors-1 ]
  #println("multipliers: ",multipliers)
  for i = p.numinteriors:-1:2
    result[i] = ch_int % multipliers[i-1]
    ch_int = div(ch_int,multipliers[i-1])
  end
  result[1] = ch_int
  result
end

function int_to_gate( g_int::Int64, i::Int64, p::Parameters, maxarity::Int64, funcs::Vector{Func}=default_funcs(p.numinputs) )
  il = inputsList( maxarity, p.numinputs+i )
  ni = (p.numinputs+i)^maxarity
  InteriorNode( funcs[div(g_int,ni)+1], il[ g_int % ni ] )
end

#= Outdated:  replaced by the newer version below
# Returns an unique integer for i^{th} interior node (gate) of chromosome ch
# Ignores any limitations due to the levelsback parameter
function gate_int( i::Int64, ch::Chromosome, funcs::Vector{Func} )
  #println("funcs: ",funcs)
  result = (findall(x->x==ch.interiors[i].func,funcs)[1]-1)*(ch.params.numinputs+1)
  #println("gi result: ",result)
  for j = 1:ch.interiors[i].func.arity
    @assert ch.interiors[i].inputs[j] <= ch.params.numinputs+i+1
    result += ch.interiors[i].inputs[j]-1
    result = j < length(ch.interiors[i].inputs) ?  result*(ch.params.numinputs+i+1) : result
    #println("j: ",j,"  result: ",result)
  end
  #println("gate_int result: ",result)
  result
end
=#

# Outdated 9/13/21
# Returns a vector of integers which is a characterization of the chromosome for the parameters c.params
# The length of the result should be p.numinteriors*(1+p.nodearity).
# Not correct in that codes for gate inputs are sometimes too large.
function circuit_code( c::Chromosome )
  result = Int64[]
  funcs = default_funcs(c.params.numinputs)
  if length(c.interiors) == 0
    return Int64[]
  end
  for i = 1:c.params.numinteriors
    # Determine index of gate function
    j = 0
    while j < length(funcs) && funcs[j+1].func != c.interiors[i].func.func
      j += 1
    end
    if j > length(funcs)
      error("illegal function in function circuit_code()")
    end
    #println("i: ",i,"  function code: ",j)
    push!(result,j)
    # Add interior node inputs to result
    index = i + c.params.numinputs  # node index in chromosome
    for k = 1:length(c.interiors[i].inputs)
      code = max(c.interiors[i].inputs[k]-(index-c.params.numlevelsback),0)
      #println("i: ",i,"  k: ",k,"  code: ",code)
      push!(result, code)
    end
  end
  # Note:  assumes that the input fields of output nodes are set by default and thus don't need to be recorded
  result
end

# Outdated 9/13/21 Not correct
# Converts a circuit code to a Chromosome.   Inverse function to circuit_code().
function code_to_circuit( code::Vector{Int64}, p::Parameters )
  funcs = default_funcs( p.numinputs )
  interior_nodes = InteriorNode[]
  for i = 1:p.numinteriors
    j = 3*(i-1)+1   # index of code component
    func = funcs[code[j] + 1]
    index = i + p.numinputs  # node index in Chromosome
    input1 = code[j+1] + max(index-p.numlevelsback,0)
    input2 = code[j+2] + max(index-p.numlevelsback,0)
    #println("i: ",i,"  j: ",j,"  func: ",func,"  input1: ",input1,"  input2: ",input2)
    push!( interior_nodes, InteriorNode(func,[input1,input2]) )
  end
  input_nodes = [ InputNode(j) for j = 1:p.numinputs ]
  output_nodes = [ OutputNode(j) for j = (p.numinputs+p.numinteriors-p.numoutputs+1):(p.numinputs+p.numinteriors) ]
  c = Chromosome(p, input_nodes, interior_nodes, output_nodes, 0.0, 0.0 )
end

# Outdated  9/13/21 by chromosome_to_int()
# An integer that characterizes the chromosome
# Doesn't overflow for 11 gates, 8 numlevelsback
function circuit_int( c::Chromosome )
  circuit_int( circuit_code(c), c.params )
end

# Outdated  9/13/21 by chromosome_to_int()
# Not correct  See diary12_25.txt for example.
function circuit_int( c_code::Vector{Int64}, p::Parameters )
  result = Int128(0)
  multiplier = 1
  funcs = default_funcs( p.numinputs )
  k = 1
  for i = 1:p.numinteriors
    multiplier = length(funcs)
    result = c_code[k] + result*multiplier
    #print("i: ",i,"  k: ",k,"  multiplier: ",multiplier,"  c_code[k]: ",c_code[k])
    #println("  result: ",result)
    k += 1
    multiplier = min(p.numlevelsback,i-1+p.numinputs)
    for j = 1:p.nodearity
      #println("result: ",result,"  multiplier: ",multiplier,"  c_code[k]: ",c_code[k])
      result = c_code[k] + result*multiplier
      #print("i: ",i,"  j: ",j,"  k: ",k,"  multiplier: ",multiplier,"  c_code[k]: ",c_code[k])
      #println("  result: ",result)
      k += 1
    end
  end
  result
end

# Outdated:  replaced in functionality by int_to_chromosome()
function int_to_circuit_code( c_int::Integer, p::Parameters )
  c_int = Int128(c_int)
  c_code = zeros(Int64,3*p.numinteriors)
  k = 3*p.numinteriors
  for i = p.numinteriors:-1:1
    multiplier = min(p.numlevelsback,i-1+p.numinputs)
    println("i: ",i,"  multiplier: ",multiplier)
    for j = p.nodearity:-1:1
      c_int_mod= c_int % multiplier
      #prev_c_int =c_int
      c_int ÷= multiplier
      #c_int_mod = prev_c_int - c_int
      println("i: ",i,"  j: ",j,"  c_int: ",c_int,"  c_int_mod: ",c_int_mod)
      c_code[k] = c_int_mod
      k -= 1
    end
    multiplier = length(funcs)
    println("i: ",i,"  multiplier: ",multiplier)
    c_int_mod= c_int % multiplier
    c_int ÷= multiplier
    println("i: ",i,"  c_int: ",c_int,"  c_int_mod: ",c_int_mod)
    c_code[k] = c_int_mod
    k -= 1
  end
  c_code
end

# Hamming distance between chromosomes normalized to be between 0.0 and 1.0
function circuit_distance( c1::Chromosome, c2::Chromosome )
  code1 = circuit_code(c1)
  code2 = circuit_code(c2)
  @assert length(code1) == length(code2)
  diff_count = 0
  for i = 1:length(code1)
    diff_count += code1[i] == code2[i] ? 0 : 1
  end
  diff_count/length(code1)
end

# Return the number of circuits for parameters p and number of funcs nfuncs if nfuncs>0
# If nfuncs==0, nfuncs is reset to be length(default_funcs(p.numinputs))
function count_circuits( p::Parameters; nfuncs::Int64=0 )
  @assert p.numoutputs == 1   # Not tested for more than 1 output, but probably works in this case.
  nfuncs = nfuncs==0 ? length(default_funcs(p.numinputs)) : nfuncs
  multiplier = Int128(1)
  mij = 0
  for i = 1:p.numinteriors
    mf = nfuncs
    multiplier *= mf
    for j = 1:p.nodearity
      mij = min(p.numlevelsback,i-1+p.numinputs)
      multiplier *= mij
    end
    #print("i: ",i,"  mf: ",mf,"  mij: ",mij,"  log multiplier: ",log10(multiplier))
    exp = trunc(log10(multiplier))
    fract = 10^(log10(multiplier)-exp)
    #println("  exp: ",exp,"  fract: ",fract)
    #@printf("  multiplier: %4.2f",fract)
    #@printf("e+%2i\n",exp)
    #=
    try
      @printf("  multiplier:  %8.2e\n",multiplier)
    catch
    end
    =#
  end
  multiplier
end

# Mutates chromosome c by inserting a gate or deleting a gate
function mutate_num_gates!( c::Chromosome; prob_insert::Float64=0.5 )
  if rand() <= prob_insert
    c = insert_gate!(c)
  else
    c = insert_gate!(c)
  end
  c
end

# Insert a gate into c so that output_values(c) is unchanged.
# The new gate is inserted at position new_gate_index.
# The following gates are shifted right by 1.
# The new gate replaces an input connection to the existing gate a position new_gate_index (before the shift)
# The output of the new gate is the input to the existing gate.
# One input to the new gate is the source of the input to the exisiting gate.
# The other is chosen randomly according to the levelsback constraint.
function insert_gate!( c::Chromosome )
  funcs = default_funcs(c.params.numinputs)
  p = c.params
  # Chose a random interior node as the "new gate".
  # Don't insert a gate to replace an output gate:
  new_gate_index = rand(1:(p.numinteriors-p.numoutputs+1))  # Position of the new gate
  #println("new_gate_index: ",new_gate_index)
  p = c.params = Parameters( p.numinputs, p.numoutputs, p.numinteriors+1, p.numlevelsback+1 )
  new_gate = deepcopy(c.interiors[new_gate_index])  # will be modified below
  input_index = rand(1:2)
  new_gate.func = input_index == 1 ? IN1 : IN2  # IN1 and IN2 are gates that return input1 and input 2 respectively.
  c.interiors = vcat(c.interiors[1:(new_gate_index-1)],[new_gate],c.interiors[new_gate_index:end])
  for i = 1:p.numoutputs
    index = p.numinputs + p.numinteriors + i - p.numoutputs # use the last numoutputs interiors
    c.outputs[i] = OutputNode(index)
  end
  for i = (new_gate_index+1):p.numinteriors
    for j = 1:p.nodearity
      if c.interiors[i].inputs[j] >= p.numinputs + new_gate_index
        c.interiors[i].inputs[j] += 1
      end
    end
  end
  # Modify one of the inputs of interiors[new_gate_index+1] to point to the new gate
  c.interiors[new_gate_index+1].inputs[input_index] = p.numinputs + new_gate_index
  set_active_to_false(c)
  #println("insert_gate!  new interiors: ",c.params.numinteriors)
  c
end

function delete_gate!( c::Chromosome, interior_to_delete::Int64 )
  #println("delete_gate! interior_to_delete: ",interior_to_delete,"  active: ",c.interiors[interior_to_delete].active)
  p = c.params
  @assert interior_to_delete <= p.numinteriors - p.numoutputs
  gate_to_delete = c.interiors[interior_to_delete]
  new_levelsback = (p.numlevelsback > p.numinputs) ? p.numlevelsback-1 : p.numlevelsback  # Don't set numlevelsback to less than p.numinputs.
  p = c.params = Parameters( p.numinputs, p.numoutputs, p.numinteriors-1, new_levelsback )
  c.interiors = vcat(c.interiors[1:(interior_to_delete-1)],c.interiors[(interior_to_delete+1):end])
  for i = 1:p.numoutputs
    index = p.numinputs + p.numinteriors + i - p.numoutputs # use the last numoutputs interiors
    c.outputs[i] = OutputNode(index)
  end
  for i = interior_to_delete:p.numinteriors
    for j = 1:p.nodearity
      if c.interiors[i].inputs[j] >= interior_to_delete + p.numinputs
        c.interiors[i].inputs[j] -= 1
      end
    end
  end
  #println("delete_gate!  new interiors: ",c.params.numinteriors)
  c
end

function delete_gate!( c::Chromosome )
  output_values(c)  # needed to set active gates to active
  if number_active_gates(c) == c.params.numinteriors
    dg = rand(1:(c.params.numinteriors-c.params.numoutputs))
    return delete_gate!( c, dg )
  end
  inactive_interior_list = filter!(i->!c.interiors[i].active, collect(1:c.params.numinteriors))
  dg = rand(inactive_interior_list)
  return delete_gate!( c, dg )
end


text:
#!/usr/bin/env julia

using Luxor, Test

using Random
Random.seed!(42)

function findanddrawsometangentialcircles(c1center, c1radius, c2center, c2radius)
    for requiredradius in 1:150
        res = circletangent2circles(requiredradius, c1center, c1radius, c2center, c2radius)
        setopacity(0.5)
        if first(res) == 1
            circle(res[2], requiredradius, :stroke)
        elseif first(res) == 2
            circle(res[2], requiredradius, :stroke)
            circle(res[3], requiredradius, :stroke)
        elseif first(res) == 0
            sethue("red")
            circle(c1center, 5, :fill)
            circle(c2center, 5, :fill)
        end
    end
end

function test_circletangents(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setline(0.5)
    tiles = Tiler(1200, 1400, 5, 5)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            circle1 = (O + Point(rand(-20:20), rand(-20:20)), rand(20:40))
            circle2 = (O + Point(rand(-20:20), rand(-20:20)), rand(20:40))
            sethue("black")
            setopacity(0.5)
            circle(circle1..., :fill)
            circle(circle2..., :fill)
            setopacity(1.0)
            circle(circle1..., :stroke)
            circle(circle2..., :stroke)
            sethue("purple")
            findanddrawsometangentialcircles(circle1[1], circle1[2], circle2[1], circle2[2])
        end
    end

    @test finish() == true
    println("...finished circle tangent test, saved in $(fname)")
end

fname = "circletangent-test.pdf"

test_circletangents(fname)

ext: jl


text:
__precompile__(true)

module SimpleImage

function load_png(filepath::String)::Array{UInt32, 2}
	imgptr = ccall((:loadPng, "simplepng"), Ptr{UInt32}, (Cstring,), filepath)

	w, h = unsafe_wrap(Array{UInt32, 1}, imgptr, (2,))

	img = unsafe_wrap(Array{UInt32, 2}, imgptr + 8, (h, w))

	finalizer(free_img, img)

	return img
end

function free_img(img::Union{Array{UInt32}, Array{UInt8}})
	ccall((:freeImg, "simplepng"), Cvoid, (Ptr{UInt8},), pointer(img) - 8)
end

function write_png(filepath::String, img::Array{UInt32, 2})
	h, w = size(img)
	ccall((:writePng, "simplepng"), Cvoid, (Cstring, Ptr{Int32}, Int32, Int32), filepath, pointer(img), w, h)
end

function load_png_bytes(filepath::String)::Array{UInt8, 3}
	imgptr = ccall((:loadPng, "simplepng"), Ptr{UInt8}, (Cstring,), filepath)
	w, h = unsafe_wrap(Array{UInt32, 1}, Ptr{UInt32}(imgptr), (2,))
	img = unsafe_wrap(Array{UInt8, 3}, imgptr + 8, (UInt32(4), w, h))

	finalizer(free_img, img)

	return img
end

function write_png(filepath::String, img::Array{UInt8, 3})
	c, w, h = size(img)
	ccall((:writePng, "simplepng"), Cvoid, (Cstring, Ptr{Int8}, Int32, Int32), filepath, pointer(img), w, h)
end

load_png_rgb(filepath::String)::Array{UInt8, 3} = load_png_bytes(filepath)[1:3, :, :]

function write_png_rgb(filepath::String, img::Array{UInt8, 3})
	c, w, h = size(img)
	a::Array{UInt8, 3} = Array{UInt8, 3}(undef, 4, w, h)
	a[4, :, :] .= 0xff
	a[1:3, :, :] .= img
	write_png(filepath, a)
end

end


text:
struct Indexes{T<:AbstractVector}
    B::Vector{SparseMatrixCSC}
    assignments::Matrix{UInt16}
    IVF::T
end

struct QuantizerData{T<:AbstractVector}
    n_codebooks::Int
    n_centers::Int
    n_dp::Int
    n_dims::Int
    n_dims_center::Int
    I::Indexes
    C::T
end

"""
    function GenerateQuantizerData(data::AbstractMatrix; n_codebooks::Int, n_centers::Int)

Initialises an empty codebook & indexes data structure with the right shape for the n_codebooks and n_centers given.
"""
function GenerateQuantizerData(data::AbstractMatrix, n_codebooks::Int, n_centers::Int)
    n_dim, n_dp = size(data)
    n_dims_center = n_dim ÷ n_codebooks
    C = zeros(n_dims_center*n_centers*n_codebooks)
    I = Indexes(Array{SparseMatrixCSC}(undef, n_dp), Matrix{UInt16}(undef, n_codebooks, n_dp),
                    [[Vector{Int}() for i in 1:n_centers] for j in 1:n_codebooks])
    return QuantizerData(n_codebooks, n_centers, n_dp, n_dim, n_dims_center, I, C)
end

"""
    function get_indexes(qd::QuantizerData, i::Int, j::Int)

Helper function for the initialization step.\n
Converts the dimension and cluster numbers to the flattened format of the codebook.
"""
function get_indexes(qd::QuantizerData, i::Int, j::Int)
    d_i1, d_i2 = (i-1)*qd.n_dims_center+1, i*qd.n_dims_center
    c_i = (i-1)*qd.n_centers*qd.n_dims_center + (j-1)*qd.n_dims_center + 1
    return d_i1, d_i2, c_i
end

"""
    function initialization!(data::AbstractMatrix, qd::QuantizerData, MT::MultiThreaded, codebook::AbstractArray)

Initialization step of the Product Quantization algorithm. Assigns codeowords to the codebook
as random selections from the data points if codebook input is '0'. \n
If an approximate codebook is given (from incremental or Euclidean estimation), it simply takes this data
and updates the empty `qd.C` with the input values.
"""
function initialization!(data::AbstractMatrix, qd::QuantizerData, MT::MultiThreaded, codebook::Int)
    Threads.@threads for i in 1:qd.n_codebooks
        random_indexes = shuffle!(collect(1:qd.n_dp))[1:qd.n_centers] # Pick random indexes
        for (j, idx) in enumerate(random_indexes)
            d_i1, d_i2, c_i = get_indexes(qd, i, j)
            qd.C[c_i:(c_i+qd.n_dims_center-1)] = deepcopy(data[d_i1:d_i2,idx])
        end
    end
end

function initialization!(data::AbstractMatrix, qd::QuantizerData, ST::SingleThreaded, codebook::Int)
    for i in 1:qd.n_codebooks
        random_indexes = shuffle!(collect(1:qd.n_dp))[1:qd.n_centers] # Pick random indexes
        for (j, idx) in enumerate(random_indexes)
            d_i1, d_i2, c_i = get_indexes(qd, i, j)
            qd.C[c_i:(c_i+qd.n_dims_center-1)] = deepcopy(data[d_i1:d_i2,idx])
        end
    end
end

function initialization!(data::AbstractMatrix, qd::QuantizerData, thread::Processing, codebook::AbstractArray)
    qd.C[:].=codebook
end


text:
using Test, Documenter, ArchGDAL
DocMeta.setdocmeta!(
    ArchGDAL,
    :DocTestSetup,
    :(using ArchGDAL, GDAL);
    recursive = true,
)
doctest(ArchGDAL)


text:
### A Nbdev script file ###
### Autogenerated file. Don't modify. ###

module Nbdev
#export
begin
include(normpath(joinpath(@__FILE__,"..","..", "src/Export.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Documenter.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Makedocs.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/CodeRunner.jl")))
include(normpath(joinpath(@__FILE__,"..","..", "src/Make.jl")))
end

#export
begin
notebook2script = Export.notebook2script
showdoc = Documenter.showdoc
build = Make.build
new = Make.new
export notebook2script, showdoc, build, new
end

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module TridiagTest
using Base.Test
debug = false

# basic tridiagonal operations
n = 5

srand(123)

d = 1 .+ rand(n)
dl = -rand(n-1)
du = -rand(n-1)
v = randn(n)
B = randn(n,2)

for elty in (Float32, Float64, Complex64, Complex128, Int)
    if elty == Int
        srand(61516384)
        d = rand(1:100, n)
        dl = -rand(0:10, n-1)
        du = -rand(0:10, n-1)
        v = rand(1:100, n)
        B = rand(1:100, n, 2)
    else
        d = convert(Vector{elty}, d)
        dl = convert(Vector{elty}, dl)
        du = convert(Vector{elty}, du)
        v = convert(Vector{elty}, v)
        B = convert(Matrix{elty}, B)
    end
    ε = eps(abs2(float(one(elty))))
    T = Tridiagonal(dl, d, du)
    Ts = SymTridiagonal(d, dl)
    @test_throws ArgumentError size(Ts,0)
    @test size(Ts,3) == 1
    @test size(T, 1) == n
    @test size(T) == (n, n)
    F = diagm(d)
    for i = 1:n-1
        F[i,i+1] = du[i]
        F[i+1,i] = dl[i]
    end
    @test Array(T) == F

    # elementary operations on tridiagonals
    @test conj(T) == Tridiagonal(conj(dl), conj(d), conj(du))
    @test transpose(T) == Tridiagonal(du, d, dl)
    @test ctranspose(T) == Tridiagonal(conj(du), conj(d), conj(dl))

    @test abs.(T) == Tridiagonal(abs.(dl),abs.(d),abs.(du))
    if elty <: Real
        @test round.(T) == Tridiagonal(round.(dl),round.(d),round.(du))
        @test isa(round.(T), Tridiagonal)
        @test trunc.(T) == Tridiagonal(trunc.(dl),trunc.(d),trunc.(du))
        @test isa(trunc.(T), Tridiagonal)
        @test floor.(T) == Tridiagonal(floor.(dl),floor.(d),floor.(du))
        @test isa(floor.(T), Tridiagonal)
        @test ceil.(T) == Tridiagonal(ceil.(dl),ceil.(d),ceil.(du))
        @test isa(ceil.(T), Tridiagonal)
    end
    @test real(T) == Tridiagonal(real(dl),real(d),real(du))
    @test imag(T) == Tridiagonal(imag(dl),imag(d),imag(du))
    @test abs.(Ts) == SymTridiagonal(abs.(d),abs.(dl))
    if elty <: Real
        @test round.(Ts) == SymTridiagonal(round.(d),round.(dl))
        @test isa(round.(Ts), SymTridiagonal)
        @test trunc.(Ts) == SymTridiagonal(trunc.(d),trunc.(dl))
        @test isa(trunc.(Ts), SymTridiagonal)
        @test floor.(Ts) == SymTridiagonal(floor.(d),floor.(dl))
        @test isa(floor.(Ts), SymTridiagonal)
        @test ceil.(Ts) == SymTridiagonal(ceil.(d),ceil.(dl))
        @test isa(ceil.(Ts), SymTridiagonal)
    end
    @test real(Ts) == SymTridiagonal(real(d),real(dl))
    @test imag(Ts) == SymTridiagonal(imag(d),imag(dl))

    # test interconversion of Tridiagonal and SymTridiagonal
    @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl)
    @test SymTridiagonal(d, dl) == Tridiagonal(dl, d, dl)
    @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du)
    @test SymTridiagonal(d, dl) + Tridiagonal(dl, d, du) == Tridiagonal(dl + dl, d+d, dl+du)
    @test convert(SymTridiagonal,Tridiagonal(Ts)) == Ts
    @test Array(convert(SymTridiagonal{Complex64},Tridiagonal(Ts))) == convert(Matrix{Complex64}, Ts)
    if elty == Int
        vv = rand(1:100, n)
        BB = rand(1:100, n, 2)
    else
        vv = convert(Vector{elty}, v)
        BB = convert(Matrix{elty}, B)
    end
    let Bs = BB, vs = vv
        for atype in ("Array", "SubArray")
            if atype == "Array"
                BB = Bs
                vv = vs
            else
                BB = view(Bs, 1:n, 1)
                vv = view(vs, 1:n)
            end
        end

        # tridiagonal linear algebra
        @test T*vv ≈ F*vv
        invFv = F\vv
        @test T\vv ≈ invFv
        # @test Base.solve(T,v) ≈ invFv
        # @test Base.solve(T, B) ≈ F\B
        Tlu = factorize(T)
        x = Tlu\vv
        @test x ≈ invFv
    end
    @test det(T) ≈ det(F)

    @test T*Base.LinAlg.UnitUpperTriangular(eye(n)) ≈ F*eye(n)
    @test T*Base.LinAlg.UnitLowerTriangular(eye(n)) ≈ F*eye(n)
    @test T*UpperTriangular(eye(n)) ≈ F*eye(n)
    @test T*LowerTriangular(eye(n)) ≈ F*eye(n)

    # symmetric tridiagonal
    if elty <: Real
        Ts = SymTridiagonal(d, dl)
        Fs = Array(Ts)
        Tldlt = factorize(Ts)
        @test_throws DimensionMismatch Tldlt\rand(elty,n+1)
        @test size(Tldlt) == size(Ts)
        if elty <: AbstractFloat
            @test typeof(convert(Base.LinAlg.LDLt{Float32},Tldlt)) == Base.LinAlg.LDLt{Float32,SymTridiagonal{elty}}
        end
        let vs = vv
            for atype in ("Array", "SubArray")
                if atype == "Array"
                    vv = vs
                else
                    vv = view(vs, 1:n)
                end
            end
            invFsv = Fs\vv
            x = Ts\vv
            @test x ≈ invFsv
            @test Array(AbstractArray(Tldlt)) ≈ Fs
        end

        # similar
        @test isa(similar(Ts), SymTridiagonal{elty})
        @test isa(similar(Ts, Int), SymTridiagonal{Int})
        @test isa(similar(Ts, Int, (3,2)), Matrix{Int})
    end

    # eigenvalues/eigenvectors of symmetric tridiagonal
    if elty === Float32 || elty === Float64
        DT, VT = @inferred eig(Ts)
        @inferred eig(Ts, 2:4)
        @inferred eig(Ts, 1.0, 2.0)
        D, Vecs = eig(Fs)
        @test DT ≈ D
        @test abs.(VT'Vecs) ≈ eye(elty, n)
        @test eigvecs(Ts) == eigvecs(Fs)
        #call to LAPACK.stein here
        Test.test_approx_eq_modphase(eigvecs(Ts,eigvals(Ts)),eigvecs(Fs))
    elseif elty != Int
        # check that undef is determined accurately even if type inference
        # bails out due to the number of try/catch blocks in this code.
        @test_throws UndefVarError Fs
    end

    # Test det(A::Matrix)
    # In the long run, these tests should step through Strang's
    #  axiomatic definition of determinants.
    # If all axioms are satisfied and all the composition rules work,
    #  all determinants will be correct except for floating point errors.

    # The determinant of the identity matrix should always be 1.
    for i = 1:10
        A = eye(elty, i)
        @test det(A) ≈ one(elty)
    end

    # The determinant of a Householder reflection matrix should always be -1.
    for i = 1:10
        A = eye(elty, 10)
        A[i, i] = -one(elty)
        @test det(A) ≈ -one(elty)
    end

    # The determinant of a rotation matrix should always be 1.
    if elty != Int
        for theta = convert(Vector{elty}, pi ./ [1:4;])
            R = [cos(theta) -sin(theta);
                 sin(theta) cos(theta)]
            @test convert(elty, det(R)) ≈ one(elty)
        end

    # issue #1490
    @test det(ones(elty,3,3)) ≈ zero(elty) atol=3*eps(real(one(elty)))

    @test det(SymTridiagonal(elty[],elty[])) == one(elty)

    #tril/triu
    @test_throws ArgumentError tril!(SymTridiagonal(d,dl),n+1)
    @test_throws ArgumentError tril!(Tridiagonal(dl,d,du),n+1)
    @test tril(SymTridiagonal(d,dl))    == Tridiagonal(dl,d,zeros(dl))
    @test tril(SymTridiagonal(d,dl),1)  == Tridiagonal(dl,d,dl)
    @test tril(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,zeros(d),zeros(dl))
    @test tril(SymTridiagonal(d,dl),-2) == Tridiagonal(zeros(dl),zeros(d),zeros(dl))
    @test tril(Tridiagonal(dl,d,du))    == Tridiagonal(dl,d,zeros(du))
    @test tril(Tridiagonal(dl,d,du),1)  == Tridiagonal(dl,d,du)
    @test tril(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,zeros(d),zeros(du))
    @test tril(Tridiagonal(dl,d,du),-2) == Tridiagonal(zeros(dl),zeros(d),zeros(du))

    @test_throws ArgumentError triu!(SymTridiagonal(d,dl),n+1)
    @test_throws ArgumentError triu!(Tridiagonal(dl,d,du),n+1)
    @test triu(SymTridiagonal(d,dl))    == Tridiagonal(zeros(dl),d,dl)
    @test triu(SymTridiagonal(d,dl),-1) == Tridiagonal(dl,d,dl)
    @test triu(SymTridiagonal(d,dl),1)  == Tridiagonal(zeros(dl),zeros(d),dl)
    @test triu(SymTridiagonal(d,dl),2)  == Tridiagonal(zeros(dl),zeros(d),zeros(dl))
    @test triu(Tridiagonal(dl,d,du))    == Tridiagonal(zeros(dl),d,du)
    @test triu(Tridiagonal(dl,d,du),-1) == Tridiagonal(dl,d,du)
    @test triu(Tridiagonal(dl,d,du),1)  == Tridiagonal(zeros(dl),zeros(d),du)
    @test triu(Tridiagonal(dl,d,du),2)  == Tridiagonal(zeros(dl),zeros(d),zeros(du))

    @test !istril(SymTridiagonal(d,dl))
    @test !istriu(SymTridiagonal(d,dl))
    @test istriu(Tridiagonal(zeros(dl),d,du))
    @test istril(Tridiagonal(dl,d,zeros(du)))
    end
end

#Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences
function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing)
    n = size(a, 1)
    @test n==size(b,1) && size(a,2)==size(b,2)
    error===nothing && (error=n^3*(eps(S)+eps(T)))
    for i=1:n
        ev1, ev2 = a[:,i], b[:,i]
        deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2)))
        if !isnan(deviation)
            @test deviation ≈ 0.0 atol=error
        end
    end
end

let n = 12 #Size of matrix problem to test
    srand(123)
    debug && println("SymTridiagonal (symmetric tridiagonal) matrices")
    for relty in (Float32, Float64), elty in (relty, Complex{relty})
        debug && println("elty is $(elty), relty is $(relty)")
        a = convert(Vector{elty}, randn(n))
        b = convert(Vector{elty}, randn(n-1))
        if elty <: Complex
            a += im*convert(Vector{elty}, randn(n))
            b += im*convert(Vector{elty}, randn(n-1))
        end

        @test_throws DimensionMismatch SymTridiagonal(a, ones(n+1))
        @test_throws ArgumentError SymTridiagonal(rand(n,n))

        A = SymTridiagonal(a, b)
        fA = map(elty <: Complex ? Complex128 : Float64, Array(A))

        debug && println("getindex")
        @test_throws BoundsError A[n+1,1]
        @test_throws BoundsError A[1,n+1]
        @test A[1,n] == convert(elty,0.0)
        @test A[1,1] == a[1]

        debug && println("setindex!")
        @test_throws ArgumentError A[n,1] = 1
        @test_throws ArgumentError A[1,n] = 1
        A[3,3] = A[3,3]
        A[2,3] = A[2,3]
        A[3,2] = A[3,2]
        @test A == fA

        debug && println("Diagonal extraction")
        @test diag(A,1) == b
        @test diag(A,-1) == b
        @test diag(A,0) == a
        @test diag(A,n-1) == zeros(elty,1)
        @test_throws ArgumentError diag(A,n+1)

        debug && println("Idempotent tests")
        for func in (conj, transpose, ctranspose)
            @test func(func(A)) == A
        end

        debug && println("Simple unary functions")
        for func in (det, inv)
            @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(relty))
        end

        debug && println("Rounding to Ints")
        if elty <: Real
            @test round.(Int,A) == round.(Int,fA)
            @test isa(round.(Int,A), SymTridiagonal)
            @test trunc.(Int,A) == trunc.(Int,fA)
            @test isa(trunc.(Int,A), SymTridiagonal)
            @test ceil.(Int,A) == ceil.(Int,fA)
            @test isa(ceil.(Int,A), SymTridiagonal)
            @test floor.(Int,A) == floor.(Int,fA)
            @test isa(floor.(Int,A), SymTridiagonal)
        end

        debug && println("Tridiagonal/SymTridiagonal mixing ops")
        B = convert(Tridiagonal{elty},A)
        @test B == A
        @test B + A == A + B
        @test B - A == A - B

        debug && println("Multiplication with strided vector")
        @test A*ones(n) ≈ Array(A)*ones(n)

        debug && println("Multiplication with strided matrix")
        @test A*ones(n, 2) ≈ Array(A)*ones(n, 2)

        debug && println("Eigensystems")
        if elty <: Real
            zero, infinity = convert(elty, 0), convert(elty, Inf)
            debug && println("This tests eigenvalue and eigenvector computations using stebz! and stein!")
            w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, a, b)
            evecs = LAPACK.stein!(a, b, w)

            (e, v) = eig(SymTridiagonal(a, b))
            @test e ≈ w
            test_approx_eq_vecs(v, evecs)

            debug && println("stein! call using iblock and isplit")
            w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, a, b)
            evecs = LAPACK.stein!(a, b, w, iblock, isplit)
            test_approx_eq_vecs(v, evecs)

            debug && println("stegr! call with index range")
            F = eigfact(SymTridiagonal(a, b),1:2)
            fF = eigfact(Symmetric(Array(SymTridiagonal(a, b))),1:2)
            Test.test_approx_eq_modphase(F[:vectors], fF[:vectors])
            @test F[:values] ≈ fF[:values]

            debug && println("stegr! call with value range")
            F = eigfact(SymTridiagonal(a, b),0.0,1.0)
            fF = eigfact(Symmetric(Array(SymTridiagonal(a, b))),0.0,1.0)
            Test.test_approx_eq_modphase(F[:vectors], fF[:vectors])
            @test F[:values] ≈ fF[:values]
        end

        debug && println("Binary operations")
        a = convert(Vector{elty}, randn(n))
        b = convert(Vector{elty}, randn(n - 1))
        if elty <: Complex
            a += im*convert(Vector{elty}, randn(n))
            b += im*convert(Vector{elty}, randn(n - 1))
        end

        B = SymTridiagonal(a, b)
        fB = map(elty <: Complex ? Complex128 : Float64, Array(B))

        for op in (+, -, *)
            @test Array(op(A, B)) ≈ op(fA, fB)
        end
        α = rand(elty)
        @test Array(α*A) ≈ α*Array(A)
        @test Array(A*α) ≈ Array(A)*α
        @test Array(A/α) ≈ Array(A)/α

        debug && println("A_mul_B!")
        @test_throws DimensionMismatch A_mul_B!(zeros(elty,n,n),B,ones(elty,n+1,n))
        @test_throws DimensionMismatch A_mul_B!(zeros(elty,n+1,n),B,ones(elty,n,n))
        @test_throws DimensionMismatch A_mul_B!(zeros(elty,n,n+1),B,ones(elty,n,n))

    end

    debug && println("Tridiagonal matrices")
    for relty in (Float32, Float64), elty in (relty, Complex{relty})
        debug && println("relty is $(relty), elty is $(elty)")
        a = convert(Vector{elty}, randn(n - 1))
        b = convert(Vector{elty}, randn(n))
        c = convert(Vector{elty}, randn(n - 1))
        if elty <: Complex
            a += im*convert(Vector{elty}, randn(n - 1))
            b += im*convert(Vector{elty}, randn(n))
            c += im*convert(Vector{elty}, randn(n - 1))
        end

        @test_throws ArgumentError Tridiagonal(a,a,a)
        A = Tridiagonal(a, b, c)
        fA = map(elty <: Complex ? Complex128 : Float64, Array(A))

        debug && println("Similar, size, and copy!")
        B = similar(A)
        @test size(B) == size(A)
        copy!(B,A)
        @test B == A
        @test isa(similar(A), Tridiagonal{elty})
        @test isa(similar(A, Int), Tridiagonal{Int})
        @test isa(similar(A, Int, (3,2)), Matrix{Int})
        @test size(A,3) == 1
        @test_throws ArgumentError size(A,0)

        debug && println("Diagonal extraction")
        @test diag(A,-1) == a
        @test diag(A,0) == b
        @test diag(A,1) == c
        @test diag(A,n-1) == zeros(elty,1)
        @test_throws ArgumentError diag(A,n+1)

        debug && println("Simple unary functions")
        for func in (det, inv)
            @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(relty))
        end

        debug && println("Rounding to Ints")
        if elty <: Real
            @test round.(Int,A) == round.(Int,fA)
            @test isa(round.(Int,A), Tridiagonal)
            @test trunc.(Int,A) == trunc.(Int,fA)
            @test isa(trunc.(Int,A), Tridiagonal)
            @test ceil.(Int,A) == ceil.(Int,fA)
            @test isa(ceil.(Int,A), Tridiagonal)
            @test floor.(Int,A) == floor.(Int,fA)
            @test isa(floor.(Int,A), Tridiagonal)
        end

        debug && println("Binary operations")
        a = convert(Vector{elty}, randn(n - 1))
        b = convert(Vector{elty}, randn(n))
        c = convert(Vector{elty}, randn(n - 1))
        if elty <: Complex
            a += im*convert(Vector{elty}, randn(n - 1))
            b += im*convert(Vector{elty}, randn(n))
            c += im*convert(Vector{elty}, randn(n - 1))
        end

        debug && println("Multiplication with strided vector")
        @test A*ones(n) ≈ Array(A)*ones(n)

        debug && println("Multiplication with strided matrix")
        @test A*ones(n, 2) ≈ Array(A)*ones(n, 2)

        B = Tridiagonal(a, b, c)
        fB = map(elty <: Complex ? Complex128 : Float64, Array(B))

        for op in (+, -, *)
            @test Array(op(A, B)) ≈ op(fA, fB)
        end
        α = rand(elty)
        @test Array(α*A) ≈ α*Array(A)
        @test Array(A*α) ≈ Array(A)*α
        @test Array(A/α) ≈ Array(A)/α

        @test_throws ArgumentError convert(SymTridiagonal{elty},A)

        debug && println("A_mul_B!")
        @test_throws DimensionMismatch Base.LinAlg.A_mul_B!(zeros(fA),A,ones(elty,n,n+1))
        @test_throws DimensionMismatch Base.LinAlg.A_mul_B!(zeros(fA),A,ones(elty,n+1,n))

        debug && println("getindex")
        @test_throws BoundsError A[n+1,1]
        @test_throws BoundsError A[1,n+1]

        debug && println("setindex!")
        @test_throws ArgumentError A[n,1] = 1
        @test_throws ArgumentError A[1,n] = 1
        A[3,3] = A[3,3]
        A[2,3] = A[2,3]
        A[3,2] = A[3,2]
        @test A == fA
    end
end

# Issue 12068
SymTridiagonal([1, 2], [0])^3 == [1 0; 0 8]

#test convert for SymTridiagonal
@test convert(SymTridiagonal{Float64},SymTridiagonal(ones(Float32,5),ones(Float32,4))) == SymTridiagonal(ones(Float64,5),ones(Float64,4))
@test convert(AbstractMatrix{Float64},SymTridiagonal(ones(Float32,5),ones(Float32,4))) == SymTridiagonal(ones(Float64,5),ones(Float64,4))

# Test constructors from matrix
@test SymTridiagonal([1 2 3; 2 5 6; 0 6 9]) == [1 2 0; 2 5 6; 0 6 9]
@test Tridiagonal([1 2 3; 4 5 6; 7 8 9]) == [1 2 0; 4 5 6; 0 8 9]

# Test constructors with range and other abstract vectors
@test SymTridiagonal(1:3, 1:2) == [1 1 0; 1 2 2; 0 2 3]
@test Tridiagonal(4:5, 1:3, 1:2) == [1 1 0; 4 2 2; 0 5 3]
end

max_stars_repo_name: felixrehren/julia

max_issues_repo_name: felixrehren/julia

max_forks_repo_name: felixrehren/julia


text:
struct RenamingFunction{S} <: AbstractRenaming where S <: Callable
    s::S
end
(s::RenamingFunction)(x::Symbol) = Symbol(s.s(string(x)))
key_prefix(s::Union{AbstractString,Char}) = RenamingFunction(x -> string(s, x))
key_suffix(s::Union{AbstractString,Char}) = RenamingFunction(x -> string(x, s))
key_map(f::Callable) = RenamingFunction(f)
key_replace(p::Pair) = RenamingFunction(x -> replace(x, first(p) => last(p)))
Base.show(io::IO, r::RenamingFunction) = print(io, "RenamingFunction(", r.s, ")")

struct RenamingSymbols{S} <: AbstractRenaming
    s::S
end
Base.show(io::IO, r::RenamingSymbols) = print(io, "RenamingSymbols(", r.s, ")")
(r::RenamingSymbols)(x) = r.s

struct RenamingSymbol <: AbstractRenaming
    s::Symbol
end
Base.show(io::IO, r::RenamingSymbol) = print(io, "RenamingSymbol(:", r.s, ")")
(r::RenamingSymbol)(x) = r.s

apply_rename(::Nothing, colname) = colname
apply_rename(c::Composition, colname) = apply_rename(c.funcs, colname)
apply_rename(fs, colname) = apply_rename(Base.tail(fs), first(fs)(colname))
apply_rename(::Tuple{}, colname) = colname

renaming(r::R) where R <: AbstractRenaming = r
renaming(r::Pair{<:Union{AbstractString,Regex}, SubstitutionString}) = key_replace(r)
renaming(s::Symbol) = RenamingSymbol(s)
renaming(s::Union{Tuple{Vararg{Symbol}}, AbstractArray{Symbol}}) = RenamingSymbols(s)
renaming(c::Callable) = key_map(c)

"""
```
key_prefix(s::Union{AbstractString,Char})
key_suffix(s::Union{AbstractString,Char})
key_map(f::Callable)
key_replace(pat=>r)
```
Modify names of columns selected within a call to `select(tab, s...)`.
Usage:
* `Selection => key_preffix(s)` or `Selection => key_suffix(s)` -- Modify all of the selected column names by concatenating string `s` to them
* `Selection => key_map(f)` -- Modify all of the selected columns names by applying function `f` to them.
* `Selection => key_replace(pat=>r)` -- Modify all of the selected columns names by replacing pattern `pat` with `r`.

See also: [`select`](@ref), [`selection`](@ref), [`transformation`](@ref)
"""
key_prefix, key_suffix, key_map, key_replace, renaming


text:
using Plots, HDF5

@testset "HDF5_Plots" begin
    fname = "tmpplotsave.hdf5"
    hdf5()

    x = 1:10
    psrc = plot(x, x .* x) #Create some plot
    Plots.hdf5plot_write(psrc, fname)

    #Read back file:
    gr() #Choose some fast backend likely to work in test environment.
    pread = Plots.hdf5plot_read(fname)

    #Make sure data made it through:
    @test psrc.subplots[1].series_list[1][:x] == pread.subplots[1].series_list[1][:x]
    @test psrc.subplots[1].series_list[1][:y] == pread.subplots[1].series_list[1][:y]

    #display(pread) #Don't display. Regression env might not support
end #testset


text:
using HDF5, JLD, KDTrees, DataFrames, Base.Dates

function selectingTimes (start::Int64, finish::Int64, radius::Float64)
	cd("/Users/bzeng/Dropbox (MIT)/7 Coding/UROP/Data")
	start = 12
	finish = 14
	startingHour = start
	endingHour = finish
	radius = 140.0
	nodePositions = load("manhattan.jld","positions")
	numRides = zeros(Int64, length(nodePositions), length(nodePositions))
	avgTime = zeros(Float32, length(nodePositions), length(nodePositions))

	nodes = Array(Float64, (2,length(nodePositions)))
	for (i,p) in enumerate(nodePositions)
	    nodes[1,i] = p[1]
	    nodes[2,i] = p[2]
	end

	tree = KDTree(nodes)

	function coordinatesToNode(x::Float64, y::Float64)
		return inball(tree, [x, y], radius, false)
	end

	for i = 1:12
		println(i)
		df = 0
		@time df = readtable("reduced_trip_data_$(i).csv");
		for j = 1:nrow(df)
			if j % 1000 == 0
				print("\r $(j * 100 / nrow(df))                   ")
			end
			pickup = df[j, :pickup_datetime]
			dropoff = df[j, :dropoff_datetime]
			pickupTime = DateTime(pickup, "y-m-d H:M:S")
			dropoffTime = DateTime(dropoff, "y-m-d H:M:S")
			pickupFlag = 1 <= dayofweek(pickupTime) <= 5 && startingHour <= hour(pickupTime) <= endingHour
			dropoffFlag = 1 <= dayofweek(dropoffTime) <= 5 && startingHour <= hour(dropoffTime) <= endingHour
			if pickupFlag && dropoffFlag
				nodeSrcList = coordinatesToNode(df[j, :pickup_x], df[j, :pickup_y])
				nodeDstList = coordinatesToNode(df[j, :dropoff_x], df[j, :dropoff_y])
				for nodeSrc in nodeSrcList, nodeDst in nodeDstList
					numRides[nodeSrc, nodeDst] += 1
					prevTime = avgTime[nodeSrc, nodeDst]
					newTime = (dropoffTime - pickupTime).value / 1000.0
					avgTime[nodeSrc, nodeDst] = (newTime + prevTime * (numRides[nodeSrc, nodeDst] - 1)) / numRides[nodeSrc, nodeDst]
				end
			end
		end
	end

	df = DataFrame(node1 = Int64[], node2 = Int64[], numberOfRides = Int64[], averageTime = Float32[])
	for src = 1:length(nodePositions), dst = 1:length(nodePositions)
		if numRides[src, dst] > 0
			push!(df, [src dst numRides[src, dst] avgTime[src, dst]])
		end
	end

	@time writetable("trip_data_$(startingHour)_$(endingHour).csv", df)
end

@time selectingTimes(12, 14, 140.0)

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
# Autogenerated wrapper script for PDFHighlights_jll for x86_64-apple-darwin
export PDFHighlightsWrapper

using Poppler_jll
JLLWrappers.@generate_wrapper_header("PDFHighlights")
JLLWrappers.@declare_library_product(PDFHighlightsWrapper, "@rpath/PDFHighlightsWrapper.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Poppler_jll)
    JLLWrappers.@init_library_product(
        PDFHighlightsWrapper,
        "lib/PDFHighlightsWrapper.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

if !isempty(ARGS)
    ARGS[1] == "0" && exit(0)
end

# Prevent this from being put into the Main namespace
let
M = Module()
@eval M begin
if !isdefined(Base, :uv_eventloop)
    Base.reinit_stdio()
end
Base.include(@__MODULE__, joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testhelpers", "FakePTYs.jl"))
import .FakePTYs: open_fake_pty

CTRL_C = '\x03'
UP_ARROW = "\e[A"
DOWN_ARROW = "\e[B"

precompile_script = """
2+2
print("")
@time 1+1
; pwd
? reinterpret
using Ra\t$CTRL_C
\\alpha\t$CTRL_C
\e[200~paste here ;)\e[201~"$CTRL_C
$UP_ARROW$DOWN_ARROW$CTRL_C
123\b\b\b$CTRL_C
\b\b$CTRL_C
f(x) = x03
f(1,2)
[][1]
cd("complet_path\t\t$CTRL_C
"""

julia_cmd() = (julia = joinpath(Sys.BINDIR, Base.julia_exename()); `$julia`)
have_repl =  haskey(Base.loaded_modules,
                    Base.PkgId(Base.UUID("3fa0cd96-eef1-5676-8a61-b3b8758bbffb"), "REPL"))
Pkg = get(Base.loaded_modules,
          Base.PkgId(Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"), "Pkg"),
          nothing)

if Pkg !== nothing
    precompile_script *= Pkg.precompile_script
end

push!(LOAD_PATH, Sys.STDLIB)
using Sockets
Sockets.__init__()
using Libdl
empty!(LOAD_PATH)

function generate_precompile_statements()
    start_time = time()

    # Precompile a package
    mktempdir() do prec_path
        push!(DEPOT_PATH, prec_path)
        push!(LOAD_PATH, prec_path)
        pkgname = "__PackagePrecompilationStatementModule"
        mkpath(joinpath(prec_path, pkgname, "src"))
        write(joinpath(prec_path, pkgname, "src", "$pkgname.jl"),
              """
              module $pkgname
              end
              """)
        @eval using __PackagePrecompilationStatementModule
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
    end

    print("Generating precompile statements...")
    sysimg = Base.unsafe_string(Base.JLOptions().image_file)
    mktemp() do precompile_file, _
        # Run a repl process and replay our script
        repl_output_buffer = IOBuffer()
        @static if Sys.iswindows()
            # Fake being cygwin
            pipename = """\\\\?\\pipe\\cygwin-$("0"^16)-pty10-abcdef"""
            server = listen(pipename)
            pty_slave = connect(pipename)
            @assert ccall(:jl_ispty, Cint, (Ptr{Cvoid},), pty_slave.handle) == 1
            pty_master = accept(server)
        else
            pty_slave, pty_master = open_fake_pty()
        end
        done = false
        blackhole = Sys.isunix() ? "/dev/null" : "nul"
        withenv("JULIA_HISTORY" => blackhole, "JULIA_PROJECT" => nothing,
                "TERM" => "") do
            if have_repl
                p = run(`$(julia_cmd()) -O0 --trace-compile=$precompile_file --sysimage $sysimg
                        --compile=all --startup-file=no --color=yes
                        -e 'import REPL; REPL.Terminals.is_precompiling[] = true'
                        -i`,
                        pty_slave, pty_slave, pty_slave; wait=false)
                readuntil(pty_master, "julia>", keep=true)
                t = @async begin
                    while true
                        sleep(0.5)
                        s = String(readavailable(pty_master))
                        write(repl_output_buffer, s)
                        if occursin("__PRECOMPILE_END__", s)
                            break
                        end
                    end
                end
                if have_repl
                    for l in split(precompile_script, '\n'; keepempty=false)
                        write(pty_master, l, '\n')
                    end
                end
                write(pty_master, "print(\"__PRECOMPILE\", \"_END__\")", '\n')
                wait(t)

                # TODO Figure out why exit() on Windows doesn't exit the process
                if Sys.iswindows()
                    print(pty_master, "ccall(:_exit, Cvoid, (Cint,), 0)\n")
                else
                    write(pty_master, "exit()\n")
                    readuntil(pty_master, "exit()\r\e[13C\r\n")
                    # @assert bytesavailable(master) == 0
                end
                wait(p)
            else
                # Is this even needed or is this already recorded just from starting this process?
                p = run(`$(julia_cmd()) -O0 --trace-compile=$precompile_file --sysimage $sysimg
                        --compile=all --startup-file=no
                        -e0`)
            end
        end
        close(pty_master)

        # Check what the REPL displayed
        # repl_output = String(take!(repl_output_buffer))
        # println(repl_output)

        # Extract the precompile statements from stderr
        statements = Set{String}()
        for statement in split(read(precompile_file, String), '\n')
            occursin("Main.", statement) && continue
            push!(statements, statement)
        end

        if have_repl
            # Seems like a reasonable number right now, adjust as needed
            # comment out if debugging script
            @assert length(statements) > 700
        end

        # Create a staging area where all the loaded packages are available
        PrecompileStagingArea = Module()
        for (_pkgid, _mod) in Base.loaded_modules
            if !(_pkgid.name in ("Main", "Core", "Base"))
                eval(PrecompileStagingArea, :($(Symbol(_mod)) = $_mod))
            end
        end

        # Execute the collected precompile statements
        include_time = @elapsed for statement in sort(collect(statements))
            # println(statement)
            # Work around #28808
            occursin("\"YYYY-mm-dd\\THH:MM:SS\"", statement) && continue
            statement == "precompile(Tuple{typeof(Base.show), Base.IOContext{Base.TTY}, Type{Vararg{Any, N} where N}})" && continue
            try
                Base.include_string(PrecompileStagingArea, statement)
            catch
                @error "Failed to precompile $statement"
                rethrow()
            end
        end
        print(" $(length(statements)) generated in ")
        tot_time = time() - start_time
        Base.time_print(tot_time * 10^9)
        print(" (overhead "); Base.time_print((tot_time - include_time) * 10^9); println(")")
    end

    return
end

generate_precompile_statements()

end # @eval
end # let

max_stars_repo_name: syntapy/julia

max_issues_repo_name: syntapy/julia

max_forks_repo_name: syntapy/julia


text:
using MatchCore, NiLang
export alloc, @auto_alloc, @auto_expand

"""
    alloc(f, args...)

allocate function output space (the first argument), where `args` only contains the last `N-1` arguments.
"""
function alloc end

macro auto_alloc(ex)
    esc(auto_alloc(ex))
end

function auto_alloc(ex)
    @smatch ex begin
        :($f($out, $(args...))) => begin
            Expr(:block, :($out ← $alloc($f, $(args...))), ex)
        end
        :($out = $f($(args...))) => begin
            if length(args) == 0
                error("number of arguments must be >= 1.")
            else
                Expr(:block, :($out ← $alloc($f, $(args...))), :($out += $f($(args...))))
            end
        end
        _ => error("can not allocate automatically for expression: `$ex`")
    end
end

for OPM in [:PlusEq, :MinusEq]
    for OP in [:+, :-, :*, :/, :^]
        @eval alloc(::$OPM{typeof($OP)}, x::T1, ::T2) where {T1<:Number,T2<:Number} = zero(promote_type(T1, T2))
    end
    for OP in [:sin, :cos, :tan, :asin, :atan, :acos, :sinh, :cosh, :tanh, :identity, :sqrt, :exp, :log]
        @eval alloc(::$OPM{typeof($OP)}, x::T) where T<:Number = zero(T)
    end
    for OP in [:abs, :abs2]
        @eval alloc(::$OPM{typeof($OP)}, x::T) where T<:Number = zero(real(T))
    end
    @eval alloc(::$OPM{typeof(sincos)}, x::T) where T<:Number = (zero(T), zero(T))
end

function auto_expand(ex)
    res = Expr[]
    auto_expand!(copy(ex), res)
    Expr(:block, res..., NiLangCore.dual_body(@__MODULE__, res[1:end-1])...)
end

function auto_expand!(ex, exprs, sym=nothing, addnew=true)
    @smatch ex begin
        :($f($(args...))) => begin
            for (i, arg) in enumerate(args)
                @smatch arg begin
                    :($_{$(_...)}($(_...))) => begin
                        auto_expand!(arg, exprs, nothing, false)
                    end
                    :($f2($(vs...))) => begin
                        sym2 = gensym()
                        auto_expand!(:(PlusEq($f2)($sym2, $(vs...))), exprs, sym2, true)
                        args[i] = sym2
                    end
                    _ => nothing
                end
            end
            if sym !== nothing
                push!(exprs, :($sym ← $alloc($f, $(args[2:end]...))))
            end
            if addnew
                push!(exprs, :($f($(args...))))
            end
        end
        :($a += $b) || :($a -= $b) || :($a *= $b) || :($a /= $b) || :($a ⊻= $b) => begin
            auto_expand!(NiLangCore.to_standard_format(ex), exprs, sym, addnew)
        end
        _ => error("Can only expand an expression like `f(args...)`, got $(ex)!")
    end
end

macro auto_expand(ex)
    esc(auto_expand(ex))
end


text:
module BenchOverhead

using BenchmarkTools
using Transducers

suite = BenchmarkGroup()

let n = 64
    # s = suite["n=$n"] = BenchmarkGroup()
    s = suite  # there is only one `n` for now
    s["default"] = @benchmarkable reduce(
        right,
        Map(_ -> nothing),
        1:$n;
        basesize = 1,
        # stoppable = nothing,
    )
    s["stoppable=true"] = @benchmarkable reduce(
        right,
        Map(_ -> nothing),
        1:$n;
        basesize = 1,
        stoppable = true,
    )
    s["stoppable=false"] = @benchmarkable reduce(
        right,
        Map(_ -> nothing),
        1:$n;
        basesize = 1,
        stoppable = false,
    )
end

end  # module
BenchOverhead.suite


text:
"""
A struct containing the results of a simulation, including the specification used.
"""
struct SpatQSimResult{P,C,S}
    pop_state::Vector{P}
    catch_record::Vector{C}
    spec::S

    function SpatQSimResult(pop_state::Vector{P},
                            catch_record::Vector{C},
                            spec::S) where {P<:PopState,C<:Catch,S<:SpatQSimSpec}
        new{P,C,S}(pop_state, catch_record, spec)
    end
end

realization(result::SpatQSimResult) = realization(result.spec)
function save(result::SpatQSimResult; csv = false, feather = true)
    files = file_paths(result.spec)
    save_pop_hdf5(result, files[1])
    if (csv)
        save_pop_csv(result, files[2])
        save_catch_csv(result, files[3])
    end
    if (feather)
        save_pop_feather(result, files[4])
        save_catch_feather(result, files[5])
    end
    files
end

function save_pop_hdf5(result::SpatQSimResult, file_name::String)
    domain_size = size(domain(result.spec))
    h5open(file_name, "w") do fid
        pop_save = create_dataset(fid, "popstate", datatype(Float64),
                                  dataspace(domain_size...,
                                            length(result.pop_state)))
        for yr in eachindex(result.pop_state)
            pop_save[:, :, yr] = result.pop_state[yr].P
        end
    end
    file_name
end

function save_pop_csv(result::SpatQSimResult, file_name::String)
    popstate = (pop = sum.(result.pop_state), )
    CSV.write(file_name, StructArray(popstate))
    file_name
end

function save_catch_csv(result::SpatQSimResult, file_name::String)
    CSV.write(file_name, StructArray(result.catch_record))
    file_name
end

function save_pop_feather(result::SpatQSimResult, file_name::String)
    popstate = (pop = sum.(result.pop_state), )
    Arrow.write(file_name, StructArray(popstate))
    file_name
end

function save_catch_feather(result::SpatQSimResult, file_name::String)
    Arrow.write(file_name, StructArray(result.catch_record))
    file_name
end


text:
using Base.Meta: isexpr

# ------------------------------------------------------------------

function getStructTypeAndFields(expr)
  if !isexpr(expr, :type)
    error("trying to create constructor for non-type expression")
  end

  #get type name
  args = expr.args
  if isa(args[2],Symbol) structType = args[2]
  elseif isexpr(args[2],:<:) structType = args[2].args[1]
  else error("Unable to extract type name!")
  end

  # loop through all the fields and build an array of types
  fields = []
  for fieldExpr in args[3].args
    if isexpr(fieldExpr, :(::))
      push!(fields, (fieldExpr.args[1], fieldExpr.args[2]))
    end
  end

  return structType, fields
end

function makebody(structType, args)
  body = :($structType())
  append!(body.args, args)
  body
end

# create a new bitstype value by converting a UInt8 pointer to a "ftype" pointer, then load from the pointer
function loadFunction(ftype, s)
  :(unsafe_load(convert(Ptr{$ftype}, pointer(buf,pos+$s))))
end

# ------------------------------------------------------------------

# this macro takes a type definition expression and creates a custom constructor
# to read in from an IO source that is packed binary data
macro packedStruct(expr)

  structType, fields = getStructTypeAndFields(expr)
  fieldTypes = [f[2] for f in fields]

  # create constructors
  defaultConstructorBody = makebody(structType, [:(zero($ft)) for ft in fieldTypes])
  readIOConstructorBody = makebody(structType, [:(read(io, $ft)) for ft in fieldTypes])

  # create the reinterpret function to construct from a position in a byte array
  fsizes = [eval(current_module(), :(sizeof($ft))) for ft in fieldTypes]
  spos = vcat(1, cumsum(fsizes)[1:end-1] + 1)
  reinterpretBody = makebody(structType, [loadFunction(ft,spos[i]) for (i,ft) in enumerate(fieldTypes)])

  # now return the full expression
  # note: the whole thing is escaped so that this block is executed in calling scope
  # note: "esc(expression)" executes in calling scope, whereas "expression" executes in macro scope
  esc(quote
    $expr
    $structType() = $defaultConstructorBody
    $structType(io::IO) = $readIOConstructorBody
    Base.zero(::Type{$structType}) = $structType()
    Base.read(io::IO, ::Type{$structType}) = $structType(io)
    Base.read(mmio::MMapIO, ::Type{$structType}) = $structType(mmio.io)
    Base.reinterpret(::Type{$structType}, buf::AbstractVector{UInt8}, pos::Integer) = $reinterpretBody
    CTechCommon.getPackedStructSize(::Type{$structType}) = $(sum(fsizes))
  end)
end

# ------------------------------------------------------------------

# # create print/show methods, assuming string(typ) is defined
# macro createIOMethods(typ::Symbol)
#   esc(quote
#     Base.print(io::IO, o::$typ) = print(io, string(o))
#     Base.show(io::IO, o::$typ) = print(io, string(o))
#   end)
# end

# # ------------------------------------------------------------------

# macro pretty(expr::Expr)

#   structType, fields = getStructTypeAndFields(expr)

#   # create a string() expression and add "$sep$fieldname: $fieldexpr" for each field
#   strExpr = :(string())
#   for (i,f) in enumerate(fields)
#     fn = f[1]
#     sep = i>1 ? ", " : ""
#     ss = "$sep$fn="
#     push!(strExpr.args, :($ss))
#     push!(strExpr.args, :(xxx.$fn))
#   end

#   # return the original type expression, along with the new methods to be added
#   quote
#     $(esc(expr))
#     Base.string(xxx::$(esc(structType))) = string($(string(structType)), "{", $strExpr, "}")
#     Base.print(io::IO, xxx::$(esc(structType))) = print(io, string(xxx))
#     Base.show(io::IO, xxx::$(esc(structType))) = print(io, string(xxx))
#   end
# end

# recurse
function replace_syms(expr::Expr, names, obj)
    for i in eachindex(expr.args)
        expr.args[i] = replace_syms(expr.args[i], names, obj)
    end
    expr
end

# replace?
function replace_syms(s::Symbol, names, obj)
    if s in names
        :($(obj).$(s))
    else
        s
    end
end

# identity
replace_syms(x, names, obj) = x

export @with

"""
Recursively replace all occurences of a field with dot notation:

```julia
type T
    a
    b
end

@with t::T a = b + 5
# becomes:
t.a = t.b + 5
```
"""
macro with(typeexpr::Expr, body::Expr)
    @assert typeexpr.head == :(::)
    @assert length(typeexpr.args) == 2
    obj, objtype = typeexpr.args
    names = fieldnames(eval(Main, objtype))
    esc(replace_syms(body, names, obj))
end


text:
using JSON
using ChainDiag # https://github.com/yomichi/ChainDiag.jl

function foo(filename, S,T,L,Jz,Jxy,G)
    param = Dict{String,Any}("S"=>S,
                             "T"=>T,
                             "L"=>L,
                             "Jz"=>Jz,
                             "Jxy"=>Jxy,
                             "Gamma"=>G
                            )
    solver = SpinChainSolver(S,L,Jz=Jz,Jxy=Jxy,h=0.0,Guni=Guni,Gstag=Gstag)
    obs = solve(solver,1.0/T,10)
    obs2 = Dict{String,Any}()
    obs2["Energy"] = obs["Energy"]
    obs2["Energy^2"] = obs["Energy^2"]
    obs2["Specific Heat"] = obs["Specific Heat"]
    obs2["Magnetization"] = obs["Order Parameter"]
    obs2["Susceptibility"] = obs["Susceptibility"]
    result = Dict("Parameter"=>param, "Result"=>obs2)
    open(filename,"w") do io
        JSON.print(io, result, 2)
    end
end

#                           S,   T, L,   Jz,  Jxy,   G
for (id,p) in enumerate(((0.5, 1.0, 6,  1.0,  1.0, 0.0),
                         (0.5, 1.0, 6,  0.0,  1.0, 0.0),
                         (0.5, 1.0, 6,  1.0,  0.0, 0.0),
                         (0.5, 1.0, 6, -1.0,  0.0, 0.0),
                         (0.5, 1.0, 6, -1.0, -1.0, 0.0),
                         (1.0, 1.0, 6,  1.0,  1.0, 0.0),
                         (0.5, 1.0, 3,  0.0,  1.0, 0.0),
                        ))
    foo("chain_$id.json", p...)
end


text:
@deprecate autodiff(mode, block) markdiff(block)
@deprecate autodiff(mode) markdiff
@deprecate opdiff faithful_grad
@deprecate statdiff faithful_grad


text:
using DrWatson
@quickactivate "OneVision"

using OneVision.Examples.Car2DExamples
using BenchmarkTools: @benchmark
import AbstractPlotting
AbstractPlotting.__init__()
AbstractPlotting.inline!(false)

Car2DExamples.live_demo()

# Car2DExamples.tracking_example(time_end=20, noise_level = 0.003, plot_result=true)
# Car2DExamples.formation_example(time_end=20,
#     noise = 0.001, sensor_noise = 0.005, plot_result=true)

# import OneVision.Car1DExample
# Car1DExample.run_example(time_end = 20.0, freq = 20.0, noise=0.0, plot_result=true)
# @benchmark Car2DExamples.formation_example(time_end=20,
#     noise = 0.001, sensor_noise = 0.003, plot_result=false)

# @profview run_example(1:20 * t_end, freq; noise=0.0, plot_result=false)
# run_example(1:20 * t_end, freq; noise=0.0, plot_result=false) modules=[OneVision] maxdepth=3

# @benchmark run_example(time_end=20, plot_result=false)


text:
import GLM, DataFrames

data = simulate_data_correlation([[0.1], [0.3]])

# Vectors
@test isa(report([3,5,1]).text, String)
@test isa(report([3,5,1, missing]).text, String)
@test isa(report([3,5,1], median=true, dispersion=false, range=false).text, String)
@test isa(report(["a", "b", "c", "e", missing], n_strings=2).text, String)

# DataFrame
results = report(data)
@test isa(results.text, String)

# LM
model = GLM.lm(GLM.@formula(y ~ Var1 * Group), data)

results = report(model)
@test length(results.values) == 35
@test length(table(results)) == 8

results = report(model, CI=90, std_coefs=false, effect_size=nothing)
@test length(results.values) == 33

results = report(model, CI=90, std_coefs=false, effect_size="sawilowsky2009")
@test length(results.values) == 35

# GLM
data = simulate_data_logistic([[0.1], [0.3]])
model = GLM.glm(GLM.@formula(y ~ Var1 * Group), data, GLM.Binomial())

results = report(model)
@test length(results.values) == 34
@test length(table(results)) == 8

results = report(model, CI=90, std_coefs=false, effect_size=nothing)
@test length(results.values) == 32


text:
using Base.Test, StatsBase

addprocs(2)
using Distributions, ParallelGLM

n = 100_000
srand(1234321)

const Xt = convert(SharedArray,vcat(ones(n)',randn(19,n)))
const βtrue = randn(size(Xt,1))
const ηtrue = Xt'*βtrue
const μtrue = similar(ηtrue)
const ll = LogitLink()
for i in 1:n
    μtrue[i] = invlink(ll,ηtrue[i])
end
const y = convert(SharedArray,similar(μtrue))
for i in 1:n
    y[i] = rand() > μtrue[i] ? 0. : 1.
end

g = PGLM(Xt,y,Bernoulli());
fit(g; verbose=true);


text:
using DelimitedFiles
using LOWDER

include("generators.jl")

# -----------------------------------------------
# Path to file
# -----------------------------------------------
source_filename = "CUTEr_selected_problems.dat"

# -----------------------------------------------
# LOWDER parameters
# -----------------------------------------------
δinit = 2.0
Δinit = 2.0

problems = readdlm( source_filename, Int64 )

# -----------------------------------------------
# Problem info
# -----------------------------------------------
i = 28
nprob = problems[i, 1]
n = problems[i, 2]
p = problems[i, 3]
rsp = convert( Bool, problems[i, 4] )
n_points = n + 1

# Generates the problem
( x, l, u, fmin ) = problem_generator_mw( nprob, n, p, rsp; unconstrained = true )

# Solves the problem using 'lowder'
sol = LOWDER.lowder( fmin, x, l, u; δ = δinit, Δ = Δinit, m = n_points, verbose = 3, filename="../data_files/simple_runtest_01.dat")