| 1 | try: |
|---|
| 2 | from sans.models.prototypes.SimCylinderF import SimCylinderF |
|---|
| 3 | except: |
|---|
| 4 | print "This test uses the prototypes module." |
|---|
| 5 | from sans.models.CylinderModel import CylinderModel |
|---|
| 6 | import time, math, random |
|---|
| 7 | |
|---|
| 8 | def reset(seed=0): |
|---|
| 9 | sim = SimCylinderF() |
|---|
| 10 | sim.setParam('radius', 20) |
|---|
| 11 | sim.setParam('length', 200) |
|---|
| 12 | sim.setParam('phi', 1) |
|---|
| 13 | sim.setParam('theta', 1) |
|---|
| 14 | sim.setParam('seed', math.fmod(time.time()*100+seed,10000)) |
|---|
| 15 | return sim |
|---|
| 16 | |
|---|
| 17 | # number of times to repeat each point |
|---|
| 18 | N_q = 50 |
|---|
| 19 | # Relative error on the simulation |
|---|
| 20 | sig_val = 0.20 |
|---|
| 21 | |
|---|
| 22 | |
|---|
| 23 | counter = 0 |
|---|
| 24 | |
|---|
| 25 | sim = reset() |
|---|
| 26 | sph = CylinderModel() |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | sph.setParam('radius', 20) |
|---|
| 30 | sph.setParam('length', 200) |
|---|
| 31 | sph.setParam('scale', 1) |
|---|
| 32 | sph.setParam('contrast',1) |
|---|
| 33 | sph.setParam('length', 200) |
|---|
| 34 | sph.setParam('cyl_phi', 1) |
|---|
| 35 | sph.setParam('cyl_theta', 1) |
|---|
| 36 | |
|---|
| 37 | vol = math.pi*20*20*200 |
|---|
| 38 | |
|---|
| 39 | f = open('sim_err_cyl.txt', 'w') |
|---|
| 40 | f.write("<q> <npts> <sigma>\n") |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | |
|---|
| 45 | for i_npts in range(1): |
|---|
| 46 | |
|---|
| 47 | # Find which q has an error of sig_val |
|---|
| 48 | #print err |
|---|
| 49 | q_sigma = 0 |
|---|
| 50 | i_min = 0 |
|---|
| 51 | i_max = 0 |
|---|
| 52 | i_last_max = 0 |
|---|
| 53 | val_min = 0 |
|---|
| 54 | val_max = 0 |
|---|
| 55 | last_error = 0 |
|---|
| 56 | last_val = 0 |
|---|
| 57 | error_vec = [] |
|---|
| 58 | |
|---|
| 59 | found_it = False |
|---|
| 60 | i_vec = 0 |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | t_0 = time.time() |
|---|
| 64 | npts = 50000 |
|---|
| 65 | #npts = 50000*i_npts+1000 |
|---|
| 66 | #npts = 1000 * math.pow(2.0,(i_npts+1)) |
|---|
| 67 | |
|---|
| 68 | sim = reset() |
|---|
| 69 | #sim.setParam('npoints', npts) |
|---|
| 70 | |
|---|
| 71 | # Vectors of errors |
|---|
| 72 | err = [] |
|---|
| 73 | for i in range(80): |
|---|
| 74 | q = 0.40*(i) /50 |
|---|
| 75 | |
|---|
| 76 | |
|---|
| 77 | # For each Q, repeat the simulation N_q times |
|---|
| 78 | sum_q_ana = 0 |
|---|
| 79 | sum_q_sim = 0 |
|---|
| 80 | sum = 0 |
|---|
| 81 | vec = [] |
|---|
| 82 | for i_q in range(N_q): |
|---|
| 83 | counter += 1 |
|---|
| 84 | sim = reset(counter) |
|---|
| 85 | #sim.setParam('npoints', npts) |
|---|
| 86 | |
|---|
| 87 | simval = sim.run([q, 0]) |
|---|
| 88 | sum_q_ana += sph.run([q, 0]) |
|---|
| 89 | sum_q_sim += simval |
|---|
| 90 | |
|---|
| 91 | vec.append(simval) |
|---|
| 92 | |
|---|
| 93 | mean = sum_q_sim/N_q |
|---|
| 94 | for val in range(N_q): |
|---|
| 95 | diff = vec[val] - mean |
|---|
| 96 | sum += diff * diff |
|---|
| 97 | sum /= (N_q-1) |
|---|
| 98 | |
|---|
| 99 | ana_value = sum_q_ana/N_q |
|---|
| 100 | sim_value = sum_q_sim/N_q |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | try: |
|---|
| 104 | err_i = [q, (sim_value-ana_value)/ana_value, ana_value, math.sqrt(sum)/ana_value] |
|---|
| 105 | except: |
|---|
| 106 | err_i = [q,0,0,0] |
|---|
| 107 | |
|---|
| 108 | err.append(err_i) |
|---|
| 109 | |
|---|
| 110 | print err_i |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | # update extrema |
|---|
| 114 | if val_min > err_i[2]: |
|---|
| 115 | val_min = err_i[2] |
|---|
| 116 | |
|---|
| 117 | if val_max < err_i[2]: |
|---|
| 118 | val_max = err_i[2] |
|---|
| 119 | |
|---|
| 120 | # Find a minimum, where the error changes sign |
|---|
| 121 | if last_val==val_min and err_i[2] > val_min: |
|---|
| 122 | |
|---|
| 123 | |
|---|
| 124 | #print " min", err_i[0] |
|---|
| 125 | if i_min > 0: |
|---|
| 126 | i_plateau = i_last_max |
|---|
| 127 | |
|---|
| 128 | # check error for plateau |
|---|
| 129 | #if math.fabs(err[i_plateau][1])>sig_val: |
|---|
| 130 | if math.fabs(err[i_max][3])>sig_val: |
|---|
| 131 | error_vec.append(err[i_plateau]) |
|---|
| 132 | t_f = time.time() |
|---|
| 133 | pred = sph.run([err[i_plateau][0],0]) |
|---|
| 134 | print '[',t_f-t_0,']',math.pow((vol/npts),-.6666), 'n=',npts, 'q=',err[i_plateau][0], \ |
|---|
| 135 | 'err=',math.fabs(err[i_plateau][1]), i_min, 'sig=',err[i_plateau][3] |
|---|
| 136 | |
|---|
| 137 | f.write("%10g %10g %10g\n" % (err[i_plateau][0], npts, err[i_plateau][3])) |
|---|
| 138 | found_it = True |
|---|
| 139 | break |
|---|
| 140 | |
|---|
| 141 | i_min = i_vec |
|---|
| 142 | # Get ready for next maximum |
|---|
| 143 | val_max = val_min |
|---|
| 144 | |
|---|
| 145 | # check if we are at max |
|---|
| 146 | elif last_val==val_max and err_i[2]<val_max: |
|---|
| 147 | #print " max", err[i_vec-1][0], err[i_vec-1][1] |
|---|
| 148 | i_last_max = i_max |
|---|
| 149 | i_max = i_vec-1 |
|---|
| 150 | # get ready for next minimum |
|---|
| 151 | val_min = val_max |
|---|
| 152 | |
|---|
| 153 | last_error = err_i[1] |
|---|
| 154 | last_val = err_i[2] |
|---|
| 155 | i_vec += 1 |
|---|
| 156 | |
|---|
| 157 | |
|---|
| 158 | if found_it == False: |
|---|
| 159 | print "Could not complete", npts |
|---|
| 160 | #print "Npts = %g; %s" %(npts, str(len(error_vec))) |
|---|
| 161 | |
|---|
| 162 | |
|---|
| 163 | f.close() |
|---|
| 164 | |
|---|
