<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Best Practices for Scientific Computing
Based on Greg Wilson et al.'s paper from 2014.
demonstration written by Yubo "Paul" Yang on Nov. 28 2016
last edited
Step3:  1. Write programs for people, not computers
example
Step4:  counter-example
Step5:  correct code $\neq$ maintainable code !
5. Plan for mistakes
checkpoints are executable documentations!
example
Step8:  counter-example
Step9:  python
while (istep != nstep)
<ASSISTANT_TASK:>
Python Code:
# what does this function do?
def ct(t):
    test = 0.0
    newt = 0.0
    
    stuff = 0.0
    junk = 100
    g = t
    while (junk>0):
        if junk == 100:
            for i in t:
                test += i
            test /= len(t)
            for i in range(len(t)):
                newt += (t[i]-test)*(t[i]-test)
            newt /= (len(t)-1)
            newt = np.sqrt(newt)
        
        for k in range(1,len(t)):
            junk = 0
            for ik in range(len(g)-k):
                junk += (t[ik]-test)*(t[ik+k]-test)
            if junk>0:
                stuff += junk/((len(g)-k)*newt**2.)
                
    return 1+2*stuff
Explanation: Best Practices for Scientific Computing
Based on Greg Wilson et al.'s paper from 2014.
demonstration written by Yubo "Paul" Yang on Nov. 28 2016
last edited: Nov. 30 2016
End of explanation
import numpy as np
def auto_corr_func(trace,k):
     calculates the correlation of a trace of scalars and 
    a shifted copy of itself by k value 
    
    mu  = np.mean(trace)       # mean
    sig = np.std(trace,ddof=1) # standard deviation 
    
    nval = len(trace)-k   # number of overlap values
    
    # find correlation between current and shifted traces
    ac = 0.0 # accumulate correlation
    for i in range(len(trace)-k):
        ac += (trace[i]-mu)*(trace[i+k]-mu)
    # end for i
    ac /= (nval*sig**2.)
    
    return ac
# end def auto_corr_func
def auto_corr_time(trace):
     calculates the auto-correlation time of a trace of scalars 
    
    corr_term = 0.0
    #  shift trace by k, and accumulate correlation until it vanishes
    for k in range(1,len(trace)):
        ac = auto_corr_func(trace,k)
        if ac>0:
            corr_term += ac
        else:
            break
        # end if
    # end for k
    
    corr_time = 1.0 + 2*corr_term
    return corr_time
# end def auto_corr_time
data = np.array([-9.  , -8.82, -8.66, -8.55, -8.44, -8.3 , -8.31, -8.18, -8.05,
       -7.93, -8.01, -8.  , -7.87, -7.88, -7.87, -7.95, -7.91, -7.86,
       -7.87, -7.82, -7.87, -7.83, -7.73, -7.8 , -7.8 , -7.78, -7.75,
       -7.75, -7.77, -7.84, -7.87, -7.84, -7.84, -7.78, -7.79, -7.75,
       -7.77, -7.73, -7.76, -7.76, -7.76, -7.85, -7.83, -7.79, -7.77,
       -7.8 , -7.72, -7.61, -7.63, -7.6 ])
auto_corr_time(data)
Explanation: 1. Write programs for people, not computers
example: meaningful variables, self-contained designs, and documentation
End of explanation
# see function defined at the beginnig
ct(data)
Explanation: counter-example: arbitrary variables, inconsistent naming, and complex design
End of explanation
# parse input file
import os
def parse_input_with_exception(some_input_file):
    type_dict = {
        'src_file_loc':str,
        'tar_file_loc':str,
        'time_step':float,
        'nstep':int,
        'natom':int
    }
    
    input_dict = {}
    for line in some_input_file.strip('\n').split('\n'):
        if line == '':
            continue
        # end if
        
        tokens = line.split()
        if len(tokens) != 3:
            raise NotImplementedError('expect "x = y" format, found %s'%line)
        # end if
        
        name,eq,val_str = tokens
        expected_type = type_dict[name]
        try:
            value = expected_type(val_str)
        except ValueError as error:
            raise ValueError( 'while reading "%s", '%name+str(error) )
        # end if
        
        if name in ['nstep','natm','time_step']:
            if value < 0:
                raise NotImplementedError('%s must be positive'%name)
            # end if
        # end if
        
        if name.endswith('loc'):
            found = os.path.isdir(value)
            if not found:
                raise NotImplementedError('%s %s not found'%(name,value))
            # end if
        # end if
        input_dict[name] = value
    # end for line
    return input_dict
# end def
Explanation: correct code $\neq$ maintainable code !
5. Plan for mistakes
checkpoints are executable documentations!
example: frequent checks, raise exceptions at the first sight of trouble
End of explanation
def just_parse(some_input_file):
    type_dict = {
        'src_file_loc':str,
        'tar_file_loc':str,
        'time_step':float,
        'nstep':int,
        'natom':int
    }
    
    input_dict = {}
    for line in some_input_file.strip('\n').split('\n'):
        name,eq,val_str = line.split()
        value = type_dict[name](val_str)
        input_dict[name] = value
    # end for line
    
    return input_dict
# end just_parse
some_input_file = 
src_file_loc = .
tar_file_loc = /tmp
time_step = 0.1
nstep     = 100
natom     = 36
parse_input_with_exception(some_input_file)
just_parse(some_input_file)
more_input_file = 
src_file_loc = .
tar_file_loc = /tmp
time_step = 0.1
nstep     = -99999999
natom     = 36
parse_input_with_exception(more_input_file)
just_parse(more_input_file)
# good luck terminating your while loop
Explanation: counter-example: relaaaxxx, it will probaly be fine (no it won't)
End of explanation
# exact diagonalization example
#  1. matrix construction vs. matrix diagonalization
#  2. locate matrix element vs. calculate matrix element
Explanation: python
while (istep != nstep):
    istep += 1
second layer of defense is automatic testing, turn bugs into test cases
use a symbolic debugger: interactive inspection
print statments: dump -> analyze -> infer
bash
python -m pdb myscript.py
pdb&gt;
bash
gdb ./myexec
gdb&gt; run input.txt
basic debugger commands
1. next: n
2. step: s 10
3. backtrace: bt
4. break
5. list
5. break if i==0
6. tell
7. return
2. Let the computer do the work
example: use for loops for parameter scan
bash
cwd=`pwd`
for folder in `ls | grep param_name`; do
  cd $folder
  ./run_my_script
  cd $cwd
done
python subprocess is your friend!
```python
import subprocess as sp
proc = sp.Popen('ls | grep param_name',stdout=sp.PIPE,stderr=sp.PIPE,shell=True)
out,err = proc.communicate()
folders = out.strip('\n').split('\n')
if len(folders)==0:
  raise IOError('no folders found')
end if
for folder in folders:
  pass # do something crazy!
end for
```
counter-example: Excel sheet
use a workflow tool: make, IPython/Jupyter
3. Make incremental changes
example: 
GitHub, BitBucket, Source Forge, and Google Code
counter-example:
DropBox, Email
4. Don't repeat yourself or others
DRY Principle
example: modularize code into functions for reuse, use high-quality libraries
counter-example: copy and paste code, and rewrite standard libraries
Libraries:
1. Numerical Integration
2. Linear Algebra Operations and Solvers
3. Standard Optimization Routines
dependency blow-up? -> docker
6. Make working code, then, and ONLY then optimize
It is hard to predict where a code would be the slowest.
End of explanation
<END_TASK>