U
    Åmœd·  ã                   @   s‚  d Z ddlZddlZddlZddlmZ ddlmZ ddlm	Z	m
Z
mZmZmZ ddlZddlZddlmZ ddlmZmZmZ dd	lmZ dd
lmZ dJed ee	e e	e e	e e	e e	e e	e e	e e	e
e ef  edœdd„Z!dd„ Z"dd„ Z#edƒddei ƒe $g ¡e $g ¡ei ƒdddf
ee ef ee e f dœdd„Z%G dd„ dƒZ&dKej'ej'eeeeeej' f dœdd „Z(dLej'eed!œd"d#„Z)dMeeeej'ej'ej'ef d&œd'd(„Z*G d)d*„ d*ƒZ+dNd+d,„Z,e-d-kr~ddl.Z.dZ/e+j0 1¡ D ]\Z2Z3e/d.e2 d/ e3 7 Z/qòe&j0 1¡ D ]\Z2Z3e/d0e2 d/ e3 7 Z/qe.j4d1e.j5d2e/ d3Z6e6j7Z8e8d4de dd5d6Z9e8d7d8d9d: e8d;edd<d= e6 :¡ Z;ee;j<ƒZ<e< =¡ j>j?dksºe. @e9d>¡‚n"e<j? Ad?¡d ZBe Cdd@eB›¡ e< D¡ r>dAe e<ƒkr>dBe<› dCZEe eFeEƒƒdDkr(e CddE¡ e G¡  ne CddF¡ e He<¡ e CdeB¡ e Cde<¡ dGeBkrre,eBe<e;jIdH ne#eBe<dI dS )Oz¨Simulate Data

Simulate stochastic dynamic systems to model gene expression dynamics and
cause-effect data.

TODO
----
Beta Version. The code will be reorganized soon.
é    N)ÚPath)ÚMappingProxyType)ÚOptionalÚUnionÚListÚTupleÚMapping)ÚAnnDataé   )Ú_utilsÚ	readwriteÚlogging)Úsettings)ÚLiteralT)Ú
krumsiek11Útoggleswitch)ÚmodelÚparams_fileÚtmaxÚ	branchingÚnrRealizationsÚnoiseObsÚnoiseDynÚstepÚseedÚwritedirÚreturnc
                 C   sj   t ƒ }
|rRt| ƒ d¡j}ddlm} t|jƒj|› d }t 	|¡}t
 ||
¡}
tf |
Ž}d|jd< |S )a      Simulate dynamic gene expression data [Wittmann09]_ [Wolf18]_.

    Sample from a stochastic differential equation model built from
    literature-curated boolean gene regulatory networks, as suggested by
    [Wittmann09]_. The Scanpy implementation is due to [Wolf18]_.

    Parameters
    ----------
    model
        Model file in 'sim_models' directory.
    params_file
        Read default params from file.
    tmax
        Number of time steps per realization of time series.
    branching
        Only write realizations that contain new branches.
    nrRealizations
        Number of realizations.
    noiseObs
        Observatory/Measurement noise.
    noiseDyn
        Dynamic noise.
    step
        Interval for saving state of system.
    seed
        Seed for generation of random numbers.
    writedir
        Path to directory for writing output files.

    Returns
    -------
    Annotated data matrix.

    Examples
    --------
    See this `use case <https://github.com/scverse/scanpy_usage/tree/master/170430_krumsiek11>`__
    Ú r
   ©Ú
sim_modelsz_params.txtr   Ziroot)Úlocalsr   Úwith_suffixÚnamer   r   Ú__file__Úparentr   Zread_paramsr   Zupdate_paramsÚsample_dynamic_dataÚuns)r   r   r   r   r   r   r   r   r   r   ÚparamsÚ	model_keyr   Z	pfile_simZdefault_paramsÚadata© r*   úJ/home/sam/Atlas/atlas_env/lib/python3.8/site-packages/scanpy/tools/_sim.pyÚsim   s    2


r,   c                 C   s"   dddt ddœi}t | |¡} | S )zi
    Update parser with tool specific arguments.

    This overwrites was is done in utils.uns_args.
    z--opfiler   Úfz=Specify a parameter file (default: "sim/${exkey}_params.txt"))ÚdefaultÚmetavarÚtypeÚhelp)Ústrr   Úadd_args)ÚpZ	dadd_argsr*   r*   r+   r3   [   s    üÿr3   c               
   K   sü  t | d ƒ d¡j}|  d¡}|dkr6tj|d  }nt |ƒ}|jddd t |d | ¡ | d	 }| d
 }| d }| d }| d }| d }d}	d}
d}d|krât	|| d}t
 |	¡}t|	ƒD ]ð}d|kr–t
 |j¡}tdƒD ]ž}d}t|jƒD ]l}t|jƒD ]\}||krRt
j ¡ dk r&dnd|||f< ||||f dkrJdnd7 }nd|||f< qqöttj |¡d ƒdk rä q„qä|||< | |¡ d}t
j |j¡}g }t||
 ƒD ]ø}d|kròt
 dd„ t|jƒD ƒ¡dt
j |j¡  }|j|||d}d}|rt|||ƒ\}}|rX|d7 }|j|dd|… |||dkrLdnd||d d|jkr |d kr |j|t
j d!d"¡ |||dkr”dnd||d ||kr¸ qÀq¸qÀt d#| ¡ › d$|j|jd  d% › ¡ nÒd&}d'}d(}t	|||| d)}g }t|ƒD ]¦}|d*kr:| |¡}|dkr:| ¡  q
d}t||
 ƒD ]b}|d*kr~|d+t
j |¡ d,t
 |¡   }ndt
j |¡ d- }|d.kr²t
 d-d-ddddg¡}|d/krâd0t
j |¡ d1 }t
 d2¡|d d…< d|kr4t
 |¡}d3||j d4 < d3||j d5 < d3||j d6 < |d7t
j |¡ 7 }|j||||d8}d}|r`t|||ƒ\}}|rœ|d7 }|j|dd|… |||dkrdnd||d ||krJ q
qJq
d}| !d9¡D ]}qÂt "d:|› ¡ tj#|ddd;}|| |j$d<< |S )=z
    Helper function.
    r   r   r   NZ_simT©ÚparentsÚexist_okz
params.txtr   r   r   r   r   r   é   iè  r   )r   r'   Úhillé
   r   çš™™™™™Ù?çffffffæ?r   c                 S   s   g | ]}d ‘qS )çš™™™™™é?r*   ©Ú.0Úir*   r*   r+   Ú
<listcomp>¨   s     z'sample_dynamic_data.<locals>.<listcomp>g{®Gáz„?)r   ÚX0r   F)Údirr   Úappendr   r   Úzerosr
   éô  é   zmean nr of offdiagonal edges z compared to total nr g       @Úrandomé   é   )ÚdimÚinitTyper   r'   Úbranchgš™™™™™©?çš™™™™™™?g333333Ó?)rJ   é   )é   é   é	   r:   g333333ã?çš™™™™™É?é   r=   ZGata2zPu.1ZCebpagü©ñÒMbP?)rB   r   Úrestartzsim*.txtzreading simulation results )Zfirst_column_namesZsuppress_cache_warningZ
tmax_write)%r   r!   r"   Úgetr   r   Úmkdirr   Zwrite_paramsÚGRNsimÚnprE   ÚrangeÚarrayÚCouplrK   rH   ÚrandÚmaxÚspZlinalgZeigÚ	set_couplÚrandnÚ	sim_modelÚ_check_branchingÚ
write_dataÚloggÚdebugÚmeanÚbranch_init_model1ÚonesÚvarNamesÚglobÚinfoÚ_readr&   )r'   r(   r   r   r   r   r   r   r   Z	nrSamplesZmaxRestartsZmaxNrSamplesZgrnsimZnrOffEdges_listÚsampler\   ZsampleCouplZ
nrOffEdgesÚgpÚgÚrealrB   ÚXsamplesrU   ÚXÚcheckrL   rK   ÚX0meanÚfilenamer)   r*   r*   r+   r%   n   sü    




  

ÿþú	ú

"ÿ



$



ú
  ÿr%   úsim/testFr   r8   )rj   Ú	boolRulesc                 C   sŒ  |j ddd |d }| ¡ rN| d¡}t| ¡ ƒ|r<dnd }W 5 Q R X nd}| d¡}d |¡}| t|ƒ¡ W 5 Q R X | jd }|sÞ|rÞ|jdkrÞ|d	|› d
  d¡}| d¡ | d¡ | d¡ | d¡ | d¡ | d|	 d ¡ | dt|
ƒ d ¡ | d¡ | d¡ | ¡ D ] \}}| |› d|› d¡ q2| d¡ t| ¡ ƒ}t
|ƒD ]`}t
|ƒD ]P}t	 |||f ¡dkr~| || d›d|| d›d|||f d›d¡ q~qrW 5 Q R X |r|dd›d7 }| ¡ D ]}||d›d7 }qü|d|› d
  |r.dnd ¡N}t	j|t	jt	 d| jd ¡| f |r`d!n|d"gd#d$„ t
|ƒD ƒ d% W 5 Q R X d&S )'z`Write simulated data.

    Accounts for saving at the same time an ID
    and a model file.
    Tr5   zid.txtÚrr   r8   Úwz{:0>6}Zmodel_ú.txtz<# For each "variable = ", there must be a right hand side: 
z?# either an empty string or a python-style logical expression 
z4# involving variable names, "or", "and", "(", ")". 
z## The order of equations matters! 
z# 
z# modelType = Ú
z# invTimeStep = z# boolean update rules: 
ú = z# coupling list: 
ç»½×Ùß|Û=Ú10ú z10.3ú 
Úitz>2z>7Zsim_ÚabÚwbr   z%4.fc                 S   s   g | ]}d ‘qS )z%7.4fr*   r>   r*   r*   r+   rA   d  s     zwrite_data.<locals>.<listcomp>)ÚheaderÚfmtN)rW   Úis_fileÚopenÚintÚreadÚformatÚwriter2   ÚshaperY   rZ   ÚsizeÚitemsÚlistÚkeysÚabsZsavetxtZc_Úarange)rs   rC   rD   r…   rj   ÚAdjr\   rx   r   Ú	modelTypeÚinvTimeSteprv   r-   ÚidrK   ÚkÚvÚnamesro   rp   r*   r*   r+   rd     sZ    $


ÿ





*ÿ"ürd   c                	   @   s  e Zd ZdZeddddZdZddd	d
dddei ƒfdd„ZdBdd„Z	dd„ Z
dCdd„ZdDdd„ZdEdd„ZdFdd„ZdGdd„Zd d!„ ZdHd"d#„Zd$d%„ Zd&d'„ Zd(d)„ Zd*d+„ Zd,d-„ Zd.d/„ Zd0d1„ Zd2d3„ ZdId5d6„Zd7d8„ Zd9d:„ Zd;d<„ Zed=ƒd>ddd?dfd@dA„ZdS )JrX   zž
    Simlulation of stochastic dynamic systems.

    Main application: simulation of gene expression dynamics.

    Also standard models are implemented.
    zŒmyeloid progenitor network, Krumsiek et al., PLOS One 6, e22649, 
      equations from Table 1 on page 3, doi:10.1371/journal.pone.0022649 
zvector autoregressive process 
zprocess with hill kinetics 
)r   Úvarr9   TrG   Zex0r›   rH   Fr   Nc	                 C   sÜ   |dkr|n|j d | _d| _d| _|| _|| _|| _|| _|| _|dkrRt	dƒ‚|| j
 ¡ krdd}	dd	lm}
 t|
jƒj|› d
 }| ¡ sžt	d|› dƒ‚|| _| j|d tj |d ¡ d| jj d | _|| _dS )zª
        Params
        ------
        model
            either string for predefined model,
            or directory with a model file and a couple matrix files
        Nr   r8   r;   )rM   rH   z'initType must be either: branch, randomz#model not among predefined models 
r
   r   r{   zModel file z does not exist)r\   r   zmodel = r   )r   rK   ÚmaxnparÚp_indepr   r•   rL   ÚshowÚ	verbosityÚRuntimeErrorÚavailModelsr‘   r   r   r   r#   r$   r‡   r`   rY   rH   r   r"   r…   r'   )ÚselfrK   r   r•   rL   rž   rŸ   r\   r'   Úmessager   r*   r*   r+   Ú__init__}  s,    ÿzGRNsim.__init__c                 C   sÂ   || _ t || jf¡}||tj | j¡  |d< td|ƒD ]‚}| jdkr\|  ||d  ¡}n.| jdkrz|  	||d  ¡}nt
d| j›ƒ‚||d  | ||< ||  |tj | j¡ 7  < q:|S )zSimulate the model.r   r8   r9   r›   zUnknown modelType )r   rY   rE   rK   rH   ra   rZ   r•   Ú
Xdiff_hillÚ	Xdiff_varÚ
ValueError)r¢   r   rB   r   rU   rs   ÚtÚXdiffr*   r*   r+   rb   ®  s    

 zGRNsim.sim_modelc                 C   s˜  | j dko| j}d| _t | j¡}t| j ¡ ƒD ]`\}}| j| r0d}|dkrXd}nq0t| j| ƒD ]Ö\}}	d}
d}t|	ƒD ]’\}}| j	| j| |  }|| }dt 
| j||f ¡ }|
|rÌ|  ||¡n
|  ||¡9 }
|dkr~||rìdnd› d| j| | › d	|d
›d7 }q~||
7 }|dkrf||dkr2dnd| 7 }qf| j|||   ||< |dkr0|› d|› d|› d| j› d|› d|› d}t d|¡ q0|S )a&  Build Xdiff from coefficients of boolean network,
        that is, using self.boolCoeff. The employed functions
        are Hill type activation and deactivation functions.

        See Wittmann et al., BMC Syst. Biol. 3, 98 (2009),
        doi:10.1186/1752-0509-3-98 for more details.
        r   Fr   r8   çš™™™™™¹?Úar@   ú(z, z.2ú)ú+Ú_ú-r}   z*()rŸ   ÚwriteOutputOncerY   rE   rK   Ú	enumerateÚpasr‘   Ú	boolCoeffrj   r’   r\   Úhill_aÚhill_ir–   r   Úm)r¢   ÚXtrŸ   r©   ÚichildÚchildZ	Xdiff_synZXdiff_syn_strZitupleÚtupleZXdiff_syn_tupleZXdiff_syn_tuple_strZivr™   ZiparentÚxÚ	thresholdZ	Xdiff_strr*   r*   r+   r¥   Á  s@    
ÿ(ÿ(ÿzGRNsim.Xdiff_hillc                 C   s   | }|t  | j|¡7 }|S ©r   )rY   Údotr\   )r¢   r¸   rŸ   r©   r*   r*   r+   r¦   ó  s    zGRNsim.Xdiff_varrª   r
   c                 C   s$   t  ||¡}t  ||¡}|||  S )zActivating hill function.©rY   Úpower©r¢   r¼   r½   rÁ   Úx_powÚthreshold_powr*   r*   r+   rµ   û  s    zGRNsim.hill_ac                 C   s$   t  ||¡}t  ||¡}|||  S )z^Inhibiting hill function.

        Is equivalent to 1-hill_a(self,x,power,threshold).
        rÀ   rÂ   r*   r*   r+   r¶     s    zGRNsim.hill_ic                 C   s,   t  ||¡}t  ||¡}|||  d|  S )z$Normalized activating hill function.r8   rÀ   )r¢   r¼   r½   rÁ   r¹   rÃ   rÄ   r*   r*   r+   Únhill_a
  s    zGRNsim.nhill_ac                 C   s,   t  ||¡}t  ||¡}|||  d|  S )zjNormalized inhibiting hill function.

        Is equivalent to 1-nhill_a(self,x,power,threshold).
        r8   rÀ   rÂ   r*   r*   r+   Únhill_i  s    zGRNsim.nhill_ic           
      C   sÚ  | j dkrt dd| j¡ g }| j ¡ D ]Ê}| d¡rtd|krt|}d|kr`| d¡dd… \}}| d¡d	  ¡ | _| d¡rÀd
|krÀ|}d|kr¨| d¡dd… \}}t	| d¡d	  ¡ ƒ| _
| d¡sä| dd„ | d¡D ƒ¡ | d¡r( qôq(t|ƒ| _t|ƒ| _dd„ t| j ¡ ƒD ƒ| _| j}t | j| jf¡| _d}| j ¡ D ]`}| d¡r^d}|rhqJ| d¡sJ| ¡  ¡ \}}}	t	|	ƒ| jt|| ƒt|| ƒf< qJt | j¡| _t t | j¡¡| _|  ¡  dS )z5Read the model and the couplings from the model file.r   zreading modelú#zmodelType =ú|Nr
   ú=r8   zinvTimeStep =c                 S   s   g | ]}|  ¡ ‘qS r*   )Ústrip©r?   Úsr*   r*   r+   rA   +  s     z%GRNsim.read_model.<locals>.<listcomp>z# coupling list:c                 S   s   i | ]\}}||“qS r*   r*   )r?   r@   rÌ   r*   r*   r+   Ú
<dictcomp>0  s      z%GRNsim.read_model.<locals>.<dictcomp>TF)rŸ   r   r·   r   rˆ   Ú
startswithÚsplitrÊ   r•   Úfloatr–   rD   ÚlenrK   Údictrx   r²   r‘   rj   rY   rE   r\   r‰   ÚsignÚ
Adj_signedr’   r[   r”   Úbuild_boolCoeff)
r¢   rx   ÚlineÚkeyvalr0   rš   ZboolContinueZgpsÚgsÚvalr*   r*   r+   Ú
read_model  sJ    




þ&zGRNsim.read_modelc                    s  dd„ t | jƒD ƒ| _| j| j ¡ kr:|dkr:|  ¡  n¾d| jjkrü|| _dd„ | j ¡ D ƒ| _	t
| j ¡ ƒ}t | jƒD ]|}g }t | jƒD ]*}t | j||f dk¡rŒ| || ¡ qŒd |dd… d	d
„ |dd… D ƒ ¡| j	|| < t |¡| _qznü| jdkrjt | j| jf¡| _d}t ddg¡‰ t |¡| jˆ ˆ f< t ˆ ¡}t ‡ fdd
„t | jƒD ƒ¡}t|ƒdkrøtjjt dt|ƒ¡ddd}	||	 }
tjjt dt|ƒ¡ddd}|| }t d|
|¡ t d¡| j||
f< | jdkrd| j|d |d f< | jdkr@d| j|d |d f< d| j|d |d f< t ||	¡}t ||¡}t ||¡}qjnŽt | j| jf¡| _t | jƒD ]n}tj d| j¡}|dkrætj d| jd ¡}tjjt d| j¡|dd}d| j||f< nd| j||f< qˆt t | j¡¡| _dS )zmConstruct the coupling matrix (and adjacancy matrix) from predefined models
        or via sampling.
        c                 S   s   i | ]}t |ƒ|“qS r*   )r2   r>   r*   r*   r+   rÍ   J  s      z$GRNsim.set_coupl.<locals>.<dictcomp>Nr›   c                 S   s   i | ]
}|d “qS r¾   r*   rË   r*   r*   r+   rÍ   P  s      r~   r   r8   c                 S   s   g | ]}d | ‘qS )z or r*   ©r?   Úpar*   r*   r+   rA   X  s     z$GRNsim.set_coupl.<locals>.<listcomp>)Ú6Ú7Ú8Ú9r   r
   r   c                    s   g | ]}|ˆ kr|‘qS r*   r*   r>   ©Z	sinknodesr*   r+   rA   h  s      F)rŽ   ÚreplacerQ   éÿÿÿÿ)rR   r:   ) rZ   rK   rj   r   r¡   r‘   rÚ   r"   r\   rx   r   rY   r’   rD   ÚjoinrÓ   rÔ   rE   r[   ri   rÑ   rH   Úchoicer“   r   r·   Údeleter”   Zbinomialr   Úrandintrœ   )r¢   r\   rš   ro   r³   rp   Zn_sinknodesZ	leafnodesZ
availnodesZ
parent_idxr$   Zchildren_idsÚchildrenr@   ZindepÚnrZj_parr*   rá   r+   r`   F  sv     ÿ
  ÿ  ÿ
  ÿzGRNsim.set_couplc                 C   s  | j dks| j dkrdS t | j| jf¡| _t| jjd ƒD ]j}t| j| ƒD ]V\}}|dkrNtj	 
d¡}dtj	 ¡  d | j||f< |dkrN| j||f  d	9  < qNq<| j dkr¼|  ¡  n:| j d
krÐ|  ¡  n&| j dkrä|  ¡  n| j dkrö|  ¡  | jdkrt d| j¡ dS )z5Using the adjacency matrix, sample a coupling matrix.r   r›   Nr   r
   ç        rª   r8   rã   rJ   )rO   rP   )rQ   rR   r:   )r   rY   rE   rK   r\   rZ   r”   r   r²   rH   rç   r]   Úcoupl_model1Úcoupl_model5Úcoupl_model6Úcoupl_model8rŸ   r   r·   )r¢   r@   Újr«   Zco_antir*   r*   r+   Úset_coupl_old˜  s(    






zGRNsim.set_coupl_oldc                 C   sH   t  | jd ¡| jd< t  | jd ¡ | jd< t  | jd ¡| jd< dS )zŽIn model 1, we want enforce the following signs
        on the couplings. Model 2 has the same couplings
        but arbitrary signs.
        ©r   r   ©r   r8   )r8   r8   N)rY   r’   r\   ©r¢   r*   r*   r+   rë   ¸  s    zGRNsim.coupl_model1c                 C   sX   d| j  | _| jd  d9  < | jd  d9  < | jd  d9  < | jd  d9  < dS )zToggle switch.çš™™™™™É¿©r
   r   rã   )rG   r   )rT   r8   )rJ   r8   N)r”   r\   ró   r*   r*   r+   rì   Á  s
    zGRNsim.coupl_model5c                 C   s   d| j  | _dS )úVariant of toggle switch.ç      à?N©rÔ   r\   ró   r*   r*   r+   rí   É  s    zGRNsim.coupl_model6c                 C   s8   d| j  | _tj| jdgdD ]}|dk rd|d< qdS )rö   r÷   r   )Zop_flagsgíµ ÷Æ°¾rô   .N)rÔ   r\   rY   Znditer)r¢   r¼   r*   r*   r+   rî   Í  s    zGRNsim.coupl_model8c                 C   s   | j | _dS )rö   Nrø   ró   r*   r*   r+   Úcoupl_model_krumsiek11Ö  s    zGRNsim.coupl_model_krumsiek11c                 C   s   | | |   |¡ S )z?Yields zero when solved for X_t
        given X_{t+1}.
        )r©   )r¢   r¸   ZXt1r*   r*   r+   Úsim_model_back_helpÚ  s    zGRNsim.sim_model_back_helpc                 C   sf   t  || jf¡}|||d < t|d ddƒD ]4}tjj| j||d  ||d  dd}|j||< q,|S )z%Simulate the model backwards in time.r8   r
   rã   Zhybr)ÚargsÚmethod)	rY   rE   rK   rZ   r_   ÚoptimizeÚrootrú   r¼   )r¢   r   rB   rs   r¨   Zsolr*   r*   r+   Úsim_model_backwardsà  s     
 
 ÿzGRNsim.sim_model_backwardséd   c                 C   sª  t  | jd | jd  dg¡}|d dks6|d dk rFt dd¡ d S | j|d |t  d	d
g¡ d}| j|d |t  d
d
g¡ d}| j||d d}| j||d d}d|d |d   }t  |¡dk sÚt  |¡dkrêt dd¡ d S | j	r¦| j
dkr¦t ¡  t |d d …df d|d d …df d¡ t |d d …df d|d d …df d¡ t |d d …df d|d d …df d¡ t |d d …df d|d d …df d¡ |S )Nrò   rñ   r8   r   g
×£p=
ï?g¸…ëQ¸ž?zP... either no fixed point in [0,1]^2! 
    or fixed point is too close to boundsrG   g{®Gáz”?g{®Gáz”¿)r   rB   r÷   rN   g333333ï?z(... initial point is too close to boundsz.bz.gÚbrp   )rY   r[   r\   r   r·   rÿ   rb   Úminr^   rž   rŸ   ÚplZfigureÚplot)r¢   r   ZXfixZXbackUpZXbackDoZXupZXdoru   r*   r*   r+   rh   ë  s8    þ ÿ ÿ((((zGRNsim.branch_init_model1c                 C   sÐ   |  dd¡  dd¡  dd¡  dd¡  dd¡}| ¡ }|s<g S g }g }|D ]n}|| j ¡ kr¤t dd¡ t dt| j ¡ ƒ¡ d	| d
 | d d d d }t|ƒ‚||krH| |¡ qH|D ]}| 	|¡ q¼|S )zYDetermine parents based on boolean updaterule.

        Returns list of parents.
        r¬   r   r­   ÚorÚandÚnotr   zlist of available variables:zprocessing of rule "z yields an invalid parent: z) | check whether the syntax is correct: 
z1only python expressions "(",")","or","and","not" zAare allowed, variable names and expressions have to be separated zby white spaces)
râ   rÏ   rj   r‘   r   r·   r   r§   rD   Úremove)r¢   ÚruleÚrule_paZpa_oldZ	pa_deleterÜ   r£   r*   r*   r+   Úparents_from_boolRule  sZ     ÿ þ ý üÿÿþýüûúùÿ
zGRNsim.parents_from_boolRulec                    s  dd„ ˆ j  ¡ D ƒˆ _dd„ ˆ j  ¡ D ƒˆ _ˆ j ¡ D ]R}ˆ j| }ˆ  |¡ˆ j|< ‡ fdd„ˆ j| D ƒ}tˆ jƒD ]r}||kr¸t 	ˆ j
ˆ j | |f ¡dk rêtd|› d|› ƒ‚qxt 	ˆ j
ˆ j | |f ¡dkrxtd	|› d|› ƒ‚qxˆ jd
kr t dd| ¡ t d|¡ t dt¡ ttjddgtˆ j| ƒdƒD ],}ˆ  |ˆ j| |¡r@ˆ j|  |¡ q@ˆ jd
kr6t dˆ j| ¡ q6dS )z%Compute coefficients for tuple space.c                 S   s   i | ]
}|g “qS r*   r*   rË   r*   r*   r+   rÍ   8  s      z*GRNsim.build_boolCoeff.<locals>.<dictcomp>c                 S   s   i | ]
}|g “qS r*   r*   rË   r*   r*   r+   rÍ   :  s      c                    s   g | ]}ˆ j | ‘qS r*   )rj   rÛ   ró   r*   r+   rA   ?  s     z*GRNsim.build_boolCoeff.<locals>.<listcomp>r~   zspecify coupling value for z <- z&there should be no coupling value for r8   r   z...FT)ÚrepeatN)rj   r‘   r´   r³   rx   r  rZ   rK   rY   r’   r\   r§   rŸ   r   r·   r
  r   Ú	itertoolsÚproductrÑ   Úprocess_rulerD   )r¢   Úkeyr	  Z
pasIndicesrp   r»   r*   ró   r+   rÕ   5  s2    
ÿÿ
zGRNsim.build_boolCoeffc                 C   s.   t |ƒD ]\}}| || t|ƒ¡}qt|ƒS )z-Process a string that denotes a boolean rule.)r²   râ   r2   Úeval)r¢   r	  rÜ   r»   r@   r™   r*   r*   r+   r  X  s    zGRNsim.process_rulerw   rê   r8   c           
      C   sÖ   | j }t|jd ƒ}	|dt|	ƒ d 7 }|dt|ƒ d 7 }|dt|ƒ d 7 }|dt|ƒ d 7 }|dt| jƒ d 7 }|dt|ƒ d 7 }||tj |	| j¡ 7 }t	||||| j
| j| j| j| j| j| jd	 d S )
Nr   ztmax = r|   zbranching = znrRealizations = znoiseObs = znoiseDyn = zseed = )rj   r”   r\   r   r•   rx   r–   )r…   r‰   r   r2   r   rY   rH   ra   rK   rd   rj   r”   r\   r   r•   rx   r–   )
r¢   rs   rC   r   rD   r   r   r   r…   r   r*   r*   r+   rd   ^  s,    
õzGRNsim.write_data)r   r   )r   )rª   r
   )rª   r
   )rª   r
   r
   )rª   r
   )N)r   ) Ú__name__Ú
__module__Ú__qualname__Ú__doc__rÒ   r¡   r±   r   r¤   rb   r¥   r¦   rµ   r¶   rÅ   rÆ   rÚ   r`   rð   rë   rì   rí   rî   rù   rú   rÿ   rh   r  rÕ   r  r   rd   r*   r*   r*   r+   rX   h  sX   ù
÷
1
2


	

	-
R 		
"(#	ørX   ç      Ð?)rs   rr   rU   r½   r   c                 C   sž   d}t |ƒ}|dkr | | ¡ nV|D ]B}t | ddd…f |ddd…f  ¡}t |d¡d |k r$d}q$|rv| | ¡ t d|› d|rŠd	nd
› d¡ ||fS )a      Check whether time series branches.

    Parameters
    ----------
    X
        current time series data.
    Xsamples
        list of previous branching samples.
    restart
        counts number of restart trials.
    threshold
        sets threshold for attractor identification.

    Returns
    -------
    check
        true if branching realization
    Xsamples
        updated list
    Tr   rã   NéþÿÿÿFzrealization ú: r   Únoz new branch)r   rD   rY   ÚabsoluteÚ	partitionre   rf   )rs   rr   rU   r½   rt   ZXcompareZ
Xtmax_diffr*   r*   r+   rc   ‚  s    &
 rc   )r”   rŸ   r   c              
   C   s†   | j d }t|ƒD ]n}t |¡}d||< t|ƒD ]N}|  |¡}|| dkr0|dkrvt d| ¡ t dd|d d|d¡   dS q0qd	S )
zò    Checks that there are no cycles in graph described by adjacancy matrix.

    Parameters
    ----------
    Adj
        adjancancy matrix of dimension (dim, dim)

    Returns
    -------
    True if there is no cycle, False otherwise.
    r   r8   r~   r
   zcontains a cycle of lengthzstarting from nodez	-> rejectFT)r   rZ   rY   rE   r¿   r   r·   )r”   rŸ   rK   rp   r™   r@   r*   r*   r+   Úcheck_nocycles¯  s&    


úr  rG   r÷   )rK   Úconnectivityr   c                 C   sÒ   d}d}t |ƒD ] }t | | f¡}d}t | ƒD ]B}t | ƒD ]4}||krHq:tj ¡ d| k r:d|||f< |d7 }q:q.tj| | fdd}	t |¡}	t |	¡}
t|
ƒr|dkrd	} q²q|sÆtd
|› dƒ‚||
|	|fS )a      Sample coupling matrix.

    Checks that returned graphs contain no self-cycles.

    Parameters
    ----------
    dim
        dimension of coupling matrix.
    connectivity
        fraction of connectivity, fully connected means 1.,
        not-connected means 0, in the case of fully connected, one has
        dim*(dim-1)/2 edges in the graph.

    Returns
    -------
    coupl
        coupling matrix
    adj
        adjancancy matrix
    adj_signed
        signed adjacancy matrix
    n_edges
        Number of edges
    r:   Fr   r÷   r<   r8   Úint_)ZdtypeTz'did not find graph without cycles afterz trials)	rZ   rY   rE   rH   r]   rÓ   r’   r  r§   )rK   r  Z	max_trialrt   Ztrialr\   Ún_edgesro   rp   rÔ   r”   r*   r*   r+   Úsample_coupling_matrixÑ  s.    


ÿr   c                	   @   sH   e Zd ZdZedddddddd	Zd
d„ Zdejdœdd„Z	dd„ Z
dS )ÚStaticCauseEffectzB
    Simulates static data to investigate structure learning.
    u	   y = Î±x 
zy = noise 
z	y = |x| 
u   y = Î±xÂ² 
zy = x - |x| 
zy = tanh(x) 
zcombinatorial regulation 
)rÖ   ÚnoiseÚabslineÚparabolaÚsawtoothÚtanhÚcombic                 C   s0   t dd„ dd„ tjdd„ dd„ dd„ d| _d S )Nc                 S   s   | S )Nr*   ©r¼   r*   r*   r+   Ú<lambda>  ó    z,StaticCauseEffect.__init__.<locals>.<lambda>c                 S   s   dS )Nr   r*   r(  r*   r*   r+   r)    r*  c                 S   s   | d S ©Nr
   r*   r(  r*   r*   r+   r)     r*  c                 S   s   d|  t  d|  ¡ S )Nr÷   )rY   Úfloorr(  r*   r*   r+   r)  !  r*  c                 S   s   t  d|  ¡S r+  )rY   r&  r(  r*   r*   r+   r)  "  r*  )rÖ   r"  r#  r$  r%  r&  )rÒ   rY   r’   Úfuncsró   r*   r*   r+   r¤     s    úzStaticCauseEffect.__init__rÖ   ©r”   c              
   C   s"  t dddddg}d}d}d}| j| }d}|jd }	t ||	f¡}
d}tt|	ƒƒ}g }t|	ƒD ]z}||d	d	…f  ¡ |kr`|d
kr tj 	d||¡|
d	d	…|f< |dkrÆtj 
| ||¡|
d	d	…|f< | |¡ | |¡ q`g }t |	¡}t|ƒ|d< td|	ƒD ]F}|D ]:}||d	d	…f  ¡ |kr| |¡ ||  d7  < qq |d dkrŒ||d |d f dkrŒ|d }|d |d< ||d< |D ]Œ}t|	ƒD ]X}|||f dkrœ|
d	d	…|f  d||d	d	…f  ¡  ||
d	d	…|f ƒ 7  < qœ|
d	d	…|f  tj 	d||¡7  < q|
S )af          Simulate data given only an adjacancy matrix and a model.

        The model is a bivariate funtional dependence. The adjacancy matrix
        needs to be acyclic.

        Parameters
        ----------
        Adj
            adjacancy matrix of shape (dim,dim).

        Returns
        -------
        Data array of shape (n_samples,dim).
        r%  ÚuniformçÍÌÌÌÌÌü?rª   )ÚfuncZgdistÚ
sigma_globÚsigma_noiser   r;   r   NZgaussianr8   g      ð?)rÒ   r-  r   rY   rE   r   rZ   ÚsumrH   Únormalr/  rD   r  rÑ   )r¢   r”   r   ZexamplesÚ	n_samplesr2  r3  r1  Z
sourcedistrK   rs   Znrparrè   r6   ro   Zchildren_sortedZnrchildren_parr1   rp   r*   r*   r+   Úsim_givenAdj%  sX    üÿ





D(zStaticCauseEffect.sim_givenAdjc                 C   s
  d}d}t  |df¡}t j | ||¡|dd…df< t j | ||¡|dd…df< | jd }||dd…df ƒtjj |dd…df dd¡ |dd…d	f< t	j
|dd…df |dd…df |dd…d	f d
d t	 ¡  t	 |dd…df |dd…d	f d¡ t	 ¡  |S )z(Simulate data to model combi regulation.rF   r0  rG   Nr   r8   r&  rS   r
   Zface)ÚcZ	edgecolorÚ.)rY   rE   rH   r/  r-  r_   ÚstatsZnormZcdfr  Zscatterrž   r  )r¢   r6  r2  rs   r1  r*   r*   r+   Ú	sim_combi‚  s"    
>   ÿ&zStaticCauseEffect.sim_combiN)rÖ   )r  r  r  r  rÒ   r¡   r¤   rY   Úndarrayr7  r;  r*   r*   r*   r+   r!    s   ù
]r!  c                 C   sä   d}d}d}|   dd¡} tj d¡ | dkr®t |¡}t|ƒD ]Z}t||ƒ\}}	}
}|dkrtt d|¡ t d|	¡ |||< t	ƒ  
|	| ¡}t|||	d	 q>t dd
| ¡ ¡ n2t	ƒ  ¡ }t d¡}	d |	d< |	d< t|||	d	 d S )Nr=   rG   é2   zstatic-r   r   r'  r8   r.  zmean edge number:)rG   rG   rõ   )r
   r8   )râ   rY   rH   r   rE   rZ   r   r   r·   r!  r7  rd   rg   r;  )r   rC   rŸ   r  rK   Zn_Couplsr  Zicouplr\   r”   rÔ   Zn_ers   r*   r*   r+   Úsample_static_data¡  s(    


r>  Ú__main__z    static-r  z    z\Simulate stochastic discrete-time dynamical systems,
in particular gene regulatory networks.zw  MODEL: specify one of the following models, or one of 
    the filenames (without ".txt") in the directory "models" 
)ÚdescriptionÚformatter_classÚepilogz--dirzgspecify directory to store data,  must start with "sim/MODEL_...", see possible values for MODEL below )Úrequiredr0   r.   r1   z--showÚ
store_truez
show plots)Úactionr1   z--verbosityz2specify integer > 0 to get more output [default 0])r0   r.   r1   zBThe parent directory of the --dir argument needs to be named 'sim'r¯   z...model is: Útestz
directory z; already exists, remove it and continue? [y/n, press enter]Úyz    ...quit program executionz*   ...removing directory and continuing...Zstatic)r   rC   rŸ   )r   rC   )	TNNNNNNNN)r  )r
   )rG   r÷   )r   )Jr  r  ÚshutilÚsysÚpathlibr   Útypesr   Útypingr   r   r   r   r   ÚnumpyrY   Zscipyr_   Zanndatar	   r   r   r   r   re   Z	_settingsr   Z_compatr   Úboolr‰   rÐ   r2   r,   r3   r%   r[   rd   rX   r<  rc   r  r   r!  r>  r  ÚargparserB  r¡   r   r˜   r™   ÚArgumentParserÚRawDescriptionHelpFormatterr4   Úadd_argumentZaaZdir_argÚ
parse_argsrû   rC   Úresolver$   r"   ÚArgumentErrorrÏ   r   r·   Úis_dirr£   ÚinputÚexitÚrmtreerŸ   r*   r*   r*   r+   Ú<module>   sü   
         öõ? $õ

øY     ÿ   þ-#   ÿ þ: 

ÿùú
ü
þ
ÿ


