Commit a2006767 authored by Michele Valentini's avatar Michele Valentini

Merge branch 'master' of git.ligo.org:finesse/pykat

parents 31b7b797 757edaa8
name: pykat_test
name: pykatcli
channels:
- gwoptics
- defaults
dependencies:
- numpy>=1.12
- flask>=0.10.1
- scipy
- scipy>=0.14.0
- six
- h5py
- pandas
- matplotlib
- finesse=2.2
- finesse>=2.3.0
- ipython
\ No newline at end of file
......@@ -18,8 +18,11 @@ requirements:
- scipy
- matplotlib
- flask >=0.10.1
- finesse >=2.2
- finesse >=2.3
- h5py
- graphviz
- python-graphviz
- click
build:
- python
- setuptools
......@@ -28,11 +31,11 @@ requirements:
- scipy
- matplotlib
- flask >=0.10.1
- finesse >=2.2
- finesse >=2.3
test:
requires:
- finesse >=2.2
- finesse >=2.3
- six
- python
- numpy
......
......@@ -2,4 +2,7 @@
(cd .. && rm -r ./**/**.pyc)
(cd .. && python setup.py sdist --formats=zip upload)
(cd .. && python setup.py sdist --formats=gztar)
(cd .. && twine upload dist/*)
......@@ -11,7 +11,7 @@ try:
except (NameError, ImportError) as ex:
__version__ = "develop"
__min_req_finesse__ = 2.2
__min_req_finesse__ = '2.3.0'
# This flag is used to switch on the gui features in pkat at import time
USE_GUI = False
......@@ -27,7 +27,7 @@ is_container = lambda c: (not isinstance(c, six.string_types)) and hasattr(c, "_
import imp
try:
try:
imp.find_module('optivis')
HAS_OPTIVIS = True
except ImportError:
......@@ -65,11 +65,16 @@ except pkex.MissingFinesse:
warn(msg)
kat = nokat()
v = str(__min_req_finesse__)
if float(v.split('-')[0]) < __min_req_finesse__:
version_found = [int(_) for _ in v.split('-')[0].split('.')]
version_required = [int(_) for _ in __min_req_finesse__.split('-')[0].split('.')]
fulfilled = all([_ == __ for _,__ in zip(version_found, version_required)]) or any([_ > __ for _,__ in zip(version_found, version_required)])
if not fulfilled:
raise pkex.BasePyKatException("Pykat %s requires Finesse version %s or higher. You have have %s" % (__version__ ,
str(__min_req_finesse__),
v))
def info():
print("Pykat version: " + __version__)
print("Pykat loaded from: " + __file__)
......
"""Finesse command line interface
Sean Leavey
<sean.leavey@ligo.org>
"""
import sys
import os
import numpy as np
import click
from . import __version__
from .finesse import kat as katparser
@click.command(help="Python interface and tools for FINESSE")
@click.argument("file", type=click.File())
@click.option("--simulate/--no-simulate", is_flag=True, default=True,
help="Simulate FILE in Finesse. Can be set to --no-simulate if for example "
"you only want to display the node graph for the model specified in FILE.")
@click.option("--xstart", type=float,
help="Set simulation start value. If specified, this overrides the xaxis start "
"value specified in the parsed file.")
@click.option("--xstop", type=float,
help="Set simulation stop value. If specified, this overrides the xaxis stop "
"value specified in the parsed file.")
@click.option("--xsteps", type=int,
help="Set number of steps to simulate between --start and --stop. If specified, "
"this overrides the number of xaxis steps specified in the parsed file.")
@click.option("--xscale", type=click.Choice(["lin", "log"]),
help="Set scaling for the xaxis.")
@click.option("--noxaxis", is_flag=True, default=False, help="Switch off x-axis.")
@click.option("--trace", type=click.Choice(["tem", "cavity", "mismatch", "beams", "gouy",
"coupling", "modechanges", "nodes", "all"]),
multiple=True, help="Show simulation trace results. This option can be specified "
"multiple times. The following values are supported: "
"'tem': list TEM modes used, "
"'cavity': cavity eigenvalues and other parameters, "
"'mismatch': mode mismatch parameters for the initial setup, "
"'beams': beam parameters for every node, "
"'gouy': Gouy phases for all spaces, "
"'coupling': coupling coefficients for all components, "
"'modechanges': mode matching parameter changes during calculations, "
"'nodes': nodes found during cavity tracing, "
"'all': all trace results.")
@click.option("--powers", type=click.Choice(["dc", "carrier", "tem00"]),
multiple=True, help="Show powers (W) at each node in the interferometer. This "
"option can be specified multiple times. The following "
"values are supported: "
"'dc': list dc powers, "
"'carrier': list powers in the f=0 fields, "
"'tem00': list powers in the TEM00 mode.")
@click.option("--maxtem", type=str, help="Set maximum transverse electric mode. Can be either "
"an integer or 'off'.")
@click.option("--phase", type=int, help="Set Gouy phase behaviour.")
@click.option("--retrace", type=click.Choice(["force", "off"]),
help="Set retrace behaviour: 'force' recomputes the Gaussian parameters at each "
"node for every data point, and will trace a cavity even if it is unstable; 'off' "
"switches off retracing even if it normally would be done.")
@click.option("--deriv-h", type=float, help="Set step size for numerical differentiation.")
@click.option("--lambda0", type=str, help="Set reference wavelength (m). Supports SI prefixes.")
@click.option("--print-frequencies", is_flag=True, default=False,
help="Print a table listing all the frequencies used in the simulation: carriers, "
"modulation sidebands and signal/quantum sidebands.")
@click.option("--print-noises", is_flag=True, default=False,
help="Print a list of all the quantum noise sources being considered in the "
"simulation and at which components and nodes it is present.")
@click.option("--ignore-block", "ignored_blocks", multiple=True,
help="Ignore the specified block. Can be specified multiple times.")
@click.option("--plot/--no-plot", default=True, show_default=True,
help="Display results as figure.")
@click.option("--save-figure", type=click.File("wb", lazy=False),
help="Save image of figure to file.")
@click.option("--display-graph", is_flag=True, help="Generate and display model node graph "
"using default system document viewer.")
@click.option("--save-input", is_flag=True, default=False,
help="Save generated Finesse input file.")
@click.option("--save-output", is_flag=True, default=False,
help="Save generated Finesse output file.")
@click.option("--finesse-dir", type=click.Path(exists=True, file_okay=False, dir_okay=True),
envvar='FINESSE_DIR', help="Path to directory containing the Finesse 'kat' "
"executable. If not specified, the environment variable FINESSE_DIR is used.")
@click.version_option(version=__version__, prog_name="Pykat")
def cli(file, simulate, xstart, xstop, xsteps, xscale, noxaxis, trace, powers, maxtem, phase,
retrace, deriv_h, lambda0, print_frequencies, print_noises, ignored_blocks, plot,
save_figure, display_graph, save_input, save_output, finesse_dir):
"""Base CLI command group"""
if finesse_dir is None:
# Use default as required by kat object.
finesse_dir = ""
kat = katparser(katdir=finesse_dir)
kat.load(file.name)
has_xaxis = hasattr(kat, "xaxis") and not noxaxis
if ignored_blocks:
for block in ignored_blocks:
kat.removeBlock(block)
if display_graph:
from .tools.plotting.graph import NodeGraph
nodegraph = NodeGraph(kat)
nodegraph.view_pdf()
if simulate:
# Default flag for showing simulation output. Overridden below as necessary.
output_to_show = False
if xstart is not None or xstop is not None or xsteps is not None or xscale is not None:
if not has_xaxis:
click.echo("Limits can only be overridden when an xaxis is defined in FILE and "
"when --noxaxis is unset.", err=True)
sys.exit(1)
# Override xaxis.
limits = kat.xaxis.limits
set_limits = False
if xstart is not None:
limits[0] = xstart
set_limits = True
if xstop is not None:
limits[1] = xstop
set_limits = True
if xsteps is not None:
kat.xaxis.steps = xsteps
if xscale is not None:
kat.xaxis.scale = xscale
if set_limits:
kat.xaxis.limits = np.array(limits).astype(float)
if not has_xaxis:
kat.noxaxis = True
if save_figure is not None:
click.echo("Cannot plot or save figure without an xaxis defined in FILE.",
err=True)
sys.exit(1)
if maxtem is not None:
kat.maxtem = maxtem
if phase is not None:
kat.phase = phase
if retrace is not None:
kat.retrace = retrace
if deriv_h is not None:
kat.deriv_h = deriv_h
if lambda0 is not None:
kat.lambda0 = lambda0
if print_frequencies:
kat.parse("frequency")
output_to_show = True
if print_noises:
kat.parse("printnoises")
output_to_show = True
if trace:
traceval = 0
if "all" in trace:
traceval = 255
else:
traceints = {"tem": 1, "cavity": 2, "mismatch": 4, "beams": 8,
"gouy": 16, "coupling": 32, "modechanges": 64, "nodes": 128}
for tracetype in trace:
traceval |= traceints[tracetype]
kat.trace = traceval
output_to_show = True
if powers:
powerval = 0
powerints = {"dc": 1, "carrier": 2, "tem00": 4}
for powertype in powers:
powerval |= powerints[powertype]
kat.parse("powers %i" % powerval)
output_to_show = True
results = kat.run(save_output=save_output, save_kat=save_input)
if output_to_show and results.rundata:
click.echo(results.rundata)
if has_xaxis and (plot or save_figure is not None):
results.plot(show=plot, filename=save_figure)
elif not output_to_show and not display_graph:
click.echo("No output requested.")
......@@ -38,7 +38,6 @@ class Command(object):
self.__removed = False
self.__name = name.strip("*")
self._putters = []
def __deepcopy__(self, memo):
"""
......@@ -53,7 +52,7 @@ class Command(object):
for _ in result._putters:
_._updateOwner(result)
result._freeze()
return result
......@@ -203,11 +202,10 @@ class lock(Command):
def __init__(self, name, variable, gain, accuracy, offset=None, singleLock=False):
Command.__init__(self, name, False)
self.__params = []
self.__variable = variable
self.__gain = gain
self.__accuracy = accuracy
self.__offset = Param("offset", self, SIfloat(offset))
self.__offset = offset
self.singleLock = singleLock
self.enabled = True
......@@ -217,12 +215,11 @@ class lock(Command):
self._putters.append(self.output)
self._freeze()
def _register_param(self, param):
self.__params.append(param)
@property
def offset(self): return self.__offset
@offset.setter
def offset(self, value): self.__offset = SIfloat(value)
@staticmethod
def parseFinesseText(line, kat):
......@@ -245,7 +242,7 @@ class lock(Command):
gain=str(self.gain),
accuracy=str(self.accuracy))
if self.offset.value is not None:
if self.offset is not None:
cmds += " %s" % str(self.offset)
if self.singleLock:
......@@ -255,9 +252,6 @@ class lock(Command):
rtn.extend(self.output._getPutFinesseText())
for p in self.__params:
rtn.extend(p.getFinesseText())
return rtn
else:
return None
......@@ -533,6 +527,36 @@ class tf2(Command):
return rtn
class tf3(Command):
def __init__(self, name, hdf, frequency_dataset, data_dataset):
Command.__init__(self, name, False)
self.hdf = hdf
self.frequency_dataset = frequency_dataset
self.data_dataset = data_dataset
self._freeze()
@staticmethod
def parseFinesseText(text):
values = text.split()
if len(values) != 5:
print(values)
raise pkex.BasePyKatException("Transfer function Finesse code format incorrect '{0}'".format(text))
_tf = tf3(*values[1:])
return _tf
def getFinesseText(self):
rtn = "tf3 {name} {hdf} {fds} {dds} ".format(name=self.name,
hdf=self.hdf,
fds=self.frequency_dataset,
dds=self.data_dataset)
return rtn
class xaxis(Command):
"""
The xaxis object is a unique object to each pykat.finesse.kat instance. It provides
......
This diff is collapsed.
......@@ -64,6 +64,11 @@ class BasePyKatException(Exception):
def __str__(self):
return self.msg
class LockLossException(BasePyKatException):
def __init__(self, step):
self.step = step
BasePyKatException.__init__(self, "Lock was lost after {} steps.".format(step))
class FinesseParse(BasePyKatException) :
def __init__(self, msg):
BasePyKatException.__init__(self, "Error parsing Finesse input\n{0}".format(msg))
......
......@@ -104,7 +104,7 @@ def vectfit_step(f, s, poles):
H = real(H)
new_poles = sort(eigvals(H))
unstable = real(new_poles) > 0
new_poles[unstable] -= 2*real(new_poles)[unstable]
#new_poles[unstable] -= 2*real(new_poles)[unstable]
return new_poles
# Dear gods of coding style, I sincerely apologize for the following copy/paste
......@@ -141,8 +141,8 @@ def calculate_residues(f, s, poles, rcond=-1):
b = concatenate((real(b), imag(b)))
cA = np.linalg.cond(A)
if cA > 1e13:
print ('Warning!: Ill Conditioned Matrix. Consider scaling the problem down')
print ('Cond(A)', cA)
print('Warning!: Ill Conditioned Matrix. Consider scaling the problem down')
print('Cond(A)', cA)
x, residuals, rnk, s = lstsq(A, b, rcond=rcond)
# Recover complex values
......@@ -160,10 +160,10 @@ def calculate_residues(f, s, poles, rcond=-1):
def print_params(poles, residues, d, h):
cfmt = "{0.real:g} + {0.imag:g}j"
print ("poles: " + ", ".join(cfmt.format(p) for p in poles))
print ("residues: " + ", ".join(cfmt.format(r) for r in residues))
print ("offset: {:g}".format(d))
print ("slope: {:g}".format(h))
print("poles: " + ", ".join(cfmt.format(p) for p in poles))
print("residues: " + ", ".join(cfmt.format(r) for r in residues))
print("offset: {:g}".format(d))
print("slope: {:g}".format(h))
def vectfit_auto(f, s, n_poles=10, n_iter=10, show=False,
inc_real=False, loss_ratio=1e-2, rcond=-1, track_poles=False):
......@@ -185,20 +185,20 @@ def vectfit_auto(f, s, n_poles=10, n_iter=10, show=False,
if track_poles:
return poles, residues, d, h, np.array(poles_list)
print_params(poles, residues, d, h)
#print_params(poles, residues, d, h)
return poles, residues, d, h
def vectfit_auto_rescale(f, s, **kwargs):
s_scale = abs(s[-1])
f_scale = abs(f[-1])
print( 'SCALED')
#print('SCALED')
poles_s, residues_s, d_s, h_s = vectfit_auto(f / f_scale, s / s_scale, **kwargs)
poles = poles_s * s_scale
residues = residues_s * f_scale * s_scale
d = d_s * f_scale
h = h_s * f_scale / s_scale
print( 'UNSCALED')
print_params(poles, residues, d, h)
#print('UNSCALED')
#print_params(poles, residues, d, h)
return poles, residues, d, h
if __name__ == '__main__':
......@@ -241,4 +241,3 @@ if __name__ == '__main__':
plot(test_s.imag, test_f.imag)
plot(test_s.imag, fitted.real)
plot(test_s.imag, fitted.imag)
show()
\ No newline at end of file
This diff is collapsed.
......@@ -1133,7 +1133,7 @@ class IFO(object):
self.__kat = kat
self.__tuning_keys = frozenset(make_list_copy(tuning_keys_list))
self.__tuning_comps = make_list_copy(tunings_components_list[:])
@property
def kat(self): return self.__kat
......@@ -1184,6 +1184,84 @@ class IFO(object):
print("New tunings")
print(self.get_tunings())
def lock_drag(self, N, *args, update_state=True, **kwargs):
"""
Performs a "lock drag" to gradually change the state of the interferometer model
from one state to another. This is used when wanting to keep a model at a chosen
operating point but you want to explore how changes in various parameters affect
properties like power buildups, noise couplings, transfer functions, etc.
Essentially what this function does is slowly drag several parameters from some
initial value to some final one, whilst keeping the locks activated. The locks
must be setup correctly and enabled before this function is called.
Example
-------
Here we take some base model that has a mirror which we want to change the curvature
of whilst keeping the interferometer locked. In this example the radius of curvature
of ETMX will be changed from its current value to +10m in x and -10m in y.
>>> kat = base.deepcopy()
>>> kat.maxtem = 2
>>> kat.verbose = True
>>>
>>> drag = lock_drag(kat, 100,
>>> ('ETMX', 'Rcx', 10),
>>> ('ETMX', 'Rcy', -10)
>>> )
To apply the lock drag to the object you need to extract the `lock` outputs and apply
them to the model.
Parameters
----------
N : int
Number of points to do the lock drag over
update_state : bool
If True, the calling kat object's state will be updated with the
final value of the lock drag
*args : tuple(target, param, final_value)
target : Name of the component
param : Parameter to change
final_value : The relative change in value to do the lock drag over
**kwargs : dict
Keyword arguments are passed to the `kat.run(**kwargs)` call.
"""
if hasattr(self, 'xaxis'): self.xaxis.remove()
if hasattr(self, 'x2axis'): self.x2axis.remove()
kat = self.kat.deepcopy()
kat.parse(f"""
var dummy 0
xaxis dummy re lin 0 1 {N}
""")
for i, (target, attr, final) in enumerate(args):
kat.parse(f"""
func LD{i} = ({final:.15G}) * $x1
put* {target} {attr} $LD{i}
""")
drag = kat.run(**kwargs)
if update_state:
self.apply_lock_drag(drag, -1, *args)
return drag
def apply_lock_drag(self, out, idx, *args):
for i, (target, attr, final) in enumerate(args):
p = [_ for _ in self.kat.components[target]._params if _.name == attr]
if len(p) != 1:
raise Exception(f"Could not find parameter {attr} for component {target}")
else:
p[0].value += final * out.x[idx]
self.apply_lock_feedback(out, idx)
def get_tuning_comps(self):
return self.__tuning_comps
......
......@@ -193,10 +193,21 @@ class ALIGO_IFO(IFO):
This function will iterate through the main mirrors
and remove any suspension settings on them. This can be
done individuallly or for z, pitch, and yaw.
kat code for the original attr commands is returned in a
dict object separated into z, pitch and yaw.
"""
old_attr={'z':[],'pitch':[],'yaw':[]}
for mirror in ["ETMY","ETMX","ITMY","ITMX","PRM","PR2","PR3","SRM","SR2","SR3","BS"]:
mirror = self.kat.components[mirror]
if mirror.mass != None: old_attr['z'].append(mirror.mass.getFinesseText()[0])
if mirror.zmech != None: old_attr['z'].append(mirror.zmech.getFinesseText()[0])
if mirror.Iy != None: old_attr['pitch'].append(mirror.Iy.getFinesseText()[0])
if mirror.rymech != None: old_attr['pitch'].append(mirror.rymech.getFinesseText()[0])
if mirror.Ix != None: old_attr['yaw'].append(mirror.Ix.getFinesseText()[0])
if mirror.rxmech != None: old_attr['yaw'].append(mirror.rxmech.getFinesseText()[0])
if z:
mirror.mass = None
......@@ -209,6 +220,29 @@ class ALIGO_IFO(IFO):
if yaw:
mirror.Ix = None
mirror.rxmech = None
return old_attr
def restore_susp(self, old_attr, z=True, pitch=True, yaw=True):
"""
This is a pair function to fix_mirrors.
It will restore any suspension settings previously stored as kat code strings
in a dict object of format
old_attr={'z':[],'pitch':[],'yaw':[]}
This can be done individuallly or for z, pitch, and yaw.
"""
if z:
for ii in old_attr['z']:
self.kat.parse(ii)
if pitch:
for ii in old_attr['pitch']:
self.kat.parse(ii)
if yaw:
for ii in old_attr['yaw']:
self.kat.parse(ii)
def lengths_status(self):
self.compute_derived_lengths()
......@@ -256,8 +290,9 @@ class ALIGO_IFO(IFO):
spaces between the laser and HAM2 and PRC. Assumes spaces exists
with name and node:
sHAM2in and node nIMCout
sPRCin and node nHAM2out
If removing HAM2, adds a replacement dbs, `FI`, directly before the PRM,
to restore the REFL port.
This function alters the kat object directly.
"""
......@@ -272,7 +307,11 @@ class ALIGO_IFO(IFO):
if removeHAM2:
self.kat.removeBlock("HAM2")
self.kat.nodes.replaceNode(self.kat.sPRCin, 'nHAM2out', 'nLaserOut')
# self.kat.nodes.replaceNode(self.kat.sPRCin, 'nHAM2out', 'nLaserOut') #without FI restoration.
self.kat.parse("""
s sIO 0 nLaserOut nFI1
dbs FI nFI1 nFI2 nHAM2out nREFL
""",addToBlock='PRC') #dummy space needed between mod2 and FI.
def remove_FI_OMC(self, removeFI=True, removeOMC=True):
......@@ -322,6 +361,25 @@ class ALIGO_IFO(IFO):
kat.lp1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def adjust_SRC_length(self, verbose=False):
"""
Adjust SRC length so that it fulfils the requirement
lSRC = (M/5) * c/(2*f1), see [1] equation C.1
In the current design M=17 and N=3.4.
This function directly alters the lengths of the associated kat object.
"""
kat = self.kat
vprint(verbose, "-- adjusting SRC length")
ltmp = 0.5 * clight / kat.IFO.f1
delta_l = 3.4 * ltmp - kat.IFO.lSRC
vprint(verbose, " adusting kat.ls1.L by {:.4g}m".format(delta_l))
kat.ls1.L += delta_l
kat.IFO.compute_derived_lengths(kat)
def apply_lock_feedback(self, out, idx=None):
"""
......@@ -819,20 +877,13 @@ def make_kat(name="design", katfile=None, verbose = False, debug=False, use_RF_D
kat.IFO.POP_f2 = Output(kat.IFO, "POP_f2", "nPOP", "f2", phase=13)
kat.IFO.REFL_f1 = Output(kat.IFO, "REFL_f1", "nREFL", "f1", phase=101)
kat.IFO.REFL_f2 = Output(kat.IFO, "REFL_f2", "nREFL", "f2", phase=14)
# If we don't have an OMC then we need to attach
# directly to the AS node. Otherwise use OMC refl
if "OMC" in kat.getBlocks():
alt_nAS = "nAS"
nAS_RF = "nOMC_ICb"
else:
alt_nAS = None
nAS_RF = "nAS"
kat.IFO.AS_f1 = Output(kat.IFO, "AS_f1", nAS_RF, "f1", phase=101, alternate_node_name=alt_nAS)
kat.IFO.AS_f2 = Output(kat.IFO, "AS_f2", nAS_RF, "f2", phase=14, alternate_node_name=alt_nAS)
kat.IFO.AS_f36 = Output(kat.IFO, "AS_f36", nAS_RF, "f36M", phase=14, alternate_node_name=alt_nAS)
kat.IFO.AS_f1 = Output(kat.IFO, "AS_f1", "nSRM2", "f1", phase=101)
kat.IFO.AS_f2 = Output(kat.IFO, "AS_f2", "nSRM2", "f2", phase=14)