l3gui/library.l3

   1 #* Welcome to the L3 worksheet interface.  For help, copyright, or
   2 # license information, see the help menu.
   3 #
   4 if ("outline",  "L3 version 0.3.1"):
   5     #* These examples introduce features of l3.  To try one:
   6     # - left-click on it; an orange outline should highlight it
   7     # - center-click on the right canvas to copy the example to it
   8     # - in the copy, left-click once on the triangle to get started
   9     # - follow the instructions in the example.
  10     if ("outline",  "Introductory examples."):
  11
  12         # The simplest example, illustrating need-based evaluation in l3.
  13         # Expand to start (left-click the triangle).
  14         if ("prog", "hello world 1"):
  15             # 1. Right-click on the headline, and choose "run code".
  16             #
  17             # 2. Watch your terminal window; it should show
  18             #       hello world
  19             #       (None, 8)
  20             #    (or some other number)
  21             #
  22             # 3. Choose "run code" again; the output should be just
  23             #       (None, 8)
  24             #    without a second "hello world"
  25             print "hello world"
  26
  27         # The simplest example, illustrating persistence in l3.
  28         # Expand to start (left-click the triangle).
  29         if ("prog", "hello world 2"):
  30             # 1. Right-click on the headline, and choose "run code".
  31             # 2. Right-click on the "hello world" string below, and
  32             #    choose "dump values".  This should show
  33             #         32774 does not repeat; value:  hello world
  34             #    or similar.
  35             print "hello world"
  36
  37             # 1. Save the worksheet via `file -> save session to`, use
  38             #    `st.1` as file name.
  39             # 2. Open the worksheet in a new session via
  40             #    l3gui -s st.1
  41             #
  42             # 3. In the new session, right-click on the "hello world"
  43             #    string, and choose "dump values".  This should show
  44             #         32774 does not repeat; value:  hello world
  45             #    again.
  46             1
  47
  48
  49
  50     #* Simple, single-purpose L3 scripts can be found below.  These
  51     # can be run by as-is, or modified and combined to form larger
  52     # scripts.
  53
  54     if ("outline",  "Short examples"):
  55
  56         #** [requires matplotlib] Collect finance data using yahoo's service.
  57         if ("outline", "Collect finance data"):
  58             inline """if 1:
  59             from matplotlib.finance import quotes_historical_yahoo
  60             import datetime
  61             """
  62             # <span background="#ffdddd">Retrieve data for stock symbol.</span>
  63             quotes = quotes_historical_yahoo(
  64                 'INTC',
  65                 datetime.date( 1995, 1, 1 ),
  66                 datetime.date( 2004, 4, 12 ),
  67                 asobject = True)
  68             if not quotes:
  69                 error("Unable to get data")
  70             else:
  71                 # <span background="#ffdddd">The data attributes of quotes.</span>
  72                 quotes.date
  73                 quotes.open
  74                 quotes.close
  75                 quotes.high
  76                 quotes.low
  77                 quotes.volume
  78
  79
  80
  81         #
  82         #** [requires matplotlib] Date plot.
  83         if ("outline", "Date plot"):
  84             # <span background="#ffdddd">Plot details.</span>
  85
  86             inline """if 1:
  87             import datetime
  88
  89             def plot_date(fname, dates, values):
  90                 from pylab import figure, show
  91                 import datetime
  92                 from matplotlib.dates import YearLocator, MonthLocator, DateFormatter
  93
  94                 fig = figure(figsize=(7,4), dpi=100)
  95                 ax = fig.add_subplot(111)
  96                 ax.plot_date(dates, values, '-')
  97
  98                 # Plot details -- ticks.
  99                 years    = YearLocator()   # every year
 100                 months   = MonthLocator()  # every month
 101                 yearsFmt = DateFormatter('%Y')
 102
 103                 ax.xaxis.set_major_locator(years)
 104                 ax.xaxis.set_major_formatter(yearsFmt)
 105                 ax.xaxis.set_minor_locator(months)
 106                 ax.autoscale_view()
 107
 108                 for label in ax.get_xticklabels() + ax.get_yticklabels():
 109                     label.set_fontsize('8')
 110
 111                 # Plot details -- coords. message box.
 112                 def price(x):
 113                     return '$%1.2f' % x
 114                 ax.fmt_xdata = DateFormatter('%Y-%m-%d')
 115                 ax.fmt_ydata = price
 116                 ax.grid(True)
 117
 118                 fig.savefig(fname)
 119                 return fname
 120             """
 121             # <span background="#ffdddd">Plot data.</span>
 122             filename = plot_date(
 123                 # <span background="#ffdddd">File name.</span>
 124                 '_plot-%d.png' % new_id(),
 125                 # <span background="#ffdddd">List of dates.</span>
 126                 [datetime.date(1995,1,1).toordinal(),
 127                  datetime.date(1997,1,1).toordinal()],
 128                 # <span background="#ffdddd">List of values.</span>
 129                 [1, 2])
 130         #
 131         #** [requires numpy] Moving average
 132         if ("outline", "Moving average."):
 133             inline ("from l3lang.external.numpy_utils import movavg")
 134             mov_avg = movavg(
 135                       # <span background="#ffdddd">Data.</span>
 136                       numpy.array([1, 2, 3, 4.0]),
 137                       # <span background="#ffdddd">Averaging interval.</span>
 138                       2)
 139         #** [requires matplotlib] Graphing X-Y data.
 140         if ("outline",  "XY Plot"):
 141             # <span background="#ffcccc">Plot details.</span>
 142             inline '''if 1:
 143             def plot_(fname, x, y):
 144                 import matplotlib as M
 145                 import pylab as P
 146                 M.use('Agg')
 147
 148                 fig = P.figure(figsize=(5,4), dpi=75)
 149                 ax = fig.add_subplot(111)
 150                 ax.plot(x, y, '-', lw=2)
 151
 152                 ax.set_xlabel('time (s)')
 153                 ax.set_ylabel('voltage (mV)')
 154                 ax.set_title('About as simple as it gets, folks')
 155                 ax.grid(True)
 156
 157                 fig.savefig(fname)
 158                 return fname
 159             '''
 160             inline 'from math import pi; import pylab as P'
 161             # <span background="#ffdddd">Plot data.</span>
 162             x = P.arange(0.0, 1.0+0.01, 0.01)
 163             y = P.cos(2 * pi * x)
 164             # <span background="#ffdddd">Plot file.</span>
 165             filename = plot_('_plot-%d.png' % new_id(), x, y)
 166
 167         if ("outline",  "XY Polygon Plot"):
 168             # <span background="#ffcccc">Plot details.</span>
 169             inline '''if 1:
 170             def plot_(fname, x, y):
 171                 import matplotlib as M
 172                 import pylab as P
 173                 M.use('Agg')
 174
 175                 fig = P.figure(figsize=(5,4), dpi=75)
 176                 ax = fig.add_subplot(111)
 177                 ax.plot(x, y, '-', lw=2)
 178
 179                 ax.set_xlabel('')
 180                 ax.set_ylabel('')
 181                 ax.set_title('Polygon plot')
 182                 ax.grid(True)
 183
 184                 fig.savefig(fname)
 185                 return fname
 186             '''
 187             inline 'from pylab import *'
 188             for N in range(3, 12):
 189                 nn = arange(0, N + 1)
 190                 xx = sin(2* nn* pi/N)
 191                 yy = cos(2* nn* pi/N)
 192                 plot_(new_name('poly', '.png'), xx, yy)
 193
 194
 195     #* Longer examples illustrating full workflows done via L3.
 196     if ("outline",  "Full examples"):
 197
 198         #* [requires sparx] An image processing example from the sparx toolkit.
 199         if ("prog", "Image processing"):
 200             # From a volume, images are projected in known directions with
 201             # known angles; an alignment algorithm is called to find the
 202             # number of directions, and reverse the angular spread.
 203             if ("outline",  "SPARX volume projection and image alignment"):
 204                 #* Form the test data, projections from volume.
 205                 if ("outline", "Project volume."):
 206                     if ("outline", "Preparation."):
 207                         inline('''if 1:
 208                                from EMAN2  import *
 209                                from sparx  import *
 210                                from random import randint
 211                                import l3lang.external.sparx_mixins
 212                                import os
 213                                ''')
 214                         ev = os.path.expandvars
 215
 216                     if ("outline", "Load volume."):
 217                         vol = getImage(ev("$SPXROOT/sparx/test/model001.tcp"))
 218
 219                     if ("outline", "Compute projection angles"):
 220                         group_angles = even_angles(20.0,
 221                                                    70.0, 90.1,
 222                                                    0.0, 110.0, "P")
 223                         sx = 0
 224                         sy = 0
 225                         alpha = 0.0
 226
 227                     # In order to have a 'correct' K in k_means lateron, the initial
 228                     # projections (phi, theta, psi) must have multiple psi per (phi,
 229                     # theta) pair.
 230                     if ("outline", "Use multiple rotation angles (psi)."):
 231                         num_groups = len(group_angles)
 232                         angles = []
 233                         for (phi, theta, psi) in group_angles:
 234                             for proj_rot in range(0,180,10):
 235                                 angles.append([phi, theta, float(proj_rot)])
 236
 237
 238                     if ("outline", "Projection setup."):
 239                         (volft, kb) = prep_vol(vol)
 240
 241                     if ("outline", "Projection output target."):
 242                         stack = "proj.hdf"
 243
 244                     if ("outline", "Project angles..."):
 245                         for ii in range(0, len(angles)):
 246                             print  ii, angles[ii][0], angles[ii][1], angles[ii][2], -sx, -sy
 247                             if ("outline", "Generate projection."):
 248                                 proj = prgs(volft, kb,
 249                                             [angles[ii][0], angles[ii][1], angles[ii][2],
 250                                              -sx, -sy])
 251                             # For 2D alignment, must set at least
 252                             # 'alpha', 'sx', 'sy', and 'mirror'.
 253                             if ("outline", "Set projection parameters"):
 254                                 if ("outline", "three angles and two shifts"):
 255                                     proj.set(phi = angles[ii][0],
 256                                              theta = angles[ii][1],
 257                                              psi = angles[ii][2],
 258                                              alpha = alpha,
 259                                              sx = sx,
 260                                              sy = sy,
 261                                              mirror = 0.)
 262                                 if ("outline", "CTF parameters"):
 263                                     proj.set(defocus = 0.0,
 264                                              amp_contrast = 0.1,
 265                                              voltage = 200,
 266                                              Cs = 2.0,
 267                                              pixel = 2.2)
 268                                 if ("outline", "flags describing the status of the image"):
 269                                     proj.set(active = 1,
 270                                              ctf_applied = 0)
 271                             if ("outline", "Save"):
 272                                 proj.write_image(stack, ii)
 273
 274
 275                 #* First, define a function in case we want to try multiple alignments.
 276                 if ("outline", "Function to align image series"):
 277                     def ali2d_c(stack, outdir, maskfile = None, ir = 1, ou = -1, rs = 1, xr = 0, yr = 0, ts = 1, center = 1, maxit = 10, CTF = False, snr = 1.):
 278                         if ("outline", "Prepare input data"):
 279                             inline '''if 1:
 280                             from utilities    import model_circle, combine_params2, dropImage, getImage, get_arb_params, center_2D
 281                             from statistics   import aves, ave_oe_series, fsc, add_oe_series
 282                             from alignment    import Numrinit, ringwe, Applyws, ali2d_s
 283                             from filter       import filt_from_fsc_bwt, filt_table, filt_ctf
 284                             from filter       import filt_params, filt_btwl
 285                             from fundamentals import fshift
 286                             from morphology   import ctf_2
 287                             from math import exp
 288                             import os
 289
 290                             from utilities import print_begin_msg, print_end_msg, print_msg
 291                             import  types
 292                             from filter import filt_ctf
 293                             '''
 294
 295                             print_begin_msg("ali2d_c")
 296                             first_ring=int(ir); last_ring=int(ou); rstep=int(rs); xrng=int(xr); yrng=int(yr); step=int(ts); max_iter=int(maxit);
 297
 298                             print_msg("Input stack                 : %s\n"%(stack))
 299                             print_msg("Output directory            : %s\n"%(outdir))
 300                             print_msg("Inner radius                : %i\n"%(first_ring))
 301
 302                             nima = EMUtil.get_image_count(stack)
 303                             ima = EMData()
 304                             ima.read_image(stack, 0)
 305                             nx = ima.get_xsize()
 306                             # default value for the last ring
 307                             if (last_ring == -1):
 308                                 last_ring = nx//2-2
 309
 310                             print_msg("Outer radius                : %i\n"%(last_ring))
 311                             print_msg("Ring step                   : %i\n"%(rstep))
 312                             print_msg("X search range              : %i\n"%(xrng))
 313                             print_msg("Y search range              : %i\n"%(yrng))
 314                             print_msg("Translational step          : %i\n"%(step))
 315                             print_msg("Center type                 : %i\n"%(center))
 316                             print_msg("Maximum iteration           : %i\n"%(max_iter))
 317                             print_msg("data with CTF               : %s\n"%(CTF))
 318                             print_msg("Signal-to-Noise Ratio       : %f\n\n"%(snr))
 319
 320                             if os.path.exists(outdir):
 321                                 os.system('rm -rf '+outdir)
 322                             os.mkdir(outdir)
 323
 324                             if maskfile:
 325                                 if(type(maskfile) is types.StringType):
 326                                     print_msg("Maskfile                    : %s\n"%(maskfile))
 327                                     mask=getImage(maskfile)
 328                                 else:
 329                                     print_msg("Maskfile                    : user provided in-core mask")
 330                                     mask = maskfile
 331                             else :
 332                                 print_msg("Maskfile                    : default, a circle with radius %i\n"%(last_ring))
 333                                 mask = model_circle(last_ring, nx, nx)
 334                             cnx = int(nx/2)+1
 335                             cny = cnx
 336                             mode = "F"
 337                             data = []
 338                             if(CTF):
 339                                 parnames = ["Pixel_size", "defocus", "voltage", "Cs", "amp_contrast", "B_factor",  "ctf_applied"]
 340                                 #                        0                1              2     3              4                5                6
 341                                 ctf_params = get_arb_params(ima, parnames)
 342                                 data_had_ctf = ctf_params[6]
 343                                 ctm = ctf_2(nx, ctf_params[0], ctf_params[1], ctf_params[2], ctf_params[3], ctf_params[4], ctf_params[5])
 344                                 lctf = len(ctm)
 345                                 ctf2 = []
 346                                 ctf2.append([0.0]*lctf)
 347                                 ctf2.append([0.0]*lctf)
 348                                 ctfb2 = [0.0]*lctf
 349                             for im in range(nima):
 350                                 if(im > 0):
 351                                     ima = EMData()
 352                                     ima.read_image(stack, im)
 353                                 if(CTF):
 354                                     ctf_params = get_arb_params(ima, parnames)
 355                                     ctm = ctf_2(nx, ctf_params[0], ctf_params[1], ctf_params[2], ctf_params[3], ctf_params[4], ctf_params[5])
 356                                     k = im%2
 357                                     for i in range(lctf):
 358                                         # ctf2[k][i] += ctm[i]
 359                                         ctf2[k][i] = ctf2[k][i] + ctm[i]
 360                                     if(ctf_params[6] == 0):
 361                                         #      filt_ctf(image,       dz,          cs,          voltage,      pixel_size,     amp_contrast=0.1,    b_factor=0.0):
 362                                         ima = filt_ctf(ima, ctf_params[1], ctf_params[3], ctf_params[2], ctf_params[0], ctf_params[4], ctf_params[5])
 363                                         ima.set_attr('ctf_applied', 1)
 364                                 data.append(ima)
 365                             if(CTF):
 366                                 for i in range(lctf):
 367                                     ctfb2[i] = 1.0/(ctf2[0][i] + ctf2[1][i] + 1.0/snr)
 368                                     for k in range(2):
 369                                         ctf2[k][i] = 1.0/(ctf2[k][i] + 1.0/snr)
 370                         if ('outline', 'algorithm start'):
 371                             # precalculate rings
 372                             numr = Numrinit(first_ring, last_ring, rstep, mode)
 373                             wr = ringwe(numr, mode)
 374                             (av1, av2) = add_oe_series(data)
 375                             Util.add_img(av1, av2)
 376                             if(CTF):
 377                                 tavg = filt_table(av1, ctfb2)
 378                             else:
 379                                 tavg = av1/nima
 380                             a0 = tavg.cmp("dot", tavg, dict(negative = 0, mask = mask))
 381                             msg = "Initial criterion = %-20.7e\n"%(a0)
 382                             print_msg(msg)
 383
 384                             cs = [0.0]*2
 385
 386                         if ('outline', 'iterate'):
 387                             for Iter in range(max_iter):
 388                                 ali2d_s(data, numr, wr, cs, tavg, cnx, cny, xrng, yrng, step, mode, range(nima))
 389                                 (av1, av2) = add_oe_series(data)
 390                                 if(CTF):
 391                                     tavg = filt_table(Util.addn_img(av1, av2), ctfb2)
 392                                     av1  = filt_table(av1, ctf2[0])
 393                                     av2  = filt_table(av2, ctf2[1])
 394                                 else:
 395                                     tavg = (av1*(int(nima/2)+(nima%2)) + av2*int(nima/2))/nima
 396
 397                                 # Apply filter (FRC) to the reference image
 398                                 dropImage(tavg, os.path.join(outdir, "aqc_%03d.hdf"%(Iter+1)))
 399                                 frsc = fsc(av1, av2, 1.0, os.path.join(outdir, "drc%03d"%(Iter+1)))
 400                                 #  apply filtration (FRC) to reference image:
 401                                 #fil_table = filt_from_fsc_bwt(frsc, low)
 402                                 #tavg = filt_table(tavg, fil_table)
 403                                 (lowfq, highfq) = filt_params(frsc, low = 0.1)
 404                                 tavg  = filt_btwl(tavg, lowfq, highfq)
 405                                 (tavg, csx, csy) = center_2D(tavg, center)
 406                                 # print  lowfq, highfq, csx, csy
 407                                 # a0 should increase; stop algorithm when it decreases.
 408                                 a1 = tavg.cmp("dot", tavg, dict(negative = 0, mask = mask))
 409                                 msg = "ITERATION #%3d        criterion = %20.7e\n"%(Iter+1,a1)
 410                                 print_msg(msg)
 411                                 # write the current average
 412                                 dropImage(tavg, os.path.join(outdir, "aqf_%03d.hdf"%(Iter+1)))
 413                                 if(a1 < a0):
 414                                     break
 415                                 else:
 416                                     a0 = a1
 417
 418                         if ('outline', 'write out headers'):
 419                             if(CTF):
 420                                 if (data_had_ctf == 0):
 421                                     for im in range(nima):
 422                                         data[im].set_attr('ctf_applied', 0)
 423                             for im in range(nima):
 424                                 # Using True writes only headers.
 425                                 data[im].write_image(stack, im, EMUtil.ImageType.IMAGE_HDF, True)
 426                         print_end_msg("ali2d_c")
 427
 428                 #* Just do one alignment for now.
 429                 if ("outline", "Perform an alignment"):
 430                     ali2d_c(stack, "_ali2d-%d" % new_id())
 431
 432                 #* Alignment parameters are not applied to images; apply
 433                 # parameters first via transform2d.
 434                 if ('outline', 'Get rotated / shifted image series.'):
 435                     shifted_stack = "proj-trans.hdf"
 436                     transform2d(stack, shifted_stack)
 437
 438                 #* Classify images via k-means, using different guesses for K.
 439                 if ('outline', 'Classify images.'):
 440                     for K in [2, 4, num_groups, num_groups + 2, num_groups + 4]:
 441                         (coleman, harabasz, ave, var, chart) = k_means_info(
 442                             shifted_stack, "_K-%d-%d" % (K, new_id()),
 443                             trials = 10, K = K, maxit = 10)
 444
 445     #* The building blocks for program assembly.
 446     if ("outline",  "Language constructs"):
 447
 448         # Programs are the only executable construct.
 449         if ("prog", "Program template"):
 450             "content_here"
 451
 452         # Outlines provide a way to visually structure the code.  They
 453         # do not affect code execution in any way.
 454         if ("outline",  "outline template"):
 455             0
 456
 457         # The loop to use when examining a sequence one item at a
 458         # time.
 459         if ("outline", "for loop"):
 460             for n in range(0,11):
 461                 print n
 462
 463         # Repeat code while a condition holds.
 464         if ("outline", "while loop"):
 465             a = 1
 466             while a < 5:
 467                 print a
 468                 a = a + 1
 469                 if a + 4 > 10:
 470                     break
 471
 472         # The most flexible looping construct.  It avoids updating
 473         # a name's value providing much cleaner data flow than the
 474         # `while` loop, and nested loops can be restarted and exited.
 475         if ("outline", "Generic tail-call loop"):
 476             def loop(a):
 477                 if a < 3:
 478                     print a
 479                     loop(a + 1)
 480             loop(1)
 481
 482
 483         #* Basic L3 constructs
 484         #
 485         if ("outline", "More..."):
 486             if ("outline", "set -- Bind a value to a name"):
 487                 (a = 2)
 488
 489             if ("outline", "inline -- Raw python code, evaluated every time."):
 490                 inline "import numpy as N"
 491
 492             if ("outline", "subdir -- A nested namespace linked to a new directory."):
 493                 subdir:
 494                     content_here
 495
 496             if ("prog", "program -- The only executable construct."):
 497                 content_here
 498
 499             if ("outline", "def / lambda -- function definitions."):
 500                 def add_num(x,y, accum = 0.0):
 501                     return x + y + accum
 502                 # lambda -- An unnamed inline function
 503                 lambda x,y:
 504                     x - y
 505                 # Alternate syntax.
 506                 { |x,y| x - y }
 507
 508             if ("outline", "call -- Call a function with arguments."):
 509                 call_function(a, b)
 510
 511             if ("outline", "list/tuple -- Ordered lists of items."):
 512                 []
 513                 # Long form.
 514                 [
 515                     ]
 516                 (1,2)
 517                 # Long form.
 518                 (1,
 519                  2)
 520
 521             if ("outline", "map -- Unordered mapping from keys to values."):
 522                 {a = 1, b = 2}
 523                 # Long form.
 524                 {a = 1,
 525                  b = 2}
 526
 527             if ("outline", "if -- Conditional execution of code."):
 528                 if condition:
 529                     yes
 530                 else:
 531                     no
 532
 533             if ("outline", "print -- Print a value one time."):
 534                 print "hello world"
 535
 536             if ("outline", "string -- string literal and text starter."):
 537                 "hello"
 538