text:
export components_matched, goals_matched_exact, goals_matched_hamming, next_chromosome!, mut_evolve, mut_evolve_increase_numints
export mut_reduce_numactive
export random_neutral_walk, match_score, findmaxall, findminall, findmaxall_1, findminall_1, findmaxrand, findminrand
export evolve_function, mut_evolve_repeat, circuit_evolve, run_circuit_evolve
export neutral_evolution, geno_circuits, geno_properties, geno_list_properties, lambda_evolution

# 5/21: To test:
# ni = 2; nc = 4
# julia> g = randgoallist(2,ni,nc)  # For this example g[1] is interpreted as output, g[2] is goal
# julia> pmin=match_score(g[1],g[2],ni)
# Not efficient if length(output) is greater than about 6
# Returns score of the best-scoring match between components of output and components of goal.
# A match between the components of outputs and the components of goal is permutation p of 1:nc
#   where ouput[i] is matched to goal[p[i]].
# There are two options depending on the keyword parameter avgfitness.
# match_score is the number of exact component matches plus the score of the best partial component match.
#
# There are two options depending on the setting of avgfitness.
# The score of of a partial match between component x and component y is 1.0 - hamming(x,y)/2^numinputs.
# Example:  let output = [0x0e, 0x0d, 0x00, 0x04], goal = [0x06, 0x04, 0x0a, 0x01]
# There is an exact match between ouput[4] and goal[2].
# One of the best partial matches is between output[1] and goal[1] where hamming(output[1],goal[1]) = 1.
# This gives a partial match score of 0.75.
# Thus, match_score(output,goal,ni) returns 1.75
function match_score( output::Vector{MyInt}, goal::Vector{MyInt},numinputs::Int64;
      avgfitness::Bool=false)
  #println("match_score: avgfitness: ",avgfitness)
  nc = length(output)  # number of components
  #println("output: ",output,"  goal: ",goal)
  @assert nc == length(goal)
  H = [ hamming_distance(output[i],goal[j],numinputs) for i = 1:nc, j=1:nc]
  #println("H: ",H)
  P = permutations(collect(1:nc))
  if avgfitness
    mxscores = [ (sum( 1.0-H[i,p[i]] for i = 1:nc ),[p[i] for i = 1:nc]) for p in P ]
    #println("mxscores: ",mxscores)
    findmaxall(mxscores)[1][1]
  else # number of exact matches plus Hamming score for best partial match
    mxscores = [ maxscore([H[i,p[i]] for i = 1:nc]) for p in P ]
    findmaxall(mxscores)[1]
  end
end
#=
function match_score( output::Vector{MyInt}, goal::Vector{MyInt},numinputs::Int64)
  nc = length(output)  # number of components
  println("match_score: output: ",output,"  goal: ",goal)
  @assert nc == length(goal)
  H = [ hamming_distance(output[i],goal[j],numinputs) for i = 1:nc, j=1:nc]
  println("H: ",H)
  P = permutations(collect(1:nc))
  mxscores = [ maxscore([H[i,p[i]] for i = 1:nc]) for p in P ]
  println("mxscores: ",mxscores)
  findmaxall(mxscores)[1]
end
=#

# return the number of zeros plus one minus the smallest non-zero value
# Thus, if V contains 2 zeros and the smallest non-zero value is 0.2, returns 2.8.
# If the components of V are Hamming distances between components x and y, then
#   maxscore is a rating of the quality of match between x and y
function maxscore( V::Vector{Float64} )
  minind = 1
  count_zeros = 0
  for i = 1:length(V)
    if V[i] > 0
      if (V[minind] <= 0 || V[i] < V[minind])
        minind = i
      end
    else
      count_zeros += 1
    end
  end
  result = count_zeros < length(V) ? (count_zeros + 1.0 - V[minind]) : Float64(count_zeros)
  @assert result <= length(V)
  result
end

# Returns a list of triples (comp,i,j) where comp==output[i]==goal[j] and for different triples
#    (comp1,i1,j1) and (comp2,i2,j2) we must have i1!=i2 and j1!=j2 (but we may have comp1==comp2).
# So the length of the result is the degree of matching between outputs and goal
# Based on a merge sort algorithm
# EXample:  let output=[0x7,0xd], let goal = [0x1,0x7]
# Result (0x07, 1, 2)  which says that output[1] == goal[2] == 0x07
function components_matched( output::Vector{MyInt}, goal::Vector{MyInt} )
  result = Tuple{MyInt,Int64,Int64}[]
  len = length(output)
  output_pairs = [(output[i],i) for i = 1:len]
  goal_pairs = [(goal[i],i) for i = 1:len]
  spoutput = sortperm(output_pairs)   # permutation that sorts the output components
  spgoal = sortperm(goal_pairs)         # permutation that sorts the goal components
  #println("spoutput: ",spoutput)
  #println("spgoal: ",spgoal)
  #println("soutput: ",output[spoutput])
  #println("sgoal: ",goal[spgoal])
  i = 1
  j = 1
  while i <= len && j <= len
    if output[spoutput[i]][1] == goal[spgoal[j]][1]
      push!(result,(output_pairs[spoutput[i]][1],output_pairs[spoutput[i]][2],goal_pairs[spgoal[j]][2]))
      i += 1
      j += 1
    elseif output[spoutput[i]][1] < goal[spgoal[j]][1]
      #println()
      i += 1
    elseif output[spoutput[i]][1] > goal[spgoal[j]][1]
      #println()
      j += 1
    end
  end
  result
end

# returns a triple (num_components_matched, goals_matched_list, components_matched_list)
#  where
#   num_components_matched is the maximum number of components of output which match a goal
#   goals_matched_list is the list of goals that achieve this maximum
#   components_matched_list is the list of components matched for each of the goals that achieves the maximum
#       (these are the outputs of the components_matched() function for each goal that achieves the maximum)
# Example:
#  output = [0x7,0xd]
#  goallist = [[0x1,0x2],[0x1,0x7],[0x7,0x3]]
#  goals_matched(output,goallist) = (1, [2, 3], Array{Tuple{UInt8,Int64,Int64},1}[[(0x07, 1, 2)], [(0x07, 1, 1)]])
#    which says that at most 1 component of output matches with a goal, and this is achieved for goal[2] and goal[3].
#    See the comments for components_matched() for the interpreation of the third component of the result triple.
# Note that fitness does not take robustness into account
function goals_matched_exact( output::Vector{MyInt}, goallist::GoalList, numinputs::Int64;
      avgfitness::Bool=false, perm_heuristic::Bool=false )   # Optional arguments not used in this function
  #println("exact")
  matched_goals = map( g->components_matched(output,g), goallist )
  #println("matched_goals: ",matched_goals)
  fm = findmaxall( map(length, matched_goals ))
  #println("fm: ",fm)
  (Float64(fm[1]), fm[2], matched_goals[fm[2]])
end

# Note that fitness does not take robustness into account
function goals_matched_hamming( output::Vector{MyInt}, goallist::GoalList, numinputs::Int64;
    avgfitness::Bool=false, perm_heuristic::Bool=false  )
  #println("goals_matched_hamming")
  if perm_heuristic
    goallist_scores = map( g->match_score_perm_heuristic(output,g,numinputs,avgfitness=avgfitness), goallist )
  else
    goallist_scores = map( g->match_score(output,g,numinputs,avgfitness=avgfitness), goallist )
  end
  matched_goals = map( g->components_matched(output,g), goallist )
  #println("goallist_scores: ",goallist_scores)
  #println("matched_goals: ",matched_goals)
  (best_score,best_score_list) = length(goallist_scores)>0 ? findmaxall( goallist_scores ) : (0.0,[])
  (best_score, best_score_list, matched_goals[best_score_list])
end

# Removed robust_sel and active_only keyword args on 6/6/20
# Fault tolerance is the  average Hamming deviation of the output from the unperturbed output under perturbation of the output of each node.
# See equation 3.3 of Macia and Sole (2009)
function next_chromosome!(c::Chromosome, goallist::GoalList, funcs::Vector{Func}, fitness::Float64=0.0;
      use_robustness::Bool=false, hamming_sel::Bool=true, active_only::Bool=false, num_mutations::Int64=1,
      avgfitness::Bool=false, perm_heuristic=perm_heuristic, fault_tol::Bool=false, ftf_param::Float64=0.95,
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  context = construct_context(c.params.numinputs)
  #println("next_chromosome!:  fault_tol: ",fault_tol,"  ftf_param: ",ftf_param)
  #print_build_chromosome( c )
  if fault_tol
    old_ftf = fault_tolerance_fitness( c )
    #println("old_ftf: ",old_ftf)
  end
  new_c = deepcopy(c)
  for _ = 1:num_mutations
    mutate_chromosome!( new_c, funcs, insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
    #print("next chromosome: ")
    #print_circuit(new_c)
  end
  output = execute_chromosome(new_c,context)
  #println("new_c output: ",output,"  new_c.fitness: ",new_c.fitness)
  goals_matched = hamming_sel ? goals_matched_hamming : goals_matched_exact
  ( new_fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs,
      avgfitness=avgfitness, perm_heuristic=perm_heuristic )
  #println("old_fitness: ",fitness,"  c.fitness: ",c.fitness, "  new_fitness: ",new_fitness)
  if !fault_tol
    if new_fitness < fitness
      new_c = c
    else
      new_c.fitness = new_fitness
    end
  else
    new_ftf = fault_tolerance_fitness( new_c )
    #println("new_ftf: ",new_ftf,"  new_c: ",(new_fitness >= fitness) && (new_ftf >= old_ftf))
    if (new_fitness < fitness) || ((new_ftf < old_ftf) && (rand() > ftf_param))
      new_c = c
    else
      new_c.fitness = new_fitness
    end
  end
  if use_robustness
    mut_robust = mutational_robustness( new_c, funcs, active_only=active_only )
    new_c.robustness = mut_robust
    #println("new_fitness: ",new_fitness,"  n_mut_robust: ",mut_robust)
  end
  #print_build_chromosome( new_c )
  return ( new_c, matched_goals, matched_goals_list )
end

# Does single-individual neutral evolution.
# Does subgoal evolution where the number of matched components of a goal is maximized if exact=true
# Does Hamming evolution where the Hamming distance of matched components of a goal is minimized if exact=false
# Removed robust_sel and active_only keyword args on 6/6/20
function mut_evolve( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, max_steps::Integer;
      hamming_sel::Bool=true, use_robustness::Bool=false, num_mutations::Int64=1, print_improvements::Bool=false,
      print_steps::Bool=true,
      avgfitness::Bool=false, perm_heuristic=false, fault_tol::Bool=false, ftf_param::Float64=0.95,
      fit_limit::Float64=Float64(c.params.numoutputs),
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  #println("mut_evolve fit limit: ",fit_limit)
  #println("mut evolve avgfitness: ",avgfitness,"  fault_tol: ",fault_tol,"  ftf_param: ",ftf_param)
  #print_build_chromosome(c)
  #print_parameters(c.params)
  #println("length(goallist): ",length(goallist))
  #println("goallist: ",goallist)
  c.fitness = 0.0  # Starting with a non-zero fitness inhibits evolution
  output = output_values(c)   # Executes c if it has not already been executed
  #println("output: ",output)
  #println("trunc(c.fitness): ",trunc(c.fitness))
  goals_matched = hamming_sel ? goals_matched_hamming : goals_matched_exact
  ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs, avgfitness=avgfitness )
  fitness = c.fitness
  orig_c = deepcopy(c)
  #println("initial fitness: ",fitness)
  robustness = 0.0
  step = 0
  worse = 0
  same = 0
  better = 0
  prev_c = deepcopy(c)  # previous generation c
  @assert output_values(prev_c) == output_values(orig_c)
  #(c, new_robustness, matched_goals_list ) =
  #    next_chromosome!(c, goallist, funcs, hamming_sel=hamming_sel, avgfitness=avgfitness )
  (c, matched_goals, matched_goals_list ) =
      next_chromosome!(c, goallist, funcs, hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations,
      perm_heuristic=perm_heuristic, avgfitness=avgfitness, fault_tol=fault_tol, ftf_param=ftf_param, insert_gate_prob=insert_gate_prob,
      delete_gate_prob=delete_gate_prob )
  #print_circuit(c)
  output = output_values(c)   # Executes c if it has not already been executed
  #println("output: ",output)
  #while step < max_steps && new_fitness < c.params.numoutputs
  #println("c.fitness: ",c.fitness,"  fit_limit: ",fit_limit)
  while step < max_steps && c.fitness < fit_limit
    #print("step: ",step,"  output: ",output," ")
    #print_circuit(c)
    #println("trunc(c.fitness): ",trunc(c.fitness))
    if c.fitness > fitness
      if print_improvements
        matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
        println("fitness improved from ",fitness," to ",c.fitness," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal componets: ", matched_goal_components)
      end
      fitness = c.fitness
      worse = 0
      same = 0
      better += 1
    elseif c.fitness == fitness
      #println("new chromosome with fitness: ",fitness )
      same += 1
    else
      #println("discarded chromosome with new_fitness: ",c.fitness)
      c = prev_c
      #print_circuit(c)
      worse += 1
    end
    if prev_c.fitness < orig_c.fitness
      println("step: ",step,"  prev_c.fitness: ",prev_c.fitness,"  orig_c.fitness: ",orig_c.fitness)
    end
    step += 1
    prev_c = deepcopy(c)
    #println("end loop ",rand(1:100))
    if step < max_steps   # Construct c for next generation
      (c, matched_goals, matched_goals_list ) =
          next_chromosome!(c, goallist, funcs, hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations,
            perm_heuristic=perm_heuristic, avgfitness=avgfitness, fault_tol=fault_tol, ftf_param=ftf_param,
            insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
            #next_chromosome!(c, goallist, funcs, hamming_sel=hamming_sel, use_robustness=use_robustness,
            #num_mutations=num_mutations, perm_heuristic=perm_heuristic, avgfitness=avgfitness, fault_tol=fault_tol,
            #ftf_param=ftf_param )
      output = output_values(c)   # Executes c if it has not already been executed
      #println("step: ",step+1,"  output: ",output)
      #print_circuit(c)
    end
  end
  if step == max_steps   # Failed to find goal
    if print_steps
      println("mut_evolve finished at step limit ",max_steps," with fitness: ", c.fitness )
    end
  else
    matched_g = map( x->[x[i][1] for i = 1:c.params.numoutputs], matched_goals_list )
    if print_steps
      println("mut_evolve finished in ",step," steps for goal ",matched_g[1]," with fitness: ", c.fitness )
    end
    #println("matched_goals: ",matched_goals,"  matched_goals_list: ",matched_goals_list)
  end
  if orig_c.fitness > c.fitness   # this should never happen
    #println("(orig_c.fitness,c.fitness): ",(orig_c.fitness,c.fitness))
    output = output_values(orig_c)
    ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
    return(orig_c,step,0,0,0,output,matched_goals,matched_goals_list)
  end
  ##println("worse: ",worse,"  same: ",same,"  better: ",better)
  #println("matched_goals: ",matched_goals)
  if step < max_steps
    if print_improvements
      matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
      println("fitness improved from ",fitness," to ",c.fitness," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal comps: ", matched_goal_components)
    end
    #println("sort(output: ",sort(output),"  sort(goallist[matched_goals[1]]): ",sort(goallist[matched_goals[1]]))
    if Int(trunc(fit_limit)) == c.params.numoutputs
      @assert sort(output) == sort(goallist[matched_goals[1]])
    end
  end
  (c,step,worse,same,better,output,matched_goals,matched_goals_list)
end

# Does up to n_repeats tries to evolve a chromosome which maps to a goal in goallist
function mut_evolve_repeat(n_repeats::Int64, p::Parameters, goallist::GoalList, funcs::Vector{Func}, max_steps::Integer;
      hamming_sel::Bool=true, use_robustness::Bool=false, num_mutations::Int64=1, print_improvements::Bool=false,
      avgfitness::Bool=false, perm_heuristic=false, fault_tol::Bool=false, ftf_param::Float64=0.95,
      fit_limit::Float64=Float64(p.numoutputs) )
  repeat = 1
  c = random_chromosome( p, funcs )
  (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, goallist, funcs, max_steps,
        hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations, print_improvements=print_improvements,
        avgfitness=avgfitness, perm_heuristic=perm_heuristic, fault_tol=fault_tol, ftf_param=ftf_param,
        fit_limit=fit_limit )
  while repeat < n_repeats && step == max_steps
    c = random_chromosome( p, funcs )
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, goallist, funcs, max_steps,
        hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations, print_improvements=print_improvements,
        avgfitness=avgfitness, perm_heuristic=perm_heuristic, fault_tol=fault_tol, ftf_param=ftf_param,
        fit_limit=fit_limit )
    repeat += 1
  end
  if step < max_steps  # circut that outputs a goal of goallist has been found
    return (c,step,worse,same,better,output,matched_goals,matched_goals_list)
  else
    return nothing
  end
end

# If mut_evolve() fails in max_steps iterations, increase the numinteriors and numlevelsback and try again.
# Continue up to repeat_limit until it suceeds.
function mut_evolve_increase_numints( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, max_steps::Integer, n_repeats::Int64=10;
      hamming_sel::Bool=true, use_robustness::Bool=false, num_mutations::Int64=1, print_improvements::Bool=false,
      avgfitness::Bool=false, perm_heuristic=false, fault_tol::Bool=false, ftf_param::Float64=0.95,
      fit_limit::Float64=Float64(c.params.numoutputs) )
  total_steps = 0
  new_numints = c.params.numinteriors
  new_levsback = c.params.numlevelsback
  repeat_limit = 10
  numinteriors_repeat = c.params.numinteriors
  numlevelsback_repeat = c.params.numlevelsback
  (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, goallist, funcs, max_steps,
      hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations, print_improvements=print_improvements,
      avgfitness=avgfitness, perm_heuristic=perm_heuristic, fault_tol=fault_tol, ftf_param=ftf_param,
      fit_limit=fit_limit )
  repeat = 0
  while repeat < repeat_limit && step == max_steps
    #println("starting repeat loop repeat: ",repeat)
    total_steps += max_steps
    numinteriors_repeat += 1
    numlevelsback_repeat += 1
    p_repeat = Parameters( numinputs=c.params.numinputs, numoutputs=c.params.numoutputs, numinteriors=numinteriors_repeat,
        numlevelsback=numlevelsback_repeat )
    println("repeating function mut_evolve_increse numints: ",numinteriors_repeat,"  and with levsback: ",numlevelsback_repeat,
        "  goallist[1]: ",goallist[1])
    c = random_chromosome( p_repeat, funcs )
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, goallist, funcs, max_steps,
        hamming_sel=hamming_sel, use_robustness=use_robustness, num_mutations=num_mutations, print_improvements=print_improvements,
        avgfitness=avgfitness, perm_heuristic=perm_heuristic, fault_tol=fault_tol, ftf_param=ftf_param,
        fit_limit=fit_limit )
    new_numints = p_repeat.numinteriors
    new_levsback = p_repeat.numlevelsback
    repeat += 1
  end
  total_steps += step
  if repeat == repeat_limit
    error(" repeat reached repeat_limit in function mut_evolve_increase_numints()")
  end
  (c,total_steps,worse,same,better,output,matched_goals,matched_goals_list,new_numints,new_levsback)
end

# Attempts to reduce the number active of the chromosome c by mutational neutral evolution
# The chomosome c is assumed to be at optimal fitness, and this function will terminate with
#   an error if it is not optimal.
# Then mutation is done nreps times, and the mutation is accepted if number_active(c) is reduced.
function mut_reduce_numactive( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, nreps::Integer;
      hamming_sel::Bool=true, num_mutations::Int64=1, print_improvements::Bool=false )
  output = output_values(c)   # Executes c if it has not already been executed
  goals_matched = hamming_sel ? goals_matched_hamming : goals_matched_exact
  ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
  if sort(output) != sort(goallist[matched_goals[1]])
    error("input chromosome must be have optimal fitness in function mut_reduce_numactive()")
  end
  best_num_active = number_active(c)
  #println("init num_active  num_active: ",best_num_active)
  orig_c = deepcopy(c)
  #println("initial fitness: ",fitness)
  step = 0
  worse = 0
  same = 0
  better = 0
  prev_c = deepcopy(c)  # previous generation c
  best_c = c
  @assert output_values(prev_c) == output_values(orig_c)
  (c, matched_goals, matched_goals_list ) = next_chromosome!(c, goallist, funcs,
      hamming_sel=hamming_sel, use_robustness=false, num_mutations=num_mutations )
  new_num_active = number_active(c)
  #println("before while new_num_active: ",new_num_active)
  #while step < nreps && new_fitness < c.params.numoutputs
  while step < nreps
    #println("step: ",step,"  output: ",output,"   ")
    if c.fitness > fitness  # This should never happen
      error("c.fitness > fitness")
      if print_improvements
        matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
        println("fitness improved from ",fitness," to ",c.fitness," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal components: ", matched_goal_components)
      end
      fitness = c.fitness
      worse = 0
      same = 0
      better += 1
    elseif c.fitness == fitness && new_num_active < best_num_active
      #println("new chromosome with fitness: ",fitness, "  and num_active: ", new_num_active )
      if print_improvements
        matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
        println("num_active decreased from ",best_num_active," to ",new_num_active," at step: ",step,"  goals_matched: ", matched_goals, "  matched goal componets: ", matched_goal_components)
      end
      worse = 0
      same = 0
      better += 1
      best_c = deepcopy(c)
      @assert new_num_active == number_active(c)
      best_num_active = new_num_active
      #println("num_active_improved  new_num_active: ",new_num_active)
    elseif c.fitness == fitness && new_num_active == best_num_active
      #println("new chromosome with fitness: ",fitness, "  and num_active: ", new_num_active )
      same += 1
    elseif c.fitness == fitness && new_num_active > best_num_active
      #println("discarded chromosome with fitness: ",fitness, " and num_active: ", new_num_active )
      c = prev_c
      worse += 1
    elseif c.fitness < fitness
      #println("discarded chromosome with new_fitness: ",c.fitness)
      c = prev_c
      worse += 1
    end
    step += 1
    prev_c = deepcopy(c)
    if step < nreps
      (c, matched_goals, matched_goals_list ) = next_chromosome!(c, goallist, funcs,
            hamming_sel=hamming_sel, use_robustness=false, num_mutations=num_mutations )
      #output = output_values(c)   # Executes c if it has not already been executed
      new_num_active = number_active(c)
      #println("end while new_num_active: ",new_num_active)
    end
  end
  if step == nreps
    #println("mut_reduce_numactive finished at step limit ",nreps," with fitness: ", best_c.fitness, " and num_active: ", best_num_active )
  else
    error("step == nreps")
    println("mut_reduce_numactive finished in ",step," steps with fitness: ", c.fitness, " and num_active: ", new_num_active )
  end
  if orig_c.fitness > c.fitness
    error("orig_c.fitness > c.fitness")
    #println("(orig_c.fitness,c.fitness): ",(orig_c.fitness,c.fitness))
    output = output_values(orig_c)
    ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
    return(orig_c,step,0,0,0,output,matched_goals,matched_goals_list)
  end
  ##println("worse: ",worse,"  same: ",same,"  better: ",better)
  #println("matched_goals: ",matched_goals)
  if step < nreps
    if print_improvements
      matched_goal_components = [map( x->x[1], mgl) for mgl in matched_goals_list]
      println("num_active decreased from ",num_active," to ",new_num_active," at step: ",step, "  goals_matched: ", matched_goals, "  matched goal componets: ", matched_goal_components)
    end
    @assert sort(output) == sort(goallist[matched_goals[1]])
  end
  ( fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, best_c.params.numinputs )
  #println("assert  best_num_active: ",best_num_active,"  number_active(best_c): ",number_active(best_c))
  @assert best_num_active == number_active(best_c)
  (best_c,step,worse,same,better,output,matched_goals,matched_goals_list)
end

# Do a random walk through circuit space starting with chromosome c which must be optimal for the given goallist.
# This function is trying to minimize num_active.
# There is another random_neutral_walk function in Evolability.jl which is trying to measure evolability
function random_neutral_walk( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, nreps::Integer;
      hamming_sel::Bool=true, num_mutations::Int64=1, print_improvements::Bool=false, reduce_numactive_reps::Int64=0 )
  sum_num_active = 0
  min_num_active = 2^20  # a very large integer
  count_saved_reps = 0
  count_min_num_active = 0
  min_active_chromes_list = Chromosome[]
  output = output_values(c)   # Executes c if it has not already been executed
  goals_matched = hamming_sel ? goals_matched_hamming : goals_matched_exact
  ( orig_fitness, matched_goals, matched_goals_list ) = goals_matched( output, goallist, c.params.numinputs )
  if sort(output) != sort(goallist[matched_goals[1]])
    error("input chromosome must be have optimal fitness in function mut_reduce_numactive()")
  end
  println("orig_fitness: ",orig_fitness)
  prev_c = deepcopy(c)
  for i = 1:nreps
    (c, matched_goals, matched_goals_list ) = next_chromosome!(c, goallist, funcs,
          hamming_sel=hamming_sel, use_robustness=false, num_mutations=num_mutations )
    if c.fitness < orig_fitness
      #println("discarding chromosome at rep: ",i)
      c = prev_c
      continue
    end
    new_num_active = number_active(c)
    count_saved_reps += 1
    #println("rep: ",i,"  num_active: ",new_num_active)
    sum_num_active += new_num_active
    if min_num_active > new_num_active
      println("new min: ",new_num_active,"  at rep: ",i)
      min_num_active = new_num_active
      count_min_num_active = 1
      @assert sort(output_values(c)) == sort(goallist[matched_goals[1]])
      min_active_chromes_list = Chromosome[c]
      if reduce_numactive_reps > 0
        (best_c,step,worse,same,better,output,matched_goals,matched_goals_list) =
            mut_reduce_numactive( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, nreps::Integer, num_mutations=num_mutations )
        best_num_active = number_active(best_c)
        if best_num_active < new_num_active
          println("best_num_active: ",best_num_active)
        end
      end
    elseif min_num_active == new_num_active
      #println("alt min: ",new_num_active,"  at rep: ",i)
      @assert sort(output_values(c)) == sort(goallist[matched_goals[1]])
      count_min_num_active += 1
      if length(min_active_chromes_list) < 4
        Base.push!(  min_active_chromes_list, c )
      end
      if reduce_numactive_reps > 0 && count_min_num_active < 10
        (best_c,step,worse,same,better,output,matched_goals,matched_goals_list) =
            mut_reduce_numactive( c::Chromosome, goallist::GoalList, funcs::Vector{Func}, nreps::Integer, num_mutations=num_mutations )
        best_num_active = number_active(best_c)
        if best_num_active < new_num_active
          println("best_num_active: ",best_num_active)
        end
      end
    end
    prev_c = deepcopy(c)
  end
  println("nreps: ",nreps,"  nmuts: ",num_mutations,"  avg num active: ",sum_num_active/count_saved_reps,
        "  min_num_active: ",min_num_active, "  count_min_num_active: ",count_min_num_active)
  return (min_num_active,count_min_num_active,min_active_chromes_list)
end

# Given a goal, for each iteraion of numiterations,
#    evolve two chromosomes that map to it, and try to find a distance nondecreasing path from one to the other
# TODO:  What happens when mut_evolve to find c1 or c2 fails?  Retry?  How many times?
function run_circuit_evolve( goal::Goal, p::Parameters, numiterations::Int64, maxreps::Int64 )
  funcs=default_funcs(p.numinputs)
  successes = 0
  sum_steps = 0
  for iter = 1:numiterations
    c = random_chromosome(p,funcs)
    res = mut_evolve(c,[goal],funcs,maxreps); c1=res[1]
    c = random_chromosome(p,funcs)
    res = mut_evolve(c,[goal],funcs,maxreps); c2=res[1]
    steps = circuit_evolve(c1,c2, maxreps )
    if steps <= maxreps
      println("success with ",steps," steps")
      sum_steps += steps
      successes += 1
    end
  end
  (successes,sum_steps/numiterations)
end

# Do neutral evolution starting at chromosome c_src trying to get to chromosome c_dest.
# c_src and c_dest must have the same parameters and the same output values (goals).
function circuit_evolve( c_src::Chromosome, c_dest::Chromosome, maxreps::Int64 )
  funcs = default_funcs(c_src.params.numinputs)
  @assert c_src.params == c_dest.params
  outputs_src = output_values( c_src )
  outputs_dest = output_values( c_dest )
  @assert outputs_dest == outputs_src
  c = deepcopy(c_src)  # current chromosome
  dist = circuit_distance( c_src, c_dest )
  #outputs_prev = output_values( c_src )
  step = 1
  c_new = deepcopy(c)
  mutate_chromosome!( c_new, funcs )
  #print("c_new before loop: ")
  #print_build_chromosome(c_new)
  while step <= maxreps
    if outputs_src != output_values(c_new)
      step += 1
      #println("continue")
      #continue  # continue with current circuit c unchanged
    else
      println("c_new code:  ",circuit_code(c_new))
      println("c_dest_code: ",circuit_code(c_dest))
      new_dist = circuit_distance( c_new, c_dest )
      println("dist:     ",dist,"  new_dist: ",new_dist)
      if new_dist > dist
        step += 1
        #continue  # continue with current circuit c unchanged
      elseif new_dist == 0  # success
        println("circuit_evolve() found a path to c_dest!")
        return step
      elseif new_dist < dist
        println("distance improved from ",dist," to ",new_dist)
        dist = new_dist
        c = deepcopy(c_new)
        step += 1
      end
    end
    c_new = deepcopy(c)
    mutate_chromosome!(c_new,funcs)
    #print("c_new end of loop: ")
    #print_build_chromosome(c_new)
  end
  println("circuit_evolve() failed to find a path to c_dest!")
  return step
end

# Change the fields of chromosome c to be the fields of chromosom c_to_copy
# TODO:  Move to Chromosome.jl
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

#=  Commented out on 4/17/21
# Example call:  build_chromosome((1,2), ((OR,[1,2]),(AND,[2,3])),(4,),0.0)
function build_chromosome( inputs::Tuple, ints::Tuple, outs::Tuple, fitness::Float64 )
  num_in = length(inputs)
  num_ints = length(ints)
  num_outs = length(outs)
  p = Parameters( numinputs=num_in, numoutputs=num_outs, numinteriors=num_ints, numlevelsback=num_ints+num_outs )
  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
=#

# The next functions should be in a Utilities.jl file
#  Returns a double whose first component is the maximum value in X, and whose second component is
#     the list of indices that give the maximum.
#=
function findmaxall( X::AbstractVector )
  max = findmax(X)[1]
  #println("max: ",max)
  Xind = [(X[i],i) for i = 1:length(X)]
  #println("Xind:",Xind)
  Xmax = filter(x->(x[1]==max),Xind)
  #println("Xmax:",Xmax)
  indices = [x[2] for x in Xmax]
  (max,indices)
end
=#

# argument A is a vector where > and == can be used to compare elements
function findmaxall( A::AbstractVector )
  if isempty(A)
    error("Empty list in findmaxall()")
  end
  indices = [1]
  for i = 2:length(A)
    if A[i] > A[indices[1]]
      indices = [i]
    elseif A[i] == A[indices[1]]
      push!(indices,i)
    end
  end
  (A[indices[1]],indices)
end

# argument A is a vector where < and == can be used to compare elements
function findminall( A::AbstractVector )
  if isempty(A)
    error("Empty list in findminall()")
  end
  indices = [1]
  for i = 2:length(A)
    if A[i] < A[indices[1]]
      indices = [i]
    elseif A[i] == A[indices[1]]
      push!(indices,i)
    end
  end
  (A[indices[1]],indices)
end

# cmp is a two-argument comparison function:  greater-than for max, less-than for min
# eq  is a two-argument equality function
function findmaxminall( cmp, eq, A::AbstractVector )
  if isempty(A)
    error("Empty list in findmaxminall()")
  end
  indices = [1]
  for i = 2:length(A)
    if cmp(A[i], A[indices[1]])
      indices = [i]
    elseif eq(A[i], A[indices[1]])
      push!(indices,i)
    end
  end
  (A[indices[1]],indices)
end

# findmaxall where the comparison functions look at the first components of the elements of A
findmaxall_1( A::AbstractVector ) = findmaxminall( (x,y)->(x[1]>y[1]), (x,y)->(x[1]==y[1]), A )

# findmnall where the comparison functions look at the first components of the elements of A
findminall_1( A::AbstractVector ) = findmaxminall( (x,y)->(x[1]<y[1]), (x,y)->(x[1]==y[1]), A )

#  Returns a double whose first component is the minimum value in X, and whose second component is
#     the list of indices that give the minimum.
function findminall( X::AbstractVector )
  min = findmin(X)[1]
  #println("min: ",min)
  Xind = [(X[i],i) for i = 1:length(X)]
  #println("Xind:",Xind)
  Xmin = filter(x->(x[1]==min),Xind)
  #println("Xmin:",Xmin)
  indices = [x[2] for x in Xmin]
  (min,indices)
end

# Returns the maximum value in A along with a random index for that value
function findmaxrand( A::AbstractVector )
  (val,indices) = findmaxall( A )
  if length(indices) == 1
    return (val,indices[1])
  else
    ind = rand(indices)
    #println("indices: ",indices,"  ind: ",ind)
    return (val,ind)
  end
end

# Returns the minimum value in A along with a random index for that value
function findminrand( A::AbstractVector )
  (val,indices) = findmaxall( A )
  if length(indices) == 1
    return (val,indices[1])
  else
    ind = rand(indices)
    #println("indices: ",indices,"  ind: ",ind)
    return (val,ind)
  end
end
# Evolves a LinCirucit that maps to g starting with chromosome c.
# max_steps is the maximum number of evolutionary steps.
# If evolution hasn't succeeeded in max_steps, return nothing.
# insert_gate_prob is the probability of inserting a gate on a mutation of a chromosome.
# delete_gate_prob is similar for deleting a gate.
# Similar to mut_evolve except that this takes a single goal instead of a goal list as an argument.
function neutral_evolution( c::Circuit, g::Goal, max_steps::Integer; print_steps::Bool=false,
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0,
      funcs::Vector{Func}=typeof(c) == LinCircuit ? lin_funcs( c.params.numinputs ) : default_funcs(c.params.numinputs) )
  LinCirc = typeof(c) == LinCircuit ? :true : :false
  #funcs = LinCirc ? lin_funcs( c.params.numinputs ) : default_funcs(c.params.numinputs)
  #println("LinCirc: ",LinCirc,"  Ones: ",Ones,"  CGP.Ones: ",CGP.Ones)
  #println("numgates: ",c.params.numinteriors)
  step = 0
  ov = output_values( c)
  current_distance = hamming_distance( ov, g, c.params.numinputs )
  while step < max_steps && ov != g
    step += 1
    #println("A step: ",step," current_distance: ",current_distance)
    if typeof(c) == Chromosome
      (new_c,active) = mutate_chromosome!( deepcopy(c), funcs, insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
    elseif typeof(c) == LinCircuit
      #println("c: ",c)
      new_c = mutate_circuit!( deepcopy(c), funcs )
    end
    new_ov = output_values( new_c )
    new_distance = hamming_distance( new_ov, g, c.params.numinputs )
    #println("step: ",step,"  ov: ",ov,"  new_ov: ",new_ov,"  cur dis: ",current_distance,"  new_dis: ",new_distance )
    #print("new_c: ")
    #print_circuit(new_c)
    #print("    c: ")
    #print_circuit(c)
    if new_ov == ov
      c = new_c
      if print_steps
        println("step: ",step," is pheno neutral.  new_ov: ",new_ov,"  new_distance: ",new_distance)
      end
    elseif new_distance == current_distance
      c = new_c
      if print_steps
        println("step: ",step," is fitness neutral.")
      end
    elseif new_distance < current_distance   # improvement
      if print_steps
        println("step: ",step,"  new_output: ",new_ov," distance improved from ",current_distance," to ",new_distance)
      end
      c = new_c
      ov = new_ov
      current_distance = new_distance
      #println("B step: ",step," current_distance: ",current_distance)
    else
      if print_steps
        #print("worse step: ",step,"  new_output: ",new_ov,"  new circuit: ")
        println("step: ",step,"  new_output: ",new_ov," current distance: ",current_distance," new: ",new_distance)
        #=
        if LinCirc
          print_circuit( new_c, funcs )
        else
          print_circuit( new_c )
        end
        =#
      end
    end
    #println("end while ",rand(1:100))
  end
  if step == max_steps
    println("neutral evolution failed with ",step," steps for goal: ",g)
    return (c, step)
  else
    println("neutral evolution succeeded at step ",step," for goal: ",g)
    @assert output_values(c) == g
    return (c, step)
  end
end

# Computes phenotype evolvability, genotype evolvability, robustness, complexity, steps
# num_circuits is the number of circuits used to compute properties.
function geno_circuits( g::Goal, p::Parameters, num_circuits::Integer, max_steps::Integer, max_attempts::Integer )
  funcs = default_funcs( p.numinputs )
  c = random_chromosome( p, funcs )
  sum_steps = 0
  circuit_list = Chromosome[]
  n_circuits = 0
  attempt = 0
  nc = nothing
  while attempt < max_attempts && n_circuits < num_circuits
    attempt += 1
    (nc,steps) = neutral_evolution( c, g, max_steps )
    sum_steps += steps
    if nc != nothing
      n_circuits += 1
      push!( circuit_list, nc )
    end
  end
  if n_circuits < num_circuits
    println("geno_properties failed to find num_circuits circuits mapping to goal: ",g," in ", attempt," attempts.")
    return (nothing, sum_steps)
  end
  return (circuit_list, sum_steps)
end

function geno_properties( cl_ss::Tuple{Vector{Chromosome}, Int64} ) # Tuple of circuit list and sum_steps.
  (cl, sum_steps) = cl_ss
  funcs = default_funcs( cl[1].params.numinputs )
  g = output_values(cl[1])
  sum_robust = 0.0
  sum_complexity = 0.0
  sum_evolvability = 0.0  # genotype evolability
  genotypes = Goal[]
  genotype_set = Set(Goal[])
  n_circuits = length(cl)
  for c in cl
    @assert output_values(c) == g
    sum_robust += mutational_robustness( c, funcs )
    sum_complexity += complexity5(c)
    genotypes = Set(mutate_all(c,funcs,output_outputs=true,output_chromosomes=false))
    sum_evolvability += length(genotypes)-1
    #println("len(genotypes): ",length(genotypes),"  genotypes: ",genotypes)
    genotype_set = union( genotype_set, genotypes )
  end
  (g, sum_robust/n_circuits, sum_complexity/n_circuits, sum_evolvability/n_circuits, length(genotype_set), sum_steps/n_circuits)
end

# Computes the properites robustness, complexity, genotype evolability, phenotype evolvability, evolutonary steps
#    for the goals of gl.  The properties are based on the evoution of num_circuits circuits that map to the goal.
# max_attempts is the maximum number of calls to neutral_evolution() that are made in the attempt to evolve num_circuits circuits
#    that evolve to a goal
# max_steps is the maximum number of evolutionary steps done by neutral_evoluion() during 1 attempt to evolve a circuit
#    that maps to a goal.
function geno_list_properties( gl::GoalList, p::Parameters, num_circuits::Integer, max_steps::Integer, max_attempts::Integer;
    csvfile::String="" )
  df = DataFrame()
  df.goal = Goal[]
  df.robustness = Float64[]
  df.complexity = Float64[]
  df.g_evolvability = Float64[]
  df.p_evolvability = Float64[]
  df.steps = Float64[]
  gp_list = Vector{Tuple{Goal,Float64,Float64,Float64,Float64,Float64}}[]
  gp_list = pmap( g->geno_properties( geno_circuits( g, p, num_circuits, max_steps, max_attempts ) ), gl )
  for gp in gp_list
    push!(df, gp)
  end
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# num_circuits: ",num_circuits)
      println(f,"# max_steps: ",max_steps)
      println(f,"# max_attempts: ",max_attempts)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# Does epochal evolution using a 1+lambda strategy where lambda is specified in c.params.
function lambda_evolution( c::Chromosome, g::Goal, maxsteps::Integer, mutrate::Float64 )
  p = c.params
  p.mutrate = mutrate
  funcs = default_funcs(p.numinputs)
  poisson_lambda = p.mutrate*p.numinteriors
  X = Poisson( poisson_lambda )
  step = 0
  ov = output_values( c)
  current_distance = hamming_distance( ov, g, c.params.numinputs )
  #print("ov: ",ov,"  cur_dist: ",current_distance,"  "); print_circuit(c)
  while step < maxsteps && ov != g
    chrome_list = Chromosome[]
    dist_list = Float64[]
    for i = 1:p.lambda
      step += 1
      num_mutations =  rand(X)
      new_c = deepcopy(c)
      for m = 1:num_mutations
        mutate_chromosome!( new_c, funcs )
      end
      new_ov = output_values( new_c)
      new_dist = hamming_distance( new_ov, g, c.params.numinputs )
      #print("new_ov: ",new_ov,"  new_dist: ",new_dist,"  "); print_circuit(new_c)
      push!(chrome_list,new_c)
      push!(dist_list,new_dist)
    end
    (best_dist,ind) = findmin( dist_list )
    if best_dist <= current_distance
      c = chrome_list[ind]
      ov = output_values( c )
      current_distance = hamming_distance( ov, g, c.params.numinputs )
      #print("ov: ",ov,"  cur_dist: ",current_distance,"  "); print_circuit(c)
    end
  end # while
  if step == maxsteps
    println("lambda evolution failed with ",step," steps for goal: ",g)
    return (c, step)
  else
    println("lambda evolution succeeded at step ",step," for goal: ",g)
    @assert output_values(c) == g
    return (c, step)
  end
end


text:
abstract type _Function end

mutable struct Variable
    data::Array
    grad::Array
    creator::_Function
    Variable(data) = new(data)
end
function backward(v::Variable)::Nothing
    if isdefined(v, :creator)
        f::_Function = v.creator
        x::Variable = f.input
        x.grad = backward(f, v.grad)
        backward(x)
    end
end

function _function(f::_Function, forward::Function, input::Variable)::Variable
    x = input.data
    y::Array = forward(x)
    output = Variable(y)
    output.creator = f
    f.input = input
    f.output = output
    return output
end

mutable struct Square <: _Function
    input::Variable
    output::Variable
    Square() = new()
end
(s::Square)(input::Variable)::Variable = _function(s, x -> x .^ 2, input)
backward(s::Square, gy::Array) = 2 * s.input.data .* gy

mutable struct Exp <: _Function
    input::Variable
    output::Variable
    Exp() = new()
end
(s::Exp)(input::Variable)::Variable = _function(s, x -> exp.(x), input)
backward(e::Exp, gy::Array) = exp.(e.input.data) .* gy

ext: jl


text:
### Day 21
## Scrambled Letters and Hash
## Author: Thanasis Georgiou
workspace()
using Combinatorics

# Parse integer shortcut
# Plus one because strings in arrays in Julia start at 1
int(str::SubString)::Int = parse(Int, str) + 1

# Swap index
function swapIndex!(str::Array{Char}, i::Int, j::Int)
  tmp = str[i]
  str[i] = str[j]
  str[j]  = tmp
end

# Swap letter
function swapCharacter!(str::Array{Char}, a::Char, b::Char)
  indexOfa = find(c -> c == a, str)
  indexOfb = find(c -> c == b, str)

  str[indexOfa] = b
  str[indexOfb] = a
end

# Rotate
function rotate!(str::Array{Char}, shift::Int)
  str[1:end] = circshift(str, shift)[1:end]
end

# Rotate based on index
function rotateByIndex!(str::Array{Char}, letter::Char)
  # Determine how much to move
  index = findfirst(c -> c == letter, str)
  if index > 4
    index += 1
  end

  # Rotate
  str[1:end] = circshift(str, index)[1:end]
end

# Reverse part of the string
function reversePart!(str::Array{Char}, from::Int, to::Int)
  str[from:to] = str[to:-1:from]
end

# Move letter
function moveIndex!(str::Array{Char}, from::Int, to::Int)
  char = str[from]
  deleteat!(str, from)
  insert!(str, to, char)
end

# Scramble a password
function scramble!(password::Array{Char}, instructions::Array{String})
  for line in instructions
    instr = split(chomp(line), ' ')

    # Find out which function to call
    # Swaps
    if instr[1] == "swap"
      if instr[2] == "position"
        int1 = int(instr[3])
        int2 = int(instr[6])
        swapIndex!(password, int1, int2)
      else # == "letter"
        c1 = instr[3][1]
        c2 = instr[6][1]
        swapCharacter!(password, c1, c2)
      end
      continue
    end

    # Rotates
    if instr[1] == "rotate"
      if length(instr) == 4
        steps = int(instr[3]) - 1
        if instr[2] == "left"
          steps *= -1
        end
        rotate!(password, steps)
      else
        letter = instr[7][1]
        rotateByIndex!(password, letter)
      end
      continue
    end

    # Reverse
    if instr[1] == "reverse"
      from = int(instr[3])
      to = int(instr[5])
      reversePart!(password, from, to)
      continue
    end

    # Move
    if instr[1] == "move"
      from = int(instr[3])
      to = int(instr[6])
      moveIndex!(password, from, to)
      continue
    end
  end
end

# Load input
f = open("input.txt")
input = readlines(f)
close(f)

# Scramble password
password = collect("abcdefgh")
scramble!(password, input)
println("Scrambled password: $(String(password))")

# Bruteforce problem 2
target = collect("fbgdceah")
for combination in permutations(password)
  scrambled = copy(combination)
  scramble!(scrambled, input)
  if scrambled == target
    println("Original password: $(String(combination))")
    break
  end
end


text:
export clip_segment, put_together

function bandpass_filter(gui,inds,wn1,wn2)

    responsetype = Bandpass(wn1,wn2; fs=30000)
    designmethod=Butterworth(4)
    df1=digitalfilter(responsetype,designmethod)
    myfilt=DF2TFilter(df1)
    test_data=gui.y_data[inds,1];
    filt!(test_data,myfilt,test_data);
    test_data
end

function thres_ts(data,thres)

    test_data=data./maximum(data)
    spike_times=falses(length(test_data))
    spike_inds=zeros(Int64,0)
    for i=1:length(test_data)

        if abs(test_data[i])>thres
            spike_times[i]=true
            i=i+30
        end
    end
    spike_times
end

#Convolve spike timestamps with fake spike waveform of 1 ms duration
function spike_audio(spike_times)
    conv(sin.(0:pi/15:2*pi),spike_times)
end

function sort_spikes(gui,t)

    gui.s.index=0
    SpikeSorting.onlinesort!(gui.s,gui.y_data[(t-30000)+1:t,1],gui.buf,gui.spike_nums)

    gui.spikes_ts=[]
    for i=1:gui.spike_nums[1]
        push!(gui.spikes_ts,t-30000+1+gui.buf[i].inds[1])
    end
    gui.spike_nums[1]=0

    nothing
end

function get_sorted_ts(gui,t)

    gui.spikes_ts=gui.pre_sorted_ts[(gui.pre_sorted_ts.>(t-30000+1))&(gui.pre_sorted_ts.<t)]

end
function event_triggered(gui,t)

    my_events=gui.event_ts[2][((gui.event_ts[2].>(t-30000)).&(gui.event_ts[2].<t))]

    myrange=-3000:300:3000

    for k=1:length(gui.myhist)
        gui.myhist[k]=0
    end

    for i=1:length(my_events)
        for j=1:length(gui.spikes_ts)
            for k=1:length(myrange)-1
               ttt = gui.spikes_ts[j]-my_events[i]
                if (ttt > myrange[k]).&(ttt<myrange[k+1])
                   gui.myhist[k] += 1
                end
            end
        end
    end

    nothing
end

#Get the video times and analog indexes that are closest
#to t1 and t2 while also being equal to a value the slider adopts
#and being a whole second value of the mp4 video
function get_clip_times(gui,t1,t2)

    slider_step = round(Int64,gui.max_time / ((gui.max_time-100) / 1000 * 10))
    total_slider_values = 30000:slider_step:gui.max_time

    v_frame_rate = 25

    #Video Frame that is equal to a second in save video nearest to desired time
    full_frame = round(Int,gui.video_ts[t1] / v_frame_rate)

    #Corresponding index of analog channel that corresponds to video
    first_ind=findfirst(gui.video_ts./v_frame_rate .== full_frame)

    #First value of the slider that corresponds to that frame
    first_slider_value = findfirst(total_slider_values.>first_ind)

    #Video time in seconds that is closest to t2
    end_frame = round(Int,gui.video_ts[t2] / v_frame_rate)

    end_ind = findfirst(gui.video_ts./v_frame_rate .== end_frame)

    last_slider_value = findfirst(total_slider_values.>end_ind)

    clip_times(t1,t2,full_frame,end_frame,first_ind,end_ind,first_slider_value,last_slider_value)
end

function clip_segment(gui,t1,t2)

    tt = get_clip_times(gui,t1,t2)

    num_frames = tt.ts2 - tt.ts1

    #Create video clip of desired length
    run(`ffmpeg -ss $(tt.ts1) -i $(gui.vid_path) -c copy -t $(num_frames) short.mp4`)

    #Cut out corresponding audio for that second

    my_voltage = gui.y_data[tt.ind1:(tt.ind2),1]

    gui.s.index=0
    #Sort so that only spikes are left to remove crumby audio
    SpikeSorting.onlinesort!(gui.s,my_voltage,gui.buf,gui.spike_nums)

    my_audio=zeros(Float32,length(my_voltage))

    for i=1:gui.spike_nums[1]
        my_audio[gui.buf[i].inds]=my_voltage[gui.buf[1].inds]
    end

    gui.spike_nums[1]=0

    #We want the duration to be the same as the length of the video.
    #so num_frames * 100
    WAV.wavwrite(my_audio,"test.wav",Fs=div(t2-t1,num_frames))

    output_name = string(t1,"_with_audio.mp4")

    #Combine audio and video
    run(`ffmpeg -i short.mp4 -i test.wav -c:v copy -c:a aac -strict experimental $(output_name)`)

    #The above video has audio and video synchronized, but the length of the video is still a
    #Function of the length of time (t2-t1), as well as how many video frames actually represent
    #one second. to be played at 1/10 of normal speed, we want a 25 fps video to represent
    # 250 camera frames
    #If the camera is higher or lower FPS, it should be adjusted accordingly

    total_time = (t2-t1)/30000

    #How do we determine frame rate of camera?
    #Probably should calculate it roughly when we load in video data
    total_camera_time = (num_frames * 25)

    #run(`ffmpeg -i $(output_name) -filter_complex "[0:v]setpts=0.5*PTS[v];[0:a]atempo=2.0[a]" -map "[v]" -map "[a]" out2.mp4`)

    nothing
end

function record_frames(gui,t1,t2)

    tt = get_clip_times(gui,t1,t2)

    io=GLVisualize.create_video_stream("analog.mp4",gui.win)

    total_frames = tt.s2 - tt.s1 - 1
    Reactive.set_value!(gui.p_slider.parents[1],tt.s1)
    Reactive.activate!(gui.p_slider.parents[1])

    last_time = gui.t
    GLVisualize.add_frame!(io)

    i=0
    while (i<total_frames)

        yield()

        if last_time != gui.t
            GLVisualize.add_frame!(io)
            i = i + 1
            last_time = gui.t
        end
    end

    close(io.io)

    sleep(5.0)
    #Because the video is saved at 35 fps, but the video is actually being
    #played at 30 fps, we need to change the video to be 30 fps
    run(`ffmpeg -y -i analog.mp4 -vf "setpts=0.833*PTS" -r 30 analog2.mp4`)

end

function put_together(gui,t1,t2)

    mkdir(string(t1))
    cd(string(t1))

    clip_segment(gui,t1,t2)
    resize!(gui.win,640,480)
    sleep(2.0)
    record_frames(gui,t1,t2)

    out_name = string(t1,"_with_audio.mp4")
    final_name = string(t1,".mp4")

    run(`ffmpeg -i analog2.mp4 -i $(out_name) -filter_complex hstack $(final_name)`)

    run(`rm short.mp4`)
    run(`rm $(out_name)`)
    run(`rm analog2.mp4`)
    run(`rm analog.mp4`)
    run(`rm test.wav`)
    cd("..")
    nothing
end


text:
@testset "histos" begin
    x = [0.2, 6.4, 3.0, 1.6]
    bins = [0.0, 1.0, 2.5, 4.0, 10.0]
    inds = ATools.digitize(x,bins)
    @test all(inds .==[1,4,3,2])

    bins = [0,1,7] # a small and a large bin
    obs = [0.5, 1.5, 1.5, 2.5]; # one observation in the small bin and three in the large
    h  = ATools.hist1d(obs, bins, false) # not normalized
    hn = ATools.hist1d(obs, bins, true) # normalized

    he = ATools.edges(h)
    xn = [h.weights[i] / (he[i+1]- he[i]) for i in 1:length(h.weights)]
    @test all(xn .== hn.weights)

    h= ATools.hist1d(qs, 100, 50.0, 150.0)

    @test all(ATools.edges(h) .== h.edges[1])
    he = ATools.edges(h)
    hc = ATools.centers(h)
    @test hc[1] == (he[1] + he[2]) / 2.0
    @test hc[end] == (he[end-1] + he[end]) / 2.0
    ic = argmax(h.weights)
    @info "position of the maximum" ic
    @info "value of the maximum" hc[ic]
    @test hc[ic] > 95.0 && hc[ic] < 105.0

    h2 = ATools.hist1d(qs, 25, 50.0, 150.0, true)
    h2n = ATools.hist1d(qs, 25, 50.0, 150.0, false)
    ic = argmax(h2.weights)
    r = h2n.weights[ic] / h2.weights[ic]  # area of bins (bin[i+1] - bin[1]) = 4
    #println(r)
    @test r == 4.0

    h2, _ = hist2d(xs, qs, 100, "xs", "qs")
    @test sum(h2.weights) == length(xs)
    @test length(h2.edges[1]) == size(h2.weights)[1] + 1
    max_q = argmax(sum(h2.weights, dims=2))[1]
    @test isapprox(centers(h2)[max_q], mean(qs), atol=2*step(h2.edges[1]))

    profile_nbins = 25
    p1df, _ = ATools.p1df(xs, qs, profile_nbins)
    @test (mean(p1df.y_mean) > 95.0 && mean(p1df.y_mean) < 105.0)
    @test (mean(p1df.x_mean) > 0.3 && mean(p1df.x_mean) < 0.7)
    @test (mean(p1df.y_std) > 8.0 && mean(p1df.y_std) < 12.0)

    p1df_filt, _ = ATools.p1df(xs, qs, profile_nbins, min_proportion=0.1)
    @test nrow(p1df_filt) == profile_nbins

end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
"""
    LU solver
"""
mutable struct LUSolver{T} <: LinearSolver
    lu::LU{T, <:AbstractMatrix{T}}
    A::Array{T, 2}
end

function lu_solver(A)
    LUSolver(
        lu!(zero(Array(A)), check = false),
        zero(Array(A)))
end

# function factorize!(solver::LUSolver{T}, A::SparseMatrixCSC{T,Int}) where T
#     solver.A .= A # sparse -> dense
#     solver.lu = lu!(solver.A, check = false)
# end

function factorize!(solver::LUSolver{T}, A::AbstractMatrix{T}) where T
    solver.A .= A # sparse -> dense
    solver.lu = lu!(solver.A, check = false)
end

function linear_solve!(solver::LUSolver{T}, x::Vector{T}, A::Matrix{T},
        b::Vector{T}; reg::T = 0.0, fact::Bool = true) where T
    #reg is useless but used for Mehrotra
    fact && factorize!(solver, A)
    ldiv!(x, solver.lu, b) # solve
end

function linear_solve!(solver::LUSolver{T}, x::Vector{T}, A::SparseMatrixCSC{T,Int},
        b::Vector{T}; reg::T = 0.0, fact::Bool = true) where T
    #reg is useless but used for Mehrotra
    fact && factorize!(solver, A)
    ldiv!(x, solver.lu, b) # solve
end

function linear_solve!(solver::LUSolver{T}, x::Matrix{T}, A::AbstractMatrix{T},
        B::AbstractMatrix{T}; reg::T = 0.0, fact::Bool = true) where T
    #reg is useless but used for Mehrotra
    fact && factorize!(solver, A)
    ldiv!(x, solver.lu, B) # solve
end


text:
# Test CategoricalVector with a hierarchical index (indexed using Tuples)
data = reshape(1.:40., 20, 2)
# v = collect(zip([:a, :b, :c][rand(1:3,20)], [:x,:y][rand(1:2,20)], [:x,:y][rand(1:2,20)]))
v = [(:b, :x, :y), (:c, :y, :y), (:b, :x, :y), (:a, :y, :y), (:b, :y, :y),
     (:c, :y, :y), (:b, :x, :x), (:c, :x, :y), (:c, :y, :y), (:a, :y, :y),
     (:a, :y, :y), (:b, :x, :y), (:c, :x, :y), (:c, :y, :y), (:b, :x, :y),
     (:a, :x, :x), (:c, :x, :x), (:c, :y, :y), (:b, :y, :x), (:b, :y, :y)]
idx = sortperm(v)
A = AxisArray(data[idx,:], AxisArrays.CategoricalVector(v[idx]), [:a, :b])
@test A[:b, :] == A[5:12, :]
@test A[[:a,:c], :] == A[[1:4;13:end], :]
@test A[(:a,:y), :] == A[2:4, :]
@test A[(:c,:y,:y), :] == A[16:end, :]
@test AxisArrays.axistrait(AxisArrays.axes(A)[1]) <: AxisArrays.Categorical

v = AxisArrays.CategoricalVector(collect([1; 8; 10:15]))
@test size(v) == (8,)
@test size(v, 1) == 8
@test size(v, 2) == 1
@test AxisArrays.axistrait(AxisArrays.axes(A)[1]) <: AxisArrays.Categorical
A = AxisArray(reshape(1:16, 8, 2), v, [:a, :b])

@test A[Axis{:row}(AxisArrays.CategoricalVector([15]))] == AxisArray(reshape(A.data[8, :], 1, 2), AxisArrays.CategoricalVector([15]), [:a, :b])
@test A[Axis{:row}(AxisArrays.CategoricalVector([15])), 1] == AxisArray([A.data[8, 1]], AxisArrays.CategoricalVector([15]))
@test A[atvalue(15), :] == AxisArray(A.data[8, :], [:a, :b])
@test A[atvalue(15), 1] == 8
@test AxisArrays.axistrait(AxisArrays.axes(A)[1]) <: AxisArrays.Categorical

# TODO: maybe make this work? Would require removing or modifying Base.getindex(A::AxisArray, idxs::Idx...)
# @test A[AxisArrays.CategoricalVector([15]), 1] == AxisArray([A.data[8, 1]], AxisArrays.CategoricalVector([15]))


text:
"""
    mix_trans(gas::A, P, T, X, mean_MW) where {A <: Arrhenius.Solution}

Compute the tranposrt properties of a mixture using mixture average formula

> Equations Ref. https://personal.ems.psu.edu/~radovic/ChemKin_Theory_PaSR.pdf
> Equations. 5-50/51/52 for viscosity and thermal conductivity

Pure species viscosities [Pa-s]

Thermal conductivity. [W/m/K].

> Equation 5-46 for diffusion

Mixture-averaged diffusion coefficients [m^2/s] relating the mass-averaged diffusive fluxes
(with respect to the mass averaged velocity) to gradients in the species mole fractions.

Test this module in _transport_test.jl

See also implementations in ReacTorch
"""
function mix_trans(gas::A, P, T, X, mean_MW) where {A <: Arrhenius.Solution}

    logT = log(T)
    trans_T = [logT^6, logT^5, logT^4, logT^3, logT^2, logT, 1.0]

    ## species_viscosities_poly

    η = trans_T' * gas.trans.species_viscosities_poly

    Wk_over_Wj = gas.MW * (1.0 ./ gas.MW')

    Wj_over_Wk = 1 ./ Wk_over_Wj

    ηk_over_ηj = η' * (1.0 ./ η)

    Φ = @. 1.0 / sqrt(8.0) / sqrt(1.0 + Wk_over_Wj) *
        (1.0 + sqrt(ηk_over_ηj) * (Wj_over_Wk)^(0.25))^2

    η_mix = sum((X .* η') ./ (Φ * X))

    ## thermal_conductivity_func
    λ = trans_T' * gas.trans.thermal_conductivity_poly

    λ_mix = (sum(X .* λ') + 1.0 / sum(X ./ λ')) / 2.0

    ## binary_diff_coeffs_func

    X_eps = clamp.(X, 1.e-12, Inf)
    D = reshape(trans_T' * gas.trans.binary_diff_coeffs_poly,
                gas.n_species, gas.n_species)

    XjWj = X_eps' * gas.MW .- X_eps .* gas.MW

    XjDjk = sum(X_eps .* (1.0 ./ D'), dims=1) .- (X_eps ./ diag(D))'

    Dkm = XjWj ./ XjDjk' ./ mean_MW / P * one_atm

    return η_mix, λ_mix, Dkm
end
export mix_trans


text:
module PlotlyJSTest
using Test

using PlotlyJS
const M = PlotlyJS

using Blink
!Blink.AtomShell.isinstalled() && Blink.AtomShell.install()

include("blink.jl")
include("kaleido.jl")

end


text:
n = 10 # rand(6:10)

A = matrixdepot("golub", n)
@test cond(A) > 1.e8

B = matrixdepot("golub", Int, n)
@test cond(A) > 1.e8

println("'golub' passed test...")


text:
"""
This file compares the performance by using dynamic and static arrays to store solutions.
"""

using Kinetic, StaticArrays, BenchmarkTools, OffsetArrays
cd(@__DIR__)

#--- dynamic ---#
ks, ctr, face, t = Kinetic.initialize("bench_ds.txt")

@btime Kinetic.solve!(ks, ctr, face, t)
"""The result on my NUC is around 707 ms."""

#--- static ---#
ctr2 = deepcopy(ctr)
for i in eachindex(ctr2)
    if i <= ks.pSpace.nx ÷ 2
        ctr2[i] = ControlVolume1D2F(
            ks.pSpace.x[i],
            ks.pSpace.dx[i],
            MVector{length(ks.ib.wL)}(ks.ib.wL),
            MVector{length(ks.ib.wL)}(ks.ib.primL),
            MVector{length(ks.ib.hL)}(ks.ib.hL),
            MVector{length(ks.ib.hL)}(ks.ib.bL),
        )
    else
        ctr2[i] = ControlVolume1D2F(
            ks.pSpace.x[i],
            ks.pSpace.dx[i],
            MVector{length(ks.ib.wL)}(ks.ib.wR),
            MVector{length(ks.ib.wL)}(ks.ib.primR),
            MVector{length(ks.ib.hL)}(ks.ib.hR),
            MVector{length(ks.ib.hL)}(ks.ib.bR),
        )
    end
end

face2 = deepcopy(face)
for i = 1:ks.pSpace.nx+1
    face2[i] = Interface1D2F(MVector{length(ks.ib.wL)}(ks.ib.wL), MVector{length(ks.ib.hL)}(ks.ib.hL))
end

@btime Kinetic.solve!(ks, ctr2, face2, t)
"""
The result on my NUC is around 523 ms.
Further improvements are possible if we make all structs (ks, ctr, face) static.
"""