[ae3ce4e] | 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 | |
---|